[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

details Non-linear Diffusion and Total Variation VIGRA

Classes

class  DiffusivityFunctor< Value >
 Diffusivity functor for non-linear diffusion. More...
 

Functions

template<... >
void anisotropicTotalVariationFilter (...)
 Performs Anisotropic Total Variation Regularization. More...
 
template<... >
void getAnisotropy (...)
 Sets up directional data for anisotropic regularization. More...
 
template<... >
void nonlinearDiffusion (...)
 Perform edge-preserving smoothing at the given scale. More...
 
template<... >
void nonlinearDiffusionExplicit (...)
 Perform edge-preserving smoothing at the given scale using an explicit scheme. More...
 
template<... >
void secondOrderTotalVariationFilter (...)
 Performs Anisotropic Total Variation Regularization. More...
 
template<... >
void totalVariationFilter (...)
 Performs standard Total Variation Regularization. More...
 

Detailed Description

Perform edge-preserving smoothing.

Function Documentation

void vigra::nonlinearDiffusion (   ...)

Perform edge-preserving smoothing at the given scale.

The algorithm solves the non-linear diffusion equation

\[ \frac{\partial}{\partial t} u = \frac{\partial}{\partial x} \left( g(|\nabla u|) \frac{\partial}{\partial x} u \right) \]

where t is the time, x is the location vector, u( x , t) is the smoothed image at time t, and g(.) is the location dependent diffusivity. At time zero, the image u( x , 0) is simply the original image. The time is proportional to the square of the scale parameter: $t = s^2$. The diffusion equation is solved iteratively according to the Additive Operator Splitting Scheme (AOS) from

J. Weickert: "Recursive Separable Schemes for Nonlinear Diffusion Filters", in: B. ter Haar Romeny, L. Florack, J. Koenderingk, M. Viergever (eds.): 1st Intl. Conf. on Scale-Space Theory in Computer Vision 1997, Springer LNCS 1252

DiffusivityFunctor implements the gradient-dependent local diffusivity. It is passed as an argument to gradientBasedTransform(). The return value must be between 0 and 1 and determines the weight a pixel gets when its neighbors are smoothed. Weickert recommends the use of the diffusivity implemented by class DiffusivityFunctor. It's also possible to use other functors, for example one that always returns 1, in which case we obtain the solution to the linear diffusion equation, i.e. Gaussian convolution.

The source value type must be a linear space with internal addition, scalar multiplication, and NumericTraits defined. The value_type of the DiffusivityFunctor must be the scalar field over wich the source value type's linear space is defined.

In addition to nonlinearDiffusion(), there is an algorithm nonlinearDiffusionExplicit() which implements the Explicit Scheme described in the above article. Both algorithms have the same interface, but the explicit scheme gives slightly more accurate approximations of the diffusion process at the cost of much slower processing.

Declarations:

pass 2D array views:

namespace vigra {
template <class T1, class S1,
class T2, class S2,
class DiffusivityFunc>
void
nonlinearDiffusion(MultiArrayView<2, T1, S1> const & src,
MultiArrayView<2, T2, S2> dest,
DiffusivityFunc const & weight, double scale);
template <class T1, class S1,
class T2, class S2,
class DiffusivityFunc>
void
nonlinearDiffusionExplicit(MultiArrayView<2, T1, S1> const & src,
MultiArrayView<2, T2, S2> dest,
DiffusivityFunc const & weight, double scale);
}

show deprecated declarations

Usage:

#include <vigra/nonlineardiffusion.hxx>
Namespace: vigra

MultiArray<2, float> src(w,h), dest(w,h);
float edge_threshold, scale;
...
nonlinearDiffusion(src, dest,
DiffusivityFunctor<float>(edge_threshold), scale);

show deprecated examples

Precondition:

scale > 0

See Also
vigra::DiffusivityFunctor
Examples:
smooth.cxx.
void vigra::nonlinearDiffusionExplicit (   ...)

Perform edge-preserving smoothing at the given scale using an explicit scheme.

See nonlinearDiffusion().

void vigra::totalVariationFilter (   ...)

Performs standard Total Variation Regularization.

The algorithm minimizes

\[ \min_u \int_\Omega \frac{1}{2} (u-f)^2\;dx + \alpha TV(u)\qquad\qquad (1) \]

where $ f=f(x)$ are the two dimensional noisy data, $ u=u(x)$ are the smoothed data, $ \alpha \ge 0 $ is the filter parameter and $ TV(u)$ is the total variation semi-norm.

Declarations:

namespace vigra {
template <class stride1,class stride2>
void totalVariationFilter(MultiArrayView<2,double,stride1> data,
MultiArrayView<2,double,stride2> out,
double alpha,
int steps,
double eps=0);
void totalVariationFilter(MultiArrayView<2,double,stride1> data,
MultiArrayView<2,double,stride2> weight,
MultiArrayView<2,double,stride3> out,
double alpha,
int steps,
double eps=0);
}

totalVariationFilter() implements a primal-dual algorithm to solve (1).

Input:

data: input data to be smoothed.
alpha: smoothing parameter.
steps: maximal number of iteration steps.
eps: The algorithm stops, if the primal-dual gap is below the threshold eps.

Output:

out contains the filtered data.

In addition, a point-wise weight ( $ \ge 0 $) for the data term can be provided (overloaded function).

Usage:

#include <vigra/tv_filter.hxx>

MultiArray<2,double> data(Shape2(width,height)); // to be initialized
MultiArray<2,double> out(Shape2(width,height));
MultiArray<2,double> weight(Shape2(width,height))); //optional argument in overloaded function, to be initialized if used
int steps; // to be initialized
double alpha,eps; // to be initialized
totalVariationFilter(data,out,alpha,steps,eps);

or

totalVariationFilter(data,weight,out,alpha,steps,eps);
Examples:
total_variation.cxx.
void vigra::getAnisotropy (   ...)

Sets up directional data for anisotropic regularization.

This routine provides a two-dimensional normalized vector field $ v $, which is normal to edges in the given data, found as the eigenvector of the structure tensor belonging to the largest eigenvalue. $ v $ is encoded by a scalar field $ \varphi $ of angles, i.e. $ v(x)=(\cos(\varphi(x)),\sin(\varphi(x)))^\top $.

In addition, two scalar fields $ \alpha $ and $ \beta $ are generated from scalar parameters $ \alpha_{par}$ and $ \beta_{par}$, such that

$ \alpha(x)= \alpha_{par}$ at edges,
$ \alpha(x)= \beta_{par}$ in homogeneous regions,
$ \beta(x)=\beta_{par}$ .

Declarations:

namespace vigra {
void getAnisotropy(MultiArrayView<2,double,stride1> data,
MultiArrayView<2,double,stride2> phi,
MultiArrayView<2,double,stride3> alpha,
MultiArrayView<2,double,stride4> beta,
double alpha_par,
double beta_par,
double sigma_par,
double rho_par,
double K_par);
}

Output:

Three scalar fields phi, alpha and beta.

Input:

data:two-dimensional scalar field.
alpha_par,beta_par:two positive values for setting up the scalar fields alpha and beta
sigma_par:non-negative parameter for presmoothing the data.
rho_par:non-negative parameter for presmoothing the structure tensor.
K_par:positive edge sensitivity parameter.

(see anisotropicTotalVariationFilter() and secondOrderTotalVariationFilter() for usage in an application).

Examples:
total_variation.cxx.
void vigra::anisotropicTotalVariationFilter (   ...)

Performs Anisotropic Total Variation Regularization.

The algorithm minimizes

\[ \min_u \int_\Omega \frac{1}{2} (u-f)^2 + \sqrt{\nabla u^\top A \nabla u}\;dx\qquad\qquad(2) \]

where $ f=f(x)$ are the noisy data, $ u=u(x)$ are the smoothed data, $ \nabla u $ is the image gradient in the sense of Total Variation and $ A $ is a locally varying symmetric, positive definite 2x2 matrix.

Matrix $ A $ is described by providing for each data point a normalized eigenvector (via angle $ \phi $) and two eigenvalues $ \alpha>0 $ and $ \beta>0 $.

getAnisotropy() can be use to set up such data $ \phi,\alpha,\beta $ by providing a vector field normal to edges.

Declarations:

namespace vigra {
template <class stride1,class stride2,class stride3,class stride4,class stride5,class stride6>
void anisotropicTotalVariationFilter(MultiArrayView<2,double,stride1> data,
MultiArrayView<2,double,stride2> weight,
MultiArrayView<2,double,stride3> phi,
MultiArrayView<2,double,stride4> alpha,
MultiArrayView<2,double,stride5> beta,
MultiArrayView<2,double,stride6> out,
int steps);
}

anisotropicTotalVariationFilter() implements a primal-dual algorithm to solve (2).

Input:

data:input data to be filtered.
steps:iteration steps.
weight :a point-wise weight ( $ \ge 0 $ ) for the data term.
phi,alpha and beta :describe matrix $ A $, see above.

Output:

out :contains filtered data.

Usage:

E.g. with a solution-dependent adaptivity cf. [1], by updating the matrix $ A=A(u)$ in an outer loop:

#include <vigra/tv_filter.hxx>

MultiArray<2,double> data(Shape2(width,height)); //to be initialized
MultiArray<2,double> out (Shape2(width,height));
MultiArray<2,double> weight(Shape2(width,height)); //to be initialized
MultiArray<2,double> phi (Shape2(width,height));
MultiArray<2,double> alpha(Shape2(width,height));
MultiArray<2,double> beta (Shape2(width,height));
double alpha0,beta0,sigma,rho,K; //to be initialized
int outer_steps,inner_steps;//to be initialized
out=data; // data serves as initial value
for (int i=0;i<outer_steps;i++){
getAnisotropy(out,phi,alpha,beta,alpha0,beta0,sigma,rho,K); // sets phi, alpha, beta
anisotropicTotalVariationFilter(data,weight,phi,alpha,beta,out,inner_steps);
}

[1] Frank Lenzen, Florian Becker, Jan Lellmann, Stefania Petra and Christoph Schnörr, A Class of Quasi-Variational Inequalities for Adaptive Image Denoising and Decomposition, Computational Optimization and Applications, Springer, 2012.

Examples:
total_variation.cxx.
void vigra::secondOrderTotalVariationFilter (   ...)

Performs Anisotropic Total Variation Regularization.

The algorithm minimizes

\[ \min_u \int_\Omega \frac{1}{2} (u-f)^2 + \sqrt{\nabla u^\top A \nabla u} + \gamma |Hu|_F\;dx \qquad\qquad (3) \]

where $ f=f(x)$ are the noisy data, $ u=u(x)$ are the smoothed data, $ \nabla u $ is the image gradient in the sense of Total Variation, $ A $ is a locally varying symmetric, positive-definite 2x2 matrix and $ |Hu|_F $ is the Frobenius norm of the Hessian of $ u $.

Matrix $ A $ is described by providing for each data point a normalized eigenvector (via angle $ \phi $) and two eigenvalues $ \alpha>0 $ and $ \beta>0 $. getAnisotropy() can be use to set up such data $ \phi,\alpha, \beta $ by providing a vector field normal to edges.

$ \gamma>0 $ is the locally varying regularization parameter for second order.

Declarations:

namespace vigra {
template <class stride1,class stride2,...,class stride9>
void secondOrderTotalVariationFilter(MultiArrayView<2,double,stride1> data,
MultiArrayView<2,double,stride2> weight,
MultiArrayView<2,double,stride3> phi,
MultiArrayView<2,double,stride4> alpha,
MultiArrayView<2,double,stride5> beta,
MultiArrayView<2,double,stride6> gamma,
MultiArrayView<2,double,stride7> xedges,
MultiArrayView<2,double,stride8> yedges,
MultiArrayView<2,double,stride9> out,
int steps);
}

secondOrderTotalVariationFilter() implements a primal-dual algorithm to solve (3).

Input:

data: the input data to be filtered.
steps : number of iteration steps.
out : contains the filtered data.
weight : point-wise weight ( $ \ge 0$ ) for the data term.
phi,alpha,beta: describe matrix $ A$, see above.
xedges and yedges : binary arrays indicating the presence of horizontal (between (x,y) and (x+1,y)) and vertical edges (between (x,y) and (x,y+1)). These data are considered in the calculation of $ Hu$, such that finite differences across edges are artificially set to zero to avoid second order smoothing over edges.

Usage:

E.g. with a solution-dependent adaptivity (cf.[1]), by updating the matrix $ A=A(u)$ in an outer loop:

#include <vigra/tv_filter.hxx>

MultiArray<2,double> data(Shape2(width,height)); //to be initialized
MultiArray<2,double> out(Shape2(width,height));
MultiArray<2,double> weight(Shape2(width,height))); //to be initialized
MultiArray<2,double> phi(Shape2(width,height);
MultiArray<2,double> alpha(Shape2(width,height);
MultiArray<2,double> beta(Shape2(width,height));
MultiArray<2,double> gamma(Shape2(width,height)); //to be initialized
MultiArray<2,double> xedges(Shape2(width,height)); //to be initialized
MultiArray<2,double> yedges(Shape2(width,height)); //to be initialized
double alpha0,beta0,sigma,rho,K; //to be initialized
int outer_steps,inner_steps;//to be initialized
out=data; // data serves as initial value
for (int i=0;i<outer_steps;i++){
getAnisotropy(out,phi,alpha,beta,alpha0,beta0,sigma,rho,K); // sets phi, alpha, beta
secondOrderTotalVariationFilter(data,weight,phi,alpha,beta,gamma,xedges,yedges,out,inner_steps);
}

[1] Frank Lenzen, Florian Becker, Jan Lellmann, Stefania Petra and Christoph Schnörr, A Class of Quasi-Variational Inequalities for Adaptive Image Denoising and Decomposition, Computational Optimization and Applications, Springer, 2012.

Examples:
total_variation.cxx.

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.1 (Fri May 19 2017)