DGtal
1.2.beta

Part of package DEC package.
The AmbrosioTortorelli functional (AT) is a classical relaxation of the MumfordShah functional [5]. Its purpose is to construct a piecewise smooth approximation of a function \( g \) defined in an open bounded domain \( \Omega \).
In the AmbrosioTortorelli functional, one wants to find a function \( u \in SBV(\Omega) \) which is a smooth approximation of the input data \( g \) (i.e. image in image processing applications). The AmbrosioTortorelli functional is defined by
\[ \displaystyle AT_{\varepsilon}(u,v) = \int_\Omega \alpha ug^2 + v^2 \nabla u^2 + \lambda \varepsilon \nabla v^2 + \frac{\lambda}{4 \varepsilon} 1v^2 dx, \]
for functions \( u,v \in W^{1,2}(\Omega)\) with \( 0 \leq v \leq 1 \).
In AT functional, function \( v \) is a smooth approximation of the set of discontinuities, and takes value close to 0 in this set, while being close to 1 outside discontinuities. A remarkable property of this functional is that it \( \Gamma \)converges to (a relaxation of) MumfordShah functional as \( \varepsilon \) tends to 0 (see [5]). The intuition is that a large \( \varepsilon \) induces a solution with a fuzzy set of discontinuities, which is then progressively narrowed to the crisp 1dimensional set of discontinuites as \( \varepsilon \) goes to 0.
The AT functional has been extended to vector data in [48].
We present a discrete calculus model of AT that works for any 2dimensional cell complex. For more information on discrete calculus, discrete exterior calculus or discrete differential calculus, see module Discrete exterior calculus, report [42], lecture [38], or the book [50] .
Our model is suited to a grid domain representing a 2D image made of pixels or a digital surface in 3D.
\[ \displaystyle AT(u,v) = \alpha \Sigma_{i=1}^n \langle u_i  g_i , u_i  g_i \rangle_2 + \langle v \wedge \bar{\mathbf{\star}} \bar{\mathbf{d_0}} \mathbf{\star} u_i, v \wedge \bar{\mathbf{\star}} \bar{\mathbf{d_0}} \mathbf{\star} u_i \rangle_1 \\ + \lambda \varepsilon \langle \mathbf{d_0} v , \mathbf{d_0} v \rangle_1 + \frac{\lambda}{4\varepsilon} \langle 1  v , 1  v \rangle_0. \]
For more details about discrete formulations and the numerical optimization of AT, we refer the reader to [27], [47] and [46].
The AT functional is a sum of quadratic terms, independently, but not simultaneously, convex in \( u \) and \( v \). It must have null derivative at optimum. We thus propose to alternatively solve for \( u \), then \( v \). The derivatives can be given explicitly as linear systems. To simplify notations, let \( \mathbf{A} := \mathbf{d_0} \), \( \mathbf{B} := \mathbf{d_1} \), \( \mathbf{A'} := \bar{\mathbf{\star}} \bar{\mathbf{d_1}} \mathbf{\star} \) and \( \mathbf{B'} := \bar{\mathbf{\star}} \bar{\mathbf{d_0}} \mathbf{\star} \).
We define \( \mathbf{M} \) the matrix which transforms a 0form into a 1form by averaging the values on the two edge extremities, i.e. \( \mathbf{M} = \frac{1}{2} \mathbf{d_0} \). It is used in the wedge product between 0form v and 1forms \( \bar{\mathbf{\star}} \bar{\mathbf{d_0}} \mathbf{\star} u_i \).
We get at optimum, if we assume trivial metrics:
\[ \displaystyle \left\{ \begin{array}{r @{=} l c} \left[\alpha \mathbf{Id}  {\mathbf{B'}}^T diag(\mathbf{M} v)^2 \mathbf{B'} \right] u & = \quad \alpha g & (1.1) \\ \left[ \frac{\lambda}{4\varepsilon} \mathbf{Id} + \lambda \varepsilon \mathbf{A}^T \mathbf{A} + \mathbf{M}^T diag(\mathbf{B'} u)^2 \mathbf{M} \right] v & = \quad \frac{\lambda}{4 \varepsilon} \mathbf{1} & (1.2) \end{array}\right. \]
for the derivative of AT. Since all matrices are symmetric, definite and positive, we use a Cholesky factorization to solve alternatively the two equations of each system.
\( \begin{array}{l l l l} \rlap{OPTIMIZEAT01( (g_1, ..., g_N ): 0forms, (\alpha, \lambda, \varepsilon_1, \varepsilon_2, \varepsilon_r): reals );} \\ \rlap{\mathbf{Var~} (u_1, ..., u_N): 0forms, (v, v'): 1forms, \varepsilon: real ; } \\ \rlap{\mathbf{begin}} & \\ & \rlap{\mathbf{foreach~} i \in \{1, 2, 3\} \mathbf{~do~} u_i \leftarrow g_i;} \\ & \rlap{v \leftarrow 1, \varepsilon \leftarrow \varepsilon_1;} \\ & \rlap{\mathbf{while~} \varepsilon \geq \varepsilon_2 \mathbf{~do~}} \\ & & \rlap{\mathbf{repeat}} \\ & & & v' \leftarrow v; \\ & & & \hbox{Solve Eq.(1.1) for all } u_i ; \\ & & & \hbox{Solve Eq.(1.2) for } v ; \\ & & \rlap{\mathbf{until~} v  v' \hbox{~is small;}} \\ & & \rlap{\varepsilon \leftarrow \varepsilon / \varepsilon_r;} \\ & \rlap{\mathbf{end}} \\ \rlap{\mathbf{end}} \end{array} \)
The AT functional is useful when trying to fit piecewisesmooth function to some input data, possibly corrupted with noise. It is therefore used in image processing for denoising [47] or for inpainting [46] . It can also work for approximating vector data, for instance a normal vector field [27] .
We show here that the AT functional can enhance a coarse approximation of the normal vector field onto a digital surface. Here, we compute such a normal vector field using Integral Invariant methods (see Integral invariant curvature estimator 2D/3D). It has the drawback to smooth normals around sharp features. Therefore, in example exampleSurfaceATNormals.cpp, we use the AT functional to compute a piecewise smooth approximation of this normal vector. Normals are smoothed everywhere except around sharp features where they have sharp transitions. Furthermore, the set of discontinuities of AT functional (places where primal 0form v is close to 0) indicates where are located sharp features onto the digital surface.
To detail a little bit how to do that, you have first to build a digital surface. For instance, you can use Shortcuts to do so quickly (see Shortcuts (for the impatient developper)):
Then you just have to build a discrete calculus over the domain where you want to compute your piecewise smooth approximation. Here, we build the calculus over the whole surfel range surfels:
After that, the AT solver can be build on this calculus. Then it is initialized with input data you wish to approximate (here the normal vector field ii_normals). Finally it suffices to call the method ATSolver2D::solveGammaConvergence to use the coarsetofine \( \Gamma \)convergence optimization method.
You recover the piecewisesmooth approximation of the input vector field with ATSolver2D::getOutputVectorFieldU2, and the function giving the locii of discontinuities with ATSolver2D::getOutputScalarFieldV0.
See exampleSurfaceATNormals.cpp for more details.
You may also avoid using ATSolver2D directly and instead use method ShortcutsGeometry::getATVectorFieldApproximation or ShortcutsGeometry::getATScalarFieldApproximation which directly builds the calculus onto the given surface and outputs the results.