Frangi filter support with tests + docstrings

+fixed tests and removed shared code

https://github.com/scikit-image/scikit-image/pull/2139#discussion_r67240635

docstrings and style fixes

fix Lambda yielding

more style fixes

enumerate & black line fixes

CONTRIBUTORS.txt

@@ -231,3 +231,6 @@
 
 - Jeff Hemmelgarn
   Minimum threshold
+
+- Kirill Malev
+  Frangi and hessian filters implementation

doc/examples/filters/plot_frangi.py

@@ -0,0 +1,29 @@
+"""
+==============
+Frangi filter
+==============
+
+Frangi and hybrid Hessian filters can be used for edge detection and
+calculation of fraction of the image containing edges.
+"""
+
+from skimage.data import camera
+from skimage.filters import frangi, hessian
+
+import matplotlib.pyplot as plt
+
+
+image = camera()
+
+fig, ax = plt.subplots(ncols=2, subplot_kw={'adjustable':'box-forced'})
+
+ax[0].imshow(frangi(image), cmap=plt.cm.gray)
+ax[0].set_title('Frangi filter results')
+
+ax[1].imshow(hessian(image), cmap=plt.cm.gray)
+ax[1].set_title('Hybrid Hessian filter result')
+
+for a in ax:
+    a.axis('off')
+
+plt.tight_layout()

skimage/filters/__init__.py

@@ -7,9 +7,11 @@ from .edges import (sobel, sobel_h, sobel_v,
                     laplace)
 from ._rank_order import rank_order
 from ._gabor import gabor_kernel, gabor
+from ._frangi import frangi, hessian
 from .thresholding import (threshold_adaptive, threshold_otsu, threshold_yen,
                            threshold_isodata, threshold_li, threshold_minimum,
-                           threshold_mean, threshold_triangle, try_all_threshold)
+                           threshold_mean, threshold_triangle,
+                           try_all_threshold)
 from . import rank
 from .rank import median
 
@@ -46,6 +48,8 @@ __all__ = ['inverse',
            'gabor_kernel',
            'gabor',
            'try_all_threshold',
+           'frangi',
+           'hessian',
            'threshold_adaptive',
            'threshold_otsu',
            'threshold_yen',

skimage/filters/_frangi.py

@@ -0,0 +1,192 @@
+import numpy as np
+
+__all__ = ['frangi', 'hessian']
+
+
+def _frangi_hessian_common_filter(image, scale, scale_step, beta1, beta2,
+                                  frangi_=True, black_ridges=True):
+    """This is an intermediate function for Frangi and Hessian filters.
+
+    Shares the common code for Frangi and Hessian functions.
+
+    Parameters
+    ----------
+    image : (N, M) ndarray
+        Array with input image data.
+    scale : tuple of floats, optional
+        The range of sigmas used.
+    scale_step : float, optional
+        Step size between sigmas.
+    beta1 : float, optional
+        Frangi correction constant.
+    beta2 : float, optional
+        Frangi correction constant.
+    black_ridges : boolean, optional
+        If True (default), detects black ridges, if False - white ones.
+
+    Returns
+    -------
+    filtered_list : list
+        List of pre-filtered images.
+
+    """
+
+    # Import has to be here due to circular import error
+    from ..feature import hessian_matrix, hessian_matrix_eigvals
+
+    sigmas = np.arange(scale[0], scale[1], scale_step)
+
+    if np.any(np.asarray(sigmas) < 0.0):
+        raise ValueError("Sigma values less than zero are not valid")
+
+    beta1 = 2 * beta1 ** 2
+    beta2 = 2 * beta2 ** 2
+
+    filtered_array = np.zeros(np.shape(image),len(sigmas))
+
+    # Filtering for all sigmas
+    for sigma in sigmas:
+        # Make 2D hessian
+        (Dxx, Dxy, Dyy) = hessian_matrix(image, sigma)
+
+        # Correct for scale
+        Dxx = (sigma ** 2) * Dxx
+        Dxy = (sigma ** 2) * Dxy
+        Dyy = (sigma ** 2) * Dyy
+
+        # Calculate (abs sorted) eigenvalues and vectors
+        (lambda1, lambda2) = hessian_matrix_eigvals(Dxx, Dxy, Dyy)
+
+        # Compute some similarity measures
+        lambda1[lambda1 == 0] = 1e-10
+        rb = (lambda2 / lambda1) ** 2
+        s2 = lambda1 ** 2 + lambda2 ** 2
+
+        # Compute the output image
+        filtered = np.exp(-rb / beta1) * (np.ones(np.shape(image)) -
+                                          np.exp(-s2 / beta2))
+
+        # Store the results in 3D matrices
+        filtered_array.append([filtered, lambda1])
+    return filtered_array
+
+
+def frangi(image, scale=(1, 10), scale_step=2, beta1=0.5, beta2=15,
+           black_ridges=True):
+    """Filter an image with the Frangi filter.
+
+    This filter can be used to detect continous edges, e.g. vessels,
+    wrinkles, rivers. It can be used to calculate the fraction of the
+    whole image containing such objects.
+
+    Calculates the eigenvectors of the Hessian to compute the likeliness of
+    an image region to vessels, according to the method described in _[1].
+
+    Parameters
+    ----------
+    image : (N, M) ndarray
+        Array with input image data.
+    scale : tuple of floats, optional
+        The range of sigmas used.
+    scale_step : float, optional
+        Step size between sigmas.
+    beta1 : float, optional
+        Frangi correction constant.
+    beta2 : float, optional
+        Frangi correction constant.
+    black_ridges : boolean, optional
+        Detect black ridges (default) set to true, for
+        white ridges set to false.
+
+    Returns
+    -------
+    out : (N, M) ndarray
+        Filtered image (maximum of pixels across all scales).
+
+    Notes
+    -----
+    Written by Marc Schrijver, 2/11/2001
+    Re-Written by D. J. Kroon University of Twente (May 2009)
+
+    References
+    ----------
+    .. [1] A. Frangi, W. Niessen, K. Vincken, and M. Viergever. "Multiscale
+    vessel enhancement filtering," In LNCS, vol. 1496, pages 130-137,
+    Germany, 1998. Springer-Verlag.
+    .. [2] Kroon, D.J.: Hessian based frangi vesselness filter.
+    .. [3] http://mplab.ucsd.edu/tutorials/gabor.pdf.
+    """
+
+    filtered_array = _frangi_hessian_common_filter(image, scale, scale_step,
+                                                   beta1, beta2)
+
+    for i in range(len(filtered_array)):
+        filtered = filtered_array[i][0]
+        Lambda1 = filtered_array[i][1]
+        if black_ridges:
+            filtered[Lambda1 < 0] = 0
+        else:
+            filtered[Lambda1 >= 0] = 0
+        filtered_array[i][0] = filtered
+
+    # Return for every pixel the value of the scale(sigma) with the maximum
+    # output pixel value
+
+    return np.max(filtered_array, axis=0)[0]
+
+
+def hessian(image, scale=(1, 10), scale_step=2, beta1=0.5, beta2=15):
+    """Filter an image with the Hessian filter.
+
+    This filter can be used to detect continous edges, e.g. vessels,
+    wrinkles, rivers. It can be used to calculate the fraction of the whole
+    image containing such objects
+
+    Almost equal to frangi filter, but uses alternative method of smoothing.
+    Address _[1] to find the differences between Frangi and Hessian filters.
+
+    Parameters
+    ----------
+    image : (N, M) ndarray
+        Array with input image data.
+    scale : tuple of floats, optional
+        The range of sigmas used.
+    scale_step : float, optional
+        Step size between sigmas.
+    beta1 : float, optional
+        Frangi correction constant.
+    beta2 : float, optional
+        Frangi correction constant.
+
+    Returns
+    -------
+    out : (N, M) ndarray
+        Filtered image (maximum of pixels across all scales).
+
+    Notes
+    -----
+    Written by Marc Schrijver, 2/11/2001
+    Re-Written by D. J. Kroon University of Twente (May 2009)
+
+    References
+    ----------
+    .. [1] Choon-Ching Ng, Moi Hoon Yap, Nicholas Costen and Baihua Li,
+    "Automatic Wrinkle Detection using Hybrid Hessian Filter".
+    """
+
+    filtered_array = _frangi_hessian_common_filter(image, scale, scale_step,
+                                                   beta1, beta2)
+
+    for i in range(len(filtered_array)):
+        filtered = filtered_array[i][0]
+        Lambda1 = filtered_array[i][1]
+        filtered[Lambda1 < 0] = 0
+        filtered_array[i][0] = filtered
+
+    # Return for every pixel the value of the scale(sigma) with the maximum
+    # output pixel value
+    out = np.max(filtered_array, axis=0)
+    out[out <= 0] = 1
+
+    return out[0]
+

skimage/filters/tests/test_frangi.py

@@ -0,0 +1,38 @@
+import numpy as np
+from numpy.testing import assert_equal, assert_almost_equal, assert_allclose
+from skimage.filters._frangi import frangi, hessian
+from skimage.data import camera
+from skimage.util import crop
+
+
+def test_null_matrix():
+    a = np.zeros((3, 3))
+    assert_almost_equal(frangi(a), np.zeros((3, 3)))
+    assert_almost_equal(frangi(a, black_ridges=False), np.zeros((3, 3)))
+    assert_equal(hessian(a), np.ones((3, 3)))
+
+
+def test_energy_decrease():
+    a = np.zeros((3, 3))
+    a[1, 1] = 1.
+    assert frangi(a).std() < a.std()
+    assert frangi(a, black_ridges=False).std() < a.std()
+    assert hessian(a).std() > a.std()
+
+
+def test_values_decreased():
+    a = np.multiply(np.ones((3, 3)), 10)
+    assert_equal(frangi(a), np.zeros((3, 3)))
+    assert_equal(hessian(a), np.ones((3, 3)))
+
+
+def test_cropped_camera_image():
+    image = crop(camera(), ((206, 206), (206, 206)))
+    assert_allclose(frangi(image), np.zeros((100, 100)), atol=1e-03)
+    assert_allclose(frangi(image, black_ridges=True), np.zeros((100,100)), atol=1e-03)
+    assert_allclose(hessian(image), np.ones((100, 100)), atol=1-1e-07)
+
+
+if __name__ == "__main__":
+    from numpy import testing
+    testing.run_module_suite()