diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index ebed896a..55451717 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -231,3 +231,6 @@ - Jeff Hemmelgarn Minimum threshold + +- Kirill Malev + Frangi and hessian filters implementation diff --git a/doc/examples/filters/plot_frangi.py b/doc/examples/filters/plot_frangi.py new file mode 100644 index 00000000..7431dcf9 --- /dev/null +++ b/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() diff --git a/skimage/filters/__init__.py b/skimage/filters/__init__.py index 61fc850c..35b321b5 100644 --- a/skimage/filters/__init__.py +++ b/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', diff --git a/skimage/filters/_frangi.py b/skimage/filters/_frangi.py new file mode 100644 index 00000000..4b95d0fb --- /dev/null +++ b/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] + diff --git a/skimage/filters/tests/test_frangi.py b/skimage/filters/tests/test_frangi.py new file mode 100644 index 00000000..9b1722f2 --- /dev/null +++ b/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()