Nonlinear anisotropic diffusion filtering is a technique for smoothing data. Unlike the smoothing filters (with the exception of the bilateral filter) in Filter 2D and Filter 3D, nonlinear anisotropic diffusion filtering adjusts the filter for local variations in structure so that important features are better preserved during the smoothing process. This implementation of nonlinear anisotropic diffusion filtering derives from the algorithms described in Fernández, JJ. and Li, S. (2003), "An improved algorithm for anisotropic nonlinear diffusion for denoising cryo-tomograms", Journal of Structural Biology, 144: 152-161 and from source code provided by Dr. José-Jesús Fernández for a filter implementing the edge enhancing algorithm and another implementing a hybrid algorithm. The Algorithm section describes the details of the algorithm.
There are a large number of parameters that you can modify to control the operation of the filter. Important ones are:
Overview | Region processing | 2D processing | Normalize | Method | Noise averaging | # of iterations | Variance limit | Noise variance limit | Contrast threshold | Structure threshold | Structure averaging | Noise area | Pixel spacing | Discretization interval | Diffusion proportion | Regularization | Balance | Eigenvalue threshold | Filter truncation | Command line | Algorithm
Filter 2D | Filter 3D | Priism
This implementation of nonlinear anisotropic diffusion filtering can process either a series of 3D volumes (xyz) where each volume is processed independently or a series of 2D slices (xy) where each slice is processed independently. Beware that 3D processing can require a large amount of memory: enough to hold the selected volume when stored as 32-bit floating-point values plus some extra space for calculations (see the Algorithm section for more precise statements about the amount of memory needed). From the graphical user interface, turn on the "2D processing" toggle button to process each 2D slice independently; turn that button off when you want to process each 3D volume independently. By default, processing is done on volumes unless the data has a z dimension size of one. Every time you select a new input file, the setting for whether or not to perform 2D processing is reset to the default value.
The filtering process includes an optional preliminary step that linearly scales the intensity values within n standard deviations of the mean so that they fall within the range of 0 to 255. Intensities less than the mean intensity minus n standard deviations are set to zero, and intensities greater than the mean intensity plus n standard deviations are set to 255. When 2D processing is done, the normalization step treats each 2D slice independently; when 3D processing is done, the normalization step treats each 3D volume independently.
One reason to normalize the data is to make it easier to compare and set the contrast threshold and structure threshold parameters when working with different data sets. However, if your data consists of a series of sections (for 2D processing) or volumes (for 3D processing) that you want to have on a consistent intensity scale, you should not use the normalization step.
From the graphical user interface, turn on the "Normalize" toggle button to enable the normalization step or turn that toggle button off to disable the normalization step. To change the value of n, the number of standard deviations from the mean to retain without saturation, adjust the value in the field immediately to the right of the "Normalize" toggle button. By default, the data will be normalized, and values more than four standard deviations from the mean saturate.
This implementation provides three different approaches to calculating the diffusion tensor used in the nonlinear anisotropic diffusion filtering; these are the three approaches described in Fernández and Li (2003):
For more information about the calculation of the diffusion tensor with the different methods, see the Algorithm section. In the graphical user interface, use the "Method" menu to select which algorithm to use. The hybrid algorithm is the default.
Nonlinear anisotropic diffusion filtering characterizes structures by calculating structure tensors from the local intensity gradients. Because calculating the gradients involves calculating differences, those calculations are sensitive to noise present in the data. For that reason, there is an option to apply a Gaussian smoothing filter to the intensities before computing the structure tensor. You control the Gaussian filter by specifying the sigmas, in units of pixels, for each of the three coordinate dimensions (only the first two are relevant when you do 2D processing). In a given dimension, the filter will extend n pixels to either side of a point where n is the smallest integer larger than or equal to the filter truncation value times the sigma. A sigma of zero for a particular direction results in no smoothing along that direction. Appropriate values for the sigmas will depend on your data: larger sigmas should help when the amount of noise relative to signal increases but unduly large sigmas will suppress the sensitivity to local variations in structure. Also, larger sigmas incur additional computation time for computations in the smoothing step and more memory needed for auxiliary storage.
In the graphical user interface, you can set the sigmas by changing the values in the "Noise averaging" field. The first value is the sigma for the x direction, the second is the sigma for the y direction, and the third is the sigma for the z direction. By default, all three have values of 0.625.
The number of iterations of filtering controls how smooth the result is: as the number of iterations increases the result becomes smoother. The amount of computation time also increases linearly with the number of iterations. For the default value of the discretization interval (0.4), Dr. Fernández recommended five to ten iterations.
In the graphical user interface, edit the value in the "# of Iterations" field to change the number of iterations performed. By default, five iterations will be performed.
In addition to the limit on the number of iterations, you can specify that processing terminate early if the variance of the block (a slice in 2D processing or a volume in 3D processing) processed falls below a given fraction of its original variance.
In the graphical user interface, you can turn on or off early termination based on the relative variance by turning on or off the "Variance limit" toggle button. Immediately to the right of that button is a field where you can set the desired relative variance level as a fraction (zero to one). By default, early termination is not enabled, and the threshold relative variance level is 0.4.
In addition to the limit on the number of iterations, you can specify that processing terminate early if the variance of the area used to estimate the noise statistics falls below a given fraction of its original variance.
In the graphical user interface, you can turn on or off early termination based on the relative noise variance by turning on or off the "Noise variance limit" toggle button. Immediately to the right of that button is a field where you can set the desired relative noise variance level as a fraction (zero to one). By default, early termination is not enabled, and the threshold relative noise variance level is 0.1.
The contrast threshold, which is the same as the parameter K in Fernández and Li (2003), is the key parameter affecting how the edge enhancing approach and the edge enhancing component of the hybrid method interpret local structures. Areas where the length of the gradient vector is greater than the contrast threshold are regarded as edges; areas where the length of the gradient vector is less than the contrast threshold are interior regions. For normalized data, Dr. Fernández found values of the contrast threshold between one and four to be useful.
In the graphical user interface, change the value in the "Contrast threshold" field to set the contrast threshold. The default value is 2.5.
The structure threshold, which is the same as the parameter C in Fernández and Li (2003), is an important parameter affecting how the coherence enhancing approach and the coherence enhancing component of the hybrid method interpret local structures. When the difference squared between the largest and smallest eigenvalues of the averaged structure tensor is greater than the structure threshold, the structure is treated as a line-like or plane-like pattern with smoothing along the line or plane. For normalized data, Dr. Fernández found values of the structure threshold between one and eight to be useful.
In the graphical user interface, change the value in the "Structure threshold" field to set the structure threshold. The default value is 4.5.
The coherence enhancing and hybrid methods use structure tensors that have been averaged over a local area. The structure averaging parameters are the sigmas for the Gaussian weighting function which defines the neighborhood for averaging. In a given direction, the neighborhood will extend n pixels to each side where n is the smallest integer larger than or equal to the filter truncation value times sigma in that direction.
In the graphical user interface, the sigmas for the weighting function are displayed in the field labeled "Structure averaging". The first value is the sigma in the x direction, the second value is the sigma in the y direction, and the third value, only used for 3D processing, is the sigma in the z direction. By default, each of the sigmas has a value of .625.
For the hybrid method or when you specify a termination criteria based on the relative noise variance, a portion of each slice or volume processed is used for estimating the statistics of the noise. The region is box shaped and specified using a bottom lower left corner and a size in each dimension. In the graphical user interface, use the field labeled "Noise area origin" to specify the bottom lower left corner of the region: the first value is the x coordinate, the second is the y coordinate, and the third, ignored in 2D processing, is the z coordinate. These coordinates must be integers and are expected in zero-based pixel units relative to the full data set size. Use the field labeled "Noise area size" to specify the dimensions of the region: the first value is the x size, the second is the y size, and the third, ignored in 2D processing, is the z size. By default, the 10 x 10 x 10 volume in the bottom lower left corner of each volume (or, for 2D processing, the 10 x 10 box in the lower left corner of each slice) is used for estimating the noise statistics. If the area for the noise statistics does not fall within the bounds of the region selected for processing, an error message will be displayed when you start processing the data and the processing will stop.
The computation of derivatives for the anisotropic nonlinear diffusion filter takes into account the sampling of the data. By default, the filter uses the sampling as recorded in the header of the input data. If you want to use a different sampling, change the values in the "Pixel sampling" field. The first value in the field is the x sampling, the second is the y sampling, and the third is the z sampling. The field is in the "other parameters" dialog and not in the main dialog; press the "Other parameters..." button in the main dialog to open the "other parameters" dialog. Note that only the relative values for the sampling matter: the filter will divide the samplings by the minimum sampling value (if you are using 2D processing, the minimum of the x and y sampling values) when determining the scale factors for the derivatives.
The discretization interval is the time step corresponding to each iteration. When adjusting the input parameters to optimize the result with your data, you should leave this parameter as is. Its intended use is for experiments to evaluate how large of a discretization interval still allows for a stable discretization of the partial differential equation. From the graphical user interface, use the field labeled "Discretization interval" to set the discretization interval. That field is in the "other parameters" dialog and not in the main dialog; press the "Other parameters..." button in the main dialog to open the "other parameters" dialog. The default value for the discretization interval is 0.4.
Fernández and Li (2003) modified the coherence enhancing diffusion method to treat line-like and plane-like features differently. For plane-like features, the diffusion proportion parameter controls the amount of diffusion along the second in-plane direction and is only relevant when processing in 3D. Allowable values for the diffusion proportion parameter are in the range of zero to one where one is the value used by Fernández and Li (2003) and a value of zero corresponds to a traditional coherence enhancing method that treats all features as line-like. From the graphical user interface, use the field labeled "Diffusion proportion" to set the diffusion proportion. That field is in the "other parameters" dialog and not in the main dialog; press the "Other parameters..." button in the main dialog to open the "other parameters" dialog. The default value for the diffusion proportion is 0.5.
The regularization parameter is the alpha in Fernández and Li (2003) and puts a lower bound on the eigenvalues of the diffusion tensor for the coherence enhancing method and for the coherence enhancing component of the hybrid method. Allowable values for the regularization parameter range from zero to one. From the graphical user interface, use the field labeled "Regularization" to set the regularization parameter. That field is in the "other parameters" dialog and not in the main dialog; press the "Other parameters..." button in the main dialog to open the "other parameters" dialog. The default value for the regularization parameter is 0.0001.
The balance parameter affects when the hybrid method uses the edge enhancing approach or the coherence enhancing approach. At positions where the difference between the largest and smallest eigenvalues of the structure tensor is greater than the balance parameter times the largest difference between the largest and smallest eigenvalues for the structure tensors of the noise area, the hybrid method will use the coherence enhancing approach; otherwise, the hybrid method will use the edge enhancing approach. The balance parameter must be greater than or equal to zero: with a value of zero the hybrid method will always use the coherence enhancing approach. From the graphical user interface, use the field labeled "Balance" to set the balance parameter. That field is in the "other parameters" dialog and not in the main dialog; press the "Other parameters..." button in the main dialog to open the "other parameters" dialog. The default value for the balance parameter is one.
The eigenvalue threshold is only relevant for the hybrid method. If the largest eigenvalue for the structure tensor at a point is less than the eigenvalue threshold times the average largest eigenvalue for the structure tensors in the noise area, then the hybrid method will use an isotropic diffusion tensor at that point rather than the coherence enhancing or edge enhancing approaches.
The eigenvalue threshold must be greater than or equal to zero. From the graphical user interface, use the field labeled "Eigenvalue threshold" to set the eigenvalue threshold. That field is in the "other parameters" dialog and not in the main dialog; press the "Other parameters..." button in the main dialog to open the "other parameters" dialog. The default value for the eigenvalue threshold is one.
Gaussian filters are used to smooth the data before computing local structure properties and, in the coherence enhancing or hybrid methods, when averaging local structure tensors. When evaluating those filters, the spatial extent is limited to a fixed number of standard deviations. To set that limit from the graphical user interface, use the field labeled "Filter truncation". That field is in the "other parameters" dialog and not in the main dialog; press the "Other parameters..." button in the main dialog to open the "other parameters" dialog. By default, the extent of the filters is limited to two standard deviations. Because of the limits of 32-bit floating-point arithmetic, there is no value to setting the limit larger than thirteen standard deviations. Larger values for the limit increase the amount of time necessary to perform the processing and also increase the amount of memory used.
The nonlinear anisotropic diffusion application accepts the command-line
arguments described in Region.html.
The -logging=debug option described there has the additional
effect causing each iteration's statistics to be written to standard output.
In addition, the nonlinear anisotropic diffusion application has the following
options (brackets are used to indicate parts that may be omitted):
-2D
-area=fx:lx:fy:ly[:fz:lz]
-balance=b
-cthresh=k
-ethresh=e
-ftrunc=c
-interval=d
-iterations=n
-method=m
-normalize=c
-prced=p
-regfac=r
-rv=c
-rnv=c
-smooth=sx:sy[:sz]
-spacing=sx:sy[:sz]
-ssmooth=sx:sy[:sz]
-sthresh=c
The following example applies ten iterations of the edge enhancing method to in.dat to generate out.dat; the contrast threshold is set at three.
ANDiffusion in.dat out.dat -method=eed -cthresh=3
The algorithm used is based on the description in Fernández and Li (2003) and source code provided by Dr. Fernández. One design goal for the algorithm was to reduce the amount of memory necessary for typical cases without introducing undue amounts of extra computation.
To reduce the amount of memory used, the algorithm employs shifting windows. One holds the image data from the previous iteration. Another holds smoothed image data for calculating structure tensors. A third holds the flux vector, the product of the diffusion tensor and the gradient vector. For the coherence enhancing algorithm and hybrid algorithm, there is a fourth shifting window for the averaged structure tensors. The hybrid algorithm has a fifth shifting window for the structure tensors. The total memory requirements for those shifting windows, the volume or slice to process, and working space for the smoothing filters are listed below. nx, ny, and nz are the dimensions of the volume to be processed. nny is the smallest integer greater than or equal to two times the y sigma for the noise averaging Gaussian filter. nnz is the smallest integer greater than or equal to two times the z sigma for the noise averaging Gaussian filter. nsy is the smallest integer greater than or equal to two times the y sigma for the Gaussian filter used to compute the average structure tensor. nsz is the smallest integer greater than or equal to two times the z sigma for the Gaussian filter used to compute the average structure tensor. min(a,b) returns a if a is less than b and otherwise returns b.
As implemented, the shifting windows introduce at least one instance of repeated calculations: the smoothed data and structure tensors for the region used to estimate the noise statistics are calculated twice every iteration (once for the statistics and a second time in the computations for the diffusion tensor).
To approximate the derivatives, the algorithm uses the approach described in Weickert, J. and Scharr, H. (2002), "A Scheme for Coherence-Enhancing Diffusion Filtering with Optimized Rotation Invariance", Journal of Visual Communication and Image Representation, 13: 103-118 and Fernández and Li (2003) where a central difference kernel, (-0.5 / d, 0, 0.5 / d) in the direction of the derivative, is convolved with smoothing kernels, (0.174654, 0.650692, 0.174654), in the other directions. d is the sampling interval in the direction of the derivative divided by the smallest of the sampling intervals from all the coordinate directions (x and y only for 2D processing; x, y, and z for 3D processing). The right hand side of the diffusion equation, the divergence of the product of the diffusion tensor and the gradient vector, is approximated, as described in Weickert and Scharr (2002), by two applications of the derivative filter: one to compute the gradient and another to compute the divergence of the flux vector.
Reflecting boundary conditions are imposed on the image data. As a consequence of that and the way the structure tensor is defined, the structure tensors and averaged structure tensors have the following boundary condition: the components of a structure tensor at some position are equal to the components of the structure tensor at that position reflected across the boundary except that the off-diagonal components coming from derivatives perpendicular to the boundary are negated. For similar reasons, the flux vectors (the product of the diffusion tensor and gradient vector) have the following boundary condition: the components of a flux vector at some position are equal to the components of the flux vector at that position reflected across the boundary except that the component of the vector perpendicular to the boundary is negated.
The calculation of the diffusion tensor is intended to be the same as the approach in Dr. Fernández's source code for the edge enhancing algorithm and the hybrid algorithm with the following differences:
In all that follows, e1, e2, and e3 are the eigenvalues of the structure tensor where e1 is greater than or equal to e2 and e2 is greater than or equal to e3; the corresponding eigenvectors are (ev11, ev12, ev13), (ev21, ev22, ev23), and (ev31, ev32, ev33). As Dr. Fernández noted to us in a personal communication, the definition of the structure tensor implies that e2 and e3 are zero and the corresponding eigenvectors are perpendicular to the gradient. We take advantage of that property below. a1, a2, and a3 are the eigenvalues of the averaged structure tensor where a1 is greater than or equal to a2 and a2 is greater than or equal to a3, and the corresponding eigenvectors are (av11, av12, av13), (av21, av22, av23), and (av31, av32, av33). d1, d2, and d3 are the eigenvalues of the diffusion tensor, and (dv11, dv12, dv13), (dv21, dv22, dv23), (dv31, dv32, dv33) are the corresponding eigenvectors. In the three dimensional case, the diffusion tensor is then the matrix product
[ dv11 dv21 dv32 ] [ d1 0 0 ] [ dv11 dv12 dv13 ]
[ dv12 dv22 dv32 ] * [ 0 d2 0 ] * [ dv21 dv22 dv23 ]
[ dv13 dv23 dv33 ] [ 0 0 d3 ] [ dv31 dv32 dv33 ]
For the two dimensional case, the diffusion tensor is the matrix product
[ dv11 dv21 ] [ d1 0 ] [ dv11 dv12 ]
[ ] * [ ] * [ ]
[ dv12 dv22 ] [ 0 d2 ] [ dv21 dv22 ]
For the edge enhancing algorithm in two dimensions, the eigenvalues and eigenvectors of the diffusion tensor are
{ g e1 > 0
d1 = {
{ 1 otherwise
dv11 = ev11; dv12 = ev12
d2 = 1
dv21 = ev21; dv22 = ev22
where g is
g = 1 - e^(-3.31488 / (e1^0.5 / K)^8)
and K is the contrast threshold parameter. For the edge enhancing algorithm in three dimensions, the eigenvalues and eigenvectors of the diffusion tensor are
{ g e1 > 0
d1 = {
{ 1 otherwise
dv11 = ev11; dv12 = ev12; dv13 = ev13
d2 = 1
dv21 = ev21; dv22 = ev22; dv23 = ev23
d3 = 1
dv31 = ev31; dv32 = ev32; dv33 = ev33
where g is
g = 1 - e^(-3.31488 / ((e1 + e2 + e3)^0.5 / K)^8)
For the coherence enhancing algorithm in two dimensions, the eigenvalues and eigenvectors of the diffusion tensor are
d1 = a
dv11 = av11; dv12 = av12
{ a + (1 - a) * e^(-C / (a1 - a2)^2) a1 > a2
d2 = {
{ a otherwise
dv21 = av21; dv22 = av22
where a is the regularization parameter and C is the structure threshold parameter.
For the coherence enhancing algorithm in three dimensions, the eigenvalues and eigenvectors of the diffusion tensor are
d1 = a
dv11 = av11; av12 = av12; av13 = av13
{ a + (1 - a) * e^(-C / (a1 - a2)^2) a1 != 0 and
{ (a1 - a2) / a1 > (a2 - a3) / a1 and
d2 = { (a1 - a2) / a1 > a3 / a1
{
{ a otherwise
dv21 = av21; dv22 = av22; dv23 = av23
{ a + (1 - a) * e^(-C / (a1 - a3)^2) a1 != 0 and
{ (((a1 - a2) / a1 > (a2 - a3) / a1 and
{ (a1 - a2) / a1 > a3 / a1)) or
{ (a2 - a3) / a1 > (a1 - a2) / a1 and
d3 = { (a2 - a3) / a1 > a3 / a1)
{
{
{ a otherwise
dv21 = av31; dv22 = av32; dv23 = av33
The hybrid algorithm uses an isotropic diffusion tensor when e1 is less than the eigenvalue threshold times the mean largest eigenvalue for the structure tensors in the region used to estimate the noise statistics. In two dimensions, the eigenvalues and eigenvectors for the isotropic diffusion tensor are
d1 = 1
dv11 = ev11; dv12 = ev12
d2 = 1
dv21 = ev21; dv22 = ev22
In three dimensions, the eigenvalues and eigenvectors for the isotropic diffusion tensor are
d1 = 1
dv11 = ev11; dv12 = ev12; dv13 = ev13
d2 = 1
dv21 = ev21; dv22 = ev22; dv23 = ev23
d3 = 1
dv31 = ev31; dv32 = ev32; dv33 = ev33
When e1 is greater than or equal to the eigenvalue threshold times the mean largest eigenvalue for the structure tensors in the region used to estimate the noise statistics, the hybrid algorithm adopts a CED or EED approach. When the largest eigenvalue of the structure tensor minus the smallest eigenvalue of the structure tensor (e1 minus e2 in two dimensions; e1 minus e3 in three dimensions) is greater than the balance parameter times the largest difference between the largest and smallest eigenvalues for structure tensors in the region used to estimate the noise statistics, the hybrid algorithm uses a CED approach; otherwise the hybrid algorithm uses an EED approach. The CED and EED approaches for calculating the diffusion tensor are the same as described above for the stand-alone implementations.