**Note:**- Most of the information in this chapter is based on the doctoral thesis of D. Tschumperle and his 2006 paper Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE's. This page is no substitute for reading these papers. Rather, it is intended to give an overview that will prepare the reader on the topic so that it might become easier to absorb the material presented in these papers.

Here, represents the image, while represents an (image-dependent) smoothing term to be further specified.

Which, through the use of the Euler-Lagrange equation gives us the following PDE (the heat equation):

This smoothing process is *isotropic*; it smooths with equal strength in all spatial directions. This means that edges in the image will soon become blurred. In LSSEG, we can apply this smoothing scheme by calling the nonlinear_gauss_filter_XD() function and setting the parameter `p`

equal to zero.

**Example of progressive smoothing of an image using the heat equation**

We see that it is somewhat similar to the integral in our last example, except that we have not squared the norm. This functional measures the * total variation * of the image, and the PDE we obtain from the Euler-Lagrange equation describes a process we call * total variation flow*. This process has several interesting *theoretical* properties; readers are refered to section 2.1.1 of Brox' thesis for an overview of these, and to section 2.3.2 to 2.3.4 of the same thesis for an examination of the resulting PDE. As for the previous example, the limit of the smoothing process is a constant image (and we will arrive there in finite time). However, the intermediate images are quite different. The edges in the image are well preserved for a long time, regions gradually merge with each other, and the intermediate images take on a segmentation-like appearance. The PDE system is:

In LSSEG, we can apply total variation flow on an image using the function nonlinear_gauss_filter_XD() with the parameter `p`

equal to 1.

**Example of progressive smoothing of an image using total variation flow**

with a scalar function that is 1 in the first example and in the second. Without departing from a functional, we can create customized schemes by specifying this function directly. What we want is to slow down diffusion in regions with high gradients (i.e., edges) while still allowing diffusion to take place at normal speed in "flat" regions. This idea was put forward by Perona and Malik in 1990, where they proposed the following choice of :

is here a fixed threshold that differentiates homogeneous regions from regions of contours. Tschumperle gives a good theoretical explanation of this in section 2.1.2 of his thesis, where he also shows that not only does this process lead to smoothing whose strength depend on the local gradient; it also smoothes with different strength (which can be explicitely determined) along the local image gradient and isophote.

where is a tensor field defined on the domain and with values in the space of tensors described by positive-definite 2x2-matrices. (In case of three-dimensional images, i.e., volumes, we use instead of ). We call a field of diffusion tensors. Contrary to described above, depends not on the *image gradient* itself, but rather on the * position in the image *. Moreover, as is a tensor rather than a scalar, its effect on the smoothing term also varies with the local direction of (matrix-vector multiplication ). The smoothing at position will in general be strongest along the direction specified by the eigenvector with largest eigenvalue of , and weakest in the direction specified by the other eigenvector (which, due to the symmetry of the matrix, is perpendicular on the first). See the section on definition of a smoothing geometry for more on how the tensor field can be specified. For more on divergence-based PDEs, refer to section 2.1.4 of Tschumperle's thesis.

Here, is the smoothing geometry and is the local Hessian matrix of the image. For more details on this equation, and a discussion about how it differs from the divergence- based formulation, please refer to 2.1.5 and section 2.4.1 respectively of Tschumperle's PhD-thesis. In section 3.2.1 of the same thesis, Tschumperle shows that this equation has an analytical solution which is:

Where is the original image, is the image at time of the evolution process, and is a bivariate oriented Gaussian kernel:

For a given image point , the eigenvectors of the structure tensor represent the two privileged smoothing directions. Due to the smoothing of the structure tensor, we cannot directly refer to these eigenvectors as the gradient and isophote of the image at , but they are perpendicular and represent the direction along which the image varies the most (and thus were we should smooth little) and the direction along which the image varies the least (and thus where we should smooth more). The corresponding eigenvalues give information about strong the variation is. If the two eigenvalues are identical, then the image variation around this point is more or less the same in all directions, and the smoothing geometry should be *isotropic* here. If the two eigenvalues are different, we should smooth more along the direction with the smaller eigenvalue than along the direction with the bigger one, and the difference in smoothing strength (the anisotropy) should be reflected by the difference between these two values. If we name the two directions and , and choose two corresponding smoothing strenths and , the smoothing geometry matrix representation for a given image point becomes:

In LSSEG, we have a function called compute_smoothing_geometry_XD() that computes a smoothing geometry based on an image's structure tensor. Here, the two smoothing directions and are the eigenvectors of the structure tensor, and the corresponding smoothing strengths are taken to be:

where and are the two eigenvalues of the structure tensor corresponding to and respectively. This choice is discussed briefly in section 2.1 of [Tschumperle06]. Note that as long as , then the smoothing will be stronger along direction than along direction , and this difference will be more pronounced as the biggest eigenvalue (an thus the denominator of the fractions above) grows.

**Note:**- Everything said above is directly generalizable to 3D-images (volumes), and LSSEG handles 2D and 3D smoothing in the same way.

is here the timestep used. In other words he carries out 1-dimensional convolutions of the image with * along the streamlines defined by the vector fields obtained from *.

where and is a tensor field that is the "square root" of , in the sense that . If , then .

where stands for the Jacobian of the vector field . However, the system is not solved using the "standard" method of discrete differencing on a regular grid. Instead, he demonstrates how this can be interpreted as a combination of 1D heat flows along the streamlines of the various vector fields .

- First, the image is read from file. This can be a normal 2D image in any standard format, or alternatively a
*stack*representing a 3D image. Stacks can be made with the program app/generate_stack.C - Then we enter a loop on the number of timesteps used. Note that since we are using line integral convolution to compute the changes, each timestep can be arbitrarily large without affecting the stability (the nonlinearity
*will*be affected, though). - Depending on whether the image in question is 2D or 3D, the 2D/3D structure tensor is now computed. Note that the image has been blurred a little prior to this.
- Based on the structure tensor, the smoothing geometry is now computed, by a call to compute_smoothing_geometry_2D() or compute_smoothing_geometry_3D(), depending on whether the image is two or three dimensional.
- Now we carry out the splitting of the smoothing geometry tensor field into a sum of vector fields, and for each such field we carry out the actual line integral convolution.
- Finally, the result is displayed and, in case of 3D images, the result is saved as a
*stack*to file.

**Note:**- due to problems with the functions LIC_2D() and LIC_3D(), these are substituted with the less precise functions LIC_2D_FS() and LIC_3D_FS(). Read the note in the documentation of these function for an explanation.
**image smoothed using our implementation of Tschumperle's algorithm. The parameters are the default ones, and the number of timesteps are 1, 4 and 8 (original image to the left)**

Generated on Tue Nov 28 18:35:47 2006 for lsseg by 1.4.7