diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 00a76e59..140842ab 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -176,3 +176,9 @@ - François Orieux Image deconvolution http://research.orieux.fr + +- Vighnesh Birodkar + Blob Detection + +- Axel Donath + Blob Detection diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 1b1cb189..3c33ffb2 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -14,6 +14,7 @@ from .censure import CENSURE from .orb import ORB from .match import match_descriptors from .util import plot_matches +from .blob import blob_dog __all__ = ['daisy', @@ -40,4 +41,5 @@ __all__ = ['daisy', 'CENSURE', 'ORB', 'match_descriptors', - 'plot_matches'] + 'plot_matches', + 'blob_dog'] diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py new file mode 100644 index 00000000..2164755b --- /dev/null +++ b/skimage/feature/blob.py @@ -0,0 +1,202 @@ +import numpy as np +from scipy.ndimage.filters import gaussian_filter +import itertools as itt +import math +from math import sqrt, hypot, log +from numpy import arccos +from skimage.util import img_as_float +from .peak import peak_local_max + + +# 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 _blob_overlap(blob1, blob2): + """Finds the overlapping area fraction between two blobs. + + Returns a float representing fraction of overlapped area. + + Parameters + ---------- + blob1 : sequence + A sequence of ``(y,x,sigma)``, where ``x,y`` are coordinates of blob + and sigma is the standard deviation of the Gaussian kernel which + detected the blob. + blob2 : sequence + A sequence of ``(y,x,sigma)``, where ``x,y`` 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. + + """ + root2 = sqrt(2) + + # extent of the blob is given by sqrt(2)*scale + r1 = blob1[2] * root2 + r2 = blob2[2] * root2 + + d = hypot(blob1[0] - blob2[0], blob1[1] - blob2[1]) + + if d > r1 + r2: + return 0 + + # one blob is inside the other, the smaller blob must die + if d <= abs(r1 - r2): + return 1 + + acos1 = arccos((d ** 2 + r1 ** 2 - r2 ** 2) / (2 * d * r1)) + acos2 = arccos((d ** 2 + r2 ** 2 - r1 ** 2) / (2 * d * r2)) + 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 _prune_blobs(blobs_array, overlap): + """Eliminated blobs with area overlap. + + Parameters + ---------- + blobs_array : ndarray + a 2d array with each row representing 3 values, the ``(y,x,sigma)`` + where ``(y,x)`` are coordinates of the blob and sigma is the standard + deviation of the Gaussian kernel which detected the blob. + 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. + + """ + + # iterating again might eliminate more blobs, but one iteration suffices + # for most cases + for blob1, blob2 in itt.combinations(blobs_array, 2): + if _blob_overlap(blob1, blob2) > overlap: + if blob1[2] > blob2[2]: + blob2[2] = -1 + else: + blob1[2] = -1 + + # return blobs_array[blobs_array[:, 2] > 0] + return np.array([b for b in blobs_array if b[2] > 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, its coordinates and area are returned. + + Parameters + ---------- + image : 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, 3) ndarray + A 2d array with each row containing the Y-Coordinate , the + X-Coordinate and the estimated area of the blob respectively. + + 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([[ 45, 336, 1608], + [ 52, 155, 1608], + [ 52, 216, 1608], + [ 54, 42, 1608], + [ 54, 276, 628], + [ 58, 100, 628], + [ 120, 272, 1608], + [ 124, 337, 628], + [ 125, 45, 1608], + [ 125, 208, 628], + [ 127, 102, 628], + [ 128, 154, 628], + [ 185, 347, 1608], + [ 193, 213, 1608], + [ 194, 277, 1608], + [ 195, 102, 1608], + [ 196, 43, 628], + [ 198, 155, 628], + [ 260, 46, 1608], + [ 261, 173, 1608], + [ 263, 245, 1608], + [ 263, 302, 1608], + [ 267, 115, 628], + [ 267, 359, 1608]]) + + """ + + if image.ndim != 2: + raise ValueError("'image' must be a grayscale ") + + 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.dstack(dog_images) + + # local_maxima = get_local_maxima(image_cube, threshold) + local_maxima = peak_local_max(image_cube, threshold_abs=threshold, + footprint=np.ones((3, 3, 3)), + threshold_rel=0.0, + exclude_border=False) + + # Convert the last index to its corresponding scale value + local_maxima[:, 2] = sigma_list[local_maxima[:, 2]] + ret_val = _prune_blobs(local_maxima, overlap) + + if len(ret_val) > 0: + ret_val[:, 2] = math.pi * \ + ((ret_val[:, 2] * math.sqrt(2)) ** 2).astype(int) + return ret_val + else: + return [] diff --git a/skimage/feature/tests/test_blob.py b/skimage/feature/tests/test_blob.py new file mode 100644 index 00000000..54ed38ce --- /dev/null +++ b/skimage/feature/tests/test_blob.py @@ -0,0 +1,38 @@ +import numpy as np +from skimage.draw import circle +from skimage.feature import blob_dog +import math + + +def test_blob_dog(): + img = np.ones((512, 512)) + + xs, ys = circle(400, 130, 5) + img[xs, ys] = 255 + + xs, ys = circle(100, 300, 25) + img[xs, ys] = 255 + + xs, ys = circle(200, 350, 45) + img[xs, ys] = 255 + + blobs = blob_dog(img, min_sigma=5, max_sigma=50) + area = lambda x: x[2] + radius = lambda x: math.sqrt(x / math.pi) + s = sorted(blobs, key=area) + thresh = 5 + + b = s[0] + assert abs(b[0] - 400) <= thresh + assert abs(b[1] - 130) <= thresh + assert abs(radius(b[2]) - 5) <= thresh + + b = s[1] + assert abs(b[0] - 100) <= thresh + assert abs(b[1] - 300) <= thresh + assert abs(radius(b[2]) - 25) <= thresh + + b = s[2] + assert abs(b[0] - 200) <= thresh + assert abs(b[1] - 350) <= thresh + assert abs(radius(b[2]) - 45) <= thresh