Merge pull request #2153 from jni/frangi

Add Frangi and Hessian filters

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,31 @@
+"""
+=============
+Frangi filter
+=============
+
+The Frangi and hybrid Hessian filters can be used to detect continuous
+edges, such as vessels, wrinkles, and rivers.
+"""
+
+from skimage.data import camera
+from skimage.filters import frangi, hessian
+
+import matplotlib.pyplot as plt
+
+image = camera()
+
+fig, ax = plt.subplots(ncols=3, subplot_kw={'adjustable': 'box-forced'})
+
+ax[0].imshow(image, cmap=plt.cm.gray)
+ax[0].set_title('Original image')
+
+ax[1].imshow(frangi(image), cmap=plt.cm.gray)
+ax[1].set_title('Frangi filter result')
+
+ax[2].imshow(hessian(image), cmap=plt.cm.gray)
+ax[2].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,179 @@
+import numpy as np
+
+__all__ = ['frangi', 'hessian']
+
+
+def _frangi_hessian_common_filter(image, scale_range, scale_step,
+                                  beta1, beta2):
+    """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_range : 2-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
+    -------
+    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_range[0], scale_range[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(sigmas.shape + image.shape)
+    lambdas_array = np.zeros(sigmas.shape + image.shape)
+
+    # Filtering for all sigmas
+    for i, sigma in enumerate(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[i] = filtered
+        lambdas_array[i] = lambda1
+    return filtered_array, lambdas_array
+
+
+def frangi(image, scale_range=(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 continuous 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 similarity of
+    an image region to vessels, according to the method described in _[1].
+
+    Parameters
+    ----------
+    image : (N, M) ndarray
+        Array with input image data.
+    scale_range : 2-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
+        When True (the default), the filter detects black ridges; when
+        False, it detects white ridges.
+
+    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, lambdas = _frangi_hessian_common_filter(image,
+                                                      scale_range, scale_step,
+                                                      beta1, beta2)
+    if black_ridges:
+        filtered[lambdas < 0] = 0
+    else:
+        filtered[lambdas > 0] = 0
+
+    # Return for every pixel the value of the scale(sigma) with the maximum
+    # output pixel value
+    return np.max(filtered, axis=0)
+
+
+def hessian(image, scale_range=(1, 10), scale_step=2, beta1=0.5, beta2=15):
+    """Filter an image with the Hessian filter.
+
+    This filter can be used to detect continuous 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.
+    Refer to _[1] to find the differences between Frangi and Hessian filters.
+
+    Parameters
+    ----------
+    image : (N, M) ndarray
+        Array with input image data.
+    scale_range : 2-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, lambdas = _frangi_hessian_common_filter(image,
+                                                      scale_range, scale_step,
+                                                      beta1, beta2)
+    filtered[lambdas < 0] = 0
+
+    # Return for every pixel the value of the scale(sigma) with the maximum
+    # output pixel value
+    out = np.max(filtered, axis=0)
+    out[out <= 0] = 1
+    return out
+

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 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()