KEYWORD 2D Normalize Method NoiseAveraging NumberOfIterations RelativeVarianceLimit RelativeNoiseVarianceLimit ContrastThreshold StructureThreshold StructureAveraging NoiseArea PixelSpacing DiscretizationInterval DiffusionProportion Regularization Balance EigenvalueThreshold FilterTruncation CommandLine Algorithm DESCRIPTION A. N. Diffusion 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 Fernandez, 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. Jose-Jesus Fernandez 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: 2D processing Use 2D processing if you want the filtering applied separately to each section rather than processing xyz volumes as a whole. Normalize This is a preprocessing step to convert the intensities to a fixed range (0 - 255). This may make it easier to choose values for the contrast threshold and structure threshold parameters. Method Fernandez and Li (2003) outline three approaches to nonlinear anisotropic diffusion filtering: edge enhancing diffusion (EED), coherence enhancing diffusion (CED), and a hybrid approach. This parameter selects which approach to use. Noise averaging Controls the amount of averaging done before character- izing local structures. An optimal setting is one that minimizes the degradation of the local structure information from the combined effects of the noise and the averaging process. Number of iterations More iterations gives a smoother result; more iterations also require more computation time. Contrast threshold and structure threshold These parameters control the amount of smoothing as a function of the strength of the local structural properties. The edge enhancing diffusion method uses the contrast threshold. The coherence enhancing method uses the structure threshold. The hybrid method uses both. Structure averaging The coherence enhancing diffusion and hybrid methods use an averaged structure tensor calculated by averaging the structure tensors over a local area. The structure averaging parameters control the size of the neighborhood over which the structure tensors are averaged. Noise area If you use the hybrid method or allow early termination based on the variance of the noise, the filtering process uses a portion of each volume (or section if processing in 2D) to estimate the properties of the noise. The noise area parameters set the region used when calculating the estimates. 2D 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. Normalize 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. Method This implementation provides three different approaches to calculating the diffusion tensor used in the non- linear anisotropic diffusion filtering; these are the three approaches described in Fernandez and Li (2003): edge enhancing Primarily enhances or preserves edges. Has fewer user parameters than the other approaches and uses less memory. coherence enhancing Intended to enhance flow-like features and complete interrupted lines. hybrid Dynamically switches between the edge enhancing and coherence enhancing approaches based on local structure properties. Has the most user parameters (all the parameters for the edge enhancing and coherence enhancing methods and some more in addition) and requires the most memory during processing. 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. NoiseAveraging 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 two 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. NumberOfIterations 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, Dr. Fernandez 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. RelativeVarianceLimit 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. RelativeNoiseVarianceLimit 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. ContrastThreshold The contrast threshold, which is the same as the parameter K in Fernandez and Li (2003), is the key parameter affecting how the dge 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. Fernandez 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. StructureThreshold The structure threshold, which is the same as the parameter C in Fernandez 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. Fernandez 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. StructureAveraging 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 two 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. NoiseArea 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. PixelSpacing 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. DiscretizationInterval 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. DiffusionProportion Fernandez 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 Fernandez 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. Regularization The regularization parameter is the alpha in Fernandez 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. Balance The balance parameter affects when the hybrid method uses the edge enhancing approach or the coherence enhancing approach. At positions where the difference between 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 a 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. EigenvalueThreshold 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. FilterTruncation 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. CommandLine The nonlinear anisotropic diffusion application accepts the command-line arguments described in Region.hlp. The -logging=debug option described there has the additional effect of 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 Turns on 2D processing; the default is to use 3D processing. -area=fx:lx:fy:ly[:fz:lz] Sets the region used for estimating the statistics of the noise. The bottom lower left corner of the region will be (fx, fy, fz) and the top upper left corner will be (lx, ly, lz). fx, fy, fz, lx, ly, and lz must be integers and are referenced to a zero-based pixel coordinate system, i.e. the lower left corner of the data set is at (0, 0, 0) and the upper left corner of the data set is at (nx-1, ny-1, nz-1) where nx, ny, nz are the dimensions, in pixels, of the data set. If the region for noise statistics does not fit within the region selected for processing, processing will terminate with an error. You can omit fz and lz: that is equivalent to setting fz equal to zero and lz equal to nine. By default, the region for noise statistics has its bottom lower left corner at (0, 0, 0) and its top upper right corner at (9, 9, 9). -balance=b Sets the balance parameter to be b. b must be a value greater than or equal to zero. The default value is one. -cthresh=k Sets the contrast threshold to be k. k must be a value greater than zero. The contrast threshold is 2.5 by default. -ethresh=e Sets the eigenvalue threshold to be e. e must be a value greater than or equal to zero. The eigenvalue threshold is one by default. -interval=d Sets the discretization interval to be d. d must be a value greater than zero. The discretization interval is 0.4 by default. -iterations=n Sets the maximum number of iterations to perform to be n. By default, at most five iterations are used. -method=m Selects which of the three possible algorithms to use for computing the diffusion tensor. m must be one of the following: ced coherence enhancing algorithm eed edge enhancing algorithm hybrid hybrid algorithm -normalize=c Enables normalization of the input data using a cutoff level of c standard deviations. c must be greater than or equal to zero. When run from the command line, normalization of the input data is off by default. -prced=p Sets the diffusion proportion parameter to be p. p must be a value between zero and one. By default, the diffusion proportion is 0.5. -regfac=r Sets the regularization factor to be r. r must be a value between zero and one. By default, the regularization factor is 0.0001. -rv=c Enables early termination before the maximum number of iterations is reached if the var- iance of the slice (if 2D) or volume (if 3D) falls below c times its original variance. c must be greater than zero and less than or equal to one. -rnv=c Enables early termination before the maximum number of iterations is reached if the variance of region used for noise statistics falls below c times its original variance. c must be greater than zero and less than or equal to one. -smooth=sx:sy[:sz] Sets the sigmas for the Gaussian filter used during computation of the structure tensors. sx is the sigma, in pixels, for the x direction, sy is the sigma, in pixels, for the y direction, and sz is the sigma, in pixels, for the z direction. sx, sy, and sz must be greater than or equal to zero. You can omit sz which is equivalent to setting sz equal to 0.625. By default, the sigmas all have values of 0.625. -spacing=sx:sy[:sz] Sets the pixel spacing to use when computing derivatives. If you omit this option, the pixel spacings recorded in the header are used. -ssmooth=sx:sy[:sz] Sets the sigmas for the Gaussian filter used for computation of the average structure tensor in the coherence enhancing and hybrid methods. sx is the sigma, in pixels, for the x direction, sy is the sigma, in pixels, for the y direction, and sz is the sigma in pixels for the z direction. sx, sy, and sz must be greater than or equal to zero. You can omit sz which is equivalent to setting sz equal to 0.625. By default, the sigmas all have values of 0.625. -sthresh=c Sets the structure threshold to be c. c must be a value greater than or equal to zero. By default, the structure threshold is 4.5. 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 Algorithm The algorithm used is based on the description in Fernandez and Li (2003) and source code provided by Dr. Fernandez. One design goal for the algorithm was to reduce the amount of memory necessary for typical cases without introducing undue amounts of extra computation. Shifting Windows 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. 1) For 2D processing with the edge enhancing algorithm the memory requirement is nx * ny + (nx + 4) * (3 * nny + 8 + min(ny, ny + 2)) + 6 * (nx + 2) 32-bit floating-point values. 2) For 3D processing with the edge enhancing algorithm the memory requirement is nx * ny * nz + (nx + 4) * (1 + ny * (3 * nnz + 8 + min(nz, nnz + 2))) + 9 * (nx + 2) * (ny + 2) 32-bit floating-point values. 3) For 2D processing with the coherence enhancing algorithm the memory requirement is nx * ny + (nx + 4) * (3 * nny + 4 * nsy + 8 + min(ny, nny + 2)) + 6 * (nx + 2) + 3 * nx * (2 * nsy + 2) 32-bit floating-point values. 4) For 3D processing with the coherence enhancing algorithm the memory requirement is nx * ny * nz + (nx + 4) * (1 + ny * (3 * nnz + 4 * nsz + 8 + min(nz, nnz + 2))) + 9 * (nx + 2) * (ny + 2) + 6 * nx * (1 + ny * (2 * nsz + 2)) 32-bit floating-point values. 5) For 2D processing with the hybrid algorithm the memory requirement is nx * ny + (nx + 4) * (3 * nny + 4 * nsy + 8 + min(ny, nny + 2)) + 6 * (nx + 2) + 3 * nx * (3 * nsy + 2) 32-bit floating-point values. 6) For 3D processing with the hybrid algorithm the memory requirement is nx * ny * nz + (nx + 4) * (1 + ny * (3 * nnz + 4 * nsz + 8 + min(nz, nnz + 2))) + 9 * (nx + 2) * (ny + 2) + 6 * nx * (1 + ny * (3 * nsz + 2)) 32-bit floating-point values. 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). Approximating Derivatives 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 Fernandez 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. Boundary Conditions 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. Diffusion Tensor The calculation of the diffusion tensor is intended to be the same as the approach in Dr. Fernandez's source code for the edge enhancing algorithm and the hybrid algorithm with the following differences: 1) The relationship between the eigenvalues for the diffusion tensor and the length of the gradient vector (or equivalently, the square root of the sum of the diagonal elements of the structure tensor) for the edge enhancing method and the edge enhancing component of the hybrid method uses the function given in Fernandez and Li (2003). 2) There is a stand-alone coherence enhancing implementation without the additional features and parameters of the hybrid method. 3) There are 2D versions as well. 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. Fernandez 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 a3, and the corresponding eigenvectors are (av11, av12, av13), (av21, av22, av23), and (av31, av32, av33). d1, d2, and d3 are the the eigenvalues of the diffusion tensor, and (dv11, dv12, dv13), (dv21, dv22, dv23), and (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 ] Edge Enhancing Diffusion (EED) 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^0.5 / K)^8) Coherence Enhancing Diffusion (CED) 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 Hybrid 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.