Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
1 contributor

Users who have contributed to this file

494 lines (384 sloc) 16.6 KB
"""
Ridge filters.
Ridge filters can be used to detect continuous edges, such as vessels,
neurites, wrinkles, rivers, and other tube-like structures. The present
class of ridge filters relies on the eigenvalues of the Hessian matrix of
image intensities to detect tube-like structures where the intensity changes
perpendicular but not along the structure.
"""
from warnings import warn
import numpy as np
from ..util import img_as_float, invert
from .._shared.utils import assert_nD
def _divide_nonzero(array1, array2, cval=1e-10):
"""
Divides two arrays.
Denominator is set to small value where zero to avoid ZeroDivisionError and
return finite float array.
Parameters
----------
array1 : (N, ..., M) ndarray
Array 1 in the enumerator.
array2 : (N, ..., M) ndarray
Array 2 in the denominator.
cval : float, optional
Value used to replace zero entries in the denominator.
Returns
-------
array : (N, ..., M) ndarray
Quotient of the array division.
"""
# Copy denominator
denominator = np.copy(array2)
# Set zero entries of denominator to small value
denominator[denominator == 0] = cval
# Return quotient
return np.divide(array1, denominator)
def _sortbyabs(array, axis=0):
"""
Sort array along a given axis by absolute values.
Parameters
----------
array : (N, ..., M) ndarray
Array with input image data.
axis : int
Axis along which to sort.
Returns
-------
array : (N, ..., M) ndarray
Array sorted along a given axis by absolute values.
Notes
-----
Modified from: http://stackoverflow.com/a/11253931/4067734
"""
# Create auxiliary array for indexing
index = list(np.ix_(*[np.arange(i) for i in array.shape]))
# Get indices of abs sorted array
index[axis] = np.abs(array).argsort(axis)
# Return abs sorted array
return array[tuple(index)]
def compute_hessian_eigenvalues(image, sigma, sorting='none'):
"""
Compute Hessian eigenvalues of nD images.
For 2D images, the computation uses a more efficient, skimage-based
algorithm.
Parameters
----------
image : (N, ..., M) ndarray
Array with input image data.
sigma : float
Smoothing factor of image for detection of structures at different
(sigma) scales.
sorting : {'val', 'abs', 'none'}, optional
Sorting of eigenvalues by values ('val') or absolute values ('abs'),
or without sorting ('none'). Default is 'none'.
Returns
-------
eigenvalues : (D, N, ..., M) ndarray
Array with (sorted) eigenvalues of Hessian eigenvalues for each pixel
of the input image.
"""
# Import has to be here due to circular import error
from ..feature import hessian_matrix, hessian_matrix_eigvals
# Convert image to float
image = img_as_float(image)
# Make nD hessian
hessian_elements = hessian_matrix(image, sigma=sigma, order='rc')
# Correct for scale
hessian_elements = [(sigma ** 2) * e for e in hessian_elements]
# Compute Hessian eigenvalues
hessian_eigenvalues = np.array(hessian_matrix_eigvals(hessian_elements))
if sorting == 'abs':
# Sort eigenvalues by absolute values in ascending order
hessian_eigenvalues = _sortbyabs(hessian_eigenvalues, axis=0)
elif sorting == 'val':
# Sort eigenvalues by values in ascending order
hessian_eigenvalues = np.sort(hessian_eigenvalues, axis=0)
# Return Hessian eigenvalues
return hessian_eigenvalues
def meijering(image, sigmas=range(1, 10, 2), alpha=None,
black_ridges=True):
"""
Filter an image with the Meijering neuriteness filter.
This filter can be used to detect continuous ridges, e.g. neurites,
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 neurites, according to the method described in [1]_.
Parameters
----------
image : (N, M[, ..., P]) ndarray
Array with input image data.
sigmas : iterable of floats, optional
Sigmas used as scales of filter
alpha : float, optional
Frangi correction constant that adjusts the filter's
sensitivity to deviation from a plate-like structure.
black_ridges : boolean, optional
When True (the default), the filter detects black ridges; when
False, it detects white ridges.
Returns
-------
out : (N, M[, ..., P]) ndarray
Filtered image (maximum of pixels across all scales).
References
----------
.. [1] Meijering, E., Jacob, M., Sarria, J. C., Steiner, P., Hirling, H.,
Unser, M. (2004). Design and validation of a tool for neurite tracing
and analysis in fluorescence microscopy images. Cytometry Part A,
58(2), 167-176.
:DOI:`10.1002/cyto.a.20022`
"""
# Check (sigma) scales
sigmas = np.asarray(sigmas)
if np.any(sigmas < 0.0):
raise ValueError('Sigma values less than zero are not valid')
# Get image dimensions
ndim = image.ndim
# Set parameters
if alpha is None:
alpha = 1.0 / ndim
# Invert image to detect bright ridges on dark background
if not black_ridges:
image = invert(image)
# Generate empty (n+1)D arrays for storing auxiliary images filtered at
# different (sigma) scales
filtered_array = np.zeros(sigmas.shape + image.shape)
# Filtering for all (sigma) scales
for i, sigma in enumerate(sigmas):
# Calculate (sorted) eigenvalues
eigenvalues = compute_hessian_eigenvalues(image, sigma, sorting='val')
if ndim > 1:
# Set coefficients for scaling eigenvalues
coefficients = [alpha] * ndim
coefficients[0] = 1
# Compute auxiliary variables l_i = e_i + sum_{j!=i} alpha * e_j
auxiliary = [np.sum([eigenvalues[i] * np.roll(coefficients, j)[i]
for j in range(ndim)], axis=0) for i in range(ndim)]
# Compute maximum over auxiliary variables
auxiliary = np.max(auxiliary, axis=0)
# Rescale image intensity
filtered = np.abs(auxiliary) / np.abs(np.max(auxiliary))
# Remove background
filtered = np.where(filtered > 0, filtered, 0)
# Store results in (n+1)D matrices
filtered_array[i] = filtered
# Return for every pixel the maximum value over all (sigma) scales
return np.max(filtered_array, axis=0)
def sato(image, sigmas=range(1, 10, 2), black_ridges=True):
"""
Filter an image with the Sato tubeness filter.
This filter can be used to detect continuous ridges, e.g. tubes,
wrinkles, rivers. It can be used to calculate the fraction of the
whole image containing such objects.
Defined only for 2-D and 3-D images. Calculates the eigenvectors of the
Hessian to compute the similarity of an image region to tubes, according to
the method described in [1]_.
Parameters
----------
image : (N, M[, P]) ndarray
Array with input image data.
sigmas : iterable of floats, optional
Sigmas used as scales of filter.
black_ridges : boolean, optional
When True (the default), the filter detects black ridges; when
False, it detects white ridges.
Returns
-------
out : (N, M[, P]) ndarray
Filtered image (maximum of pixels across all scales).
References
----------
.. [1] Sato, Y., Nakajima, S., Shiraga, N., Atsumi, H., Yoshida, S.,
Koller, T., ..., Kikinis, R. (1998). Three-dimensional multi-scale line
filter for segmentation and visualization of curvilinear structures in
medical images. Medical image analysis, 2(2), 143-168.
:DOI:`10.1016/S1361-8415(98)80009-1`
"""
# Check image dimensions
assert_nD(image, [2, 3])
# Check (sigma) scales
sigmas = np.asarray(sigmas)
if np.any(sigmas < 0.0):
raise ValueError('Sigma values less than zero are not valid')
# Invert image to detect bright ridges on dark background
if not black_ridges:
image = invert(image)
# Generate empty (n+1)D arrays for storing auxiliary images filtered
# at different (sigma) scales
filtered_array = np.zeros(sigmas.shape + image.shape)
# Filtering for all (sigma) scales
for i, sigma in enumerate(sigmas):
# Calculate (sorted) eigenvalues
lamba1, *lambdas = compute_hessian_eigenvalues(image, sigma,
sorting='val')
# Compute tubeness, see equation (9) in reference [1]_.
# np.abs(lambda2) in 2D, np.sqrt(np.abs(lambda2 * lambda3)) in 3D
filtered = np.abs(np.multiply.reduce(lambdas)) ** (1/len(lambdas))
# Remove background and store results in (n+1)D matrices
filtered_array[i] = np.where(lambdas[-1] > 0, filtered, 0)
# Return for every pixel the maximum value over all (sigma) scales
return np.max(filtered_array, axis=0)
def frangi(image, sigmas=range(1, 10, 2), scale_range=None, scale_step=None,
beta1=None, beta2=None, alpha=0.5, beta=0.5, gamma=15,
black_ridges=True):
"""
Filter an image with the Frangi vesselness filter.
This filter can be used to detect continuous ridges, e.g. vessels,
wrinkles, rivers. It can be used to calculate the fraction of the
whole image containing such objects.
Defined only for 2-D and 3-D images. 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[, P]) ndarray
Array with input image data.
sigmas : iterable of floats, optional
Sigmas used as scales of filter, i.e.,
np.arange(scale_range[0], scale_range[1], scale_step)
scale_range : 2-tuple of floats, optional
The range of sigmas used.
scale_step : float, optional
Step size between sigmas.
alpha : float, optional
Frangi correction constant that adjusts the filter's
sensitivity to deviation from a plate-like structure.
beta = beta1 : float, optional
Frangi correction constant that adjusts the filter's
sensitivity to deviation from a blob-like structure.
gamma = beta2 : float, optional
Frangi correction constant that adjusts the filter's
sensitivity to areas of high variance/texture/structure.
black_ridges : boolean, optional
When True (the default), the filter detects black ridges; when
False, it detects white ridges.
Returns
-------
out : (N, M[, P]) ndarray
Filtered image (maximum of pixels across all scales).
Notes
-----
Written by Marc Schrijver, November 2001
Re-Written by D. J. Kroon, University of Twente, May 2009, [2]_
Adoption of 3D version from D. G. Ellis, Januar 20017, [3]_
References
----------
.. [1] Frangi, A. F., Niessen, W. J., Vincken, K. L., & Viergever, M. A.
(1998,). Multiscale vessel enhancement filtering. In International
Conference on Medical Image Computing and Computer-Assisted
Intervention (pp. 130-137). Springer Berlin Heidelberg.
:DOI:`10.1007/BFb0056195`
.. [2] Kroon, D. J.: Hessian based Frangi vesselness filter.
.. [3] Ellis, D. G.: https://github.com/ellisdg/frangi3d/tree/master/frangi
"""
# Check deprecated keyword parameters
if beta1:
warn('Use keyword parameter `beta` instead of `beta1` which '
'will be removed in version 0.17.')
beta = beta1
if beta2:
warn('Use keyword parameter `gamma` instead of `beta2` which '
'will be removed in version 0.17.')
gamma = beta2
if scale_range and scale_step:
warn('Use keyword parameter `sigmas` instead of `scale_range` and '
'`scale_range` which will be removed in version 0.17.')
sigmas = np.arange(scale_range[0], scale_range[1], scale_step)
# Check image dimensions
assert_nD(image, [2, 3])
# Check (sigma) scales
sigmas = np.asarray(sigmas)
if np.any(sigmas < 0.0):
raise ValueError('Sigma values less than zero are not valid')
# Rescale filter parameters
alpha_sq = 2 * alpha ** 2
beta_sq = 2 * beta ** 2
gamma_sq = 2 * gamma ** 2
# Get image dimensions
ndim = image.ndim
# Invert image to detect dark ridges on light background
if black_ridges:
image = invert(image)
# Generate empty (n+1)D arrays for storing auxiliary images filtered
# at different (sigma) scales
filtered_array = np.zeros(sigmas.shape + image.shape)
lambdas_array = np.zeros(sigmas.shape + image.shape)
# Filtering for all (sigma) scales
for i, sigma in enumerate(sigmas):
# Calculate (abs sorted) eigenvalues
lambda1, *lambdas = compute_hessian_eigenvalues(image, sigma,
sorting='abs')
# Compute sensitivity to deviation from a plate-like structure
# see equations (11) and (15) in reference [1]_
r_a = np.inf if ndim == 2 else _divide_nonzero(*lambdas) ** 2
# Compute sensitivity to deviation from a blob-like structure,
# see equations (10) and (15) in reference [1]_,
# np.abs(lambda2) in 2D, np.sqrt(np.abs(lambda2 * lambda3)) in 3D
filtered_raw = np.abs(np.multiply.reduce(lambdas)) ** (1/len(lambdas))
r_b = _divide_nonzero(lambda1, filtered_raw) ** 2
# Compute sensitivity to areas of high variance/texture/structure,
# see equation (12)in reference [1]_
r_g = sum([lambda1 ** 2] + [lambdai ** 2 for lambdai in lambdas])
# Compute output image for given (sigma) scale and store results in
# (n+1)D matrices, see equations (13) and (15) in reference [1]_
filtered_array[i] = ((1 - np.exp(-r_a / alpha_sq))
* np.exp(-r_b / beta_sq)
* (1 - np.exp(-r_g / gamma_sq)))
lambdas_array[i] = np.max(lambdas, axis=0)
# Remove background
filtered_array[lambdas_array > 0] = 0
# Return for every pixel the maximum value over all (sigma) scales
return np.max(filtered_array, axis=0)
def hessian(image, sigmas=range(1, 10, 2), scale_range=None, scale_step=None,
beta1=None, beta2=None, alpha=0.5, beta=0.5, gamma=15,
black_ridges=True):
"""Filter an image with the Hybrid 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.
Defined only for 2-D and 3-D images. 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[, P]) ndarray
Array with input image data.
sigmas : iterable of floats, optional
Sigmas used as scales of filter, i.e.,
np.arange(scale_range[0], scale_range[1], scale_step)
scale_range : 2-tuple of floats, optional
The range of sigmas used.
scale_step : float, optional
Step size between sigmas.
beta = beta1 : float, optional
Frangi correction constant that adjusts the filter's
sensitivity to deviation from a blob-like structure.
gamma = beta2 : float, optional
Frangi correction constant that adjusts the filter's
sensitivity to areas of high variance/texture/structure.
black_ridges : boolean, optional
When True (the default), the filter detects black ridges; when
False, it detects white ridges.
Returns
-------
out : (N, M[, P]) ndarray
Filtered image (maximum of pixels across all scales).
Notes
-----
Written by Marc Schrijver (November 2001)
Re-Written by D. J. Kroon University of Twente (May 2009) [2]_
References
----------
.. [1] Ng, C. C., Yap, M. H., Costen, N., & Li, B. (2014,). Automatic
wrinkle detection using hybrid Hessian filter. In Asian Conference on
Computer Vision (pp. 609-622). Springer International Publishing.
:DOI:`10.1007/978-3-319-16811-1_40`
.. [2] Kroon, D. J.: Hessian based Frangi vesselness filter.
"""
filtered = frangi(image, sigmas=sigmas, scale_range=scale_range,
scale_step=scale_step, beta1=beta1, beta2=beta2,
alpha=alpha, beta=beta, gamma=gamma,
black_ridges=black_ridges)
filtered[filtered <= 0] = 1
return filtered
You can’t perform that action at this time.