Numerical Implementation of the Ambrosio-Tortorelli Functional Using Discrete Calculus and Application to Image Restoration and Inpainting

Abstract

The Mumford-Shah (MS) functional is one of the most influential variational model in image segmentation, restoration, and cartooning. Difficult to solve, the Ambrosio-Tortorelli (AT) functional is of particular interest, because minimizers of AT can be shown to converge to a minimizer of MS. This paper takes an interest in a new method for numerically solving the AT model [11]. This method formulates the AT functional in a discrete calculus setting, and by this way is able to capture the set of discontinuities as a one-dimensional set. It is also shown that this model is competitive with total variation restoration methods. We present here the discrete AT models in details, and compare its merit with recent convex relaxations of AT and MS functionals. We also examine the potential of this model for inpainting, and describe its implementation in the DGtal library, an open-source project.

This work has been partly funded by CoMeDiC ANR-15-CE40-0006 research grant.

1 Introduction

This paper takes an interest in a new numerical method to solve the Ambrosio-Tortorelli (AT) functional [1]. This functional is a well known non-convex relaxation of the famous Mumford-Shah (MS) functional [17], which has many applications in image segmentation [2, 6, 18], restoration [2, 4, 5, 8, 12, 19, 23], cartooning [20] or inpainting [10]. This new method for solving AT redefines this functional into a discrete calculus setting and is presented in details in [11].

Contributions. The present paper is a companion to the original paper [11], and complements it on several points. First of all, the possible formulations of AT in discrete calculus are discussed and compared in more details. Secondly, we examine other convex relaxations of the Ambrosio-Tortorelli and Mumford-Shah functionals (Kee and Kim [15], Strekalovskiy and Cremers [20]) and we compare our approach to them. Third, we study the influence of parameters in the method, as well as the possible scale spaces induced by the model. Fourth, we test the possible application of our discrete AT model to inpainting problems and discuss its potential and drawbacks. Last, we describe its implementation in the DGtal library [22], which is an open source C++ library for image and geometry processing of digital data.

Outline. The paper is organized as follows. We first briefly introduce the MS and AT functionals as well as the many variants and discretizations proposed in the literature (Sect. 2). We then present our discrete formulations of the AT functional, which is expressed in a discrete calculus framework, and the algorithm for its numerical optimization (Sect. 3). Our approach is compared to several state-of-the-art convex relaxations of MS or AT in Sect. 4. Parameters and their associated scale spaces are discussed in Sect. 5. The potential of our discrete AT for inpainting is encompassed in Sect. 6. Section 7 concludes by describing the implementation of our AT models in the DGtal library.

2 The Mumford-Shah and Ambrosio-Tortorelli Functionals

The Mumford-Shah (MS) model was introduced in [17] to solve image segmentation and restoration problems. This model was seminal to many important contributions in image processing and analysis, such as the Rudin-Osher-Fatemi (ROF) restoration method, the total variation (TV) denoising, Chan-Vese or Boykov-Jolly segmentation methods. The main idea is to see the input grey-level image as a function $g \in L^{\infty }(\varOmega )$ where $\varOmega$ is some domain of $\mathbb {R}^2$, and to seek the solution as a piecewise smooth function u that approaches g. The loci of discontinuities K is also a variable of the problem. Another strong idea is that the “size” of K is penalized. So the MS model is the optimal solution (K, u) to the following functional:

$$\begin{aligned} \mathcal {MS}(K,u) = \alpha \displaystyle \int _{\varOmega \backslash K} |u - g|^{2} \mathop {}\mathopen {}\mathrm {\ dx}+ \displaystyle \int _{\varOmega \backslash K} |\nabla u|^{2} \mathop {}\mathopen {}\mathrm {\ dx}+ \lambda \mathop {}\mathopen {}\mathcal {H}^{1}(K \cap \varOmega ) \end{aligned}$$

(1)

where $\alpha , \lambda >0$ and $\mathop {}\mathopen {}\mathcal {H}^{1}$ denotes the 1-dimensional Hausdorff measure. We recognize in the first term the $L_2$ approximation to the input image, the second term requires the smoothness of u, while the last term penalizes the length of discontinuities K.

The MS functional is difficult to solve in this form, so many relaxations have been proposed to overcome this issue. We take an interest in the Ambrosio-Tortorelli relaxation [1], defined as:

$$\begin{aligned} \mathrm {AT}_{\varepsilon }(u,v) = \displaystyle \int _\varOmega \alpha |u-g|^2 + v^2 |\nabla u|^2 + \lambda \varepsilon |\nabla v|^2 + \dfrac{\lambda }{4 \varepsilon } |1-v|^2 \mathop {}\mathopen {}\mathrm {\ dx}, \end{aligned}$$

(2)

for functions $u,v \in W^{1,2}(\varOmega )$ with $0 \le v \le 1$.

In (2), function u is a smooth approximation of the input image g. 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 $\varGamma$-converges to (a relaxation of) the MS functional as $\varepsilon$ tends to 0 [1]. The intuition is that a large $\varepsilon$ induces a solution with a fuzzy set of discontinuities, which is then progressively narrowed to the crisp 1-dimensional set of discontinuites as $\varepsilon$ goes to 0.

However the minimization of (2) with standard discretization schemes such as finite differences is not possible. Indeed, to achieve convergence, parameters $\varepsilon$, the grid-step h and the ratio $\frac{\varepsilon }{h}$ must all tend to zero [4]. So in practice we set $\varepsilon \approx 5h$ but the resulting set of discontinuities is thick and fuzzy (sometime more than 6 pixels wide), and the restoration is very poor in these regions. To overcome this issue, Bourdin and Chambolle [4] proposed a finite element method with adaptive domain refinement around discontinuities. It will be shown later that this method is still not satisfactory, despite the added complexity. These difficulties explains why the AT model is not a popular approach in image processing, and why TV, ROF, or graph cuts methods have been heavily used instead.

For more details about other similar variational models for image processing and more background about related works, we refer the reader to [11].

3 Discrete Calculus Formulations of AT Functional

In [11], a different formulation of the AT functional was proposed, which retains its main characteristics (piecewise smoothness, penalization of the length of the discontinuities), while being able to capture thin discontinuities. This reformulation uses the setting of discrete exterior calculus (DEC), as described in the context of computer graphics in [7] and image analysis in [13]. The advantage of this setting is that one can control easily where the variables are defined, and where they are coupled. In particular, the function v does not become identically one even when $\varepsilon$ tends to zero.

Discrete exterior calculus. The idea is to decompose the image domain $\varOmega$ into a cell complex $K$. Here, the faces of $K$ are simply the pixels of the image, the edges of $K$ are the sides shared by two pixels, while the vertices of $K$ are the those shared by four pixels. For technical reasons, the dual complex $\bar{K}$ of $K$ is also needed and is defined in the usual way: a vertex of $K$ is thus associated one-to-one to a face of $\bar{K}$, etc. A discrete k-dimensional form is just a map that associates a scalar to a k-dimensional cell.¹ If the number of k-dimensional cells of $K$ is denoted by $n_k$, then a discrete k-form is simply represented by a column vector of size $n_k \times 1$.

We then define the usual linear operators between k-forms (see Fig. 1). We denote by $\mathbf {d}_{k}$ and $\bar{\mathbf {d}}_{k}$ the standard discrete exterior primal and dual derivative operators. The derivative operator $\mathbf {d}_{0}$ is the oriented vertex-edge incidence matrix of $K$, and has size $n_1 \times n_0$. Since the complex K is the regular grid with edge size 1, the dual derivative $\bar{\mathbf {d}}_{1}$ is the transpose of $\mathbf {d}_{0}$. Similarly, the primal derivative $\mathbf {d}_{1}$ is the oriented edge-face incidence matrix of $K$, of size $n_2 \times n_1$, and the dual derivative $\bar{\mathbf {d}}_{0}$ is its transpose.

The Hodge star is an operator that maps a differential form onto its complementary differential form. Discrete Hodge star operators $\mathbf {\star }$ sends k-forms of the primal complex $K$ onto $n-k$-forms of the dual complex $\bar{K}$ (see again Fig. 1). In our setting, these operators reduce to identity operations except for a change of sign in some cases, and they also change the meaning of vectors. For instance, if f is a 0-form on the vertices of $K$, $\mathbf {\star }f$ is a dual 2-form defined on the faces of $\bar{K}$. Clearly, it is also a vector of size $n_0 \times 1$ since there are $n_0$ faces in $\bar{K}$.

We define $\mathbf {M}_{01}$ the matrix which transforms a 0-form into a 1-form by averaging the values on the two edge extremities, i.e. $\mathbf {M}_{01}= \frac{1}{2} |\mathbf {d}_{0}|$. Moreover, we use the edge laplacian defined in [13] by $\bar{\mathbf {\star }} \bar{\mathbf {d}}_{0} \mathbf {\star } \mathbf {d}_{1} + \mathbf {d}_{0} \bar{\mathbf {\star }} \bar{\mathbf {d}}_{1} \mathbf {\star }$.

Fig. 1.

Primal and dual complex and common Discrete Exterior Calculus (DEC) operators.

Discrete formulations of AT. We first set u and g to live on the faces and v to live on the vertices and edges. Pixels are faces, so functions u and g are 2-forms since they represent the gray levels of each pixel. On the contrary, v is an estimation of the set of discontinuities of u, and should be of null Hausdorff-1 measure when $\varepsilon \rightarrow 0$. Thus we set v in-between cells of non null measure, so in this case on vertices as a 0-form, and on edges by averaging with $\mathbf {M}_{01}$. We call this formulation $\mathrm {AT}^{2,0}_{\varepsilon }$. Looking at (2), the DEC reformulation is straightforward, except for the second term, where v is a 0-form and $\nabla u$ a 1-form. Hence we use matrix $\mathbf {M}_{01}$ to transport the 0-form v onto edges by simple averaging:

$$\begin{aligned} \mathrm {AT}^{2,0}_{\varepsilon }(u,v) =\,&\alpha \langle {u - g}{,}{u - g}\rangle _{2} + \langle {\mathbf {M}_{01}v \odot \bar{\mathbf {\star }} \bar{\mathbf {d}}_{0} \mathbf {\star } u}{,}{\mathbf {M}_{01}v \odot \bar{\mathbf {\star }} \bar{\mathbf {d}}_{0} \mathbf {\star } u}\rangle _{1} \nonumber \\&+ \lambda \varepsilon \langle { \mathbf {d}_{0} v }{,}{ \mathbf {d}_{0} v }\rangle _{1} + \frac{\lambda }{4\varepsilon } \langle {1 - v}{,}{1 - v}\rangle _{0}. \end{aligned}$$

(3)

where $\odot$ denotes the point-wise multiplication.

A second possibility is to define u and g on the vertices and v on the edges. We denote this formulation $\mathrm {AT}^{0,1}_{\varepsilon }$. Contrary to the previous formulation, the gray levels are seen as point mass on the center of pixels, so that functions u and g are both 0-forms. An alternative interpretation is to say that pixel values are dual 2-forms, while v is a dual 1-form in between u. It follows:

$$\begin{aligned} \mathrm {AT}^{0,1}_{\varepsilon }(u,v) =\,&\alpha \langle {u - g}{,}{u - g}\rangle _{0} + \langle {v \odot \mathbf {d}_{0} u}{,}{v \odot \mathbf {d}_{0} u}\rangle _{1} \nonumber \\&+ \lambda \varepsilon \langle { (\mathbf {d}_{1} + \bar{\mathbf {\star }} \bar{\mathbf {d}}_{1} \mathbf {\star }) v }{,}{ (\mathbf {d}_{1} + \bar{\mathbf {\star }} \bar{\mathbf {d}}_{1} \mathbf {\star }) v }\rangle _{1} + \frac{\lambda }{4\varepsilon } \langle {1 - v}{,}{1 - v}\rangle _{1}. \end{aligned}$$

(4)

We further extend (3) and (4) to the vectorial case by associating to each component of the input image n forms $\{g_1,\dots , g_n\}$, and the corresponding forms $\{u_1$,..., $u_n\}$, and summing over the coordinates. In this paper, we will use $n=1$ for gray-level images and $n=3$ for RGB color images.

Optimization. The functionals $\mathrm {AT}^{2,0}_{\varepsilon }$ and $\mathrm {AT}^{0,1}_{\varepsilon }$ are both sums of quadratic terms, independently, but not simultaneously, convex in u and v. They 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 }$. It is worth to note $\mathbf {A'}$ and $\mathbf {B'}$ are the respective transpose of $\mathbf {A}$ of $\mathbf {B}$, except on the boundary of the image. We get at optimum:

$$\begin{aligned} \left\{ \begin{array}{r @{=} l} \left[ \alpha \mathbf {Id}- {\mathbf {B'}}^{\intercal } \mathrm {diag}\left( \mathbf {M}_{01}v\right) ^2 \mathbf {B'}\right] u &{}\alpha g, \\ \left[ \frac{\lambda }{4\varepsilon } \mathbf {Id}+ \lambda \varepsilon {\mathbf {A}}^{\intercal }\mathbf {A}+ {\mathbf {M}_{01}}^{\intercal } \mathrm {diag}\left( \mathbf {B'}u\right) ^2 \mathbf {M}_{01}\right] v &{}\frac{\lambda }{4 \varepsilon } \mathbf {1}. \end{array}\right. \end{aligned}$$

(5)

for the derivative of $\mathrm {AT}^{2,0}_{\varepsilon }$, and

$$\begin{aligned} \left\{ \begin{array}{r @{=} l} \left[ \alpha \mathbf {Id}- {\mathbf {A}}^{\intercal } \mathrm {diag}\left( v\right) ^2 \mathbf {A}\right] u &{}\alpha g, \\ \left[ \frac{\lambda }{4\varepsilon } \mathbf {Id}+ \lambda \varepsilon ({\mathbf {A'}}^{\intercal }\mathbf {A'}+ {\mathbf {B}}^{\intercal }\mathbf {B}) + \mathrm {diag}\left( \mathbf {A}u\right) ^2 \right] v &{}\frac{\lambda }{4 \varepsilon } \mathbf {1}. \end{array}\right. \end{aligned}$$

(6)

for the derivative of $\mathrm {AT}^{0.1}_{\varepsilon }$. Since all matrices are symmetric, definite and positive, we use a Cholesky factorization to solve alternatively the two equations of (5) (resp. (6)). Because of the $\mathbf {Id}$ additive term, the left-hand side is full rank, yielding a unique minimum at each iteration. It is a known result in convex analysis linked to block coordinate descent algorithms [3, Prop. 2.7.1], that these iterations must converge to a stationary point of $\mathrm {AT}^{2,0}_{\varepsilon }$ (resp. $\mathrm {AT}^{0,1}_{\varepsilon }$). In the optimization process, we start with a large enough $\varepsilon _1$, in order to capture the discontinuities, and after convergence for a given $\varepsilon$, we decrease it by dividing it by $\varepsilon _r$, until $\varepsilon = \varepsilon _2$. We denote this process by $\varepsilon = \varepsilon _1 \searrow \varepsilon _2$. The main loop stops when v reaches stability. Figure 2 illustrates that our AT model are able to capture discontinuities as a 1-dimensional set, since v is defined on 0- and 1-dimensional cells.

Fig. 2.

Minimization of AT with finite differences (FD), finite elements (FE) from [4] and our approach $\mathrm {AT}^{2,0}_{\varepsilon }$ on a triple point image with PSNR = 20.3 dB.

Fig. 3.

In this figure we have added a spatially-variant noise to the Mandrill image. In this case it is not possible to find an optimal $(\alpha , \lambda )$ pair that denoises and segments the image appropriately with our current model.

Limit cases. As is the case with many image restoration methods, our approach assumes stationary noise and blur. Some spatially variant noise occurs in nature, for instance Poisson noise [21]. As is illustrated on Fig. 3, when the noise is spatially variant, suboptimal results are achieved. This is due to the parameters $\alpha$ and $\lambda$ being constant for the entire image. For future work, we plan to let $\alpha$ and $\lambda$ depend on x and thus choose locally the optimal pair $(\alpha , \lambda )$, and to adapt AT to use data fidelity terms that are appropriate for some non-stationary noise models [14].

Fig. 4.

Comparison of the minimization of convex relaxations of AT (extracted from [15]), MS [20] and our approach on images extracted from [15]. In all cases, parameters were fixed to get best possible SNR.

Fig. 5.

Comparison with competing methods. Noise level is specified in the first column.

4 Comparisons with Other Relaxations of MS and AT Functionals

A frequent criticism about MS and AT functionals is that they are non convex. Current optimization techniques thus do not extract the global optimum, and results are subject to initialization. This is for instance why TV models have gained so much popularity, despite known problems (most notably staircasing instead of smoothing). Other people have followed a different path and proposed convex approximations of MS or AT. Strekalovskiy and Cremers have proposed a truncated quadratic potential to approach directly the MS model [20], which leads them to a fast algorithm for segmentation and cartooning. Kee an Kim proposed in [15] three convex relaxations of AT: a quadratic relaxation, a linear relaxation, and the last one based on a factorization theorem due to McCormick [16]. They used a gradient descend algorithm to compute the minimizers, and observed that the best results are obtained with the linear relaxation. We compare in Fig. 4 the results of Kee and Kim, the results of Strekalovskiy and Cremers, and our approach on the images used in [15].

We can see that the McCormick relaxation (Fig. 4(c)) retains many artifacts while the quadratic relaxation (Fig. 4(d)) leads to blurry reconstructions. The results obtained with the linear relaxation (Fig. 4(e)) are piecewise smooth, but we observe a poor reconstruction of angles in the vectorial case. Furthermore we see on Fig. 4(f) that they obtain diffuse and sometimes false contours.

The convex relaxation of Strekalovskiy and Cremers leads to piecewise smooth results (Fig. 4(b)), with sharp discontinuities, similar to our approaches (Fig. 4(g)–(h)).

In Fig. 5, we compare the results of TV model implemented in [9], the results of Strekalovskiy et al.  [20], and our approach on synthetic images. For all methods, we set the parameters that maximize the PSNR (see Table 1). We observe the characteristic staircasing effect with the TV method. Strekalovskiy et al. and our method are robust to noise and yield piecewise smooth results.

Table 1. Best PSNRs results are in bold.

5 Influence of Parameters and Scale Spaces

The $\varGamma$-convergence parameter $\varepsilon$ controls the thickness of the contours. Large values of $\varepsilon$ convexify the AT functional and help to detect the discontinuities. Then, as $\varepsilon$ goes to 0, the contours become thinner and thinner (see Fig. 6).

Fig. 6.

$\varepsilon$ scale space, for fixed $\lambda = 0.006$ and $\alpha = 0.1$.

On the other hand, for a large enough and fixed $\lambda$, the minimization of AT is equivalent to a diffusion with a small perturbation dependent on $\alpha$. Hence parameter $\alpha$ can be chosen for an initial blur of the input data: the smaller $\alpha$ is, the stronger the initial blur is (see Fig. 7).

Fig. 7.

$\alpha$ scale space, for fixed $\lambda = 1.0$ and $\varepsilon = 2 \searrow 0.25$.

Finally, for fixed $\alpha$ (chosen so as to remove noise), the parameter $\lambda$ controls the length of discontinuities. Hence, for a large enough $\lambda$, the set of discontinuities is empty because the minimization of AT is equivalent to a diffusion. As the parameter $\lambda$ decreases, more and more discontinuities are detected and delineated, as we can see in Fig. 8. Note that results of $\mathrm {AT}^{0,1}_{\varepsilon }$ and $\mathrm {AT}^{2,0}_{\varepsilon }$ are almost identical for this image.

Fig. 8.

$\lambda$ scale space of $\mathrm {AT}^{2,0}_{\varepsilon }$, for fixed $\alpha = 0.48$ and $\varepsilon = 2 \searrow 0.25$.

6 Potential of the Model in Image Inpainting

It is easy to adapt the AT models for image inpainting. It suffices to have a parameter $\alpha$ that varies across the image. More precisely, it is set to 0 wherever the user selected a damaged zone while it is kept to its normal value elsewhere. Furthermore, forms u are initialized with random values at these places.

We have checked if our discrete AT models can be interesting for image inpainting. As illustrated on Fig. 9, it appears that model $\mathrm {AT}^{0,1}_{\varepsilon }$ introduces too many artefacts on artificial examples. Indeed, it optimizes the L1-length of discontinuities, which is not a natural inpainting for a human eye. The problem is that the L1-length has a very flat minimum, which is composed of all the discrete paths in some rectangle joining two constraints. Model $\mathrm {AT}^{0,1}_{\varepsilon }$ picks anyone of these paths: this is governed by the order of cells in the linear systems.

However, model $\mathrm {AT}^{2,0}_{\varepsilon }$ seems much more interesting. It is able to detect the two triple points that optimize a four region contact as well as the $120{^\circ }$ angle for a three region contact. The reason is that this model averages the L1-length of discontinuities especially when $\varepsilon$ is large. Hence, the diffused L1-length is closer to the natural L2-length and the model chooses a better L1-path for large $\varepsilon$. As $\varepsilon$ decreases, discontinuities stay at the correct location while getting crispier. This explains why triple points are correctly detected by $\mathrm {AT}^{2,0}_{\varepsilon }$, if we start from a large enough $\varepsilon$, which be must proportional to the width of the damaged zone. In Fig. 9, setting $\varepsilon =4 \searrow 0.25$ was good for all examples.

Of course, variational methods for image inpainting should be combined with some non-local/texture filtering to compete with state-of-the-art approaches, but we may conclude that $\mathrm {AT}^{2,0}_{\varepsilon }$ is a candidate for such a hybrid inpainting method.

Fig. 9.

Using discrete AT functionals for image inpaiting. Top row: input images. Second row: results of $\mathrm {AT}^{0,1}_{\varepsilon }$ with $\varepsilon = 4 \searrow 0.25$. Third row: results of $\mathrm {AT}^{2,0}_{\varepsilon }$ with $\varepsilon = 4 \searrow 0.25$. Fourth row: results of $\mathrm {AT}^{2,0}_{\varepsilon }$ with $\varepsilon = 1 \searrow 0.25$. Fifth row: inpainting of real picture with $\mathrm {AT}^{0,1}_{\varepsilon }$ and $\varepsilon = 4 \searrow 0.25$.

7 Implementation in the DGtal Library

All presented material and programs have been implemented as a tool in the DGtal library [22], which is an open-source C++ library, bundled with a set of tools for digital geometry and image processing. Our discrete AT implementation relies heavily on the Discrete Exterior Calculus package. This package provides data-structures and algorithms for manipulating discrete forms on cellular complexes, linear operators such as derivatives, and wrappers for linear system solvers.

We provide two command-line tools in the imageProcessing package of DGtalTools (http://dgtal.org/doc/tools/nightly/imageProcessing), one called at-u0-v1 that implements $\mathrm {AT}^{0,1}_{\varepsilon }$, the other called at-u2-v0 that implements $\mathrm {AT}^{2,0}_{\varepsilon }$. They provide exactly the same functionalities and options.

Among functionalities, these tools can process grey-level or color images given in PBM format. You may either restore the whole image or perform inpainting by giving a mask image specifying damaged zone. SNR computation is provided if you provide ground truth. Both tools export the restored image u as well as an (possibly zoomed) image displaying both u and the set of discontinuities v. You may set parameter $\alpha$, the range for parameter $\varepsilon$, and either a value or a range of values for parameter $\lambda$. In the latter case, computations are performed in descending $\lambda$ order, delineating progressively more and more discontinuities. All options are described in pages associated to these tools. Note also that several generic classes have been created to facilitate the development of further image tools that use a discrete exterior calculus formulation.

For further works, we plan to develop an on-line demo of these image restoration methods, which uses the IPOL (Image Processing On-Line) platform. This will facilitate even more the reproducibility of our research as well as fair comparisons with other restoration methods.