Solution of Ambrosio–Tortorelli model for image segmentation by generalized relaxation method

doi:10.1016/j.cnsns.2014.06.036

Communications in Nonlinear Science and Numerical Simulation

Volume 20, Issue 3, March 2015, Pages 819-831

https://doi.org/10.1016/j.cnsns.2014.06.036 Get rights and content

Highlights

•
We propose a new numerical approach for solving an image segmentation model.
•
The Ambrosio–Tortorelli approximation of the Mumford–Shah model is considered.
•
Finite-difference discretization of the Euler–Lagrange equations is proposed.
•
Non-linear Gauss–Seidel coupled with linear geometric multigrid is used.
•
Experiments show that the method is efficient for increasing image sizes.

Abstract

Image segmentation addresses the problem to partition a given image into its constituent objects and then to identify the boundaries of the objects. This problem can be formulated in terms of a variational model aimed to find optimal approximations of a bounded function by piecewise-smooth functions, minimizing a given functional. The corresponding Euler–Lagrange equations are a set of two coupled elliptic partial differential equations with varying coefficients. Numerical solution of the above system often relies on alternating minimization techniques involving descent methods coupled with explicit or semi-implicit finite-difference discretization schemes, which are slowly convergent and poorly scalable with respect to image size. In this work we focus on generalized relaxation methods also coupled with multigrid linear solvers, when a finite-difference discretization is applied to the Euler–Lagrange equations of Ambrosio–Tortorelli model. We show that non-linear Gauss–Seidel, accelerated by inner linear iterations, is an effective method for large-scale image analysis as those arising from high-throughput screening platforms for stem cells targeted differentiation, where one of the main goal is segmentation of thousand of images to analyze cell colonies morphology.

Introduction

Image segmentation is one of the central tasks in image processing and has the aim of partitioning an image into its constituent regions or objects, leading to also identify the boundaries of these objects. Usually, two types of techniques are used for image segmentation: the discontinuity-based approach, where main aim is to identify the set of pixels in which intensity function has high gradients or jumps (the edges) or the similarity-based approach, where the aim is to identify regions with similar characteristics in intensity function. The problem can be formulated in terms of a variational model, i.e., in terms of an energy minimization criterion, which in some sense merges features of both the above approaches. The original variational model for segmentation is a free-discontinuity problem formulated by Mumford and Shah [17] that proposed to look for a piecewise smooth approximation u of the original image function f, with u discontinuous across a closed set K included in the image domain $Ω$ . In more details, let $Ω \subset R^{2}$ be a bounded open set and $f \in L^{\infty} (Ω)$ the observed gray-level image, the problem consists in the minimization of the following functional: $E (u, K) = \int_{Ω} {(u - f)}^{2} dxdy + β \int_{Ω ⧹ K} | \nabla u |^{2} dxdy + α | K |,$ where $u \in C^{1} (Ω ⧹ K), K \subset Ω$ is a closed set whose $| K |$ is the length of K, and $α$ and $β$ are positive coefficients.

This problem is not simple to deal with due to the presence of the last term. However, an extensive theory has been developed since its introduction and many efforts have been devoted to approximate the problem with a relaxed one which can be solved by classical approach of Calculus of Variations. On the other hand, since the functional is not convex, large efforts have been also devoted to relax the model in order to get convexity. It is beyond the focus of this work to discuss theoretical aspects or advantages and drawbacks of the Mumford–Shah (MS) model, however, we can observe that it is the most general way to formulate a segmentation problem and the well-based theory developed during the last 20 years motivates its use to obtain robust and accurate software modules for large-scale analysis. For classical books and recent reviews on theory, numerical approximations and applications of the MS model we refer to [1], [2], [6], [16], [22].

One of the most general approximations of the MS model was proposed by Ambrosio and Tortorelli who showed that the model can be approximated, in the sense of $Γ$ -convergence, by a sequence of elliptic functionals [3], [4]. For details on $Γ$ -convergence and its role in the study of asymptotic variational problems we refer to [9]. Ambrosio and Tortorelli introduced a new variable $0 ⩽ z ⩽ 1$ , which controls $| \nabla u |$ and gives an approximate representation of the set K in a tubular neighbourhood of radius $∊$ of the minimizer K. The possibility to approximate the measure of K by an elliptic functional depending on z leads to the following sequence of functionals depending on $∊$ : $E_{∊} (u, z) = \int_{Ω} {(u - f)}^{2} dxdy + β \int_{Ω} z^{2} | \nabla u |^{2} dxdy + α \int_{Ω} (∊ | \nabla z |^{2} + \frac{{(z - 1)}^{2}}{4 ∊}) dxdy,$ which $Γ$ -converges to $E (u, K)$ for $∊ \to 0$ . Note that $E_{∊} (u, z)$ remains not convex and a drawback of this approximation is its dependence on the choice of a good $∊$ parameter in numerical solution. In [3], [21] the reader can find a clear discussion, which we do not report here for sake of brevity, on the behaviour of the function z introduced by Ambrosio and Tortorelli in their work. However, we notice that this function is forced to tend to 1 as $∊$ tends to 0 because the term ${(z - 1)}^{2}$ is positive and vanishes only for $z = 1$ ; furthermore, the term $z^{2}$ must go to zero near the discontinuities of u, where $| \nabla u |$ is large, to keep bounded the second term in (2). Therefore, z makes a sharp transition to zero inside a neighbourhood of the set K of thickness of radius $∊$ and an accurate numerical solution should require a suitable spatial discretization of the model to well resolve the above tubular region.

Some efforts have been devoted to minimize the functional in (2) by finite-element approximations, see for example [7], [8]. On the other hand, the form of the functional in (2) allows us to apply the classical approach of Calculus of Variations, i.e. writing Euler–Lagrange equations to obtain stationary solutions. Numerical solutions of the Euler–Lagrange equations associated to (2) are usually based on finite-difference discretizations of the equations and alternating minimization schemes are obtained by applying gradient descent iterations [6], [19], [22]. In [21] the authors proposed to use the non-linear Gauss–Seidel method for solving discrete equations arising from a finite-difference approximation of (2). In this paper we show that non-linear Gauss–Seidel relaxation, coupled with a fixed-point approach producing inner linear systems at each application of a basic step of the non-linear method to a discrete form of the equations, is very effective for large-scale image segmentation. Inner iterations largely reduce the number of non-linear iterations already for very low accuracy requests on inner solutions but also improves robustness when those accuracy requests are increased, leading to reliable, flexible and efficient solver. We discuss convergence results and computational cost of the method in the analysis of 2D gray-scale images of embryonic stem cells colonies, since our final aim is to develop an efficient segmentation module for automatic screening of colonies morphology useful to discriminate different phenotypic transitions.

The paper is organized as follows. In Section 2 we introduce a finite-difference discretization of the Euler–Lagrange equations for model (2). In Section 3 we briefly describe the non-linear Gauss–Seidel method when applied to a discrete form of (2), then we describe our fixed-point approach coupled with the non-linear relaxation and outline the linear multigrid solver which can be used to accelerate inner convergence. In Section 4 we present results obtained on real images of cultures of mouse embryonic stem cells. An extensive discussion on the performance of our method, in terms of robustness, accuracy, computational costs and comparison with standard gradient descent methods, varying model parameters, is reported. Concluding remarks and future work are included in Section 5.

Section snippets

Finite-difference discretization of Euler–Lagrange equations

Euler–Lagrange equations for Ambrosio–Tortorelli model are the following system of two coupled elliptic PDEs, associated with Neumann boundary conditions: $\{\begin{matrix} 2 (u - f) - 2 β \nabla \cdot (z^{2} \nabla u) = 0 \\ 2 β z | \nabla u |^{2} - 2 α ∊ \nabla^{2} z + \frac{α}{2 ∊} (z - 1) = 0 \end{matrix} (x, y) \in Ω,$ $\{\begin{matrix} \nabla u \cdot \vec{n} = 0 \\ \nabla z \cdot \vec{n} = 0 \end{matrix} (x, y) \in \partial Ω,$ where $\vec{n}$ is the exterior normal to $\partial Ω$ .

We note that the first equation in (3) is a form of anisotropic diffusion where function $z^{2}$ represents anisotropy and observe that each one of the equations is linear in one of the unknowns if the other is fixed. Main approaches

Numerical algorithms

The non-linear system (5) can be solved by a generalized linear method such as that arising from basic point-wise iterative methods, which have a simple interpretation if we consider the system in the form (6). For example, the basic step of non-linear Gauss–Seidel, starting from a vector $(u^{k}, z^{k})$ , requires the application of a basic step of the linear Gauss–Seidel relaxation to the system: $A (z^{k}) u = f_{1},$ in order to obtain a new approximation $u^{k + 1}$ , and then a basic step of the linear Gauss–Seidel

Numerical results

In the following we discuss results related to the application of the algorithms described in the above section to gray-scale real images of plates of cultured mouse embryonic stem cells, whose original size is $1040 \times 1392$ .

First of all we present simulation results on an image with size ${nsize}^{2} = 1024^{2}$ , for varying model coefficients, in order to show their impact on the results’ quality. Since in case of convergence, as expected, there is no significative impact of the chosen numerical method on

Concluding remarks

In this work we focus on efficient and reliable numerical solution of the Ambrosio–Tortorelli model for image segmentation. The main motivation was the need to develop an effective software tool for large-scale analysis of images coming from high-throughput screening platforms for stem cells targeted differentiation. Due to its generality and well-based theory, we choose the phase-field approximation of the original Mumford–Shah variational model, developed by Ambrosio and Tortorelli, and we

Acknowledgements

We would like to thank Dr. Laura Casalino of the Stem Cell Fate Laboratory at the Institute of Genetics and Biophysics (IGB) of the National Research Council (CNR) of Italy for images used in our analysis. We also thank Dr. Riccardo March for useful discussions on the theoretical aspects of the Ambrosio–Tortorelli model and the practical choice of the model parameters. Finally, we would like to thank the anonymous reviewers for their useful comments, which helped us to improve the paper.

References (23)

R. March et al.
A variational method for the recovery of smooth boundaries
Image Vision Comput
(1997)
A. Vitti
The Mumford–Shah variational model for image segmentation: an overview of the theory, implementation and use
ISPRS J Photogrammetry Remote Sens
(2012)
L. Ambrosio et al.
Functions of bounded variations and free discontinuity problems
(2000)
G. Aubert et al.
Mathematical problems in image processing. Partial differential equations and the calculus of variations
(2002)
L. Ambrosio et al.
Approximation of functional depending on jumps by elliptic functionals via $γ$ -convergence
Commun Pure Appl Math
(1990)
L. Ambrosio et al.
On the approximation of free discontinuity problem by elliptic functionals
Boll Un Mat Ital B
(1992)
L. Bar et al.
Generalized Newton-type methods for energy formulations in image processing
SIAM J Imaging Sci
(2009)
L. Bar
Mumford and Shah model and its applications to image segmentation and image restoration
G. Bellettini et al.
Discrete approximation of a free discontinuity problem
Numer Funct Anal Optim
(1994)
B. Bourdin
Image segmentation with a finite element method
Math Model Numer Anal
(1999)

A. Braides

$Γ$ -convergence for beginners

(2002)

Cited by (12)

Mesh adaptation-aided image segmentation
2019, Communications in Nonlinear Science and Numerical Simulation
Citation Excerpt :
For this reason, we focus on mesh adaptation, as it both increases the resolution of the edges and greatly reduces the computational cost, in the spirit of a mesh adaptation-aided image segmentation. In the context of finite differences, in [16], the authors consider generalized relaxation methods coupled with multigrid linear solvers applied to the Euler-Lagrange equations of the Ambrosio-Tortorelli model. To the best of our knowledge, no other paper in the literature has addressed anisotropic mesh adaptation in the context of image segmentation based on the Ambrosio-Tortorelli variational framework.
We focus on a variational approach to image segmentation based on the Ambrosio-Tortorelli functional. To make the procedure more effective with respect to standard algorithms, we combine the functional minimization with the employment of an optimal discretization. More precisely, we perform a finite element approximation of the Ambrosio-Tortorelli functional on a triangular adapted mesh able to follow exactly the contours present in the images, in the spirit of a mesh adaptation-aided image segmentation. This challenging goal is reached via a rigorous a posteriorierror analysis enriched with anisotropic information. The benefits due to the proposed algorithm are evident both in terms of increased resolution in the edge detection and in a considerable reduction of the computational costs, as confirmed by an extensive numerical investigation.
Serial and parallel approaches for image segmentation by numerical minimization of a second-order functional
2018, Applied Mathematics and Computation
Because of its attractive features, second order segmentation has shown to be a promising tool in remote sensing. A known drawback about its implementation is computational complexity, above all for large set of data. Recently in Zanetti et al. [1], an efficient version of the block-coordinate descent algorithm (BCDA) has been proposed for the minimization of a second order elliptic approximation of the Blake–Zissermann functional. Although the parallelization of linear algebra operations is expected to increase the performance of BCDA when addressing the segmentation of large-size gridded data (e.g., full-scene images or Digital Surface Models (DSMs)), numerical evidence shows that this is not sufficient to get significant reduction of computational time. Therefore a novel approach is proposed which exploits a decomposition technique of the image domain into tiles. The solution can be computed by applying BCDA on each tile in parallel way and combining the partial results corresponding to the different blocks of variables through a proper interconnection rule. We prove that this parallel method (OPARBCDA) generates a sequence of iterates which converges to a critical point of the functional on the level set devised by the starting point. Furthermore, we show that the parallel method can be efficiently implemented even in a commodity multicore CPU. Numerical results are provided to evaluate the efficiency of the parallel scheme on large images in terms of computational cost and its effectiveness with respect to the behavior on the tile junctions.
A parallel generalized relaxation method for high-performance image segmentation on GPUS
2016, Journal of Computational and Applied Mathematics
Citation Excerpt :
In Section 4 we discuss performance results and, finally, in Section 5 we give some remarks and future plans. In [5] we propose to solve system (3) by using a generalized relaxation method, such as the Gauss–Seidel method, accelerated by inner linear iterations. We proved that this approach is more efficient and reliable when compared both with standard non-linear Gauss–Seidel iterations and with other first-order alternating minimization schemes, such as gradient descent.
Fast and scalable software modules for image segmentation are needed for modern high-throughput screening platforms in Computational Biology. Indeed, accurate segmentation is one of the main steps to be applied in a basic software pipeline aimed to extract accurate measurements from a large amount of images. Image segmentation is often formulated through a variational principle, where the solution is the minimum of a suitable functional, as in the case of the Ambrosio–Tortorelli model. Euler–Lagrange equations associated with the above model are a system of two coupled elliptic partial differential equations whose finite-difference discretization can be efficiently solved by a generalized relaxation method, such as Jacobi or Gauss–Seidel, corresponding to a first-order alternating minimization scheme. In this work we present a parallel software module for image segmentation based on the Parallel Sparse Basic Linear Algebra Subprograms (PSBLAS), a general-purpose library for parallel sparse matrix computations, using its Graphics Processing Unit (GPU) extensions that allow us to exploit in a simple and transparent way the performance capabilities of both multi-core CPUs and of many-core GPUs. We discuss performance results in terms of execution times and speed-up of the segmentation module running on GPU as well as on multi-core CPUs, in the analysis of 2D gray-scale images of mouse embryonic stem cells colonies coming from biological experiments.
Numerical minimization of a second-order functional for image segmentation
2016, Communications in Nonlinear Science and Numerical Simulation
Citation Excerpt :
In their approximation the discontinuity set is replaced by an auxiliary function that plays the role of indicator function. Numerical solutions based on the Ambrosio–Tortorelli approximation are given in the framework of Finite Element Method (FEM) in [13], and via finite-difference discretization of Euler–Lagrange equations in [14]. In [15], a Γ-convergence approximation using local integral functionals defined on a discrete space is given.
In this paper we address the numerical minimization of a variational approximation of the Blake–Zisserman functional given by Ambrosio, Faina and March. Our approach exploits a compact matricial formulation of the objective functional and its decomposition into quadratic sparse convex sub-problems. This structure is well suited for using a block-coordinate descent method that cyclically determines a descent direction with respect to a block of variables by few iterations of a preconditioned conjugate gradient algorithm. We prove that the computed search directions are gradient related and, with convenient step-sizes, we obtain that any limit point of the generated sequence is a stationary point of the objective functional. An extensive experimentation on different datasets including real and synthetic images and digital surface models, enables us to conclude that: (1) the numerical method has satisfying performance in terms of accuracy and computational time; (2) a minimizer of the proposed discrete functional preserves the expected good geometrical properties of the Blake–Zisserman functional, i.e., it is able to detect first and second order edge-boundaries in images and (3) the method allows the segmentation of large images.
Preface to the special issue NUMTA 2013
2015, Communications in Nonlinear Science and Numerical Simulation
On an effective multigrid solver for solving a class of variational problems with application to image segmentation
2020, International Journal of Computer Mathematics

View all citing articles on Scopus

^☆: This work has been partially supported by Italian Flagship Project InterOmics and INdAM-GNCS Project on Numerical Methods and Software for Large-Scale Optimization with Applications to Image Processing.

View full text

Solution of Ambrosio–Tortorelli model for image segmentation by generalized relaxation method☆

Highlights

Abstract

Introduction

Section snippets

Finite-difference discretization of Euler–Lagrange equations

Numerical algorithms

Numerical results

Concluding remarks

Acknowledgements

Image Vision Comput

ISPRS J Photogrammetry Remote Sens

Functions of bounded variations and free discontinuity problems

Mathematical problems in image processing. Partial differential equations and the calculus of variations

Approximation of functional depending on jumps by elliptic functionals via γ-convergence

Commun Pure Appl Math

On the approximation of free discontinuity problem by elliptic functionals

Boll Un Mat Ital B

Generalized Newton-type methods for energy formulations in image processing

SIAM J Imaging Sci

Mumford and Shah model and its applications to image segmentation and image restoration

Discrete approximation of a free discontinuity problem

Numer Funct Anal Optim

Image segmentation with a finite element method

Math Model Numer Anal

Γ-convergence for beginners

Approximation of functional depending on jumps by elliptic functionals via $γ$ -convergence

$Γ$ -convergence for beginners