From f456245d57aa13caeba4356b0f195e1a6be864c4 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Fri, 14 Mar 2014 20:57:17 +0530
Subject: [PATCH 01/22] added cython file

---
 skimage/feature/_hessian_det_appx.pyx | 69 +++++++++++++++++++++++++++
 skimage/feature/setup.py              |  3 ++
 2 files changed, 72 insertions(+)
 create mode 100644 skimage/feature/_hessian_det_appx.pyx

diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
new file mode 100644
index 00000000..fa5f4cef
--- /dev/null
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -0,0 +1,69 @@
+import numpy as np
+cimport numpy as np
+
+from skimage.transform import integral_image, integrate
+from skimage import util
+
+
+cdef inline int  clip(np.int_t x, np.int_t low, np.int_t high):
+    if(x > high):
+        return high
+    if(x < low):
+        return low
+    return x
+
+
+cdef inline int integ(np.int_t[:, :] img, np.int_t r1, np.int_t c1, np.int_t rl, np.int_t cl):
+
+    r1 = clip(r1, 0, img.shape[0] - 1)
+    c1 = clip(c1, 0, img.shape[1] - 1)
+
+    r2 = clip(r1 + rl, 0, img.shape[0] - 1)
+    c2 = clip(c1 + cl, 0, img.shape[1] - 1)
+
+    cdef np.int_t r = img[r2, c2] + img[r1, c1] - img[r1, c2] - img[r2, c1]
+
+    if (r < 0):
+        return 0
+    return r
+
+
+def hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
+
+    cdef np.int_t[:, :] img = image
+    cdef int size = int(3 * sigma)
+    cdef np.ndarray[np.float_t, ndim = 2] out = np.zeros_like(img).astype(np.float)
+
+    cdef int height = img.shape[0]
+    cdef int width = img.shape[1]
+    cdef int r, c
+    cdef int s2 = (size - 1) / 2
+    cdef int s3 = size / 3
+    cdef int l = size / 3
+    cdef int w = size
+
+    cdef float dxx, dyy, dxy
+
+    if not size % 2:
+        size += 1
+
+    for r in range(height):
+        for c in range(width):
+        
+            dxy =  integ(img, r - s3, c + 1, s3, s3) + \
+                   integ(img, r + 1, c - s3, s3, s3) - \
+                   integ(img, r - s3, c - s3, s3, s3) - \
+                   integ(img, r + 1, c + 1, s3, s3)
+            dxy = -dxy / w / w
+
+            dxx = integ(img, r - s3 + 1, c - s2, 2 * s3 - 1,w) - \
+                  integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3) * 3
+            dxx = -dxx / w / w
+
+            dyy = integ(img, r - s2, c - s2 + 1, w, 2 * s3 - 1) - \
+                  integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1) * 3
+            dyy = -dyy / w / w
+
+            out[r, c] = (dxx * dyy - 0.81 * (dxy * dxy))
+
+    return out
diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py
index 9d8b4fbb..d7be3dcf 100644
--- a/skimage/feature/setup.py
+++ b/skimage/feature/setup.py
@@ -17,6 +17,7 @@ def configuration(parent_package='', top_path=None):
     cython(['orb_cy.pyx'], working_path=base_path)
     cython(['brief_cy.pyx'], working_path=base_path)
     cython(['_texture.pyx'], working_path=base_path)
+    cython(['_hessian_det_appx.pyx'], working_path=base_path)
 
     config.add_extension('corner_cy', sources=['corner_cy.c'],
                          include_dirs=[get_numpy_include_dirs()])
@@ -28,6 +29,8 @@ def configuration(parent_package='', top_path=None):
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_texture', sources=['_texture.c'],
                          include_dirs=[get_numpy_include_dirs(), '../_shared'])
+    config.add_extension('_hessian_det_appx', sources=['_hessian_det_appx.c'],
+                         include_dirs=[get_numpy_include_dirs()])
 
     return config
 

From 448c4114aefcd8865bac75585937ecad91b21a9c Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Fri, 14 Mar 2014 20:59:46 +0530
Subject: [PATCH 02/22] made methods private

---
 skimage/feature/_hessian_det_appx.pyx | 28 +++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index fa5f4cef..92c4c08e 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -5,7 +5,7 @@ from skimage.transform import integral_image, integrate
 from skimage import util
 
 
-cdef inline int  clip(np.int_t x, np.int_t low, np.int_t high):
+cdef inline int  _clip(np.int_t x, np.int_t low, np.int_t high):
     if(x > high):
         return high
     if(x < low):
@@ -13,13 +13,13 @@ cdef inline int  clip(np.int_t x, np.int_t low, np.int_t high):
     return x
 
 
-cdef inline int integ(np.int_t[:, :] img, np.int_t r1, np.int_t c1, np.int_t rl, np.int_t cl):
+cdef inline int _integ(np.int_t[:, :] img, np.int_t r1, np.int_t c1, np.int_t rl, np.int_t cl):
 
-    r1 = clip(r1, 0, img.shape[0] - 1)
-    c1 = clip(c1, 0, img.shape[1] - 1)
+    r1 = _clip(r1, 0, img.shape[0] - 1)
+    c1 = _clip(c1, 0, img.shape[1] - 1)
 
-    r2 = clip(r1 + rl, 0, img.shape[0] - 1)
-    c2 = clip(c1 + cl, 0, img.shape[1] - 1)
+    r2 = _clip(r1 + rl, 0, img.shape[0] - 1)
+    c2 = _clip(c1 + cl, 0, img.shape[1] - 1)
 
     cdef np.int_t r = img[r2, c2] + img[r1, c1] - img[r1, c2] - img[r2, c1]
 
@@ -50,18 +50,18 @@ def hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
     for r in range(height):
         for c in range(width):
         
-            dxy =  integ(img, r - s3, c + 1, s3, s3) + \
-                   integ(img, r + 1, c - s3, s3, s3) - \
-                   integ(img, r - s3, c - s3, s3, s3) - \
-                   integ(img, r + 1, c + 1, s3, s3)
+            dxy =  _integ(img, r - s3, c + 1, s3, s3) + \
+                   _integ(img, r + 1, c - s3, s3, s3) - \
+                   _integ(img, r - s3, c - s3, s3, s3) - \
+                   _integ(img, r + 1, c + 1, s3, s3)
             dxy = -dxy / w / w
 
-            dxx = integ(img, r - s3 + 1, c - s2, 2 * s3 - 1,w) - \
-                  integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3) * 3
+            dxx = _integ(img, r - s3 + 1, c - s2, 2 * s3 - 1,w) - \
+                  _integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3) * 3
             dxx = -dxx / w / w
 
-            dyy = integ(img, r - s2, c - s2 + 1, w, 2 * s3 - 1) - \
-                  integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1) * 3
+            dyy = _integ(img, r - s2, c - s2 + 1, w, 2 * s3 - 1) - \
+                  _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1) * 3
             dyy = -dyy / w / w
 
             out[r, c] = (dxx * dyy - 0.81 * (dxy * dxy))

From 2781d30bb736358fa6a2b1540bb67c0d6c8e405a Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Fri, 14 Mar 2014 21:09:31 +0530
Subject: [PATCH 03/22] doc strings

---
 skimage/feature/_hessian_det_appx.pyx | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index 92c4c08e..f5d0fc6e 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -6,6 +6,27 @@ from skimage import util
 
 
 cdef inline int  _clip(np.int_t x, np.int_t low, np.int_t high):
+    """Clips coordinate between high and low.
+    
+    This method was created so that `hessian_det_appx` does not have to make 
+    a Python call.
+
+    Parameters
+    ----------
+    x : int
+        Coordinate to be clipped.
+    low : int
+        The lower bound.
+    high : int
+        The higher bound.
+
+    Returns
+    -------
+    x : int
+        `x` clipped between 'high' and `low`.
+
+    """
+
     if(x > high):
         return high
     if(x < low):

From d58cab146fe31fdb7bed43ea5a3106b88ae14a51 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Fri, 14 Mar 2014 22:27:13 +0530
Subject: [PATCH 04/22] pep8 cython file

---
 skimage/feature/_hessian_det_appx.pyx | 106 ++++++++++++++++++++------
 1 file changed, 84 insertions(+), 22 deletions(-)

diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index f5d0fc6e..5f776ac6 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -7,8 +7,8 @@ from skimage import util
 
 cdef inline int  _clip(np.int_t x, np.int_t low, np.int_t high):
     """Clips coordinate between high and low.
-    
-    This method was created so that `hessian_det_appx` does not have to make 
+
+    This method was created so that `hessian_det_appx` does not have to make
     a Python call.
 
     Parameters
@@ -34,27 +34,80 @@ cdef inline int  _clip(np.int_t x, np.int_t low, np.int_t high):
     return x
 
 
-cdef inline int _integ(np.int_t[:, :] img, np.int_t r1, np.int_t c1, np.int_t rl, np.int_t cl):
+cdef inline int _integ(np.int_t[:, :] img, np.int_t r, np.int_t c,
+                       np.int_t rl, np.int_t cl):
+    """Integrate over the integral image in the given window
 
-    r1 = _clip(r1, 0, img.shape[0] - 1)
-    c1 = _clip(c1, 0, img.shape[1] - 1)
+    This method was created so that `hessian_det_appx` does not have to make
+    a Python call.
 
-    r2 = _clip(r1 + rl, 0, img.shape[0] - 1)
-    c2 = _clip(c1 + cl, 0, img.shape[1] - 1)
+    Parameters
+    ----------
+    img : array
+        The integral image over which to integrate.
+    r : int
+        The row number of the top left corner.
+    c : int
+        The column number of the top left corner.
+    rl : int
+        The number of rows over which to integrate.
+    cl : int
+        The number of columns over which to integrate.
 
-    cdef np.int_t r = img[r2, c2] + img[r1, c1] - img[r1, c2] - img[r2, c1]
+    Returns
+    -------
+    ans : int
+        The integral over the given window.
 
-    if (r < 0):
+    """
+
+    r = _clip(r, 0, img.shape[0] - 1)
+    c = _clip(c, 0, img.shape[1] - 1)
+
+    r2 = _clip(r + rl, 0, img.shape[0] - 1)
+    c2 = _clip(c + cl, 0, img.shape[1] - 1)
+
+    cdef np.int_t ans = img[r, c] + img[r2, c2] - img[r, c2] - img[r2, c]
+
+    if (ans < 0):
         return 0
-    return r
+    return ans
 
 
-def hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
+def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
+    """Computes the approximate Hessian Determinant over an image.
+
+    This method uses box filters over integral images to compute the
+    approximate Hessian Determinant as described in [1].
+
+    Parameters
+    ----------
+    image : array
+        The integral image over which to compute Hessian Determinant.
+    sigma : float
+        Standard deviation used for the Gaussian kernel, used for the Hessian
+        matrix
+
+    Returns
+    -------
+    out : array
+        The array of the Determinant of Hessians.
+
+    References
+    ----------
+    .. [1] ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf
+
+    Notes
+    -----
+    The running time of this method only depends on size of the image. It is
+    independent of `sigma` as one would expect. The downside is that the
+    result for `sigma` less than `3` is not accurate, i.e., not similar to
+    the result obtained if someone computed the Hessian and took it's
+    determinant.
+    """
 
     cdef np.int_t[:, :] img = image
     cdef int size = int(3 * sigma)
-    cdef np.ndarray[np.float_t, ndim = 2] out = np.zeros_like(img).astype(np.float)
-
     cdef int height = img.shape[0]
     cdef int width = img.shape[1]
     cdef int r, c
@@ -62,6 +115,9 @@ def hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
     cdef int s3 = size / 3
     cdef int l = size / 3
     cdef int w = size
+    cdef int mid, side
+    zeros = np.zeros_like(img)
+    cdef np.ndarray[np.float_t, ndim = 2] out = zeros.astype(np.float)
 
     cdef float dxx, dyy, dxy
 
@@ -70,19 +126,25 @@ def hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
 
     for r in range(height):
         for c in range(width):
-        
-            dxy =  _integ(img, r - s3, c + 1, s3, s3) + \
-                   _integ(img, r + 1, c - s3, s3, s3) - \
-                   _integ(img, r - s3, c - s3, s3, s3) - \
-                   _integ(img, r + 1, c + 1, s3, s3)
+
+            tl = _integ(img, r - s3, c - s3, s3, s3)  # top left
+            br = _integ(img, r + 1, c + 1, s3, s3)  # bottom right
+            bl = _integ(img, r - s3, c + 1, s3, s3)  # bottom left
+            tr = _integ(img, r + 1, c - s3, s3, s3)  # top right
+
+            dxy = bl + tr - tl - br
             dxy = -dxy / w / w
 
-            dxx = _integ(img, r - s3 + 1, c - s2, 2 * s3 - 1,w) - \
-                  _integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3) * 3
+            mid = _integ(img, r - s3 + 1, c - s2, 2 * s3 - 1, w)  # middle box
+            side = _integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3)  # sides
+
+            dxx = mid - 3 * side
             dxx = -dxx / w / w
 
-            dyy = _integ(img, r - s2, c - s2 + 1, w, 2 * s3 - 1) - \
-                  _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1) * 3
+            mid = _integ(img, r - s2, c - s2 + 1, w, 2 * s3 - 1)
+            side = _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1) * 3
+
+            dyy = mid - 3 * side
             dyy = -dyy / w / w
 
             out[r, c] = (dxx * dyy - 0.81 * (dxy * dxy))

From cca0c74cf7e2035f263a760ccf61161cc525b627 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Fri, 14 Mar 2014 23:46:17 +0530
Subject: [PATCH 05/22] added python call and doc strings

---
 skimage/feature/__init__.py           |   3 +-
 skimage/feature/_hessian_det_appx.pyx |   8 +-
 skimage/feature/blob.py               | 110 +++++++++++++++++++++++++-
 3 files changed, 115 insertions(+), 6 deletions(-)

diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py
index 54b7bac3..49167a36 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, blob_log
+from .blob import blob_dog, blob_log, blob_doh
 
 
 __all__ = ['daisy',
@@ -43,4 +43,5 @@ __all__ = ['daisy',
            'match_descriptors',
            'plot_matches',
            'blob_dog',
+           'blob_doh',
            'blob_log']
diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index 5f776ac6..f4bd213d 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -1,6 +1,5 @@
 import numpy as np
 cimport numpy as np
-
 from skimage.transform import integral_image, integrate
 from skimage import util
 
@@ -113,8 +112,9 @@ def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
     cdef int r, c
     cdef int s2 = (size - 1) / 2
     cdef int s3 = size / 3
-    cdef int l = size / 3
+    cdef int l = size/3
     cdef int w = size
+    cdef int b = (size - 1)/2
     cdef int mid, side
     zeros = np.zeros_like(img)
     cdef np.ndarray[np.float_t, ndim = 2] out = zeros.astype(np.float)
@@ -141,8 +141,8 @@ def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
             dxx = mid - 3 * side
             dxx = -dxx / w / w
 
-            mid = _integ(img, r - s2, c - s2 + 1, w, 2 * s3 - 1)
-            side = _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1) * 3
+            mid = _integ(img, r - s2, c - s3 + 1, w, 2 * s3 - 1)
+            side = _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1)
 
             dyy = mid - 3 * side
             dyy = -dyy / w / w
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 1b465c5f..aa373c86 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -4,8 +4,11 @@ import itertools as itt
 import math
 from math import sqrt, hypot, log
 from numpy import arccos
-from skimage.util import img_as_float
+from skimage.util import img_as_float, img_as_ubyte
 from .peak import peak_local_max
+from ._hessian_det_appx import _hessian_det_appx
+from skimage.transform import integral_image
+
 
 
 # This basic blob detection algorithm is based on:
@@ -298,3 +301,108 @@ def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2,
     # Convert the last index to its corresponding scale value
     local_maxima[:, 2] = sigma_list[local_maxima[:, 2]]
     return _prune_blobs(local_maxima, overlap)
+
+
+def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
+             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 : 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 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_Laplacian_of_Gaussian
+    .. [2] ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf
+
+    Examples
+    --------
+    >>> from skimage import data, feature, exposure
+    >>> img = data.coins()
+    >>> img = exposure.equalize_hist(img)  # improves detection
+    >>> feature.blob_log(img, threshold = .3)
+    array([[113, 323,   1],
+           [121, 272,  17],
+           [124, 336,  11],
+           [126,  46,  11],
+           [126, 208,  11],
+           [127, 102,  11],
+           [128, 154,  11],
+           [185, 344,  17],
+           [194, 213,  17],
+           [194, 276,  17],
+           [197,  44,  11],
+           [198, 103,  11],
+           [198, 155,  11],
+           [260, 174,  17],
+           [263, 244,  17],
+           [263, 302,  17],
+           [266, 115,  11]])
+
+    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
+    mathods line :py:meth:`blob_dog` and :py:math:`blob_log` the computation
+    of Gaussians for larger `sigma` takes more time.
+    """
+    if image.ndim != 2:
+        raise ValueError("'image' must be a grayscale ")
+
+    image = img_as_ubyte(image).astype(np.uint8)
+    image = integral_image(image).astype(np.int)
+    print 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_det_appx(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, 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]]
+    return _prune_blobs(local_maxima, overlap)

From a4029c0e71914845248495905f96332ff47699f4 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sat, 15 Mar 2014 00:02:26 +0530
Subject: [PATCH 06/22] removed print

---
 skimage/feature/blob.py            | 22 +++++++--------
 skimage/feature/tests/test_blob.py | 44 +++++++++++++++++++++++++++++-
 2 files changed, 53 insertions(+), 13 deletions(-)

diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index aa373c86..0be10f39 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -10,7 +10,6 @@ from ._hessian_det_appx import _hessian_det_appx
 from skimage.transform import integral_image
 
 
-
 # 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)
@@ -288,8 +287,8 @@ def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2,
     else:
         sigma_list = np.linspace(min_sigma, max_sigma, num_sigma)
 
-    #computing gaussian laplace
-    #s**2 provides scale invariance
+    # computing gaussian laplace
+    # s**2 provides scale invariance
     gl_images = [-gaussian_laplace(image, s) * s ** 2 for s in sigma_list]
     image_cube = np.dstack(gl_images)
 
@@ -307,10 +306,10 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
              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 
+    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]
+    detected the blob. Determinant of Hessians is approximated using [2]_
 
     Parameters
     ----------
@@ -318,10 +317,10 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
         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 used to compute 
+        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 
+        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
@@ -347,7 +346,7 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
 
     References
     ----------
-    .. [1] http://en.wikipedia.org/wiki/Blob_detection#The_Laplacian_of_Gaussian
+    .. [1] http://en.wikipedia.org/wiki/Blob_detection#The_determinant_of_the_Hessian
     .. [2] ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf
 
     Examples
@@ -377,9 +376,9 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
     Notes
     -----
     The radius of each blob is approximately `sigma`.
-    Computation of Determinant of Hessians is independent of the standard 
+    Computation of Determinant of Hessians is independent of the standard
     deviation. Therefore detecting larger blobs won't take more time. In
-    mathods line :py:meth:`blob_dog` and :py:math:`blob_log` the computation
+    methods line :py:meth:`blob_dog` and :py:meth:`blob_log` the computation
     of Gaussians for larger `sigma` takes more time.
     """
     if image.ndim != 2:
@@ -387,14 +386,13 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
 
     image = img_as_ubyte(image).astype(np.uint8)
     image = integral_image(image).astype(np.int)
-    print 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_det_appx(image, s) for s in sigma_list]
     image_cube = np.dstack(hessian_images)
 
diff --git a/skimage/feature/tests/test_blob.py b/skimage/feature/tests/test_blob.py
index 3ae949c2..b68f144d 100644
--- a/skimage/feature/tests/test_blob.py
+++ b/skimage/feature/tests/test_blob.py
@@ -1,6 +1,6 @@
 import numpy as np
 from skimage.draw import circle
-from skimage.feature import blob_dog, blob_log
+from skimage.feature import blob_dog, blob_log, blob_doh
 import math
 
 
@@ -79,3 +79,45 @@ def test_blob_log():
     assert abs(b[0] - 200) <= thresh
     assert abs(b[1] - 350) <= thresh
     assert abs(radius(b) - 30) <= thresh
+    
+def test_blob_doh():
+    r2 = math.sqrt(2)
+    img = np.ones((512, 512), dtype = np.uint8)
+
+    xs, ys = circle(400, 130, 20)
+    img[xs, ys] = 255
+
+    xs, ys = circle(160, 50, 30)
+    img[xs, ys] = 255
+
+    xs, ys = circle(100, 300, 40)
+    img[xs, ys] = 255
+
+    xs, ys = circle(200, 350, 50)
+    img[xs, ys] = 255
+
+    blobs = blob_doh(img, min_sigma=1, max_sigma=60, num_sigma=10)
+
+    radius = lambda x: r2*x[2]
+    s = sorted(blobs, key=radius)
+    thresh = 3
+
+    b = s[0]
+    assert abs(b[0] - 400) <= thresh
+    assert abs(b[1] - 130) <= thresh
+    assert abs(radius(b) - 20) <= thresh
+
+    b = s[1]
+    assert abs(b[0] - 160) <= thresh
+    assert abs(b[1] - 50) <= thresh
+    assert abs(radius(b) - 30) <= thresh
+
+    b = s[2]
+    assert abs(b[0] - 100) <= thresh
+    assert abs(b[1] - 300) <= thresh
+    assert abs(radius(b) - 40) <= thresh
+
+    b = s[3]
+    assert abs(b[0] - 200) <= thresh
+    assert abs(b[1] - 350) <= thresh
+    assert abs(radius(b) - 50) <= thresh

From 2e441b2e5b4e6fbfe2a08ea38291f352ff973b73 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sat, 15 Mar 2014 00:16:39 +0530
Subject: [PATCH 07/22] tests and minor corrections

---
 skimage/feature/blob.py            | 39 +++++++++++++++---------------
 skimage/feature/tests/test_blob.py | 23 +++++++++++-------
 2 files changed, 33 insertions(+), 29 deletions(-)

diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 0be10f39..e47c60c9 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -351,27 +351,25 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
 
     Examples
     --------
-    >>> from skimage import data, feature, exposure
+    >>> from skimage import data, feature
     >>> img = data.coins()
-    >>> img = exposure.equalize_hist(img)  # improves detection
-    >>> feature.blob_log(img, threshold = .3)
-    array([[113, 323,   1],
-           [121, 272,  17],
-           [124, 336,  11],
-           [126,  46,  11],
-           [126, 208,  11],
-           [127, 102,  11],
-           [128, 154,  11],
-           [185, 344,  17],
-           [194, 213,  17],
-           [194, 276,  17],
-           [197,  44,  11],
-           [198, 103,  11],
-           [198, 155,  11],
-           [260, 174,  17],
-           [263, 244,  17],
-           [263, 302,  17],
-           [266, 115,  11]])
+    >>> feature.blob_doh(img,threshold = 700)
+    array([[121, 271,  30],
+           [123,  44,  23],
+           [123, 205,  20],
+           [124, 336,  20],
+           [126, 101,  20],
+           [126, 153,  20],
+           [156, 302,  30],
+           [185, 348,  30],
+           [192, 212,  23],
+           [193, 275,  23],
+           [195, 100,  23],
+           [197, 153,  20],
+           [260, 173,  30],
+           [262, 243,  23],
+           [265, 113,  23],
+           [270, 363,  30]])
 
     Notes
     -----
@@ -381,6 +379,7 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
     methods line :py:meth:`blob_dog` and :py:meth:`blob_log` the computation
     of Gaussians for larger `sigma` takes more time.
     """
+    
     if image.ndim != 2:
         raise ValueError("'image' must be a grayscale ")
 
diff --git a/skimage/feature/tests/test_blob.py b/skimage/feature/tests/test_blob.py
index b68f144d..4166ff2a 100644
--- a/skimage/feature/tests/test_blob.py
+++ b/skimage/feature/tests/test_blob.py
@@ -18,7 +18,7 @@ def test_blob_dog():
     img[xs, ys] = 255
 
     blobs = blob_dog(img, min_sigma=5, max_sigma=50)
-    radius = lambda x: r2*x[2]
+    radius = lambda x: r2 * x[2]
     s = sorted(blobs, key=radius)
     thresh = 5
 
@@ -56,7 +56,7 @@ def test_blob_log():
 
     blobs = blob_log(img, min_sigma=5, max_sigma=20, threshold=1)
 
-    radius = lambda x: r2*x[2]
+    radius = lambda x: r2 * x[2]
     s = sorted(blobs, key=radius)
     thresh = 3
 
@@ -79,15 +79,15 @@ def test_blob_log():
     assert abs(b[0] - 200) <= thresh
     assert abs(b[1] - 350) <= thresh
     assert abs(radius(b) - 30) <= thresh
-    
+
+
 def test_blob_doh():
-    r2 = math.sqrt(2)
-    img = np.ones((512, 512), dtype = np.uint8)
+    img = np.ones((512, 512), dtype=np.uint8)
 
     xs, ys = circle(400, 130, 20)
     img[xs, ys] = 255
 
-    xs, ys = circle(160, 50, 30)
+    xs, ys = circle(460, 50, 30)
     img[xs, ys] = 255
 
     xs, ys = circle(100, 300, 40)
@@ -96,9 +96,14 @@ def test_blob_doh():
     xs, ys = circle(200, 350, 50)
     img[xs, ys] = 255
 
-    blobs = blob_doh(img, min_sigma=1, max_sigma=60, num_sigma=10)
+    blobs = blob_doh(
+        img,
+        min_sigma=1,
+        max_sigma=60,
+        num_sigma=10,
+        threshold=1000)
 
-    radius = lambda x: r2*x[2]
+    radius = lambda x: x[2]
     s = sorted(blobs, key=radius)
     thresh = 3
 
@@ -108,7 +113,7 @@ def test_blob_doh():
     assert abs(radius(b) - 20) <= thresh
 
     b = s[1]
-    assert abs(b[0] - 160) <= thresh
+    assert abs(b[0] - 460) <= thresh
     assert abs(b[1] - 50) <= thresh
     assert abs(radius(b) - 30) <= thresh
 

From 507a2159b94b337ddd57022a11d70dd51f4ecc5e Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sat, 15 Mar 2014 00:22:35 +0530
Subject: [PATCH 08/22] notes change

---
 skimage/feature/blob.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index e47c60c9..2bb03db8 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -377,7 +377,9 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
     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.
+    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.
     """
     
     if image.ndim != 2:

From 3cf6bffe84cd5bacae8e0c81701961c4dcbb31bf Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sat, 15 Mar 2014 00:44:49 +0530
Subject: [PATCH 09/22] bento change

---
 bento.info | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/bento.info b/bento.info
index 45f57a36..f1cbabc1 100644
--- a/bento.info
+++ b/bento.info
@@ -151,6 +151,9 @@ Library:
     Extension: skimage.restoration._denoise_cy
         Sources:
             skimage/restoration/_denoise_cy.pyx
+    Extension: skimage.feature._hessian_det_appx
+        Sources:
+            skimage/exposure/_hessian_det_appx.pyx
 
 Executable: skivi
     Module: skimage.scripts.skivi

From d5f9ebfc8f28484c865dd2e52e9eec2ad6d499ce Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sat, 15 Mar 2014 01:10:58 +0530
Subject: [PATCH 10/22] minor doc changes

---
 skimage/feature/_hessian_det_appx.pyx | 4 ++--
 skimage/feature/blob.py               | 5 ++---
 2 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index f4bd213d..225f340c 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -22,7 +22,7 @@ cdef inline int  _clip(np.int_t x, np.int_t low, np.int_t high):
     Returns
     -------
     x : int
-        `x` clipped between 'high' and `low`.
+        `x` clipped between `high` and `low`.
 
     """
 
@@ -77,7 +77,7 @@ def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
     """Computes the approximate Hessian Determinant over an image.
 
     This method uses box filters over integral images to compute the
-    approximate Hessian Determinant as described in [1].
+    approximate Hessian Determinant as described in [1]_.
 
     Parameters
     ----------
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 2bb03db8..c0126a52 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -309,13 +309,12 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
     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]_
+    detected the blob. Determinant of Hessians is approximated using [2]_.
 
     Parameters
     ----------
     image : ndarray
-        Input grayscale image, blobs are assumed to be light on dark
-        background (white on black).
+        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.

From 630136ae91f40ea5c0d9833caf1e26aec89a5753 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sat, 15 Mar 2014 02:18:33 +0530
Subject: [PATCH 11/22] added full paper desc.

---
 skimage/feature/_hessian_det_appx.pyx | 4 +++-
 skimage/feature/blob.py               | 9 ++++++---
 2 files changed, 9 insertions(+), 4 deletions(-)

diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index 225f340c..f620c35b 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -94,7 +94,9 @@ def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
 
     References
     ----------
-    .. [1] ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf
+    .. [1] 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
 
     Notes
     -----
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index c0126a52..3d602b3a 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -346,7 +346,10 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
     References
     ----------
     .. [1] http://en.wikipedia.org/wiki/Blob_detection#The_determinant_of_the_Hessian
-    .. [2] ftp://ftp.vision.ee.ethz.ch/publications/articles/eth_biwi_00517.pdf
+
+    .. [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
     --------
@@ -377,10 +380,10 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
     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` 
+    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.
     """
-    
+
     if image.ndim != 2:
         raise ValueError("'image' must be a grayscale ")
 

From fbf40b0ae7d55ad2479ccecc61e657bab213cbb4 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sat, 15 Mar 2014 13:43:10 +0530
Subject: [PATCH 12/22] corrections to cython file, switched dtype to double

---
 skimage/feature/_hessian_det_appx.pyx | 50 +++++++++++++++------------
 skimage/feature/blob.py               | 10 +++---
 skimage/feature/tests/test_blob.py    |  2 +-
 3 files changed, 34 insertions(+), 28 deletions(-)

diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index f620c35b..358a4973 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -1,10 +1,14 @@
+# cython: cdivision=True
+# cython: boundscheck=False
+# cython: nonecheck=False
+# cython: wraparound=False
 import numpy as np
-cimport numpy as np
+cimport numpy as cnp
 from skimage.transform import integral_image, integrate
 from skimage import util
 
 
-cdef inline int  _clip(np.int_t x, np.int_t low, np.int_t high):
+cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low, Py_ssize_t high):
     """Clips coordinate between high and low.
 
     This method was created so that `hessian_det_appx` does not have to make
@@ -33,8 +37,9 @@ cdef inline int  _clip(np.int_t x, np.int_t low, np.int_t high):
     return x
 
 
-cdef inline int _integ(np.int_t[:, :] img, np.int_t r, np.int_t c,
-                       np.int_t rl, np.int_t cl):
+cdef inline cnp.double_t _integ(
+    cnp.double_t[:, :] img, Py_ssize_t r, Py_ssize_t c,
+        Py_ssize_t rl, Py_ssize_t cl):
     """Integrate over the integral image in the given window
 
     This method was created so that `hessian_det_appx` does not have to make
@@ -66,14 +71,14 @@ cdef inline int _integ(np.int_t[:, :] img, np.int_t r, np.int_t c,
     r2 = _clip(r + rl, 0, img.shape[0] - 1)
     c2 = _clip(c + cl, 0, img.shape[1] - 1)
 
-    cdef np.int_t ans = img[r, c] + img[r2, c2] - img[r, c2] - img[r2, c]
+    cdef cnp.double_t ans = img[r, c] + img[r2, c2] - img[r, c2] - img[r2, c]
 
     if (ans < 0):
         return 0
     return ans
 
 
-def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
+def _hessian_det_appx(cnp.double_t[:, :] img, float sigma):
     """Computes the approximate Hessian Determinant over an image.
 
     This method uses box filters over integral images to compute the
@@ -81,7 +86,7 @@ def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
 
     Parameters
     ----------
-    image : array
+    img : array
         The integral image over which to compute Hessian Determinant.
     sigma : float
         Standard deviation used for the Gaussian kernel, used for the Hessian
@@ -107,19 +112,18 @@ def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
     determinant.
     """
 
-    cdef np.int_t[:, :] img = image
-    cdef int size = int(3 * sigma)
-    cdef int height = img.shape[0]
-    cdef int width = img.shape[1]
-    cdef int r, c
-    cdef int s2 = (size - 1) / 2
-    cdef int s3 = size / 3
-    cdef int l = size/3
-    cdef int w = size
-    cdef int b = (size - 1)/2
-    cdef int mid, side
-    zeros = np.zeros_like(img)
-    cdef np.ndarray[np.float_t, ndim = 2] out = zeros.astype(np.float)
+    cdef Py_ssize_t size = int(3 * sigma)
+    cdef Py_ssize_t height = img.shape[0]
+    cdef Py_ssize_t width = img.shape[1]
+    cdef Py_ssize_t r, c
+    cdef Py_ssize_t s2 = (size - 1) / 2
+    cdef Py_ssize_t s3 = size / 3
+    cdef Py_ssize_t l = size / 3
+    cdef Py_ssize_t w = size
+    cdef Py_ssize_t b = (size - 1) / 2
+    cdef cnp.double_t mid, side, tl, tr, bl, br
+    cdef cnp.double_t[:, ::1] out = np.zeros_like(img, dtype=np.double)
+    cdef cnp.double_t w_i = 1.0 / size / size
 
     cdef float dxx, dyy, dxy
 
@@ -135,19 +139,19 @@ def _hessian_det_appx(np.ndarray[np.int_t, ndim=2] image, float sigma):
             tr = _integ(img, r + 1, c - s3, s3, s3)  # top right
 
             dxy = bl + tr - tl - br
-            dxy = -dxy / w / w
+            dxy = -dxy * w_i
 
             mid = _integ(img, r - s3 + 1, c - s2, 2 * s3 - 1, w)  # middle box
             side = _integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3)  # sides
 
             dxx = mid - 3 * side
-            dxx = -dxx / w / w
+            dxx = -dxx * w_i
 
             mid = _integ(img, r - s2, c - s3 + 1, w, 2 * s3 - 1)
             side = _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1)
 
             dyy = mid - 3 * side
-            dyy = -dyy / w / w
+            dyy = -dyy * w_i
 
             out[r, c] = (dxx * dyy - 0.81 * (dxy * dxy))
 
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 3d602b3a..7c6a5649 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -4,7 +4,7 @@ import itertools as itt
 import math
 from math import sqrt, hypot, log
 from numpy import arccos
-from skimage.util import img_as_float, img_as_ubyte
+from skimage.util import img_as_float
 from .peak import peak_local_max
 from ._hessian_det_appx import _hessian_det_appx
 from skimage.transform import integral_image
@@ -302,7 +302,7 @@ def blob_log(image, min_sigma=1, max_sigma=50, num_sigma=10, threshold=.2,
     return _prune_blobs(local_maxima, overlap)
 
 
-def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
+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.
 
@@ -367,12 +367,14 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
            [192, 212,  23],
            [193, 275,  23],
            [195, 100,  23],
+           [197,  44,  20],
            [197, 153,  20],
            [260, 173,  30],
            [262, 243,  23],
            [265, 113,  23],
            [270, 363,  30]])
 
+
     Notes
     -----
     The radius of each blob is approximately `sigma`.
@@ -387,8 +389,8 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=500,
     if image.ndim != 2:
         raise ValueError("'image' must be a grayscale ")
 
-    image = img_as_ubyte(image).astype(np.uint8)
-    image = integral_image(image).astype(np.int)
+    image = img_as_float(image)
+    image = integral_image(image)
 
     if log_scale:
         start, stop = log(min_sigma, 10), log(max_sigma, 10)
diff --git a/skimage/feature/tests/test_blob.py b/skimage/feature/tests/test_blob.py
index 4166ff2a..10167214 100644
--- a/skimage/feature/tests/test_blob.py
+++ b/skimage/feature/tests/test_blob.py
@@ -101,7 +101,7 @@ def test_blob_doh():
         min_sigma=1,
         max_sigma=60,
         num_sigma=10,
-        threshold=1000)
+        threshold=.05)
 
     radius = lambda x: x[2]
     s = sorted(blobs, key=radius)

From 1b2ca1145b8810cc203461463a437cd25b5d6dcd Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sat, 15 Mar 2014 13:45:38 +0530
Subject: [PATCH 13/22] doc test

---
 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 7c6a5649..96de1d65 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -355,7 +355,7 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=0.01,
     --------
     >>> from skimage import data, feature
     >>> img = data.coins()
-    >>> feature.blob_doh(img,threshold = 700)
+    >>> feature.blob_doh(img)
     array([[121, 271,  30],
            [123,  44,  23],
            [123, 205,  20],

From a15c5fac87cc68177b0d854791aa01e99ab00ef5 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sat, 15 Mar 2014 19:58:01 +0530
Subject: [PATCH 14/22] indexing

---
 skimage/feature/_hessian_det_appx.pyx | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index 358a4973..a8ca2734 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -38,7 +38,7 @@ cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low, Py_ssize_t high):
 
 
 cdef inline cnp.double_t _integ(
-    cnp.double_t[:, :] img, Py_ssize_t r, Py_ssize_t c,
+    cnp.double_t[:, ::1] img, Py_ssize_t r, Py_ssize_t c,
         Py_ssize_t rl, Py_ssize_t cl):
     """Integrate over the integral image in the given window
 
@@ -78,7 +78,7 @@ cdef inline cnp.double_t _integ(
     return ans
 
 
-def _hessian_det_appx(cnp.double_t[:, :] img, float sigma):
+def _hessian_det_appx(cnp.double_t[:, ::1] img, float sigma):
     """Computes the approximate Hessian Determinant over an image.
 
     This method uses box filters over integral images to compute the

From 61a8a78657bd6fd996747fea9664bc7e9a80c8ea Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Sat, 15 Mar 2014 20:02:09 +0530
Subject: [PATCH 15/22] made hessian_det_appx public

---
 skimage/feature/__init__.py           | 2 ++
 skimage/feature/_hessian_det_appx.pyx | 2 +-
 skimage/feature/blob.py               | 4 ++--
 3 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py
index 49167a36..7d1d9add 100644
--- a/skimage/feature/__init__.py
+++ b/skimage/feature/__init__.py
@@ -15,6 +15,7 @@ from .orb import ORB
 from .match import match_descriptors
 from .util import plot_matches
 from .blob import blob_dog, blob_log, blob_doh
+from ._hessian_det_appx import hessian_det_appx
 
 
 __all__ = ['daisy',
@@ -44,4 +45,5 @@ __all__ = ['daisy',
            'plot_matches',
            'blob_dog',
            'blob_doh',
+           'hessian_det_appx',
            'blob_log']
diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index a8ca2734..5d3f7511 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -78,7 +78,7 @@ cdef inline cnp.double_t _integ(
     return ans
 
 
-def _hessian_det_appx(cnp.double_t[:, ::1] img, float sigma):
+def hessian_det_appx(cnp.double_t[:, ::1] img, float sigma):
     """Computes the approximate Hessian Determinant over an image.
 
     This method uses box filters over integral images to compute the
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 96de1d65..57818a98 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -6,7 +6,7 @@ from math import sqrt, hypot, log
 from numpy import arccos
 from skimage.util import img_as_float
 from .peak import peak_local_max
-from ._hessian_det_appx import _hessian_det_appx
+from ._hessian_det_appx import hessian_det_appx
 from skimage.transform import integral_image
 
 
@@ -398,7 +398,7 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=0.01,
     else:
         sigma_list = np.linspace(min_sigma, max_sigma, num_sigma)
 
-    hessian_images = [_hessian_det_appx(image, s) for s in sigma_list]
+    hessian_images = [hessian_det_appx(image, s) for s in sigma_list]
     image_cube = np.dstack(hessian_images)
 
     local_maxima = peak_local_max(image_cube, threshold_abs=threshold,

From e818a42fba281fef30cfb75f2fcd620adc38b79b Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Fri, 21 Mar 2014 03:49:31 +0530
Subject: [PATCH 16/22] moved public interface of hessian determinant to
 corner.py

---
 skimage/feature/__init__.py           |  8 ++---
 skimage/feature/_hessian_det_appx.pyx |  2 +-
 skimage/feature/blob.py               |  5 +--
 skimage/feature/corner.py             | 50 ++++++++++++++++++++++++++-
 skimage/feature/tests/test_corner.py  |  9 ++++-
 5 files changed, 65 insertions(+), 9 deletions(-)

diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py
index 7d1d9add..6a933056 100644
--- a/skimage/feature/__init__.py
+++ b/skimage/feature/__init__.py
@@ -6,7 +6,7 @@ from .corner import (corner_kitchen_rosenfeld, corner_harris,
                      corner_shi_tomasi, corner_foerstner, corner_subpix,
                      corner_peaks, corner_fast, structure_tensor,
                      structure_tensor_eigvals, hessian_matrix,
-                     hessian_matrix_eigvals)
+                     hessian_matrix_eigvals, hessian_matrix_det)
 from .corner_cy import corner_moravec, corner_orientations
 from .template import match_template
 from .brief import BRIEF
@@ -15,7 +15,6 @@ from .orb import ORB
 from .match import match_descriptors
 from .util import plot_matches
 from .blob import blob_dog, blob_log, blob_doh
-from ._hessian_det_appx import hessian_det_appx
 
 
 __all__ = ['daisy',
@@ -27,6 +26,7 @@ __all__ = ['daisy',
            'structure_tensor',
            'structure_tensor_eigvals',
            'hessian_matrix',
+           'hessian_matrx_det',
            'hessian_matrix_eigvals',
            'corner_kitchen_rosenfeld',
            'corner_harris',
@@ -45,5 +45,5 @@ __all__ = ['daisy',
            'plot_matches',
            'blob_dog',
            'blob_doh',
-           'hessian_det_appx',
-           'blob_log']
+           'blob_log',
+           'hessian_matrix_det']
diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index 5d3f7511..e66c0d53 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -78,7 +78,7 @@ cdef inline cnp.double_t _integ(
     return ans
 
 
-def hessian_det_appx(cnp.double_t[:, ::1] img, float sigma):
+def _hessian_matrix_det(cnp.double_t[:, ::1] img, float sigma):
     """Computes the approximate Hessian Determinant over an image.
 
     This method uses box filters over integral images to compute the
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 57818a98..9c35d2a5 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -1,3 +1,4 @@
+
 import numpy as np
 from scipy.ndimage.filters import gaussian_filter, gaussian_laplace
 import itertools as itt
@@ -6,7 +7,7 @@ from math import sqrt, hypot, log
 from numpy import arccos
 from skimage.util import img_as_float
 from .peak import peak_local_max
-from ._hessian_det_appx import hessian_det_appx
+from ._hessian_det_appx import _hessian_matrix_det
 from skimage.transform import integral_image
 
 
@@ -398,7 +399,7 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=0.01,
     else:
         sigma_list = np.linspace(min_sigma, max_sigma, num_sigma)
 
-    hessian_images = [hessian_det_appx(image, s) for s in sigma_list]
+    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,
diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py
index f4894aa5..67bed0a1 100644
--- a/skimage/feature/corner.py
+++ b/skimage/feature/corner.py
@@ -7,6 +7,8 @@ from skimage.util import img_as_float, pad
 from skimage.feature import peak_local_max
 from skimage.feature.util import _prepare_grayscale_input_2D
 from skimage.feature.corner_cy import _corner_fast
+from ._hessian_det_appx import _hessian_matrix_det
+from ..transform import integral_image
 
 
 def _compute_derivatives(image, mode='constant', cval=0):
@@ -29,7 +31,7 @@ def _compute_derivatives(image, mode='constant', cval=0):
     imy : ndarray
         Derivative in y-direction.
 
-    """
+v    """
 
     imy = ndimage.sobel(image, axis=0, mode=mode, cval=cval)
     imx = ndimage.sobel(image, axis=1, mode=mode, cval=cval)
@@ -170,6 +172,52 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0):
     return Hxx, Hxy, Hyy
 
 
+def hessian_matrix_det(image, sigma, integral=True):
+    """Computes the approximate Hessian Determinant over an image.
+
+    This method uses box filters over integral images to compute the
+    approximate Hessian Determinant as described in [1]_.
+
+    Parameters
+    ----------
+    image : array
+        The integral image over which to compute Hessian Determinant.
+    sigma : float
+        Standard deviation used for the Gaussian kernel, used for the Hessian
+        matrix
+    integral : bool
+        If `False`, `image` is assumed to be integral and intergral image is
+        not computed. If `True` the integral image is computed for `image`
+        and used for finding the Hessian Determinant.
+
+
+    Returns
+    -------
+    out : array
+        The array of the Determinant of Hessians.
+
+    References
+    ----------
+    .. [1] 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
+
+    Notes
+    -----
+    The running time of this method only depends on size of the image. It is
+    independent of `sigma` as one would expect. The downside is that the
+    result for `sigma` less than `3` is not accurate, i.e., not similar to
+    the result obtained if someone computed the Hessian and took it's
+    determinant.
+    """
+
+    image = img_as_float(image)
+    if(integral):
+        image = integral_image(image)
+
+    return np.array(_hessian_matrix_det(image, sigma))
+
+
 def _image_orthogonal_matrix22_eigvals(M00, M01, M11):
     l1 = (M00 + M11) / 2 + np.sqrt(4 * M01 ** 2 + (M00 - M11) ** 2) / 2
     l2 = (M00 + M11) / 2 - np.sqrt(4 * M01 ** 2 + (M00 - M11) ** 2) / 2
diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py
index 21070779..dd9f7703 100644
--- a/skimage/feature/tests/test_corner.py
+++ b/skimage/feature/tests/test_corner.py
@@ -12,7 +12,8 @@ from skimage.feature import (corner_moravec, corner_harris, corner_shi_tomasi,
                              corner_kitchen_rosenfeld, corner_foerstner,
                              corner_fast, corner_orientations,
                              structure_tensor, structure_tensor_eigvals,
-                             hessian_matrix, hessian_matrix_eigvals)
+                             hessian_matrix, hessian_matrix_eigvals,
+                             hessian_matrix_det)
 
 
 def test_structure_tensor():
@@ -91,6 +92,12 @@ def test_hessian_matrix_eigvals():
                                      [0, 0, 0, 0, 0]]))
 
 
+def test_hessian_matrix_det():
+    image = np.ones((5, 5))
+    det = hessian_matrix_det(image, 3, False)
+    assert_array_equal(det, 0)
+
+
 def test_square_image():
     im = np.zeros((50, 50)).astype(float)
     im[:25, :25] = 1.

From 7a4401e262f0ca79324727a6cb98189d1de38c0f Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Fri, 21 Mar 2014 03:54:52 +0530
Subject: [PATCH 17/22] doc string change

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

diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py
index 67bed0a1..a7b7c5b7 100644
--- a/skimage/feature/corner.py
+++ b/skimage/feature/corner.py
@@ -181,7 +181,7 @@ def hessian_matrix_det(image, sigma, integral=True):
     Parameters
     ----------
     image : array
-        The integral image over which to compute Hessian Determinant.
+        The image over which to compute Hessian Determinant.
     sigma : float
         Standard deviation used for the Gaussian kernel, used for the Hessian
         matrix

From a66f2a0b373881231f2d76b1f446362742365470 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Wed, 9 Apr 2014 20:30:33 +0530
Subject: [PATCH 18/22] extra line and spelling mistake

---
 skimage/feature/blob.py   | 2 +-
 skimage/feature/corner.py | 1 -
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 9c35d2a5..343863a9 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -388,7 +388,7 @@ def blob_doh(image, min_sigma=1, max_sigma=30, num_sigma=10, threshold=0.01,
     """
 
     if image.ndim != 2:
-        raise ValueError("'image' must be a grayscale ")
+        raise ValueError("'image' must be grayscale ")
 
     image = img_as_float(image)
     image = integral_image(image)
diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py
index a7b7c5b7..1889d799 100644
--- a/skimage/feature/corner.py
+++ b/skimage/feature/corner.py
@@ -190,7 +190,6 @@ def hessian_matrix_det(image, sigma, integral=True):
         not computed. If `True` the integral image is computed for `image`
         and used for finding the Hessian Determinant.
 
-
     Returns
     -------
     out : array

From a4a939e74bd684e9e94add262f4d0f474de410de Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Wed, 9 Apr 2014 23:19:25 +0530
Subject: [PATCH 19/22] improved test coverage

---
 skimage/feature/tests/test_blob.py | 86 ++++++++++++++++++++++++++++++
 1 file changed, 86 insertions(+)

diff --git a/skimage/feature/tests/test_blob.py b/skimage/feature/tests/test_blob.py
index 10167214..a46c442b 100644
--- a/skimage/feature/tests/test_blob.py
+++ b/skimage/feature/tests/test_blob.py
@@ -2,11 +2,13 @@ import numpy as np
 from skimage.draw import circle
 from skimage.feature import blob_dog, blob_log, blob_doh
 import math
+from numpy.testing import assert_raises
 
 
 def test_blob_dog():
     r2 = math.sqrt(2)
     img = np.ones((512, 512))
+    img3 = np.ones((5, 5, 5))
 
     xs, ys = circle(400, 130, 5)
     img[xs, ys] = 255
@@ -37,10 +39,13 @@ def test_blob_dog():
     assert abs(b[1] - 350) <= thresh
     assert abs(radius(b) - 45) <= thresh
 
+    assert_raises(ValueError, blob_dog, img3)
+
 
 def test_blob_log():
     r2 = math.sqrt(2)
     img = np.ones((512, 512))
+    img3 = np.ones((5, 5, 5))
 
     xs, ys = circle(400, 130, 5)
     img[xs, ys] = 255
@@ -80,9 +85,40 @@ def test_blob_log():
     assert abs(b[1] - 350) <= thresh
     assert abs(radius(b) - 30) <= thresh
 
+    # Testing log scale
+    blobs = blob_log(
+        img,
+        min_sigma=5,
+        max_sigma=20,
+        threshold=1,
+        log_scale=True)
+
+    b = s[0]
+    assert abs(b[0] - 400) <= thresh
+    assert abs(b[1] - 130) <= thresh
+    assert abs(radius(b) - 5) <= thresh
+
+    b = s[1]
+    assert abs(b[0] - 160) <= thresh
+    assert abs(b[1] - 50) <= thresh
+    assert abs(radius(b) - 15) <= thresh
+
+    b = s[2]
+    assert abs(b[0] - 100) <= thresh
+    assert abs(b[1] - 300) <= thresh
+    assert abs(radius(b) - 25) <= thresh
+
+    b = s[3]
+    assert abs(b[0] - 200) <= thresh
+    assert abs(b[1] - 350) <= thresh
+    assert abs(radius(b) - 30) <= thresh
+
+    assert_raises(ValueError, blob_log, img3)
+
 
 def test_blob_doh():
     img = np.ones((512, 512), dtype=np.uint8)
+    img3 = np.ones((5, 5, 5))
 
     xs, ys = circle(400, 130, 20)
     img[xs, ys] = 255
@@ -126,3 +162,53 @@ def test_blob_doh():
     assert abs(b[0] - 200) <= thresh
     assert abs(b[1] - 350) <= thresh
     assert abs(radius(b) - 50) <= thresh
+
+    # Testing log scale
+    blobs = blob_doh(
+        img,
+        min_sigma=1,
+        max_sigma=60,
+        num_sigma=10,
+        log_scale=True,
+        threshold=.05)
+
+    b = s[0]
+    assert abs(b[0] - 400) <= thresh
+    assert abs(b[1] - 130) <= thresh
+    assert abs(radius(b) - 20) <= thresh
+
+    b = s[1]
+    assert abs(b[0] - 460) <= thresh
+    assert abs(b[1] - 50) <= thresh
+    assert abs(radius(b) - 30) <= thresh
+
+    b = s[2]
+    assert abs(b[0] - 100) <= thresh
+    assert abs(b[1] - 300) <= thresh
+    assert abs(radius(b) - 40) <= thresh
+
+    b = s[3]
+    assert abs(b[0] - 200) <= thresh
+    assert abs(b[1] - 350) <= thresh
+    assert abs(radius(b) - 50) <= thresh
+
+    assert_raises(ValueError, blob_doh, img3)
+
+
+def test_blob_overlap():
+    img = np.ones((512, 512), dtype=np.uint8)
+
+    xs, ys = circle(100, 100, 20)
+    img[xs, ys] = 255
+
+    xs, ys = circle(120, 100, 30)
+    img[xs, ys] = 255
+
+    blobs = blob_doh(
+        img,
+        min_sigma=1,
+        max_sigma=60,
+        num_sigma=10,
+        threshold=.05)
+
+    assert len(blobs) == 1

From 5658704a355e8e2b8a216da42bc5f9cf9176f165 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Thu, 10 Apr 2014 19:54:49 +0530
Subject: [PATCH 20/22] changed parameter of hessian_matrix_det

---
 skimage/feature/__init__.py           |  5 ++---
 skimage/feature/_hessian_det_appx.pyx |  4 +---
 skimage/feature/corner.py             | 14 ++++----------
 skimage/feature/tests/test_corner.py  |  9 ++++++---
 4 files changed, 13 insertions(+), 19 deletions(-)

diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py
index 6a933056..c46fde01 100644
--- a/skimage/feature/__init__.py
+++ b/skimage/feature/__init__.py
@@ -26,7 +26,7 @@ __all__ = ['daisy',
            'structure_tensor',
            'structure_tensor_eigvals',
            'hessian_matrix',
-           'hessian_matrx_det',
+           'hessian_matrix_det',
            'hessian_matrix_eigvals',
            'corner_kitchen_rosenfeld',
            'corner_harris',
@@ -45,5 +45,4 @@ __all__ = ['daisy',
            'plot_matches',
            'blob_dog',
            'blob_doh',
-           'blob_log',
-           'hessian_matrix_det']
+           'blob_log']
diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index e66c0d53..7f5a403e 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -4,8 +4,6 @@
 # cython: wraparound=False
 import numpy as np
 cimport numpy as cnp
-from skimage.transform import integral_image, integrate
-from skimage import util
 
 
 cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low, Py_ssize_t high):
@@ -78,7 +76,7 @@ cdef inline cnp.double_t _integ(
     return ans
 
 
-def _hessian_matrix_det(cnp.double_t[:, ::1] img, float sigma):
+def _hessian_matrix_det(cnp.double_t[:, ::1] img, double sigma):
     """Computes the approximate Hessian Determinant over an image.
 
     This method uses box filters over integral images to compute the
diff --git a/skimage/feature/corner.py b/skimage/feature/corner.py
index 1889d799..e1dd519b 100644
--- a/skimage/feature/corner.py
+++ b/skimage/feature/corner.py
@@ -31,7 +31,7 @@ def _compute_derivatives(image, mode='constant', cval=0):
     imy : ndarray
         Derivative in y-direction.
 
-v    """
+    """
 
     imy = ndimage.sobel(image, axis=0, mode=mode, cval=cval)
     imx = ndimage.sobel(image, axis=1, mode=mode, cval=cval)
@@ -172,7 +172,7 @@ def hessian_matrix(image, sigma=1, mode='constant', cval=0):
     return Hxx, Hxy, Hyy
 
 
-def hessian_matrix_det(image, sigma, integral=True):
+def hessian_matrix_det(image, sigma):
     """Computes the approximate Hessian Determinant over an image.
 
     This method uses box filters over integral images to compute the
@@ -184,11 +184,7 @@ def hessian_matrix_det(image, sigma, integral=True):
         The image over which to compute Hessian Determinant.
     sigma : float
         Standard deviation used for the Gaussian kernel, used for the Hessian
-        matrix
-    integral : bool
-        If `False`, `image` is assumed to be integral and intergral image is
-        not computed. If `True` the integral image is computed for `image`
-        and used for finding the Hessian Determinant.
+        matrix.
 
     Returns
     -------
@@ -211,9 +207,7 @@ def hessian_matrix_det(image, sigma, integral=True):
     """
 
     image = img_as_float(image)
-    if(integral):
-        image = integral_image(image)
-
+    image = integral_image(image)
     return np.array(_hessian_matrix_det(image, sigma))
 
 
diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py
index dd9f7703..3d59a495 100644
--- a/skimage/feature/tests/test_corner.py
+++ b/skimage/feature/tests/test_corner.py
@@ -93,9 +93,12 @@ def test_hessian_matrix_eigvals():
 
 
 def test_hessian_matrix_det():
-    image = np.ones((5, 5))
-    det = hessian_matrix_det(image, 3, False)
-    assert_array_equal(det, 0)
+    image = np.zeros((5, 5))
+    image[2, 2] = 1
+    det = hessian_matrix_det(image, 5)
+    
+
+    assert_almost_equal(det, 0, decimal = 3)
 
 
 def test_square_image():

From 848a0f5fbe91fdce0bd9b57d70f1f04edd40a54a Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Thu, 10 Apr 2014 21:03:37 +0530
Subject: [PATCH 21/22] removed whitespaces

---
 skimage/feature/tests/test_corner.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/skimage/feature/tests/test_corner.py b/skimage/feature/tests/test_corner.py
index 3d59a495..7d2af0df 100644
--- a/skimage/feature/tests/test_corner.py
+++ b/skimage/feature/tests/test_corner.py
@@ -96,8 +96,6 @@ def test_hessian_matrix_det():
     image = np.zeros((5, 5))
     image[2, 2] = 1
     det = hessian_matrix_det(image, 5)
-    
-
     assert_almost_equal(det, 0, decimal = 3)
 
 

From 4d11868c76037c26a3e4fbf56627b94775da83f4 Mon Sep 17 00:00:00 2001
From: Vighnesh Birodkar <vighneshbirodkar@gmail.com>
Date: Fri, 11 Apr 2014 13:08:25 +0530
Subject: [PATCH 22/22] removed blanklines in doc strings

---
 skimage/feature/_hessian_det_appx.pyx | 3 ---
 skimage/feature/blob.py               | 2 --
 2 files changed, 5 deletions(-)

diff --git a/skimage/feature/_hessian_det_appx.pyx b/skimage/feature/_hessian_det_appx.pyx
index 7f5a403e..0c54cba7 100644
--- a/skimage/feature/_hessian_det_appx.pyx
+++ b/skimage/feature/_hessian_det_appx.pyx
@@ -25,7 +25,6 @@ cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low, Py_ssize_t high):
     -------
     x : int
         `x` clipped between `high` and `low`.
-
     """
 
     if(x > high):
@@ -60,7 +59,6 @@ cdef inline cnp.double_t _integ(
     -------
     ans : int
         The integral over the given window.
-
     """
 
     r = _clip(r, 0, img.shape[0] - 1)
@@ -130,7 +128,6 @@ def _hessian_matrix_det(cnp.double_t[:, ::1] img, double sigma):
 
     for r in range(height):
         for c in range(width):
-
             tl = _integ(img, r - s3, c - s3, s3, s3)  # top left
             br = _integ(img, r + 1, c + 1, s3, s3)  # bottom right
             bl = _integ(img, r - s3, c + 1, s3, s3)  # bottom left
diff --git a/skimage/feature/blob.py b/skimage/feature/blob.py
index 343863a9..134e7026 100644
--- a/skimage/feature/blob.py
+++ b/skimage/feature/blob.py
@@ -36,7 +36,6 @@ def _blob_overlap(blob1, blob2):
     -------
     f : float
         Fraction of overlapped area.
-
     """
     root2 = sqrt(2)
 
@@ -81,7 +80,6 @@ def _prune_blobs(blobs_array, overlap):
     -------
     A : ndarray
         `array` with overlapping blobs removed.
-
     """
 
     # iterating again might eliminate more blobs, but one iteration suffices