How Anisotropic Smoothing in LSSEG Works

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.

Using partial differential equations for smoothing digital images

The reasons why we would want to smooth a digital image might be to remove noise, simplify the image or remove local variations that might become obstacles for further processing (like segmentation). When using a partial differential equation (PDE) to smooth an image, we consider the image a discretized function over its domain, and apply a time-dependent PDE that makes the image evolve timestep by timestep towards progressively simplified and smooth versions of itself. The limit of such a process is generally a image with a constant value everywhere; not a very interesting solution, but we usually stop the process long before this happens. All PDE's that we discuss here are on the form:

$\frac{\partial I}{\partial t} = R$

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

How do we obtain the equations

The simplest smoothing schemes can be derived directly from a functional using the Euler-Lagrange equation. The functional in question gives a value that in some sense represents the amount of variation in the image, and the process of minimizing this functional leads to a gradual transformation of the original image into simplified versions of itself. We give here two examples:

Example 1: isotropic smoothing

This is perhaps the simplest case one can think of, and it leads to a smoothing process that is specified by the heat equation, and leads to linear diffusion. The functional that we aim to minimize is the squared image gradient norm integrated over the image domain $\Omega$ :

$E(I) = \int_\Omega ||\nabla I||^2 d\Omega$

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

$\begin{eqnarray*} I_{(t=0)} & = & I_0 \\ \frac{\partial I}{\partial t} & = & \Delta I \end{eqnarray*}$

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

Example 2: total variation flow

This time, we depart from the following functional:

$E(I) = \int_\Omega ||\nabla I|| d\Omega$

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:

$\begin{eqnarray*} I_{(t=0)} & = & I_0 \\ \frac{\partial I}{\partial t} & = & div\left(\frac{\nabla I}{||\nabla I||}\right)\end{eqnarray*}$

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

However, smoothing PDEs are not limited to those obtained by the Euler-Lagrange equations on some functional. Although this kind of model has the teoretical benefit of reducing some explicitely defined criterion, it does not always have the effect on the image that we would like. In general, what we seek are schemes that, in a predictable way, smooth the image regions with different strength depending on the amount of local variation, and with different strength along the local image gradient and isophote directions. For this, we can make modified PDEs that no longer represent the gradient of some functional, and whose evolution does not lead to a minimization of some defined quantity, but which may in turn produce results that are tailored to our requirements.

Anisotropic smoothing formulations

If we have a look at the time-dependent evolution equation from the two examples above, we see that they both can be written on the general form:

$\frac{\partial I}{\partial t} = div\left(g(||\nabla I||) \nabla I\right)$

with $g : R \rightarrow R$ a scalar function that is 1 in the first example and $\frac{1}{||\nabla I||}$ in the second. Without departing from a functional, we can create customized schemes by specifying this function $g$ 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 $g$ :

$g(||\nabla I||) = \exp\left(-\frac{||\nabla I||^2}{K^2}\right)$

$K$ 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.

Divergence-based PDEs

We can go further, however. Above, $g$

is a scalar function that only depends on the image gradient norm. A more general formulation would be:

$\frac{\partial I}{\partial t} = div\left(D \nabla I\right)$

where $D: \Omega \rightarrow P(2)$ is a tensor field defined on the domain $\Omega$ and with values in the space $P(2)$ of tensors described by positive-definite 2x2-matrices. (In case of three-dimensional images, i.e., volumes, we use $P(3)$ instead of $P(2)$ ). We call $D$ a field of diffusion tensors. Contrary to $g$ described above, $D$ depends not on the image gradient itself, but rather on the position in the image . Moreover, as $D(x,y)$ is a tensor rather than a scalar, its effect on the smoothing term also varies with the local direction of $\nabla I$ (matrix-vector multiplication $D \nabla I$ ). The smoothing at position $(x, y)$ will in general be strongest along the direction specified by the eigenvector with largest eigenvalue of $D(x,y)$ , 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 $D$ can be specified. For more on divergence-based PDEs, refer to section 2.1.4 of Tschumperle's thesis.

Smoothing based on oriented 1D Laplacians

When expanding the expression on the right hand side of the equation above, the placement of $D$

inside the $div$

term leads to a term which depend on the derivative of $D$

. For this reason, we cannot in general say that the smoothing process locally follows the directions specified by the tensor field $D$

. D. Tschumperle shows that a more "predictable" process can be obtained by using the smoothing geometry in a trace-based formulation:

$\frac{\partial I}{\partial t} = trace(TH)$

Here, $T$ is the smoothing geometry and $H = H(x, y)$ 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:

$I(t) = I_0 * G^{(T, t)}$

Where $I_0$ is the original image, $I(t)$ is the image at time $t$ of the evolution process, and $G^{(T, t)}$ is a bivariate oriented Gaussian kernel:

$G^{(T, t)}(x) = \frac{1}{4\pi t}\exp\left(-\frac{[x,y]T^{-1}[x,y]^T}{4t}\right)$

Defining a smoothing geometry

The purpose of the tensor field referred to as the "smoothing geometry" is to define local smothing directions at each point $(x, y)$

in the image. Typically, in regular areas of the image, we want to smooth equally in all directions, whereas in regions with important gradients, it is important to smooth the image in a way that does not destroy the gradient, i.e., smooth much less along the gradient than in the direction perpendicular to it (the isophote direction). Thus, in order to define the smoothing geometry tensor field, we should, for each point in the image, define two directions (2D-vectors) and the smoothing strength along these who directions (scalars). In order to define these, we base ourselves on a smoothed version of the image's structure tensor. (The reason we smooth the structure tensor is to determine a gradient for each pixel that is somewhat less "local", thus representing structures in the image that go beyond the size of one pixel and also obtaining gradient estimates that are more robust against image noise).

For a given image point $(x, y)$ , 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 $u$ and $v$ , and choose two corresponding smoothing strenths $c_1$ and $c_2$ , the smoothing geometry matrix representation for a given image point becomes:

$T = c_1 u u^T + c_2 v v^T$

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 $u$ and $v$ are the eigenvectors of the structure tensor, and the corresponding smoothing strengths are taken to be:

$c_1 = \frac{1}{(1 + \lambda_u + \lambda_v)^{p_1}}\;,\;\;\;\;c_2 = \frac{1}{(1 + \lambda_u + \lambda_v)^{p_2}}$

where $\lambda_u$ and $\lambda_v$ are the two eigenvalues of the structure tensor corresponding to $u$ and $v$ respectively. This choice is discussed briefly in section 2.1 of [Tschumperle06]. Note that as long as $p_1 < p_2$ , then the smoothing will be stronger along direction $u$ than along direction $v$ , 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.

Images with multiple channels

When smoothing images containing multiple channels, it is usually beneficial to smooth all image channels using the same smoothing geometry. If we smooth each image channel separately using a smoothing geometry computed independently of the other image channels, we risk ending up with an image where image edges are positioned differently in each channel, giving a undersirable blurred effect. Instead, we compute a common smoothing geometry based on the multi-channel definition of the image's structure tensor. This is covered in detail in section 2.2 of Tschumperle's PhD thesis.

Smoothing that preserves curvatures

We saw earlier that the trace-based formulation behaves locally as an oriented gaussian smoothing. The strength and orientation of this smoothing is directly reflected by the given structure tensor $T$

. This does a good job in preserving edges, but on curved structures, especially sharp corners, this leads to an undesirable rounding effect. The problem is that the gaussian kernel does not take curvatures into account. As a remedy, Tschumperle proposes an approach which can be seen as local filtering of the image with curved gaussian kernels . He develops the theory and algorithm in section 3 and 4 of his paper Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE's. In short, he splits the smoothing geometry tensor field $T$

into a sum of vector fields, and then performs line integral convolutions of the image with each of these vector fields and a univariate gaussian kernel function:

$G_t(s) = \frac{1}{4t}\exp\left(-\frac{s^2}{4t}\right)$

$t$ is here the timestep used. In other words he carries out 1-dimensional convolutions of the image with $G_t(s)$ along the streamlines defined by the vector fields obtained from .

decomposing the tensor field

In order to define vector fields along which to carry out line integral convolution, Tschumperle decompose the smoothing geometry $T$

in the following way:

$T = \frac{2}{\pi} \int_{\alpha = 0}^\pi (\sqrt{T}a_\alpha)(\sqrt{T}a_\alpha)^T$

where $a_\alpha = [\cos \alpha, \sin \alpha]^T$ and $\sqrt{T}$ is a tensor field that is the "square root" of $T$ , in the sense that $(\sqrt{T})^2=T$ . If $T = c_1 uu^T + c_2 vv^T$ , then $\sqrt{T} = \sqrt{c_1} uu^T + \sqrt{c_2} vv^T$ .

The resulting partial differential equation

Formally, the PDE that Tschumperle solves in the scheme proposed is:

$\frac{\partial I}{\partial t} = trace(TH) + \frac{2}{\pi}\nabla I^T \int_{\alpha=0}^\pi J_{\sqrt{T}a_\alpha}\sqrt{T}a_\alpha d\alpha$

where $J_{\sqrt{T}a_\alpha}$ stands for the Jacobian of the vector field $\Omega \rightarrow \sqrt{T}a_\alpha$ . 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 $\sqrt{T}a_\alpha$ .

The LSSEG implementation of Tschumperles curvature-preserving smoothing algorithm

We have implemented a 2D/3D version of the algorithm outlined above, which can be found in the example program app/smoothcurvpres.C that comes with LSSEG. Here, we will briefly present what is going on inside this implementation.

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