from __future__ import division import numpy as np from scipy.ndimage import gaussian_filter, gaussian_laplace import math from math import sqrt, log from scipy import spatial from ..util import img_as_float from .peak import peak_local_max from ._hessian_det_appx import _hessian_matrix_det from ..transform import integral_image from .._shared.utils import assert_nD # This basic blob detection algorithm is based on: # http://www.cs.utah.edu/~jfishbau/advimproc/project1/ (04.04.2013) # Theory behind: http://en.wikipedia.org/wiki/Blob_detection (04.04.2013) def _compute_disk_overlap(d, r1, r2): """ Compute surface overlap between two disks of radii ``r1`` and ``r2``, with centers separated by a distance ``d``. Parameters ---------- d : float Distance between centers. r1 : float Radius of the first disk. r2 : float Radius of the second disk. Returns ------- vol: float Volume of the overlap between the two disks. """ ratio1 = (d ** 2 + r1 ** 2 - r2 ** 2) / (2 * d * r1) ratio1 = np.clip(ratio1, -1, 1) acos1 = math.acos(ratio1) ratio2 = (d ** 2 + r2 ** 2 - r1 ** 2) / (2 * d * r2) ratio2 = np.clip(ratio2, -1, 1) acos2 = math.acos(ratio2) a = -d + r2 + r1 b = d - r2 + r1 c = d + r2 - r1 d = d + r2 + r1 area = (r1 ** 2 * acos1 + r2 ** 2 * acos2 - 0.5 * sqrt(abs(a * b * c * d))) return area / (math.pi * (min(r1, r2) ** 2)) def _compute_sphere_overlap(d, r1, r2): """ Compute volume overlap between two spheres of radii ``r1`` and ``r2``, with centers separated by a distance ``d``. Parameters ---------- d : float Distance between centers. r1 : float Radius of the first sphere. r2 : float Radius of the second sphere. Returns ------- vol: float Volume of the overlap between the two spheres. Notes ----- See for example http://mathworld.wolfram.com/Sphere-SphereIntersection.html for more details. """ vol = (math.pi / (12 * d) * (r1 + r2 - d)**2 * (d**2 + 2 * d * (r1 + r2) - 3 * (r1**2 + r2**2) + 6 * r1 * r2)) return vol / (4./3 * math.pi * min(r1, r2) ** 3) def _blob_overlap(blob1, blob2): """Finds the overlapping area fraction between two blobs. Returns a float representing fraction of overlapped area. Parameters ---------- blob1 : sequence of arrays A sequence of ``(row, col, sigma)`` or ``(pln, row, col, sigma)``, where ``row, col`` (or ``(pln, row, col)``) are coordinates of blob and ``sigma`` is the standard deviation of the Gaussian kernel which detected the blob. blob2 : sequence of arrays A sequence of ``(row, col, sigma)`` or ``(pln, row, col, sigma)``, where ``row, col`` (or ``(pln, row, col)``) are coordinates of blob and ``sigma`` is the standard deviation of the Gaussian kernel which detected the blob. Returns ------- f : float Fraction of overlapped area (or volume in 3D). """ n_dim = len(blob1) - 1 root_ndim = sqrt(n_dim) # extent of the blob is given by sqrt(2)*scale r1 = blob1[-1] * root_ndim r2 = blob2[-1] * root_ndim d = sqrt(np.sum((blob1[:-1] - blob2[:-1])**2)) if d > r1 + r2: return 0 # one blob is inside the other, the smaller blob must die if d <= abs(r1 - r2): return 1 if n_dim == 2: return _compute_disk_overlap(d, r1, r2) else: # http://mathworld.wolfram.com/Sphere-SphereIntersection.html return _compute_sphere_overlap(d, r1, r2) def _prune_blobs(blobs_array, overlap): """Eliminated blobs with area overlap. Parameters ---------- blobs_array : ndarray A 2d array with each row representing 3 (or 4) values, ``(row, col, sigma)`` or ``(pln, row, col, sigma)`` in 3D, where ``(row, col)`` (``(pln, row, col)``) are coordinates of the blob and ``sigma`` is the standard deviation of the Gaussian kernel which detected the blob. This array must not have a dimension of size 0. overlap : float A value between 0 and 1. If the fraction of area overlapping for 2 blobs is greater than `overlap` the smaller blob is eliminated. Returns ------- A : ndarray `array` with overlapping blobs removed. """ sigma = blobs_array[:, -1].max() distance = 2 * sigma * sqrt(blobs_array.shape[1] - 1) tree = spatial.cKDTree(blobs_array[:, :-1]) pairs = np.array(list(tree.query_pairs(distance))) if len(pairs) == 0: return blobs_array else: for (i, j) in pairs: blob1, blob2 = blobs_array[i], blobs_array[j] if _blob_overlap(blob1, blob2) > overlap: if blob1[-1] > blob2[-1]: blob2[-1] = 0 else: blob1[-1] = 0 return np.array([b for b in blobs_array if b[-1] > 0]) def blob_dog(image, min_sigma=1, max_sigma=50, sigma_ratio=1.6, threshold=2.0, overlap=.5,): """Finds blobs in the given grayscale image. Blobs are found using the Difference of Gaussian (DoG) method [1]_. For each blob found, the method returns its coordinates and the standard deviation of the Gaussian kernel that detected the blob. Parameters ---------- image : 2D or 3D ndarray Input grayscale image, blobs are assumed to be light on dark background (white on black). min_sigma : float, optional The minimum standard deviation for Gaussian Kernel. Keep this low to detect smaller blobs. max_sigma : float, optional The maximum standard deviation for Gaussian Kernel. Keep this high to detect larger blobs. sigma_ratio : float, optional The ratio between the standard deviation of Gaussian Kernels used for computing the Difference of Gaussians threshold : float, optional. The absolute lower bound for scale space maxima. Local maxima smaller than thresh are ignored. Reduce this to detect blobs with less intensities. overlap : float, optional A value between 0 and 1. If the area of two blobs overlaps by a fraction greater than `threshold`, the smaller blob is eliminated. Returns ------- A : (n, image.ndim + 1) ndarray A 2d array with each row representing 3 values for a 2D image, and 4 values for a 3D image: ``(r, c, sigma)`` or ``(p, r, c, sigma)`` where ``(r, c)`` or ``(p, r, c)`` are coordinates of the blob and ``sigma`` is the standard deviation of the Gaussian kernel which detected the blob. References ---------- .. [1] http://en.wikipedia.org/wiki/Blob_detection#The_difference_of_Gaussians_approach Examples -------- >>> from skimage import data, feature >>> feature.blob_dog(data.coins(), threshold=.5, max_sigma=40) array([[ 267. , 359. , 16.777216], [ 267. , 115. , 10.48576 ], [ 263. , 302. , 16.777216], [ 263. , 245. , 16.777216], [ 261. , 173. , 16.777216], [ 260. , 46. , 16.777216], [ 198. , 155. , 10.48576 ], [ 196. , 43. , 10.48576 ], [ 195. , 102. , 16.777216], [ 194. , 277. , 16.777216], [ 193. , 213. , 16.777216], [ 185. , 347. , 16.777216], [ 128. , 154. , 10.48576 ], [ 127. , 102. , 10.48576 ], [ 125. , 208. , 10.48576 ], [ 125. , 45. , 16.777216], [ 124. , 337. , 10.48576 ], [ 120. , 272. , 16.777216], [ 58. , 100. , 10.48576 ], [ 54. , 276. , 10.48576 ], [ 54. , 42. , 16.777216], [ 52. , 216. , 16.777216], [ 52. , 155. , 16.777216], [ 45. , 336. , 16.777216]]) Notes ----- The radius of each blob is approximately :math:`\sqrt{2}\sigma` for a 2-D image and :math:`\sqrt{3}\sigma` for a 3-D image. """ image = img_as_float(image) # k such that min_sigma*(sigma_ratio**k) > max_sigma k = int(log(float(max_sigma) / min_sigma, sigma_ratio)) + 1 # a geometric progression of standard deviations for gaussian kernels sigma_list = np.array([min_sigma * (sigma_ratio ** i) for i in range(k + 1)]) gaussian_images = [gaussian_filter(image, s) for s in sigma_list] # computing difference between two successive Gaussian blurred images # multiplying with standard deviation provides scale invariance dog_images = [(gaussian_images[i] - gaussian_images[i + 1]) * sigma_list[i] for i in range(k)] image_cube = np.stack(dog_images, axis=-1) # local_maxima = get_local_maxima(image_cube, threshold) local_maxima = peak_local_max(image_cube, threshold_abs=threshold, footprint=np.ones((3,) * (image.ndim + 1)), threshold_rel=0.0, exclude_border=False) # Catch no peaks if local_maxima.size == 0: return np.empty((0, 3)) # Convert local_maxima to float64 lm = local_maxima.astype(np.float64) # Convert the last index to its corresponding scale value lm[:, -1] = sigma_list[local_maxima[:, -1]] return _prune_blobs(lm, overlap) def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2, overlap=.5, log_scale=False): """Finds blobs in the given grayscale image. Blobs are found using the Laplacian of Gaussian (LoG) method [1]_. For each blob found, the method returns its coordinates and the standard deviation of the Gaussian kernel that detected the blob. Parameters ---------- image : 2D or 3D ndarray Input grayscale image, blobs are assumed to be light on dark background (white on black). min_sigma : float, optional The minimum standard deviation for Gaussian Kernel. Keep this low to detect smaller blobs. max_sigma : float, optional The maximum standard deviation for Gaussian Kernel. Keep this high to detect larger blobs. num_sigma : int, optional The number of intermediate values of standard deviations to consider between `min_sigma` and `max_sigma`. threshold : float, optional. The absolute lower bound for scale space maxima. Local maxima smaller than thresh are ignored. Reduce this to detect blobs with less intensities. overlap : float, optional A value between 0 and 1. If the area of two blobs overlaps by a fraction greater than `threshold`, the smaller blob is eliminated. log_scale : bool, optional If set intermediate values of standard deviations are interpolated using a logarithmic scale to the base `10`. If not, linear interpolation is used. Returns ------- A : (n, image.ndim + 1) ndarray A 2d array with each row representing 3 values for a 2D image, and 4 values for a 3D image: ``(r, c, sigma)`` or ``(p, r, c, sigma)`` where ``(r, c)`` or ``(p, r, c)`` are coordinates of the blob and ``sigma`` is the standard deviation of the Gaussian kernel which detected the blob. References ---------- .. [1] http://en.wikipedia.org/wiki/Blob_detection#The_Laplacian_of_Gaussian Examples -------- >>> from skimage import data, feature, exposure >>> img = data.coins() >>> img = exposure.equalize_hist(img) # improves detection >>> feature.blob_log(img, threshold = .3) array([[ 266. , 115. , 11.88888889], [ 263. , 302. , 17.33333333], [ 263. , 244. , 17.33333333], [ 260. , 174. , 17.33333333], [ 198. , 155. , 11.88888889], [ 198. , 103. , 11.88888889], [ 197. , 44. , 11.88888889], [ 194. , 276. , 17.33333333], [ 194. , 213. , 17.33333333], [ 185. , 344. , 17.33333333], [ 128. , 154. , 11.88888889], [ 127. , 102. , 11.88888889], [ 126. , 208. , 11.88888889], [ 126. , 46. , 11.88888889], [ 124. , 336. , 11.88888889], [ 121. , 272. , 17.33333333], [ 113. , 323. , 1. ]]) Notes ----- The radius of each blob is approximately :math:`\sqrt{2}\sigma` for a 2-D image and :math:`\sqrt{3}\sigma` for a 3-D image. """ image = img_as_float(image) if log_scale: start, stop = log(min_sigma, 10), log(max_sigma, 10) sigma_list = np.logspace(start, stop, num_sigma) else: sigma_list = np.linspace(min_sigma, max_sigma, num_sigma) # computing gaussian laplace # s**2 provides scale invariance gl_images = [-gaussian_laplace(image, s) * s ** 2 for s in sigma_list] image_cube = np.stack(gl_images, axis=-1) local_maxima = peak_local_max(image_cube, threshold_abs=threshold, footprint=np.ones((3,) * (image.ndim + 1)), threshold_rel=0.0, exclude_border=False) # Catch no peaks if local_maxima.size == 0: return np.empty((0, 3)) # Convert local_maxima to float64 lm = local_maxima.astype(np.float64) # Convert the last index to its corresponding scale value lm[:, -1] = sigma_list[local_maxima[:, -1]] return _prune_blobs(lm, overlap) def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=0.01, overlap=.5, log_scale=False): """Finds blobs in the given grayscale image. Blobs are found using the Determinant of Hessian method [1]_. For each blob found, the method returns its coordinates and the standard deviation of the Gaussian Kernel used for the Hessian matrix whose determinant detected the blob. Determinant of Hessians is approximated using [2]_. Parameters ---------- image : 2D ndarray Input grayscale image.Blobs can either be light on dark or vice versa. min_sigma : float, optional The minimum standard deviation for Gaussian Kernel used to compute Hessian matrix. Keep this low to detect smaller blobs. max_sigma : float, optional The maximum standard deviation for Gaussian Kernel used to compute Hessian matrix. Keep this high to detect larger blobs. num_sigma : int, optional The number of intermediate values of standard deviations to consider between `min_sigma` and `max_sigma`. threshold : float, optional. The absolute lower bound for scale space maxima. Local maxima smaller than thresh are ignored. Reduce this to detect less prominent blobs. overlap : float, optional A value between 0 and 1. If the area of two blobs overlaps by a fraction greater than `threshold`, the smaller blob is eliminated. log_scale : bool, optional If set intermediate values of standard deviations are interpolated using a logarithmic scale to the base `10`. If not, linear interpolation is used. Returns ------- A : (n, 3) ndarray A 2d array with each row representing 3 values, ``(y,x,sigma)`` where ``(y,x)`` are coordinates of the blob and ``sigma`` is the standard deviation of the Gaussian kernel of the Hessian Matrix whose determinant detected the blob. References ---------- .. [1] http://en.wikipedia.org/wiki/Blob_detection#The_determinant_of_the_Hessian .. [2] Herbert Bay, Andreas Ess, Tinne Tuytelaars, Luc Van Gool, "SURF: Speeded Up Robust Features" ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf Examples -------- >>> from skimage import data, feature >>> img = data.coins() >>> feature.blob_doh(img) array([[ 270. , 363. , 30. ], [ 265. , 113. , 23.55555556], [ 262. , 243. , 23.55555556], [ 260. , 173. , 30. ], [ 197. , 153. , 20.33333333], [ 197. , 44. , 20.33333333], [ 195. , 100. , 23.55555556], [ 193. , 275. , 23.55555556], [ 192. , 212. , 23.55555556], [ 185. , 348. , 30. ], [ 156. , 302. , 30. ], [ 126. , 153. , 20.33333333], [ 126. , 101. , 20.33333333], [ 124. , 336. , 20.33333333], [ 123. , 205. , 20.33333333], [ 123. , 44. , 23.55555556], [ 121. , 271. , 30. ]]) Notes ----- The radius of each blob is approximately `sigma`. Computation of Determinant of Hessians is independent of the standard deviation. Therefore detecting larger blobs won't take more time. In methods line :py:meth:`blob_dog` and :py:meth:`blob_log` the computation of Gaussians for larger `sigma` takes more time. The downside is that this method can't be used for detecting blobs of radius less than `3px` due to the box filters used in the approximation of Hessian Determinant. """ assert_nD(image, 2) image = img_as_float(image) image = integral_image(image) if log_scale: start, stop = log(min_sigma, 10), log(max_sigma, 10) sigma_list = np.logspace(start, stop, num_sigma) else: sigma_list = np.linspace(min_sigma, max_sigma, num_sigma) hessian_images = [_hessian_matrix_det(image, s) for s in sigma_list] image_cube = np.dstack(hessian_images) local_maxima = peak_local_max(image_cube, threshold_abs=threshold, footprint=np.ones((3,) * image_cube.ndim), threshold_rel=0.0, exclude_border=False) # Catch no peaks if local_maxima.size == 0: return np.empty((0, 3)) # Convert local_maxima to float64 lm = local_maxima.astype(np.float64) # Convert the last index to its corresponding scale value lm[:, -1] = sigma_list[local_maxima[:, -1]] return _prune_blobs(lm, overlap)