From 5f070397fd5627312caef29cf97e2901fc384c90 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Tue, 4 Mar 2014 21:59:54 +0530
Subject: [PATCH 01/13] added blob_dog

---
 skimage/feature/__init__.py |   4 +-
 skimage/feature/blob.py     | 233 ++++++++++++++++++++++++++++++++++++
 2 files changed, 236 insertions(+), 1 deletion(-)
 create mode 100644 skimage/feature/blob.py

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..6ed2f385
--- /dev/null
+++ b/skimage/feature/blob.py
@@ -0,0 +1,233 @@
+import numpy as np
+from scipy.ndimage.filters import gaussian_filter, maximum_filter
+import itertools as itt
+import math
+from math import sqrt, hypot, log
+from numpy import arccos
+from skimage.util import img_as_float
+
+
+# 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)
+
+# A lot of this code is borrowed from here
+# https://github.com/adonath/blob_detection/tree/master/blob_detection
+
+
+def _get_local_maxima_3d(array, threshold):
+    """Finds local maxima in a 3d array.
+
+    A pixel is considered to be a maximum if it is greater than or equal to all
+    its 28 neighbors in the 3d cube.
+
+    Parameters
+    ----------
+    array : ndarray
+        The 3d array whose local maximas are sought.
+    thresh : float
+        Local maximas lesser than thresh are ignored.
+
+    Returns
+    -------
+    A : (n, 3) ndarray
+        A 2d array in which each row contains 3 values, the indices of local
+        maxima.
+
+    """
+    # computing max filter using all neighbors in cube
+    fp = np.ones((3, 3, 3))
+    max_array = maximum_filter(array, footprint=fp)
+    peaks = (max_array == array) & (array > threshold)
+    return np.argwhere(peaks)
+
+
+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=25, 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
+        `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 `thresh`, the smaller blob is eliminated.
+    log_scale : boolean, optional
+        If set to True, the standard deviations of Gaussian Kernels are
+        interpolated using a logarithmic scale. This is useful when finding
+        blobs with a large variation in size. If set, scales are interpolated
+        with log to the base 10.
+
+    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())
+    array([[  46,  336, 2513],
+           [  53,  156, 2035],
+           [  53,  217, 1608],
+           [  54,  276, 1231],
+           [  55,   42, 1608],
+           [  57,  100, 1231],
+           [ 121,  272, 2035],
+           [ 124,  337, 1413],
+           [ 125,   45, 1815],
+           [ 125,  207, 1608],
+           [ 126,  102, 1231],
+           [ 128,  154, 1231],
+           [ 185,  347, 2513],
+           [ 194,  213, 1815],
+           [ 194,  277, 1608],
+           [ 196,   42, 1231],
+           [ 196,  101, 1608],
+           [ 197,  155, 1231],
+           [ 260,   46, 2513],
+           [ 261,  174, 2035],
+           [ 263,  245, 2035],
+           [ 263,  302, 2035],
+           [ 266,  114, 1608],
+           [ 268,  358, 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 succesive gaussian blurred images
+    # multipying with square of standard deviation provides scale invariance
+    dog_images = [(gaussian_images[i] - gaussian_images[i + 1])
+                  * sigma_list[i] ** 2 for i in range(k)]
+    image_cube = np.dstack(dog_images)
+
+    local_maxima = _get_local_maxima_3d(image_cube, threshold)
+
+    # 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 []

From 6b0b9db1d0729ab882a02c0a3f0bb9296db5ab73 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Tue, 4 Mar 2014 22:09:14 +0530
Subject: [PATCH 02/13] doctest

---
 skimage/feature/blob.py | 53 +++++++++++++++++++++--------------------
 1 file changed, 27 insertions(+), 26 deletions(-)

diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 6ed2f385..86d77f8f 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -170,32 +170,33 @@ def blob_dog(image, min_sigma=1, max_sigma=25, sigma_ratio=1.6, threshold=2.0,
 
     Examples
     --------
-    >>> from skimage import data,feature
+    >>> from skimage import data, feature
     >>> feature.blob_dog(data.coins())
-    array([[  46,  336, 2513],
-           [  53,  156, 2035],
-           [  53,  217, 1608],
-           [  54,  276, 1231],
-           [  55,   42, 1608],
-           [  57,  100, 1231],
-           [ 121,  272, 2035],
-           [ 124,  337, 1413],
-           [ 125,   45, 1815],
-           [ 125,  207, 1608],
-           [ 126,  102, 1231],
-           [ 128,  154, 1231],
-           [ 185,  347, 2513],
-           [ 194,  213, 1815],
+    array([[  45,  336, 1608],
+           [  51,  277, 1608],
+           [  52,  155, 1608],
+           [  52,  216, 1608],
+           [  54,   42, 1608],
+           [  56,  101, 1608],
+           [ 120,  272, 1608],
+           [ 124,  206, 1608],
+           [ 124,  339, 1608],
+           [ 125,   45, 1608],
+           [ 125,  102, 1608],
+           [ 127,  154, 1608],
+           [ 185,  347, 1608],
+           [ 193,  213, 1608],
            [ 194,  277, 1608],
-           [ 196,   42, 1231],
-           [ 196,  101, 1608],
-           [ 197,  155, 1231],
-           [ 260,   46, 2513],
-           [ 261,  174, 2035],
-           [ 263,  245, 2035],
-           [ 263,  302, 2035],
-           [ 266,  114, 1608],
-           [ 268,  358, 1608]])
+           [ 195,  102, 1608],
+           [ 196,   41, 1608],
+           [ 197,  154, 1608],
+           [ 260,   46, 1608],
+           [ 261,  173, 1608],
+           [ 263,  245, 1608],
+           [ 263,  302, 1608],
+           [ 267,  115, 1608],
+           [ 267,  359, 1608]])
+
 
     """
 
@@ -213,8 +214,8 @@ def blob_dog(image, min_sigma=1, max_sigma=25, sigma_ratio=1.6, threshold=2.0,
 
     gaussian_images = [gaussian_filter(image, s) for s in sigma_list]
 
-    # computing difference between two succesive gaussian blurred images
-    # multipying with square of standard deviation provides scale invariance
+    # computing difference between two successive Gaussian blurred images
+    # multiplying with square of standard deviation provides scale invariance
     dog_images = [(gaussian_images[i] - gaussian_images[i + 1])
                   * sigma_list[i] ** 2 for i in range(k)]
     image_cube = np.dstack(dog_images)

From d9adb82b145d88962493d618ade4e68d7b10f861 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Tue, 4 Mar 2014 22:30:42 +0530
Subject: [PATCH 03/13] added test_blob

---
 skimage/feature/tests/test_blob.py | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)
 create mode 100644 skimage/feature/tests/test_blob.py

diff --git a/skimage/feature/tests/test_blob.py b/skimage/feature/tests/test_blob.py
new file mode 100644
index 00000000..e1d2baf9
--- /dev/null
+++ b/skimage/feature/tests/test_blob.py
@@ -0,0 +1,28 @@
+import numpy as np
+from skimage.draw import circle
+from skimage.feature import blob_dog
+
+
+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, 30)
+    img[xs, ys] = 255
+
+    blobs = blob_dog(img)
+    coords = blobs[:, 0:2]
+
+    if coords[coords == [400, 130]].shape != (2,):
+        assert False
+
+    if coords[coords == [100, 300]].shape != (2,):
+        assert False
+
+    if coords[coords == [200, 350]].shape != (2,):
+        assert False

From c5a6f67654f1a68c29ec1eddbb1d8560c2670a78 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Tue, 4 Mar 2014 22:41:04 +0530
Subject: [PATCH 04/13] doc string

---
 skimage/feature/blob.py | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 86d77f8f..c702ecad 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -26,7 +26,7 @@ def _get_local_maxima_3d(array, threshold):
     array : ndarray
         The 3d array whose local maximas are sought.
     thresh : float
-        Local maximas lesser than thresh are ignored.
+        Local maximas lesser than `thresh` are ignored.
 
     Returns
     -------
@@ -144,7 +144,6 @@ def blob_dog(image, min_sigma=1, max_sigma=25, sigma_ratio=1.6, threshold=2.0,
     sigma_ratio : float, optional
         The ratio between the standard deviation of Gaussian Kernels used for
         computing the Difference of Gaussians
-        `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
@@ -152,11 +151,6 @@ def blob_dog(image, min_sigma=1, max_sigma=25, sigma_ratio=1.6, threshold=2.0,
     overlap : float, optional
         A value between 0 and 1. If the area of two blobs overlaps by a
         fraction greater than `thresh`, the smaller blob is eliminated.
-    log_scale : boolean, optional
-        If set to True, the standard deviations of Gaussian Kernels are
-        interpolated using a logarithmic scale. This is useful when finding
-        blobs with a large variation in size. If set, scales are interpolated
-        with log to the base 10.
 
     Returns
     -------

From 66bd0c4818a17ecfd1644c7cd05281f1110167ce Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Tue, 4 Mar 2014 22:44:33 +0530
Subject: [PATCH 05/13] typo

---
 skimage/feature/blob.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index c702ecad..f00b1785 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -191,7 +191,6 @@ def blob_dog(image, min_sigma=1, max_sigma=25, sigma_ratio=1.6, threshold=2.0,
            [ 267,  115, 1608],
            [ 267,  359, 1608]])
 
-
     """
 
     if image.ndim != 2:

From 59fa47cb5ef852c3c719f0aed1c5cacb42d79d7a Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Tue, 4 Mar 2014 22:56:11 +0530
Subject: [PATCH 06/13] threshold

---
 skimage/feature/blob.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index f00b1785..11d4c717 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -150,7 +150,7 @@ def blob_dog(image, min_sigma=1, max_sigma=25, sigma_ratio=1.6, threshold=2.0,
         intensities.
     overlap : float, optional
         A value between 0 and 1. If the area of two blobs overlaps by a
-        fraction greater than `thresh`, the smaller blob is eliminated.
+        fraction greater than `threshold`, the smaller blob is eliminated.
 
     Returns
     -------

From 5a988fdae79c66c0b8ffec6a1328be42bd8b0c32 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Wed, 5 Mar 2014 00:46:58 +0530
Subject: [PATCH 07/13] Changed formula, this improves detected size accuracy

---
 skimage/feature/blob.py            | 26 +++++++++++++-------------
 skimage/feature/tests/test_blob.py | 28 +++++++++++++++++++---------
 2 files changed, 32 insertions(+), 22 deletions(-)

diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 11d4c717..4ba18b33 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -123,7 +123,7 @@ def _prune_blobs(blobs_array, overlap):
     return np.array([b for b in blobs_array if b[2] > 0])
 
 
-def blob_dog(image, min_sigma=1, max_sigma=25, sigma_ratio=1.6, threshold=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.
 
@@ -165,30 +165,30 @@ def blob_dog(image, min_sigma=1, max_sigma=25, sigma_ratio=1.6, threshold=2.0,
     Examples
     --------
     >>> from skimage import data, feature
-    >>> feature.blob_dog(data.coins())
+    >>> feature.blob_dog(data.coins(),threshold=.5,max_sigma=40)
     array([[  45,  336, 1608],
-           [  51,  277, 1608],
            [  52,  155, 1608],
            [  52,  216, 1608],
            [  54,   42, 1608],
-           [  56,  101, 1608],
+           [  54,  276,  628],
+           [  58,  100,  628],
            [ 120,  272, 1608],
-           [ 124,  206, 1608],
-           [ 124,  339, 1608],
+           [ 124,  337,  628],
            [ 125,   45, 1608],
-           [ 125,  102, 1608],
-           [ 127,  154, 1608],
+           [ 125,  208,  628],
+           [ 127,  102,  628],
+           [ 128,  154,  628],
            [ 185,  347, 1608],
            [ 193,  213, 1608],
            [ 194,  277, 1608],
            [ 195,  102, 1608],
-           [ 196,   41, 1608],
-           [ 197,  154, 1608],
+           [ 196,   43,  628],
+           [ 198,  155,  628],
            [ 260,   46, 1608],
            [ 261,  173, 1608],
            [ 263,  245, 1608],
            [ 263,  302, 1608],
-           [ 267,  115, 1608],
+           [ 267,  115,  628],
            [ 267,  359, 1608]])
 
     """
@@ -208,9 +208,9 @@ def blob_dog(image, min_sigma=1, max_sigma=25, sigma_ratio=1.6, threshold=2.0,
     gaussian_images = [gaussian_filter(image, s) for s in sigma_list]
 
     # computing difference between two successive Gaussian blurred images
-    # multiplying with square of standard deviation provides scale invariance
+    # multiplying with standard deviation provides scale invariance
     dog_images = [(gaussian_images[i] - gaussian_images[i + 1])
-                  * sigma_list[i] ** 2 for i in range(k)]
+                  * sigma_list[i] for i in range(k)]
     image_cube = np.dstack(dog_images)
 
     local_maxima = _get_local_maxima_3d(image_cube, threshold)
diff --git a/skimage/feature/tests/test_blob.py b/skimage/feature/tests/test_blob.py
index e1d2baf9..54ed38ce 100644
--- a/skimage/feature/tests/test_blob.py
+++ b/skimage/feature/tests/test_blob.py
@@ -1,6 +1,7 @@
 import numpy as np
 from skimage.draw import circle
 from skimage.feature import blob_dog
+import math
 
 
 def test_blob_dog():
@@ -12,17 +13,26 @@ def test_blob_dog():
     xs, ys = circle(100, 300, 25)
     img[xs, ys] = 255
 
-    xs, ys = circle(200, 350, 30)
+    xs, ys = circle(200, 350, 45)
     img[xs, ys] = 255
 
-    blobs = blob_dog(img)
-    coords = blobs[:, 0:2]
+    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
 
-    if coords[coords == [400, 130]].shape != (2,):
-        assert False
+    b = s[0]
+    assert abs(b[0] - 400) <= thresh
+    assert abs(b[1] - 130) <= thresh
+    assert abs(radius(b[2]) - 5) <= thresh
 
-    if coords[coords == [100, 300]].shape != (2,):
-        assert False
+    b = s[1]
+    assert abs(b[0] - 100) <= thresh
+    assert abs(b[1] - 300) <= thresh
+    assert abs(radius(b[2]) - 25) <= thresh
 
-    if coords[coords == [200, 350]].shape != (2,):
-        assert False
+    b = s[2]
+    assert abs(b[0] - 200) <= thresh
+    assert abs(b[1] - 350) <= thresh
+    assert abs(radius(b[2]) - 45) <= thresh

From e449ce98701591fb19f8012bbfd35b6aa483288d Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Wed, 5 Mar 2014 15:37:10 +0530
Subject: [PATCH 08/13] Made get_local_maxima_3d public It also takes
 connectivity as a paramater now

---
 skimage/feature/blob.py | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 4ba18b33..c5cca3bc 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -1,5 +1,6 @@
 import numpy as np
 from scipy.ndimage.filters import gaussian_filter, maximum_filter
+from scipy.ndimage.morphology import generate_binary_structure
 import itertools as itt
 import math
 from math import sqrt, hypot, log
@@ -15,11 +16,11 @@ from skimage.util import img_as_float
 # https://github.com/adonath/blob_detection/tree/master/blob_detection
 
 
-def _get_local_maxima_3d(array, threshold):
+def get_local_maxima_3d(array, threshold, connectivity=3):
     """Finds local maxima in a 3d array.
 
     A pixel is considered to be a maximum if it is greater than or equal to all
-    its 28 neighbors in the 3d cube.
+    its neighbors in the 3d cube.
 
     Parameters
     ----------
@@ -27,6 +28,11 @@ def _get_local_maxima_3d(array, threshold):
         The 3d array whose local maximas are sought.
     thresh : float
         Local maximas lesser than `thresh` are ignored.
+    connectivity : float, optional
+        Elements up to a squared distance of `connectivity` from a point are
+        considered neighbors. If `connectivity` is 1, 6 neighbors are
+        considered, if `connectivity` is 2, 18 neighbors are considered and if
+        `connectivity` is 3, all 26 neighbors are considered.
 
     Returns
     -------
@@ -36,7 +42,7 @@ def _get_local_maxima_3d(array, threshold):
 
     """
     # computing max filter using all neighbors in cube
-    fp = np.ones((3, 3, 3))
+    fp = generate_binary_structure(3, connectivity)
     max_array = maximum_filter(array, footprint=fp)
     peaks = (max_array == array) & (array > threshold)
     return np.argwhere(peaks)
@@ -213,7 +219,7 @@ def blob_dog(image, min_sigma=1, max_sigma=50, sigma_ratio=1.6, threshold=2.0,
                   * sigma_list[i] for i in range(k)]
     image_cube = np.dstack(dog_images)
 
-    local_maxima = _get_local_maxima_3d(image_cube, threshold)
+    local_maxima = get_local_maxima_3d(image_cube, threshold)
 
     # Convert the last index to its corresponding scale value
     local_maxima[:, 2] = sigma_list[local_maxima[:, 2]]

From 7681598380b7426cabe4fd7404b0f3a54627169a Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Wed, 5 Mar 2014 21:08:37 +0530
Subject: [PATCH 09/13] changed params and functionality of get_local_maxima

---
 skimage/feature/blob.py | 27 ++++++++++++++-------------
 1 file changed, 14 insertions(+), 13 deletions(-)

diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index c5cca3bc..d1ab9ecc 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -16,23 +16,24 @@ from skimage.util import img_as_float
 # https://github.com/adonath/blob_detection/tree/master/blob_detection
 
 
-def get_local_maxima_3d(array, threshold, connectivity=3):
-    """Finds local maxima in a 3d array.
+def get_local_maxima(ar, threshold, connectivity=3):
+    """Finds local maxima in an array.
 
-    A pixel is considered to be a maximum if it is greater than or equal to all
-    its neighbors in the 3d cube.
+    A point is considered to be a maximum if it is greater than or equal to all
+    its neighbors.
 
     Parameters
     ----------
-    array : ndarray
-        The 3d array whose local maximas are sought.
+    ar : ndarray
+        The array whose local maximas are sought.
     thresh : float
         Local maximas lesser than `thresh` are ignored.
     connectivity : float, optional
         Elements up to a squared distance of `connectivity` from a point are
-        considered neighbors. If `connectivity` is 1, 6 neighbors are
-        considered, if `connectivity` is 2, 18 neighbors are considered and if
-        `connectivity` is 3, all 26 neighbors are considered.
+        considered neighbors. For example in a 3 Dimensional array, if
+        `connectivity` is 1, 6 neighbors are considered, if `connectivity` is
+        2, 18 neighbors are considered and if `connectivity` is 3, all 26
+        neighbors are considered.
 
     Returns
     -------
@@ -42,9 +43,9 @@ def get_local_maxima_3d(array, threshold, connectivity=3):
 
     """
     # computing max filter using all neighbors in cube
-    fp = generate_binary_structure(3, connectivity)
-    max_array = maximum_filter(array, footprint=fp)
-    peaks = (max_array == array) & (array > threshold)
+    fp = generate_binary_structure(ar.ndim, connectivity)
+    max_ar = maximum_filter(ar, footprint=fp)
+    peaks = (max_ar == ar) & (ar > threshold)
     return np.argwhere(peaks)
 
 
@@ -219,7 +220,7 @@ def blob_dog(image, min_sigma=1, max_sigma=50, sigma_ratio=1.6, threshold=2.0,
                   * sigma_list[i] for i in range(k)]
     image_cube = np.dstack(dog_images)
 
-    local_maxima = get_local_maxima_3d(image_cube, threshold)
+    local_maxima = get_local_maxima(image_cube, threshold)
 
     # Convert the last index to its corresponding scale value
     local_maxima[:, 2] = sigma_list[local_maxima[:, 2]]

From 579c3d5aa3e2c0438772ec11c669772d5f9c3467 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Thu, 6 Mar 2014 13:04:09 +0530
Subject: [PATCH 10/13] added get_local_maxima to feature

---
 skimage/feature/__init__.py |  5 +++--
 skimage/feature/blob.py     | 10 ++++++++++
 2 files changed, 13 insertions(+), 2 deletions(-)

diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py
index 3c33ffb2..ef4aa113 100644
--- a/skimage/feature/__init__.py
+++ b/skimage/feature/__init__.py
@@ -14,7 +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
+from .blob import blob_dog, get_local_maxima
 
 
 __all__ = ['daisy',
@@ -42,4 +42,5 @@ __all__ = ['daisy',
            'ORB',
            'match_descriptors',
            'plot_matches',
-           'blob_dog']
+           'blob_dog',
+           'get_local_maxima']
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index d1ab9ecc..23ce7452 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -41,6 +41,16 @@ def get_local_maxima(ar, threshold, connectivity=3):
         A 2d array in which each row contains 3 values, the indices of local
         maxima.
 
+    Examples
+    --------
+    >>> a = np.array([[ 0 , 0 , 0 , 0 , 0],
+    ...               [ 0 , 0 , 3 , 0 , 0],
+    ...               [ 0 , 0 , 1 , 0 , 0],
+    ...               [ 0 , 1 , 0 , 0 , 0],
+    ...               [ 0 , 0 , 0 , 0 , 0]])
+    >>> get_local_maxima(a, threshold = 1, connectivity = 2)
+    array([[1, 2]])
+
     """
     # computing max filter using all neighbors in cube
     fp = generate_binary_structure(ar.ndim, connectivity)

From d6910ba9bffde657265d29cb0ea95066f675397f Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sun, 9 Mar 2014 21:56:45 +0530
Subject: [PATCH 11/13] removed get_local_maxima, now using peak_local_max

---
 skimage/feature/__init__.py |  2 +-
 skimage/feature/blob.py     | 53 +++++--------------------------------
 2 files changed, 8 insertions(+), 47 deletions(-)

diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py
index ef4aa113..37f67794 100644
--- a/skimage/feature/__init__.py
+++ b/skimage/feature/__init__.py
@@ -14,7 +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, get_local_maxima
+from .blob import blob_dog
 
 
 __all__ = ['daisy',
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 23ce7452..5389954c 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -1,11 +1,11 @@
 import numpy as np
-from scipy.ndimage.filters import gaussian_filter, maximum_filter
-from scipy.ndimage.morphology import generate_binary_structure
+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:
@@ -16,49 +16,6 @@ from skimage.util import img_as_float
 # https://github.com/adonath/blob_detection/tree/master/blob_detection
 
 
-def get_local_maxima(ar, threshold, connectivity=3):
-    """Finds local maxima in an array.
-
-    A point is considered to be a maximum if it is greater than or equal to all
-    its neighbors.
-
-    Parameters
-    ----------
-    ar : ndarray
-        The array whose local maximas are sought.
-    thresh : float
-        Local maximas lesser than `thresh` are ignored.
-    connectivity : float, optional
-        Elements up to a squared distance of `connectivity` from a point are
-        considered neighbors. For example in a 3 Dimensional array, if
-        `connectivity` is 1, 6 neighbors are considered, if `connectivity` is
-        2, 18 neighbors are considered and if `connectivity` is 3, all 26
-        neighbors are considered.
-
-    Returns
-    -------
-    A : (n, 3) ndarray
-        A 2d array in which each row contains 3 values, the indices of local
-        maxima.
-
-    Examples
-    --------
-    >>> a = np.array([[ 0 , 0 , 0 , 0 , 0],
-    ...               [ 0 , 0 , 3 , 0 , 0],
-    ...               [ 0 , 0 , 1 , 0 , 0],
-    ...               [ 0 , 1 , 0 , 0 , 0],
-    ...               [ 0 , 0 , 0 , 0 , 0]])
-    >>> get_local_maxima(a, threshold = 1, connectivity = 2)
-    array([[1, 2]])
-
-    """
-    # computing max filter using all neighbors in cube
-    fp = generate_binary_structure(ar.ndim, connectivity)
-    max_ar = maximum_filter(ar, footprint=fp)
-    peaks = (max_ar == ar) & (ar > threshold)
-    return np.argwhere(peaks)
-
-
 def _blob_overlap(blob1, blob2):
     """Finds the overlapping area fraction between two blobs.
 
@@ -230,7 +187,11 @@ def blob_dog(image, min_sigma=1, max_sigma=50, sigma_ratio=1.6, threshold=2.0,
                   * sigma_list[i] for i in range(k)]
     image_cube = np.dstack(dog_images)
 
-    local_maxima = get_local_maxima(image_cube, threshold)
+    # 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]]

From 7986714c5d302b1937c16b6ec19a2529f30ad0b0 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sun, 9 Mar 2014 22:00:26 +0530
Subject: [PATCH 12/13] rectified __all__

---
 skimage/feature/__init__.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py
index 37f67794..3c33ffb2 100644
--- a/skimage/feature/__init__.py
+++ b/skimage/feature/__init__.py
@@ -42,5 +42,4 @@ __all__ = ['daisy',
            'ORB',
            'match_descriptors',
            'plot_matches',
-           'blob_dog',
-           'get_local_maxima']
+           'blob_dog']

From 0f44e6f1ba35d7a4ecf505597c1d784bd74954f2 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Mon, 10 Mar 2014 15:39:15 +0530
Subject: [PATCH 13/13] added names in CONTRIBUTIRS.txt

---
 CONTRIBUTORS.txt        | 6 ++++++
 skimage/feature/blob.py | 3 ---
 2 files changed, 6 insertions(+), 3 deletions(-)

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/blob.py b/skimage/feature/blob.py
index 5389954c..2164755b 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -12,9 +12,6 @@ from .peak import peak_local_max
 # http://www.cs.utah.edu/~jfishbau/advimproc/project1/ (04.04.2013)
 # Theory behind: http://en.wikipedia.org/wiki/Blob_detection (04.04.2013)
 
-# A lot of this code is borrowed from here
-# https://github.com/adonath/blob_detection/tree/master/blob_detection
-
 
 def _blob_overlap(blob1, blob2):
     """Finds the overlapping area fraction between two blobs.