i3d::ImplicitAOSScheme< VOXELIN, VOXELOUT > Class Template Reference

This class is root abstract class for wide family of nonlinear diffusion filters. This class implements the skeleton of numerical semi-implicit AOS scheme. Instances of this class will not produce any usable output. To obtain fully functional nonlinear diffusion filter try to instantiate some offspring of this class. More...

#include <LevelSet.h>

Inheritance diagram for i3d::ImplicitAOSScheme< VOXELIN, VOXELOUT >:

List of all members.


Public Member Functions
virtual bool	Convergence ()
	If returns true the iteration proccess stops.
virtual void	ComputeChange ()
virtual void	ApplyChange ()
	ImplicitAOSScheme (Image3d< VOXELIN > &in, Image3d< VOXELOUT > &out)
virtual void	PrepareDiagonals (size_t start, size_t incr, Vector3d< size_t > &size, size_t axe)
	Computes the coefficients in the diagonals of the matrix.
virtual void	ComputeDiffusivity ()
	This function computes the diffusivity function g for each voxel of the image.
void	SetSymetric (bool s=true)
	If Set to true, the threediagonal matrix A is considered to be symetric and the faster routines to solve this systems will be used.
bool	GetSymetric ()
	Get the symetric flag.
virtual void	SetSpacing (float x, float y, float z)
	Set spacing of the cartesian grid. Typical values are (1.0,1.0,1.0 * z).
virtual void	SetSpacing (Vector3d< float > vec)
	Set spacing of the cartesian grid. Typical values are (1.0,1.0,1.0 * z).
void	AllocateDiagonals (size_t Max)
Public Attributes
GradientMagnitudeEstimator< VOXELOUT > *	GradMag
valarray< VOXELOUT > *	GradMagData
	Pointer to valarray with values of gradient estimation.
Protected Attributes
valarray< VOXELOUT >	BufferArray
	The solution of the AOS scheme is saved in this array.
valarray< VOXELOUT >	Diag
	The coefficients of the main diagonal are saved in this valraay.
valarray< VOXELOUT >	UDiag
	The coefficients of the upper diagonal are saved in this valraay.
valarray< VOXELOUT >	LDiag
	The coefficients of the lower diagonal are saved in this valraay.
valarray< VOXELOUT >	row
	The right hand side of the linear system is saved in this valarray.

Detailed Description

template<class VOXELIN, class VOXELOUT>
class i3d::ImplicitAOSScheme< VOXELIN, VOXELOUT >

The main idea of nonlinear diffusion filtering is to smooth the image inside relatively continuous regions and do not smooth the image near edges. Nonlinear diffusion filtering is ussualy performed with explicit schemes (see ExplicitScheme and PMFilter class for examples). This schemes are stable only for small time steps (i.e. $\tau < 0.25$ for 2D images, see documentation of MCFFilter class). This limits their practical use and leads to poor efficiency. This abstract class implements the skeleton of numerical semi-implicit AOS (additive oerator splitting). The AOS scheme is absolute stable fo all time steps, but there is one limitation to the size of time step $\tau$ . The scheme tends to produce larger and larger errors for time step $\tau > 5.0$ (considering that the distance between grid points is $ h = 1$ for all dimensions). The offsprings of this abstract class modificate the input image by solving following partial differential equation:

$u_t = {\mathrm {div}} (g(|\nabla u|) \nabla u),$

where $ g(s) $ is decreasing function of the image gradient, and $ u_0$ is the input image. The equation is discretised on cartesian grid, the time is discretized by forward Euler operator. The image gradient is discretised by central differences. The divergence operator is discretised also with central differences. The semi-implicit scheme leads to following discretisation of the equation:

$\frac {u_{i,j}^{n+1} - u_{i,j}^n}{\tau} = {\mathrm {div}} (g(|\nabla u_{i,j}^n|) \nabla u_{i,j}^{n+1}) \\$

The AOS schemes solves this discretisation dimension by dimension, so following equations will be shown only for the direction of x axe.

$\begin{array}[b]{l} \frac {u_{i}^{n+1} - u_{i}^n}{\tau} = {\mathrm {div}} (g(|\nabla u_{i}^n|) \nabla u_{i}^{n+1}) \\ \\ \frac {u_{i}^{n+1} - u_{i}^n}{\tau} =\frac{g_{i+1}^n + g_{i}^{n}}{2h^2}(u^{n+1}_{i+1} - u^{n+1}_{i}) - \frac{g_{i-1}^n + g_{i}^{n}}{2h^2}(u^{n+1}_{i} - u^{n+1}_{i-1}) \\ \end{array}$

This can be written in more compact form:

$\begin{array}[b]{l} \frac {u_{i}^{n+1} - u_{i}^n}{\tau} =\frac{g_{i+1}^n + g_{i}^{n}}{2h^2}(u^{n+1}_{i+1} - u^{n+1}_{i}) + \frac{g_{i-1}^n + g_{i}^{n}}{2h^2}(u^{n+1}_{i-1} - u^{n+1}_{i}) \\ \\ \frac {u_{i}^{n+1} - u_{i}^n}{\tau} =\sum \limits_{j \in N(i)} \frac{g_{j}^n + g_{i}^{n}}{2h^2}(u^{n+1}_{j} - u^{n+1}_{i}) \\ \end{array}$

which in matrix-vector notation becomes to

$\begin{array}[b]{l} \frac {u^{n+1} - u^n}{\tau} =A(u^n)u^{n+1} \\ \\ (I - \tau A(u^n)) u^{n+1} = u^n \\ \end{array}$

with $A(u^n) = (a_{ij}(u^n))$ and

$\begin{array}[b]{l} a_{ij}(u^n) = \frac{g_{j}^n + g_{i}^{n}}{2h^2} \quad (j \in N(i)) \\ \\ a_{ij}(u^n) = - \sum \limits_{k \in N(i)} \frac{g_{k}^n + g_{i}^{n}}{2h^2} \quad (j = i) \\ \\ a_{ij}(u^n) = 0 \quad {\mathrm {else}}\\ \\ \\ g_i^n = g(\sqrt {u_x^2 + u_y^2}) \quad {\mathrm {discretised \, with \, central \, differences}}\\ \end{array}$

The AOS scheme solves this equation dimension by dimension by adding solutions of the folowing equation in particular directions (x,y,z) to the aggregate final solution.

$u^{n+1} = \frac{1}{m} (I - m \tau A(u^n))^{-1} u^n$

where m is number of dimensions. The solution of above equation in each direction leads to solving the linear system of equations in one iteration step. If we order the voxels extra for each direction, the matrix A becomes threediagonal and there exist fast methods for finding solutions of such matrices. This class uses LAPACK routines to find solutions of three diagonal linear system. The i3d library must be linked with lapack library to obtain fully functionality of this class.

The implementation of this class is based on the article
Joachim Weickert, Bart M. ter Haar Romeny, Max A. Viergever, Efficient and Reliable Schemes for Nonlinear Diffusion Filtering , IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 6, NO. 3, PP. 398-410, MARCH 1998.

The typical code for using the anisotropic diffusion filters, which are offsprings of this class, is shown bellow. The example code uses the AOS implementation of CLMC diffusion filter (class AOSCLMCFilter).

    .....
   
   Image3d<GRAY8> i;
   Image3d<float> fimg;
   i.ReadImage('input.tif'); // read input image
   AOSCLMCFilter<GRAY8,float> CLMC(i,fimg); // instantiate  the filter  (the two parameters are the input and the output image)
   CLMC.SetMaximumNumberOfIterations(2); // Set the number of iterations to 2
   CLMC.SetSpacing(1.0,1.0,2.0); // Set the spacing of spatial discretisation to 1,1,2
   CLMC.SetSigma(1.0); // Set the sigma of the gaussian blur to 1.0
   CLMC.SetConductance(3.0); // Set the conductance parameter to 3.0
   CLMC.Execute(); // Perform the work of the filter
   FloatToGray8NoWeight(fimg,i); //Convert the output float image to gray8 image
   i.SaveImage('output.tif'); // save the output of the filter

   ......

Author:: Jan Hubený xhubeny@fi.muni.cz

Date:: 2005

Example input and output:

Left: Input to anisotropic diffusion, Right: Output of Total variation filter (AOSTVFilter class), 10 iterations with delta t = 5.0