A Convex Geodesic Selective Model for Image Segmentation

Michael Roberts^†, Ke Chen^† and Klaus L. Irion^‡
^†Centre for Mathematical Imaging Techniques and
Department of Mathematical Sciences,
University of Liverpool, United Kingdom
^‡Department of Radiology, Manchester University NHS Foundation Trust, Manchester, UK
Email k.chen@liverpool.ac.uk, Web: www.liv.ac.uk/cmit (corresponding author). Work supported by UK EPSRC grant EP/N014499/1.

Abstract

Selective segmentation is an important application of image processing. In contrast to global segmentation in which all objects are segmented, selective segmentation is used to isolate specific objects in an image and is of particular interest in medical imaging – permitting segmentation and review of a single organ. An important consideration is to minimise the amount of user input to obtain the segmentation; this differs from interactive segmentation in which more user input is allowed than selective segmentation. To achieve selection, we propose a selective segmentation model which uses the edge-weighted geodesic distance from a marker set as a penalty term. It is demonstrated that this edge-weighted geodesic penalty term improves on previous selective penalty terms. A convex formulation of the model is also presented, allowing arbitrary initialisation. It is shown that the proposed model is less parameter dependent and requires less user input than previous models. Further modifications are made to the edge-weighted geodesic distance term to ensure segmentation robustness to noise and blur. We can show that the overall Euler-Lagrange equation admits a unique viscosity solution. Numerical results show that the result is robust to user input and permits selective segmentations that are not possible with other models.

Keywords. Variational model, partial differential equations, image segmentation, additive operator splitting, viscosity solution, geodesic.

1. Introduction

Segmentation of an image into its individual objects is one incredibly important application of image processing techniques. Segmentation can take two forms; firstly global segmentation for isolation of all foreground objects in an image from the background and secondly, selective segmentation for isolation of a subset of the objects in an image from the background. A comprehensive review of selective segmentation can be found in [7, 19] and in [45] for medical image segmentation where selection refers to extraction of single organs.

Approaches to image segmentation broadly fall into two classes; region-based and edge-based. Some region-based approaches are region growing [1], watershed algorithms [40], Mumford-Shah [29] and Chan-Vese [11]. The final two of these are partial differential equations (PDEs)-based variational approaches to the problem of segmentation. There are also models which mix the two classes to use the benefits of the region-based and edge-based approaches and will incorporate features of each. Edge-based methods aim to encourage an evolving contour towards the edges in an image and normally require an edge detector function [8]. The first edge-based variational approach was devised by Kass et al. [22] with the famous snakes model, this was further developed by Casselles et al. [8] who introduced the Geodesic Active Contour (GAC) model. Region-based global segmentation models include the well known works of Mumford-Shah [29] and Chan-Vese [11]. Importantly they are non-convex and hence a minimiser of these models may only be a local, not the global, minimum. Further work by Chan et al. [10] gave rise to a method to find the global minimiser for the Chan-Vese model under certain conditions.

This paper is mainly concerned with selective segmentation of objects in an image, given a set of points near the object or objects to be segmented. It builds in such user input to a model using a set $\mathcal{M}=\{(x_{i},y_{i})\in\Omega,1\leq i\leq k\}$ where $\Omega\subset\mathbb{R}^{2}$ is the image domain [4, 5, 17]. Nguyen et al. [30] considered marker sets $\mathcal{M}$ and $\mathcal{A}$ which consist of points inside and outside, respectively, the object or objects to be segmented. Gout et al. [17] combined the GAC approach with the geometrical constraint that the contour passes through the points of $\mathcal{M}$ . This was enforced with a distance function which is zero at $\mathcal{M}$ and non-zero elsewhere. Badshah and Chen [4] then combined the Gout et al. model with [11] to incorporate a constraint on the intensity in the selected region, thereby encouraging the contour to segment homogeneouss regions. Rada and Chen [36] introduced a selective segmentation method based on two-level sets which was shown to be more robust than the Badshah-Chen model. We also refer to [5, 23] for selective segmentation models which include different fitting constraints, using coefficient of variation and the centroid of $\mathcal{M}$ respectively. None of these models have a restriction on the size of the object or objects to be detected and depending on the initialisation these methods have the potential to detect more or fewer objects than the user desired. To address this and to improve on [36], Rada and Chen [37] introduced a model combining the Badshah-Chen [4] model with a constraint on the area of the objects to be segmented. The reference area used to constrain the area within the contour is that of the polygon formed by the markers in $\mathcal{M}$ . Spencer and Chen [39] introduced a model with the distance fitting penalty as a standalone term in the energy functional, unbounding it from the edge detector term of the Gout et al. model.

All of the above selective segmentation models discussed are non-convex and hence the final result depends on the initialisation. Spencer and Chen [39], in the same paper, reformulated the model they introduced to a convex form using convex relaxation and an exact penalty term as in [10]. Their model uses Euclidean distance from the marker set $\mathcal{M}$ as a distance penalty term, however we propose replacing this with the edge-weighted geodesic distance from $\mathcal{M}$ (we call this simply the geodesic distance). This distance increases at edges in the image and is more intuitive for selective segmentation. The proposed model is given as a convex relaxed model with exact penalty term and we give a general existence and uniqueness proof for the viscosity solution to the PDE given by its Euler-Lagrange equation, which is also applicable to a whole class of PDEs arising in image segmentation. We note that the use of geodesic distance for segmentation has been considered before [6, 34], however the models only use geodesic distance as the fitting term within the regulariser, so are liable to make segmentation errors for poor initialisation or complex images. Here we take a different approach, by including geodesic distance as a standalone fitting term, separate from the regulariser, and using intensity fitting terms to ensure robustness.

In this paper we only consider 2D images, however for completion we remark that 3D segmentation models do exist [25, 44] and it is simple to extend the proposed model to 3D. The contributions of this paper can be summarised as follows:

•

We incorporate the geodesic distance as a distance penalty term within the variational framework.
•

We propose a convex selective segmentation model using this penalty term and demonstrate how it can achieve results which cannot be achieved by other models.
•

We improve the geodesic penalty term, focussing on improving robustness to noise and improving segmentation when object edges are blurred.
•

We give an existence and uniqueness proof for the viscosity solution for the PDEs associated with a whole class of segmentation models (both global and selective).

We find that the proposed model gives accurate segmentation results for a wide range of parameters and, in particular, when segmenting the same objects from the same modality images, i.e. segmenting lungs from CT scans, the parameters are very similar from one image to the next to obtain accurate results. Therefore, this model may be used to assist the preparation of large training sets for deep learning studies [32, 41, 42] that concern segmentation of particular objects from images.

The paper is structured as follows; in §2 we review some global and selective segmentation models. In §3 we discuss the geodesic distance penalty term, propose a new convex model and also address weaknesses in the naïve implementation of the geodesic distance term. In §4 we discuss the non-standard AOS scheme, introduced in [39], which we use to solve the model. In §5 we give an existence and uniqueness proof for a general class of PDEs arising in image segmentation, thereby showing that for a given initialisation the solution to our model is unique. In §6 we compare the results of the proposed model to other selective segmentation models, show that the proposed model is less parameter dependent than other models and is more robust to user input. Finally, in §7 we provide some concluding remarks.

2. Review of Variational Segmentation Models

Although we focus on selective segmentation, it is illuminating to introduce some global segmentation models first. Throughout this paper we denote the original image by $z(x,y)$ with image domain $\Omega\subset\mathbb{R}^{2}$ .

2.1. Global Segmentation

The model of Mumford and Shah [29] is one of the most famous and important variational models in image segmentation. We will review its two-dimensional piecewise constant variant, commonly known as the Chan-Vese model [11], which takes the form

F_{CV}(\Gamma,c_{1},c_{2})=\mu\cdot length(\Gamma)+\lambda_{1}\int_{\Omega_{1}}|z(x,y)-c_{1}|^{2}\,\mathrm{d}\Omega+\lambda_{2}\int_{\Omega_{2}}|z(x,y)-c_{2}|^{2}\,\mathrm{d}\Omega

(1)

where the foreground $\Omega_{1}$ is the subdomain to be segmented, the background $\Omega_{2}=\Omega\backslash\Omega_{1}$ and $\mu,\lambda_{1},\lambda_{2}$ are fixed non-negative parameters. The values $c_{1}$ and $c_{2}$ are the average intensities of $z(x,y)$ inside $\Omega_{1}$ and $\Omega_{2}$ respectively. We use a level set function

\phi(x,y)=\begin{cases}>0,&(x,y)\in\Omega_{1},\\ 0,&(x,y)\in\Gamma,\\ <0,&otherwise,\\ \end{cases}

to track $\Gamma=\{(x,y)\in\Omega\,|\,\phi(x,y)=0\}$ (an idea popularised by Osher and Sethian [31]) and reformulate (1) as

\begin{split}F_{CV}(\phi,c_{1},c_{2})=&\mu\int_{\Omega}|\nabla H_{\varepsilon}(\phi)|\,\mathrm{d}\Omega+\lambda_{1}\int_{\Omega}(z(x,y)-c_{1})^{2}H_{\varepsilon}(\phi)\,\mathrm{d}\Omega\\ &\ \hskip 76.82234pt+\lambda_{2}\int_{\Omega}(z(x,y)-c_{2})^{2}(1-H_{\varepsilon}(\phi))\,\mathrm{d}\Omega,\\ \end{split}

(2)

with $H_{\varepsilon}(\phi)$ a smoothed Heaviside function such as $H_{\varepsilon}(\phi)=\frac{1}{2}+\frac{1}{\pi}\arctan(\frac{\phi}{\varepsilon})$ for some $\varepsilon$ , we set $\varepsilon=1$ throughout. We solve this in two stages, first with $\phi$ fixed we minimise $F_{CV}$ with respect to $c_{1}$ and $c_{2}$ , obtaining

c_{1}=\frac{\int_{\Omega}H_{\varepsilon}(\phi)\cdot z(x,y)\,\mathrm{d}\Omega}{\int_{\Omega}H_{\varepsilon}(\phi)\,\mathrm{d}\Omega},\hskip 36.135ptc_{2}=\frac{\int_{\Omega}(1-H_{\varepsilon}(\phi))\cdot z(x,y)\,\mathrm{d}\Omega}{\int_{\Omega}(1-H_{\varepsilon}(\phi))\,\mathrm{d}\Omega},

(3)

and secondly, with $c_{1}$ and $c_{2}$ fixed we minimise (2) with respect to $\phi$ . This requires the calculation of the associated Euler-Lagrange equations. A drawback of the Chan-Vese energy functional (2) is that it is non-convex. Therefore a minimiser may only be a local minimum and not the global minimum and the final segmentation result is dependent on the initialisation. Chan et al. [10] reformulated (2) using an exact penalty term to obtain an equivalent convex model – we use this same technique in §2.2 for the Geodesic Model.

2.2. Selective Segmentation

Selective segmentation models make use of user input, i.e. a marker set $\mathcal{M}$ of points near the object or objects to be segmented. Let $\mathcal{M}=\{(x_{i},y_{i})\in\Omega,1\leq i\leq k\}$ be such a marker set. The aim of selective segmentation is to design an energy functional where the segmentation contour $\Gamma$ is close to the points of $\mathcal{M}$ .

Early work. An early model by Caselles et al. [8], commonly known as the Geodesic Active Contour (GAC) model, uses an edge detector function to ensure the contour follows edges, the functional to minimise is given by

\bigintsss_{\Gamma}g(|\nabla z(x,y)|)d\Gamma.

The term $g(|\nabla z(x,y)|)$ is an edge detector, one example is $g(s)=1/(1+\beta s^{2})$ with $\beta$ a tuning parameter. It is common to smooth the image with a Gaussian filter $G_{\sigma}$ where $\sigma$ is the kernel size, i.e. use $g(|\nabla\left(G_{\sigma}*z(x,y)\right)|)$ as the edge detector. This mitigates the effect of noise in the image, giving a more accurate edge detector. Gout et al. [25] built upon the GAC model by incorporating a distance term $\mathcal{D}(x,y)$ into this integral, i.e. the integrand is $\mathcal{D}(x,y)g(|\nabla z|)$ . The distance term is a penalty on the distance from $\mathcal{M}$ , this model encourages the contour to be near to the set $\mathcal{M}$ whilst also lying on edges. However this model struggles when boundaries between objects and their background are fuzzy or blurred. To address this, Badshah and Chen [4] introduced a new model which adds the intensity fitting terms from the Chan-Vese model (1) to the Gout et al. model. However, their model has poor robustness [36]. To improve on this, Rada and Chen [37] introduced a model which adds an area fitting term into the Badshah-Chen model and is far more robust.

The Rada-Chen model [37]. We first briefly introduce this model, defined by

\begin{split}F_{RC}(\phi,c_{1},c_{2})=&\mu\int_{\Omega}\mathcal{D}(x,y)g(|\nabla z(x,y)|)|\nabla H_{\varepsilon}(\phi)|\,\mathrm{d}\Omega\\ &+\lambda_{1}\int_{\Omega}(z(x,y)-c_{1})^{2}H_{\varepsilon}(\phi)\,\mathrm{d}\Omega+\lambda_{2}\int_{\Omega}(z(x,y)-c_{2})^{2}(1-H_{\varepsilon}(\phi))\,\mathrm{d}\Omega\\ &+\gamma\bigg{[}\left(\int_{\Omega}H_{\varepsilon}(\phi)\,\mathrm{d}\Omega-A_{1}\right)^{2}+\left(\int_{\Omega}(1-H_{\varepsilon}(\phi))\,\mathrm{d}\Omega-A_{2}\right)^{2}\bigg{]},\end{split}

(4)

where $\mu,\lambda_{1},\lambda_{2},\gamma$ are fixed non-negative parameters. There is freedom in choosing the distance term $\mathcal{D}(x,y)$ , see [37] for some examples. $A_{1}$ is the area of the polygon formed from the points of $\mathcal{M}$ and $A_{2}=|\Omega|-A_{1}$ . The final term of this functional puts a penalty on the area inside a contour being very different to $A_{1}$ . One drawback of the Rada-Chen model is that the selective fitting term uses no location information from the marker set $\mathcal{M}$ . Therefore the result can be a contour which is separated over the domain into small parts, whose sum area totals the area fitting term.

Nguyen et al. [30]. This model is based on the GAC model and uses likelihood functions as fitting terms, it has the energy functional

\begin{split}F_{NG}(\phi)=&\mu\int_{\Omega}g(|\nabla z(x,y)|)|\nabla H_{\varepsilon}(\phi)|\,\mathrm{d}\Omega\\ &+\lambda\int_{\Omega}\alpha\left(P_{B}(x,y)-P_{F}(x,y)\right)+\left(1-\alpha\right)\left(1-2P(x,y)\right)\phi\,\mathrm{d}\Omega\end{split}

where $P_{B}(x,y)$ and $P_{F}(x,y)$ are the normalised log-likelihoods that the pixel $(x,y)$ is in the foreground and background respectively. $P(x,y)$ is the probability that pixel $(x,y)$ belongs to the foreground, $\alpha\in[0,1]$ and minimisation is constrained, requiring $\phi\in[0,1]$ , so $F_{NG}(\phi)$ is convex. This model is good for many examples, see [30], however fails when the boundary of the object to segment is non-smooth or has fine structures. Also, the final result is sometimes sensitive to the marker sets used.

The Spencer-Chen model [39]. The authors introduced the following model

\begin{split}F_{SC}(\phi,c_{1},c_{2})=&\mu\int_{\Omega}g(|\nabla z(x,y)|)|\nabla H_{\varepsilon}(\phi)|\,\mathrm{d}\Omega+\lambda_{1}\int_{\Omega}(z(x,y)-c_{1})^{2}H_{\varepsilon}(\phi)\,\mathrm{d}\Omega\\ &+\lambda_{2}\int_{\Omega}(z(x,y)-c_{2})^{2}(1-H_{\varepsilon}(\phi))\,\mathrm{d}\Omega+\theta\int_{\Omega}\mathcal{D}_{E}(x,y)H_{\varepsilon}(\phi)\,\mathrm{d}\Omega,\end{split}

(5)

where $\mu,\lambda_{1},\lambda_{2},\theta$ are fixed non-negative parameters. Note that the regulariser of this model differs from the Rada-Chen model (4) as the distance function $\mathcal{D}(x,y)$ has been separated from the edge detector term and is now a standalone penalty term $\mathcal{D}_{E}(x,y)$ . The authors use normalised Euclidean distance $\mathcal{D}_{E}(x,y)$ from the marker set $\mathcal{M}$ as their distance penalty term. We will discuss this later in §3 as it is one of the key improvements we make to the Spencer-Chen model, replacing the Euclidean distance term with a geodesic distance term.

Convex Spencer-Chen model [39]. Spencer and Chen use the ideas of [10] to reformulate (5) into a convex minimisation problem. It can be shown that the Euler-Lagrange equations for $F_{SC}(\phi,c_{1},c_{2})$ have the same stationary solutions as for

\begin{split}F_{SC1}(u,c_{1},c_{2})=&\mu\int_{\Omega}g(|\nabla z(x,y)|)|\nabla u|\,\mathrm{d}\Omega+\int_{\Omega}\left[\lambda_{1}(z(x,y)-c_{1})^{2}-\lambda_{2}(z(x,y)-c_{2})^{2}\right]u\,\mathrm{d}\Omega\\ &+\theta\int_{\Omega}\mathcal{D}_{E}(x,y)u\,\mathrm{d}\Omega,\end{split}

(6)

with the minimisation constrained to $u\in[0,1]$ . This is a constrained convex minimisation which can be reformulated to an unconstrained minimisation using an exact penalty term $\nu(u):=\max\{0,2|u-\frac{1}{2}|-1\}$ in the functional, which encourages the minimiser to be in the range $[0,1]$ . In [39] the authors use a smooth approximation $\nu_{\varepsilon}(u)$ to $\nu(u)$ given by

\nu_{\varepsilon}(u)=H_{\varepsilon}\left(\sqrt{(2u-1)^{2}+\varepsilon}-1\right)\left[\sqrt{(2u-1)^{2}+\varepsilon}-1\right],

(7)

and perform the unconstrained minimisation of

\begin{split}F_{SC2}(u,c_{1},c_{2})=&\mu\int_{\Omega}g(|\nabla z(x,y)|)|\nabla u|\,\mathrm{d}\Omega+\int_{\Omega}\left[\lambda_{1}(z(x,y)-c_{1})^{2}-\lambda_{2}(z(x,y)-c_{2})^{2}\right]u\,\mathrm{d}\Omega\\ &+\theta\int_{\Omega}\mathcal{D}_{E}(x,y)u\,\mathrm{d}\Omega+\alpha\int_{\Omega}\nu_{\varepsilon}(u)\,\mathrm{d}\Omega.\end{split}

(8)

When $\alpha>\frac{1}{2}\left|\left|\left[\lambda_{1}(z(x,y)-c_{1})^{2}-\lambda_{2}(z(x,y)-c_{2})^{2}\right]+\theta\mathcal{D}_{E}(x,y)\right|\right|_{L^{\infty}}$ , the above functional has the same set of stationary solutions as $F_{SC1}(u,c_{1},c_{2})$ . It permits us to choose arbitrary $u$ initialisation to obtain the desired selective segmentation result due to its complexity.

Convex Liu et al. model [26]. Recently, a convex model was introduced by Liu et al. which applies a weighting to the data fitting terms, the functional to minimise is given by

\begin{split}F_{LIU}(u)=&\mu\int_{\Omega}|\nabla u|\,\mathrm{d}\Omega+\mu_{2}\int_{\Omega}|\nabla u|^{2}\,\mathrm{d}\Omega+\lambda\int_{\Omega}\omega^{2}(x,y)\left|z-u\right|^{2}\,\mathrm{d}\Omega,\end{split}

(9)

where $\mu,\mu_{2},\lambda$ are non-negative parameters and $\omega(x,y)=1-\mathcal{D}(x,y)g(|\nabla z|)$ where $\mathcal{D}(x,y)$ is a distance function from marker set $\mathcal{M}$ (see [26] for examples).

3. Proposed Convex Geodesic Selective Model

We propose an improved selective model, based on the Spencer-Chen model, which uses geodesic distance from the marker set $\mathcal{M}$ as the distance term, rather than the Euclidean distance. Increasing the distance when edges in the image are encountered gives a more accurate reflection of the true similarity of pixels in an image from the marker set. We propose minimising the convex functional

\begin{split}F_{CG}(u,c_{1},c_{2})=&\mu\int_{\Omega}g(|\nabla z(x,y)|)|\nabla u|\,\mathrm{d}\Omega+\int_{\Omega}\left[\lambda_{1}(z(x,y)-c_{1})^{2}-\lambda_{2}(z(x,y)-c_{2})^{2}\right]u\,\mathrm{d}\Omega\\ &+\theta\int_{\Omega}\mathcal{D}_{M}(x,y)u\,\mathrm{d}\Omega+\alpha\int_{\Omega}\nu_{\varepsilon}(u)\,\mathrm{d}\Omega,\end{split}

(10)

where $\mathcal{D}_{M}(x,y)$ is the edge-weighted geodesic distance from the marker set. In Figure 1, we compare the normalised geodesic distance and the Euclidean distance from the same marker point (i.e. set $\mathcal{M}$ has one point in it); clearly the former gives a more intuitively correct distance penalty than the latter. We will refer to this proposed model as the Geodesic Model.

Refer to caption — Figure 1: Comparison of distance measures. (i) Simple binary image with marker point; (ii) normalised Euclidean distance from marker point; (iii) edge map function $f(x)$ for the image; (iv) normalised geodesic distance from marker point.

3.1. Computing the Geodesic Distance Term $\mathcal{D}_{M}(x,y)$

The geodesic distance from the marker set $\mathcal{M}$ is given by $\mathcal{D}_{M}(x,y)=0$ for $(x,y)\in\mathcal{M}$ and $\mathcal{D}_{M}(x,y)=\frac{\mathcal{D}_{M}^{0}(x,y)}{||\mathcal{D}_{M}^{0}(x,y)||_{L^{\infty}}}$ for $(x,y)\not\in\mathcal{M}$ , where $\mathcal{D}_{M}^{0}(x,y)$ is the solution of the following PDE

|\nabla\mathcal{D}_{M}^{0}(x,y)|=f(x,y),\qquad\mathcal{D}^{0}_{M}(x_{0},y_{0})=0,\,(x_{0},y_{0})\in\mathcal{M}.

(11)

where $f(x,y)$ is defined later on with respect to the image contents.

If $f(x,y)\equiv 1$ (i.e. $|\nabla\mathcal{D}_{M}^{0}(x,y)|=1$ ) then the distance penalty $D_{M}(x,y)$ is simply the normalised Euclidean distance $\mathcal{D}_{E}(x,y)$ as used in the Spencer-Chen model (5). We have free rein to design $f(x,y)$ as we wish. Looking at the PDE in (11), we see that when $f(x,y)$ small this results in a small gradient in our distance function and it is almost flat. When $f(x,y)$ is large, we have a large gradient in our distance map. In the case of selective image segmentation, we want small gradients in homogeneous areas of the image and large gradients at edges. If we set

f(x,y)=\varepsilon_{\mathcal{D}}+\beta_{G}|\nabla z(x,y)|^{2}

(12)

this gives us the desired property that in areas where $|\nabla z(x,y)|\approx 0$ , the distance function increases by some small $\varepsilon_{\mathcal{D}}$ ; here image $z(x,y)$ is scaled to $[0,1]$ . At edges, $|\nabla z(x,y)|$ is large and the geodesic distance increases here. We set value of $\beta_{G}=1000$ and $\varepsilon_{\mathcal{D}}=10^{-3}$ throughout. In Figure 1, we see that the geodesic distance plot gives a low distance penalty on the triangle, which the marker indicates we would like segmented. There is a reasonable penalty on the background, and all other objects in the image have a very high distance penalty (as the geodesic to these points must cross two edges). This contrasts with the Euclidean distance, which gives a low distance penalty to some background pixels and maximum penalty to the pixels furthest away.

3.2. Comparing Euclidean and Geodesic Distance Terms

We briefly give some advantages of using the geodesic distance as a penalty term rather than Euclidean distance and a remark on the computational complexity for both distances.

1.

Parameter Robustness. The Geodesic Model is more robust to the choice of the fitting parameter $\theta$ , as the penalty on the inside of the shape we want segmented is consistently small. It is only outside the shape where the penalty is large. Whereas with the Euclidean distance term we always have a penalty inside the shape we actually want to segment. This is due to the nature of the Euclidean distance which does not discriminate on intensity – this penalty can also be quite high if our marker set is small and doesn’t cover the whole object.
2.

Robust to Marker Set Selection. The geodesic distance term is far more robust to point selection, for example we can choose just one point inside the object we want to segment and this will give a nearly identical geodesic distance compared to choosing many more points. This is not true of the Euclidean distance term which is very sensitive to point selection and requires markers to be spread in all areas of the object you want to segment (especially at extrema of the object).

Remark 1 (Computational Complexity.).

The main concern of using the geodesic penalty term, which we obtain by solving PDE (11), would be that it takes a significant amount of time to compute compared to the Euclidean distance. However, using the fast marching algorithm of Sethian [38], the complexity of computing $\mathcal{D}_{M}(x,y)$ is $\mathcal{O}(N\log(N))$ for an image with $N$ pixels. This is is only marginally more complex than computing the Euclidean distance which has $\mathcal{O}(N)$ complexity [28].

3.3. Improvements to Geodesic Distance Term

We now propose some modifications to the geodesic distance. Although the geodesic distance presents many advantages for selective image segmentation, we have three key disadvantages of this fitting term, which the Euclidean fitting term does not suffer.

1.

Not robust to noise. The computation of the geodesic distance depends on $|\nabla z(x,y)|^{2}$ in $f(x,y)$ (see (11)). So, if an image contains a lot of noise, each noisy pixel appears as an edge and we get a misleading distance term.
2.

Objects far from $\mathcal{M}$ with low penalty. As the geodesic distance only uses marker set $\mathcal{M}$ for its initial condition (see (11)), this can result in objects far from $\mathcal{M}$ having a low distance penalty, which is clearly not desired.
3.

Blurred edges. If we have two objects separated by a blurry edge and we have marker points only in one object, the geodesic distance will be low to the other object, as the edge penalty is weakly enforced for a blurry edge. We would desire low penalty inside the object with markers and a reasonable penalty in the joined object.

In Figure 2, each column shows an example for each of the problems listed above. We now propose solutions to each of these problems.

Problem 1: Noise Robustness. A naïve solution to the problem of noisy images would be to apply a Gaussian blur to $z(x,y)$ to remove the effect of the noise, so we change $f(x,y)$ to

\tilde{f}(x,y)=\varepsilon_{\mathcal{D}}+\beta_{G}|\nabla G_{\sigma}*z(x,y)|^{2}

(13)

where $G_{\sigma}$ is a Gaussian convolution with standard deviation $\sigma$ . However, the effect of Gaussian convolution is that it also blurs edges in the image. This then gives us the same issues described in Problem 3. We see in Figure 3 column 3, that the Gaussian convolution reduces the sharpness of edges and this results in the geodesic distance being very similar in adjacent objects – therefore we see more pixels with high geodesic distance.

Our alternative to Gaussian blur is to consider anisotropic TV denoising. We refer the reader to [9, 33] for information on the model, here we just give the PDE which results from its minimisation:

\begin{gathered}\begin{aligned} \tilde{\mu}\nabla\cdot\bigg{(}g(|\nabla z(x,y)|)\frac{\nabla u}{|\nabla u|_{\varepsilon_{2}}}\bigg{)}+\iota(z(x,y)-u)=0,\end{aligned}\end{gathered}

(14)

where $\tilde{\mu},\iota$ are non-negative parameters (we fix throughout $\tilde{\mu}=10^{-3},\iota=5\times 10^{-4}$ ). It is proposed to apply a relatively small number of cheap fixed point Gauss-Seidel iterations (between 100 and 200) to the discretised PDE. We cycle through all pixels $(i,j)$ and update $u_{i,j}$ as follows

\begin{gathered}\begin{aligned} u_{i,j}=\frac{A_{i,j}u_{i+1,j}+B_{i,j}u_{i-1,j}+C_{i,j}u_{i,j+1}+D_{i,j}u_{i,j-1}}{A_{i,j}+B_{i,j}+C_{i,j}+D_{i,j}+\iota}\end{aligned}\end{gathered}

(15)

where $A_{i,j}=\frac{\tilde{\mu}}{h_{x}^{2}}g(|\nabla z(x,y)|)_{i+1/2,j}$ , $B_{i,j}=\frac{\tilde{\mu}}{h_{x}^{2}}g(|\nabla z(x,y)|)_{i-1/2,j}$ , $C_{i,j}=\frac{\tilde{\mu}}{h_{y}^{2}}g(|\nabla z(x,y)|)_{i,j+1/2}$ and $D_{i,j}=\frac{\tilde{\mu}}{h_{y}^{2}}g(|\nabla z(x,y)|)_{i,j-1/2}$ . We update all pixels once per iteration and solve the PDE in (11) with $f(x,y)$ replaced by

{f}_{1}(x,y)=\varepsilon_{\mathcal{D}}+\beta_{G}|\nabla S^{k}(z(x,y))|^{2}

(16)

where $S$ represents the Gauss-Seidel iterative scheme and $k$ is the number of iterations performed (we choose $k=100$ in our tests). In the final column of Figure 3 we see that the geodesic distance map more closely resembles that of the clean image than the Gaussian blurred map in column 3 and in Figure 4 we see that the segmentation results are qualitatively and quantitatively better using the anisotropic smoothing technique.

Problem 2: Objects far from $\mathcal{M}$ with low penalty.

In Figure 2 column 2 we see that the geodesic distance to the outside of the patient is lower than to their ribs. This is due to the fact that the region outside the body is homogeneous and there is almost zero distance penalty in this region. Similarly for Figure 3 column $4$ , the distances from the marker set to many surrounding objects is low, even though their Euclidean distance from the marker set is high. We wish to have the Euclidean distance $\mathcal{D}_{E}(x,y)$ incorporated somehow. Our solution is to modify the term ${f}_{1}(x,y)$ from (16) to

{f}_{2}(x,y)=\varepsilon_{\mathcal{D}}+\beta_{G}|\nabla S^{k}(z(x,y))|^{2}+\vartheta\mathcal{D}_{E}(x,y).

(17)

In Figure 5 the effect of this is clear, as $\vartheta$ increases, the distance function resembles the Euclidean distance more. We use value $\vartheta=10^{-1}$ in all experiments as it adds a reasonable penalty to pixels far from the marker set.

Problem 3: Blurred edges.

If there are blurred edges between objects in an image, the geodesic distance will not increase significantly at this edge. Therefore the final segmentation result is liable to include unwanted objects. We look to address this problem through the use of anti-markers. These are markers which indicate objects that we do not want to segment, i.e. the opposite of marker points, we denote the set of anti-marker points by $\mathcal{AM}$ .

We propose to use a geodesic distance map from the set $\mathcal{AM}$ denoted by $\mathcal{D}_{AM}(x,y)$ which penalises pixels near to the set $\mathcal{AM}$ and doesn’t add any penalty to those far away. We could naïvely choose $\mathcal{D}_{AM}(x,y)=1-\tilde{\mathcal{D}}_{GAM}(x,y)$ where $\tilde{\mathcal{D}}_{GAM}(x,y)$ is the normalised geodesic distance from $\mathcal{AM}$ . However this puts a large penalty on those pixels inside the object we actually want to segment (as $\tilde{\mathcal{D}}_{GAM}(x,y)$ to those pixels is small). To avoid this problem, we propose the following anti-marker distance term

\mathcal{D}_{AM}(x,y)=\frac{\exp\left(-\tilde{\alpha}\tilde{\mathcal{D}}_{GAM}(x,y)\right)-\exp\left(-\tilde{\alpha}\right)}{1-\exp\left(-\tilde{\alpha}\right)}

where $\tilde{\alpha}$ is a tuning parameter. We choose $\tilde{\alpha}=200$ throughout. This distance term ensures rapid decay of the penalty away from the set $\mathcal{AM}$ but still enforces high penalty around the anti-marker set itself. See Figure 6 where a segmentation result with and without anti-markers is shown. As $\mathcal{D}_{AM}(x,y)$ decays rapidly from $\mathcal{AM}$ , we do require that the anti-marker set be close to the blurred edge and away from the object we desire to segment.

3.4. The new model and its Euler-Lagrange equation

The Proposed Geodesic Model. Putting the above $3$ ingredients together, we propose the model

$\begin{split}\min_{u,c_{1},c_{2}}\Big{\{}F_{GEO}(u,c_{1},c_{2})=&\int_{\Omega}\left[\lambda_{1}(z(x,y)-c_{1})^{2}-\lambda_{2}(z(x,y)-c_{2})^{2}\right]u\,\mathrm{d}\Omega\\ &+\mu\int_{\Omega}g(|\nabla z(x,y)|)|\nabla u|\,\mathrm{d}\Omega+\theta\int_{\Omega}\mathcal{D}_{G}(x,y)u\,\mathrm{d}\Omega+\alpha\int_{\Omega}\nu_{\varepsilon}(u)\,\mathrm{d}\Omega\Big{\}},\end{split}$

(18)

where $\mathcal{D}_{G}(x,y)=\left(\mathcal{D}_{M}(x,y)+\mathcal{D}_{AM}(x,y)\right)/2$ and $\mathcal{D}_{M}(x,y)$ is the geodesic distance from the marker set $\mathcal{M}$ . We compute $\mathcal{D}_{M}(x,y)$ using (11) where $f(x,y)=f_{2}(x,y)$ defined in (17). Using Calculus of Variations, solving (18) with respect to $c_{1},\ c_{2}$ , with $u$ fixed, leads to

c_{1}(u)=\frac{\int_{\Omega}u\cdot z(x,y)\,\mathrm{d}\Omega}{\int_{\Omega}u\,\mathrm{d}\Omega},\hskip 36.135ptc_{2}(u)=\frac{\int_{\Omega}(1-u)\cdot z(x,y)\,\mathrm{d}\Omega}{\int_{\Omega}(1-u)\,\mathrm{d}\Omega},

(19)

and the minimisation with respect to $u$ (with $c_{1}$ and $c_{2}$ fixed) gives the PDE

\begin{gathered}\begin{aligned} \mu\nabla\cdot\left(g(|\nabla z(x,y)|)\frac{\nabla u}{|\nabla u|_{\varepsilon_{2}}}\right)-\Big{[}\lambda_{1}(z(x,y)-c_{1})^{2}&-\lambda_{2}(z(x,y)-c_{2})^{2}\Big{]}\\ &-\theta\mathcal{D}_{G}(x,y)-\alpha\nu^{\prime}_{\varepsilon}(u)=0\end{aligned}\end{gathered}

(20)

in $\Omega$ , where we replace $|\nabla u|$ with $|\nabla u|_{\varepsilon_{2}}=\sqrt{u_{x}^{2}+u_{y}^{2}+\varepsilon_{2}}$ to avoid zero denominator; we choose $\varepsilon_{2}=10^{-6}$ throughout. We also have Neumann boundary conditions $\frac{\partial u}{\partial{\bm{n}}}=0$ on $\partial\Omega$ where ${\bm{n}}$ is the outward unit normal vector.

Next we discuss a numerical scheme for solving this PDE (20). However it should be remarked that updating $c_{1}(u),c_{2}(u)$ should be done as soon as $u$ is updated; practically $c_{1},c_{2}$ converge very quickly since the object intensity $c_{1}$ does not change much.

4. An additive operator splitting algorithm

Additive Operator Splitting (AOS) is a widely used method [14, 27, 43] as seen from more recent works [2, 3, 4, 5, 37, 39] on the diffusion type equation such as

\frac{\partial u}{\partial t}=\mu\nabla\cdot(G(u)\nabla u)-f.

(21)

AOS allows us to split the two dimensional problem into two one-dimensional problems, which we solve and then combine. Each one dimensional problem gives rise to a tridiagonal system of equations which can be solved efficiently, hence AOS is a very efficient method for solving diffusion-like equations. AOS is a semi-implicit method and permits far larger time-steps than the corresponding explicit schemes would. Hence AOS is more stable than an explicit method [43]. We rewrite the above equation as

\frac{\partial u}{\partial t}=\mu\bigg{(}\partial_{x}\left(G(u)\partial_{x}u\right)+\partial_{y}\left(G(u)\partial_{y}u\right)\bigg{)}-f.

and after discretisation, we can rewrite this as [43]

u^{k+1}=\frac{1}{2}\sum_{\ell=1}^{2}\bigg{(}I-2\tau\mu A_{\ell}(u^{k})\bigg{)}^{-1}\left(u^{k}+\tau f\right)

where $\tau$ is the time-step, $A_{1}(u)=\partial_{x}(G(u)\partial_{x})$ and $A_{2}(u)=\partial_{y}(G(u)\partial_{y})$ . For notational convenience we write $G=G(u)$ . The matrix $A_{1}(u)$ can be obtained as follows

$\begin{gathered}\begin{aligned} \bigg{(}A_{1}(u^{k})u^{k+1}\bigg{)}_{i,j}=\bigg{(}\partial_{x}\Big{(}G\partial_{x}u^{k+1}\Big{)}\bigg{)}_{i,j}=&\Bigg{(}\frac{G_{i+\frac{1}{2},j}}{h_{x}^{2}}\Bigg{)}u^{k+1}_{i+1,j}+\Bigg{(}\frac{G_{i-\frac{1}{2},j}}{h_{x}^{2}}\Bigg{)}u^{k+1}_{i-1,j}-\Bigg{(}\frac{G_{i+\frac{1}{2},j}+G_{i-\frac{1}{2},j}}{h_{x}^{2}}\Bigg{)}u^{k+1}_{i,j}\end{aligned}\end{gathered}$

and similarly to [37, 39], for the half points in $G$ we take the average of the surrounding pixels, e.g. $G_{i+\frac{1}{2},j}=\frac{G_{i+1,j}+G_{i,j}}{2}$ . Therefore we must solve two tridiagonal systems to obtain $u^{k+1}$ , the Thomas algorithm allows us to solve each of these efficiently [43]. The AOS method described here assumes $f$ does not depend on $u$ , however in our case it depends on $\nu^{\prime}_{\varepsilon}(u)$ (see (20)) which has jumps around 0 and 1, so the algorithm has stability issues. This was noted in [39] and the authors adapted the formulation of (20) to offset the changes in $f$ . Here we repeat their arguments for adapting AOS when the exact penalty term $\nu^{\prime}_{\varepsilon}(u)$ is present (we refer to Figures 7 and 8 for plots of the penalty function and its derivative, respectively).

The main consideration is to extract a linear part out of the nonlinearity in $f=f(u)$ . If we evaluate the Taylor expansion of $\nu^{\prime}_{\varepsilon}(u)$ around $u=0$ and $u=1$ and group the terms into the constant and linear components in $u$ , we can respectively write $\nu^{\prime}_{\varepsilon}(u)=a_{0}(\varepsilon)+b_{0}(\varepsilon)u+\mathcal{O}(u^{2})$ and $\nu^{\prime}_{\varepsilon}(u)=a_{1}(\varepsilon)+b_{1}(\varepsilon)u+\mathcal{O}(u^{2})$ . We actually find that $b_{0}(\varepsilon)=b_{1}(\varepsilon)$ and denote the linear term as $b$ from now on. Therefore, for a change in $u$ of $\delta u$ around $u=0$ and $u=1$ , we can approximate the change in $\nu^{\prime}_{\varepsilon}(u)$ by $b\cdot\delta u$ .

To focus on the jumps, define the interval in which $\nu^{\prime}_{\varepsilon}(u)$ jumps as

I_{\zeta}:=[0-\zeta,0+\zeta]\cup[1-\zeta,1+\zeta]

and refine the linear function by

\tilde{b}_{i,j}^{k}=\begin{cases}b,&u_{i,j}^{k}\in I_{\zeta}\\ 0,&else.\end{cases}

Using these we can now offset the change in $\nu^{\prime}_{\varepsilon}(u^{k})$ by changing the formulation (21) to

\frac{\partial u}{\partial t}=\mu\nabla\cdot(G(u)\nabla u)-\alpha\tilde{b}^{k}u+\big{[}\alpha\tilde{b}^{k}u-f\big{]}

or in AOS form $u^{k+1}=u^{k}+\tau\mu\nabla\cdot(G(u^{k})\nabla u^{k+1})-\tau\alpha\tilde{b}^{k}u^{k+1}+\big{[}\tau\alpha\tilde{b}^{k}u^{k}-f^{k}\big{]}$ which, following the derivation in [39], can be reformulated as

u^{k+1}=\frac{1}{2}\sum_{\ell=1}^{2}{\underbrace{\bigg{(}I+\tilde{B}^{k}-2\tau\mu A_{\ell}\left(u^{k}\right)\bigg{)}}_{Q_{1}}}^{-1}\left(\left(I+\tilde{B}^{k}\right)u^{k}+\tau f^{k}\right)

where $\tilde{B}^{k}=\text{diag}(\tau\alpha\tilde{b}^{k})$ . We note that $Q_{1}$ is invertible as it is strictly diagonally dominant. This scheme improves on (21) as now, changes in $f^{k}$ are damped. However, it is found in [39] that although this scheme does satisfy most of the discrete scale space conditions of Weickert [43] (which guarantee convergence of the scheme), it does not satisfy all of them. In particular the matrix $Q_{1}$ doesn’t have unit row sum and is not symmetrical. The authors adapt the scheme above to the equivalent

u^{k+1}=\frac{1}{2}\sum_{\ell=1}^{2}{\underbrace{\bigg{(}I-2\tau\mu\left(I+\tilde{B}^{k}\right)^{-1}A_{\ell}\left(u^{k}\right)\bigg{)}}_{Q_{2}}}^{-1}\left(u^{k}+\tau\left(I+\tilde{B}^{k}\right)^{-1}f^{k}\right),

(22)

where the matrix $Q_{2}$ does have unit row sum, however the matrix is not always symmetrical. We can guarantee convergence for $\zeta=0.5$ (in which case $Q_{2}$ must be symmetrical) but we desire to use a small $\zeta$ to give a small interval $I_{\zeta}$ . We find experimentally that convergence is achieved for any small value of $\zeta$ , this is due to the fact that at convergence the solution $u$ is almost binary [10]. Therefore, although initially $Q_{2}$ is asymmetrical at some pixels, at convergence all pixels have values which fall within $I_{\zeta}$ and $I+\tilde{B}^{k}$ is a matrix with all diagonal entries $1+\tau\alpha b$ . Therefore we find that at convergence $Q_{2}$ is symmetrical and the discrete scale space conditions are all satisfied. In all of our tests we fix $\zeta=0.01$ .

Set

\mu,\lambda,\theta

. Compute

g(|\nabla z(x,y)|)=\frac{1}{1+\beta_{G}|\nabla z(x,y)|^{2}}

and

\mathcal{D}_{G}(x,y)=\frac{\mathcal{D}_{G}^{0}(x,y)}{||\mathcal{D}_{G}^{0}(x,y)||_{L^{\infty}}}

with

\mathcal{D}_{G}^{0}(x,y)

the solution of (11). Initialise

u^{(0)}

arbitrarily.

for

iter=1

max\_iterations

Calculate

c_{1}

and

c_{2}

using (19).

Calculate

r=\lambda_{1}(z-c_{1})^{2}-\lambda_{2}(z-c_{2})^{2}+\theta\mathcal{D}_{G}

Set

\alpha=||r||_{L^{\infty}}

Calculate

f^{k}=r+\alpha\nu^{\prime}_{\varepsilon}(u^{k})

Update

u^{k}

u^{k+1}

using the AOS scheme (22).

end for

u^{*}\leftarrow u^{k}

Algorithm 1 Solution of the Geodesic Model

5. Existence and Uniqueness of the Viscosity Solution

In this section we use the viscosity solution framework and the work of Ishii and Sato [20] to prove that, for a class of PDEs in image segmentation, the solution exists and is unique. In particular, we prove the existence and uniqueness of the viscosity solution for the PDE which is determined by the Euler-Lagrange equation for the Geodesic Model. Throughout, we will assume $\Omega$ is a bounded domain with $C^{1}$ boundary.

From the work of [12, 20], we have the following Theorem for analysing the solution of a partial differential equation of the form $F(\bm{x},u,Du,D^{2}u)=0$ where $F:\mathbb{R}^{n}\times\mathbb{R}\times\mathbb{R}^{n}\times\mathscr{M}^{n}\rightarrow\mathbb{R}$ , $\mathscr{M}^{n}$ is the set of $n\times n$ symmetric matrices, $Du$ is the gradient of $u$ and $D^{2}u$ is the Hessian of $u$ . For simplicity, and in a slight abuse of notation, we use $x:=\bm{x}$ for the vector of a general point in $\mathbb{R}^{n}$ .

Theorem 2 (Theorem 3.1 [20]).

Assume that the following conditions (C1)–(C2) and (I1)–(I7) hold. Then for each $u_{0}\in C(\overline{\Omega})$ there is a unique viscosity solution $u\in C([0,T)\times\overline{\Omega})$ of (23) and (24) satisfying (25).

\frac{\partial u}{\partial t}+F(t,x,u,Du,D^{2}u)=0\hskip 14.45377pt\text{in}\hskip 14.45377ptQ=(0,T)\times\Omega,

(23)

B(x,Du)=0\hskip 14.45377pt\text{in}\hskip 14.45377ptS=(0,T)\times\partial\Omega,

(24)

u(0,x)=u_{0}(x)\hskip 14.45377pt\text{for}\hskip 14.45377ptx\in\overline{\Omega}.

(25)

Conditions (C1)–(C2).

(C1)

$F(t,x,u,p,X)\leq F(t,x,v,p,X)$ for $u\leq v$ .
(C2)

$F(t,x,u,p,X)\leq F(t,x,u,p,Y)$ for $X,Y\in\mathscr{M}^{n}$ and $Y\leq X$ .

Conditions (I1)–(I7). Assume $\Omega$ is a bounded domain in $\mathbb{R}^{n}$ with $C^{1}$ boundary.

(I1)

$F\in C\left([0,T]\times\overline{\Omega}\times\mathbb{R}\times\left(\mathbb{R}^{n}\backslash\{0\}\right)\times\mathscr{M}^{n}\right)$ .
(I2)

There exists a constant $\gamma\in\mathbb{R}$ such that for each $(t,x,p,X)\in[0,T]\times\overline{\Omega}\times\left(\mathbb{R}^{n}\backslash\{0\}\right)\times\mathscr{M}^{n}$ the function $u\mapsto F(t,x,u,p,X)-\gamma u$ is non-decreasing on $\mathbb{R}$ .
(I3)

$F$ is continuous at $(t,x,u,0,0)$ for any $(t,x,u)\in[0,T]\times\overline{\Omega}\times\mathbb{R}$ in the sense that

$-\infty<F_{*}(t,x,u,0,0)=F^{*}(t,x,u,0,0)<\infty$

holds. Here $F^{*}$ and $F_{*}$ denote, respectively, the upper and lower semi-continuous envelopes of $F$ , which are defined on $[0,T]\times\overline{\Omega}\times\mathbb{R}\times\mathbb{R}^{n}\times\mathscr{M}^{n}$ .
(I4)

$B\in C\left(\mathbb{R}^{n}\times\mathbb{R}^{n}\right)\cap C^{1,1}\left(\mathbb{R}^{n}\times\left(\mathbb{R}^{n}\backslash\{0\}\right)\right)$ , where $C^{1,1}$ is the Hölder functional space.
(I5)

For each $x\in\mathbb{R}^{n}$ the function $p\mapsto B(x,p)$ is positively homogeneous of degree one in $p$ , i.e. $B(x,\lambda p)=\lambda B(x,p)$ for all $\lambda\geq 0$ and $p\in\mathbb{R}^{n}\backslash\{0\}$ .
(I6)

There exists a positive constant $\Theta$ such that $\langle{\bm{n}}(x),D_{p}B(x,p)\rangle\geq\Theta$ for all $x\in\partial\Omega$ and $p\in\mathbb{R}^{n}\backslash\{0\}$ . Here ${\bm{n}}(x)$ denotes the unit outward normal vector of $\Omega$ at $x\in\partial\Omega$ .

(I7)

For each $R>0$ there exists a non-decreasing continuous function $\omega_{R}:[0,\infty)\rightarrow[0,\infty)$ satisfying $\omega_{R}(0)=0$ such that if $X,Y\in\mathscr{M}^{n}$ and $\mu_{1},\mu_{2}\in[0,\infty)$ satisfy

\begin{bmatrix}X&0\\ 0&Y\\ \end{bmatrix}\leq\mu_{1}\begin{bmatrix}I&-I\\ -I&I\\ \end{bmatrix}+\mu_{2}\begin{bmatrix}I&0\\ 0&I\\ \end{bmatrix}

(26)

then

\begin{gathered}\begin{aligned} F(t,x,u,p,X)-F(t,y,u,q,-Y)\geq&-\omega_{R}\Big{(}\mu_{1}\left(|x-y|^{2}+\rho(p,q)^{2}\right)+\mu_{2}+|p-q|\\ &+|x-y|\left(\max(|p|,|q|)+1\right)\Big{)}\\ \end{aligned}\end{gathered}

for all $t\in[0,T],x,y\in\overline{\Omega},u\in\mathbb{R}$ , with $|u|\leq R$ , $p,q\in\mathbb{R}^{n}\backslash\{0\}$ and $\rho(p,q)=\min\left(\frac{|p-q|}{\min(|p|,|q|)},1\right)$ .

5.1. Existence and uniqueness for the Geodesic Model

We now prove that there exists a unique solution for the PDE (20) resulting from the minimisation of the functional for the Geodesic Model (18).

Remark 3.

It is important to note that although the values of $c_{1}$ and $c_{2}$ depend on $u$ , they are fixed when we solve the PDE for $u$ and therefore the probem is a local one and Theorem 2 can be applied. Once we update $c_{1}$ and $c_{2}$ , using the updated $u$ , then we fix them again and apply Theorem 2. In practice, as we near convergence, we find $c_{1}$ and $c_{2}$ stabilise so we typically stop updating $c_{1}$ and $c_{2}$ once the change in both values is below a tolerance.

To apply the above theorem to the proposed model (20), the key step will be to verify the nine conditions. First, we multiply (20) by the factor $|\nabla u|_{\varepsilon_{2}}$ , obtaining the nonlinear PDE

	$\displaystyle-\mu\|\nabla u\|_{\varepsilon_{2}}\nabla\cdot\bigg{(}G(x,\nabla z)\frac{\nabla u}{\|\nabla u\|_{\varepsilon_{2}}}\bigg{)}+\|\nabla u\|_{\varepsilon_{2}}\bigg{[}\lambda_{1}(z(x,y)-c_{1})^{2}-\lambda_{2}(z(x,y)-c_{2})^{2}\bigg{.}$			(27)
	$\displaystyle\bigg{.}+\theta\mathcal{D}_{G}(x,y)+\alpha\nu^{\prime}_{\varepsilon}(u)$	$\displaystyle\bigg{]}=0$		(27)

where $G(x,\nabla z)=g(|\nabla z(x,y)|)$ . We can rewrite this as

F(x,u,p,X)=-\mu\,\text{trace}\left(A(x,p)X\right)-\mu\langle\nabla G(x,\nabla z),p\rangle+|p|k(u)+|p|f(x)=0

(28)

where $f(x)=\lambda_{1}(z(x)-c_{1})^{2}-\lambda_{2}(z(x)-c_{2})^{2}$ , $k(u)=\alpha\nu^{\prime}_{\varepsilon}(u)$ , $p=(p_{1},p_{2})=|\nabla u|_{\varepsilon_{2}}$ , $X$ is the Hessian of $u$ and

A(x,p)=\begin{bmatrix}G(x,\nabla z)\frac{p_{2}^{2}}{|p|^{2}}&-G(x,\nabla z)\frac{p_{1}p_{2}}{|p|^{2}}\\ -G(x,\nabla z)\frac{p_{1}p_{2}}{|p|^{2}}&G(x,\nabla z)\frac{p_{1}^{2}}{|p|^{2}}\\ \end{bmatrix}

(29)

Theorem 4 (Theory for the Geodesic Model).

The parabolic PDE $\frac{\partial u}{\partial t}+F(t,x,u,Du,D^{2}u)=0$ with $u_{0}=u(0,x)\in C(\overline{\Omega})$ , $F$ as defined in (28) and Neumann boundary conditions has a unique solution $u=u(t,x)$ in $C([0,T)\times\overline{\Omega})$ .

Proof. By Theorem 2, it remains to verify that $F$ satisfies (C1)–(C2) and (I1)–(I7). We will show that each of the conditions is satisfied. Most are simple to show, the exception being (I7) which is non-trivial.

(C1): Equation (28) only has dependence on $u$ in the term $k(u)$ , we therefore have a restriction on the choice of $k$ , requiring $k(v)\geq k(u)$ for $v\geq u$ . This is satisfied for $k(u)=\alpha\nu^{\prime}_{\varepsilon}(u)$ with $\nu^{\prime}_{\varepsilon}(u)$ defined as in (7).

(C2): We find for arbitrary $s=(s_{1},s_{2})\in\mathbb{R}^{2}$ that $s^{T}A(x,p)s\geq 0$ and so $A(x,p)\geq 0$ . It follows that $-\text{trace}(A(x,p)X)\leq-\text{trace}(A(x,p)Y)$ , therefore this condition is satisfied.

(I1): $A(x,p)$ is only singular at $p=0$ , however it is continuous elsewhere and satisfies this condition.

(I2): In $F$ the only term which depends on $u$ is $k(u)=\alpha\nu^{\prime}_{\varepsilon}(u)$ . With $\nu^{\prime}_{\varepsilon}(u)$ defined as in (7), in the limit $\varepsilon\rightarrow 0$ this function is a step function from $-2$ on $(\infty,0)$ , $0$ on $[0,1]$ and $2$ on $(0,\infty)$ . So we can choose any constant $\varepsilon<-2$ . With $\varepsilon\neq 0$ there is smoothing at the end of the intervals, however there is still a lower bound on $L$ for $\nu^{\prime}_{\varepsilon}(u)$ and we can choose any constant $\gamma<L$ .

(I3): $F$ is continuous at $(x,0,0)$ for any $x\in\Omega$ because $F^{*}(x,0,0)=F_{*}(x,0,0)=0.$ Hence this condition is satisfied.

(I4): The Euler-Lagrange equations give Neumann boundary conditions

B(x,\nabla u)=\frac{\partial u}{\partial{\bm{n}}}={\bm{n}}\cdot\nabla~u=\langle{\bm{n}},\nabla u\rangle=0

on $\partial\Omega$ , where ${\bm{n}}$ is the outward unit normal vector, and we see that $B(x,\nabla u)\in C^{1,1}\left(\mathbb{R}^{n}\times\mathbb{R}^{n}\backslash\{0\}\right)$ and therefore this condition is satisfied.

(I5): By the definition above, $B(x,\lambda\nabla u)=\langle{\bm{n}},\lambda\nabla u\rangle=\lambda\langle{\bm{n}},\nabla u\rangle=\lambda B(x,\nabla u)$ . So this condition is satisfied.

(I6): As before we can use the definition, $\langle\bm{n}(x),D_{p}B(x,p)\rangle=\langle\bm{n}(x),\bm{n}(x)\rangle=|\bm{n}(x)|^{2}$ . So we can choose $\Theta=1$ and the condition is satisfied.

(I7): This is the most involved condition to prove and uses many other results. For clarity of the overall paper, we postpone the proof to Appendix A. $\Box$

5.2. Generalisation to other related models

Theorems 2 and 4 can be generalised to a few other models. This amounts to writing each model as a PDE of the form (28) where $k(u)$ is monotone and $f(x),k(u)$ are bounded. This is summarised in the following Corollary:

Corollary 5.

Assume that $c_{1}$ and $c_{2}$ are fixed, with the terms $f(x)$ and $k(u)$ respectively defined as follows for a few related models:

•

Chan-Vese [11]: $f(x)=f_{CV}(x):=\lambda_{1}(z(x)-c_{1})^{2}-\lambda_{2}(z(x)-c_{2})^{2}$ , $k(u)=0$ .
•

Chan-Vese (Convex) [10]: $f(x)=f_{CV}(x)$ , $k(u)=\alpha\nu^{\prime}_{\varepsilon}(u)$ .
•

Geodesic Active Contours [8] and Gout et al. [25]: $f(x)=0$ , $k(u)=0$ .
•

Nguyen et al. [30]: $f(x)=\alpha\left(P_{B}(x,y)-P_{F}(x,y)\right)+\left(1-\alpha\right)\left(1-2P(x,y)\right)$ , $k(u)=0$ .
•

Spencer-Chen (Convex) [39]: $f(x)=f_{CV}(x)+\theta\mathcal{D}_{E}(x)$ , $k(u)=\alpha\nu^{\prime}_{\varepsilon}(u)$ .

Then if we define a PDE of the general form

-\mu\nabla\cdot\left(G(x)\frac{\nabla u}{|\nabla u|_{\varepsilon_{2}}}\right)+k(u)+f(x)=0

with

(i)

Neumann boundary conditions $\frac{\partial u}{\partial{\bm{n}}}=0$ ( ${\bm{n}}$ the outward normal unit vector)
(ii)

$k(u)$ satisfies $k(u)\geq k(v)$ if $u\geq v$
(iii)

$k(u)$ and $f(x)$ are bounded; and
(iv)

$G(x)=Id$ or $G(x)=f(|\nabla z(x)|)=\frac{1}{1+|\nabla z(x)|^{2}}$ ,

we have a unique solution $u\in C([0,T)\times\overline{\Omega})$ for a given initialisation. Consequently we conclude that all above models admit a unique solution.

Proof. The conditions (i)–(iv) are hold for all of these models. All of these models require Neumann boundary conditions and use the permitted G(x). The monotonicity of $\nu^{\prime}_{\varepsilon}(u)$ is discussed in the proof of (C1) for Theorem 4 and the boundedness of $f(x)$ and $k(u)$ is clear in all cases. $\Box$

Remark 6.

Theorem 4 and Corollary 5 also generalise to cases where $G(x)=\frac{1}{1+\beta|\nabla z|^{2}}$ and to $G(x)=\mathcal{D}(x)g(|\nabla z|)$ where $\mathcal{D}(x)$ is a distance function such as in [15, 16, 17, 39]. The proof is very similar to that shown in §5.1, relying on Lipschitz continuity of the function $G(x)$ .

Remark 7.

We cannot apply the classical viscosity solution framework to the Rada-Chen model [37] as this is a non-local problem with $k(u)=2\nu\left(\int_{\Omega}H_{\varepsilon}(u)\,\mathrm{d}\Omega-A_{1}\right)$ .

6. Numerical Results

In this section we will demonstrate the advantages of the Geodesic Model for selective image segmentation over related and previous models. Specifically we shall compare

•

M1 — the Nguyen et al. (2012) model [30];
•

M2 — the Rada-Chen (2013) model [37];
•

M3 — the convex Spencer-Chen (2015) model [39];
•

M4 — the convex Liu et al. (2017) model [26];
•

M5 — the reformulated Rada-Chen model with geodesic distance penalty (see Remark 8);
•

M6 — the reformulated Liu et al. model with geodesic distance penalty (see Remark 8);
•

M7 — the proposed convex Geodesic Model (Algorithm 1).

Remark 8 (A note on M5 and M6).

We include M5 – M6 to test how the geodesic distance penalty term can improve M2 [37] and M4 [26]. These were obtained as follows:

•

we extend M2 to M5 simply by including the geodesic distance function $\mathcal{D}_{G}(x,u)$ in the functional.

•

we extend M4 to M6 with a minor reformulation to include data fitting terms. Specifically, the model M6 is

\begin{split}\min_{u,c_{1},c_{2}}\Big{\{}F_{CV\omega}(u,c_{1},c_{2})=&\int_{\Omega}\omega^{2}(x,y)\left[\lambda_{1}(z(x,y)-c_{1})^{2}-\lambda_{2}(z(x,y)-c_{2})^{2}\right]u\,\mathrm{d}\Omega\\ +&\mu\int_{\Omega}g(|\nabla z|))|\nabla u|\,\mathrm{d}\Omega+\theta\int_{\Omega}\mathcal{D}_{G}(x,y)u\,\mathrm{d}\Omega+\alpha\int_{\Omega}\nu_{\varepsilon}(u)\,\mathrm{d}\Omega\Big{\}}\end{split}

(30)

for $\mu,\lambda_{1},\lambda_{2}$ non-negative fixed parameters, $\alpha$ and $\nu_{\varepsilon}(u)$ as defined in (7) and $\omega$ as defined for the convex Liu et al. model. This is a convex model and is the same as the proposed Geodesic Model M7 but with weighted intensity fitting terms.

Four sets of test results are shown below. In Test 1 we compare models M1 – M6 to the proposed model M7 for two images which are hard to segment. The first is a CT scan from which we would like to segment the lower portion of the heart, the second is an MRI scan of a knee and we would like to segment the top of the Tibia. See Figure 9 for the test images and the marker sets used in the experiments. In Test 2 we will review the sensitivity of the proposed model to the main parameters. In Test 3 we will give several results achieved by the model using marker and anti-marker sets. In Test 4 we show the initialisation independence and marker independence of the Geodesic Model on real images.

For M7, we denote by $\tilde{u}$ the thresholded $u>\tilde{\gamma}$ at some value $\tilde{\gamma}\in(0,1)$ to define the segmented region. Although the threshold can be chosen arbitrarily in $(0,1)$ from the work by [10, Thm 1] and [39], we usually take $\tilde{\gamma}=0.5$ .

Quantitative Comparisons. To measure the quality of a segmentation, we use the Tanimoto Coefficient (TC) (or Jaccard Coefficient [21]) defined by

TC(\tilde{u},GT)=\frac{|\tilde{u}\cap GT|}{|\tilde{u}\cup GT|}

where GT is the ‘ground truth’ segmentation and $\tilde{u}$ is the result from a particular model. This measure takes value one for a segmentation which coincides perfectly with the ground truth and reduces to zero as the quality of the segmentation gets worse. In the other tests, where a ground truth is not available, we use visual plots.

Parameter Choices and Implementation. We set $\mu=1$ , $\tau=10^{-2}$ and vary $\lambda=\lambda_{1}=\lambda_{2}$ and $\theta$ . Following [10] we let $\alpha=||\lambda_{1}(z-c_{1})^{2}-\lambda_{2}(z-c_{2})^{2}+\theta\mathcal{D}_{G}(x,y)||_{L^{\infty}}$ . To implement the marker points in MATLAB we use roipoly for choosing a small number of points by clicking and also freedraw which allows the user to draw a path of marker points. The stopping criteria used is the dynamic residual falling below a given threshold, i.e. once $||u^{k+1}-u^{k}||/||u^{k}||<tol$ the iterations stop (we use $tol=10^{-6}$ in the tests shown).

Test 1 – Comparison of models M1 – M7.

In this test we give the segmentation results for models M1 – M7 for the two challenging test images shown in Figure 9. The marker and anti-marker sets used in the experiments are also shown in this figure. After extensive parameter tuning, the best final segmentation results for each of the models are shown in Figures 10 and 11. For M1 – M4 we obtain incorrect segmentations in both cases. In particular, the results of M2 and M4 are interesting as the former gives poor results for both images, and the latter gives a reasonable result for Test Image 1 and a poor result for Test Image 2. In the case of M2, the regularisation term includes the edge detector and the distance penalty term (see (4)). It is precisely this which permits the poor result in Figures 10(b) and 11(b) as the edge detector is zero along the contour and the fitting terms are satisfied there (both intensity and area constraints) – the distance term is not large enough to counteract the effect of these. In the case of M4, the distance term and edge detector are separated from the regulariser and are used to weight the Chan-Vese fitting terms (see (9)). The poor segmentation in Figure 11(b) is due to the Chan-Vese terms encouraging segmentation of bright objects (in this case), weighting $\omega$ enforces these terms at all edges in the image and near $\mathcal{M}$ . In experiments, we find that M4 performs well when the object to segment is of approximately the highest or lowest intensity in the image, however when this is not the case, results tend to be poor. We see that, in both cases, models M5 and M6 give much improved results to M2 and M4 (obtained by incorporating the geodesic distance penalty into each). The proposed Geodesic Model M7 gives an accurate segmentation in both cases. It remains to compare M5, M6 and M7. We see that M5 is a non-convex model (and cannot be made convex [39]), therefore results are initialisation dependent. It also requires one more parameter than M6 and M7, and an accurate set $\mathcal{M}$ to give a reasonable area constraint in (4). These limitations lead us to conclude M6 and M7 are better choices than M5. In the case of M6, it has the same number of parameters as M7 and gives good results. M6 can be viewed as the model M7 with weighted intensity fitting terms (compare (18) and (30)). Experimentally, we find that the same quality of segmentation result can be achieved with both models generally, however M6 is more parameter sensitive than M7. This can be seen in the parameter map in Figure 12 with M7 giving an accurate result for a wider range of parameters than M6. To show the improvement of M7 over previous models, we also give an image in Figure 13 which can be accurately segmented with M7 but the correct result is never achieved with M6 (or M3). Therefore we find that M7 outperforms all other models tested M1 – M6.

Remark 9.

Models M2 – M7 are coded in MATLAB and use exactly the same marker/anti-marker set. For model M1, the software of Nguyen et al. requires marker and anti-marker sets to be input to an interface. These have been drawn as close as possible to match those used in the MATLAB code.

Test 2 – Test of M7’s sensitivity to changes in its main parameters. In this test we demonstrate that the proposed Geodesic Model is robust to changes in the main parameters. The main parameters in (20) are $\mu,\lambda_{1},\lambda_{2},\theta$ and $\varepsilon_{2}$ . In all tests we set $\mu=1$ , which is simply a rescaling of the other parameters, and we set $\lambda=\lambda_{1}=\lambda_{2}$ . In the first example, in Figure 12, we compare the TC value for various $\lambda$ and $\theta$ values for segmentation of a bone in a knee scan. We see that the segmentation is very good for a larger range of $\theta$ and $\lambda$ values. For the second example, in Figure 13, we show an image and marker set for which the Spencer-Chen model (M3) and modified Liu et al. model M6 cannot achieve the desired segmentation for any parameter range, but which can be attained for the Geodesic Model for a vast range of parameters. The final example, in Table 1, compares the TC values for various $\varepsilon_{2}$ values with fixed parameters $\lambda=2$ and $\theta=2$ . We use the images and ground truth as shown in Figures 12 and 13: on the synthetic circles image we obtain a perfect segmentation for all values of $\varepsilon_{2}$ tested, and in the case of the knee segmentation the results are almost identical for any $\varepsilon_{2}<10^{-6}$ , above which the quality slowly deteriorates.

$\varepsilon_{2}$	Knee Segmentation (Figure 12)	Circle Segmentation (Figure 13)
$10^{-10}$	0.97287	1.00000
$10^{-8}$	0.97287	1.00000
$10^{-6}$	0.97235	1.00000
$10^{-4}$	0.96562	1.00000
$10^{-2}$	0.94463	1.00000
$10^{0}$	0.90660	1.00000
$10^{2}$	0.89573	1.00000
$10^{4}$	0.89159	1.00000

Table 1: The Tanimoto Coeffcient for various

\varepsilon_{2}

values, segmenting the images in Figures 12 and 13.

Test 3 – Further Results from the Geodesic Model M7. In this test we give some medical segmentation results obtained using the Geodesic Model M7. The results are shown in Figure 14. In the final two columns we use anti-markers to demonstrate how to overcome blurred edges and low contrast edges in an image. These are challenging and it is pleasing to see the correctly segmented results.

Test 4 – Initialisation and Marker Set Independence. In the first example, in Figure 16, we see how the convex Geodesic Model M7 gives the same segmentation result regardless of initialisation, as expected of a convex model. Hence the model is flexible in implementation. From many experiments it is found that using the polygon formed by the marker points as the initialisation converges to the final solution faster than using an arbitrary initialisation. In the second example, in Figure 16, we show intuitively how Model M7 is robust to the number of markers and the location of the markers within the object to be segmented. The Euclidean distance term, used in the Spencer-Chen model M3, is sensitive to the position and number of marker points, however, regardless of where the markers are chosen, and how many are chosen, the geodesic distance map will be almost identical.

7. Conclusions

In this paper a new convex selective segmentation model has been proposed, using geodesic distance as a penalty term. This model gives results that are unachievable by alternative selective segmentation models and is also more robust to the parameter choices. Adaptations to the penalty term have been discussed which make it robust to noisy images and blurry edges whilst also penalising objects far from the marker set (in a Euclidean distance sense). A proof for the existence and uniqueness of the viscosity solution to the PDE given by the Euler-Lagrange equation for the model has been given (which applies to an entire class of image segmentation PDEs). Finally we have confirmed the advantages of using the geodesic distance with some experimental results. Future works will look for further extension of selective segmentation to other frameworks such as using high order regularizers [46, 13] where only incomplete theories exist.

Acknowledgements

The first author wishes to thank the UK EPSRC, the Smith Institute for Industrial Mathematics, and the Liverpool Heart and Chest Hospital for supporting the work through an Industrial CASE award. Thanks also must go to Professor Joachim Weickert (Saarland, Germany) for fruitful discussions at the early stages of this work. The second author is grateful to the UK EPSRC for the grant EP/N014499/1.

Appendix A — Proof that Condition (I7) Holds in Theorem 4

Using the assumption in (26), we write

\begin{gathered}\begin{aligned} (Xr,r)+(Ys,s)=r^{T}Xr+s^{T}Ys&\leq\mu_{1}\begin{bmatrix}r^{T}&s^{T}\end{bmatrix}\begin{bmatrix}I&-I\\ -I&I\\ \end{bmatrix}\begin{bmatrix}r\\ s\\ \end{bmatrix}+\mu_{2}\begin{bmatrix}r^{T}&s^{T}\end{bmatrix}\begin{bmatrix}I&0\\ 0&I\\ \end{bmatrix}\begin{bmatrix}r\\ s\\ \end{bmatrix}\\ &=\mu_{1}|r-s|^{2}+\mu_{2}\left(|r|^{2}+|s|^{2}\right).\\ \end{aligned}\end{gathered}

Note that matrix $A$ from (29) is a real symmetric matrix and decomposes as $A=QDQ^{T}=QD^{1/2}D^{1/2}Q^{T}=BB^{T}$ with $Q$ orthonormal and $B=QD^{1/2}$ . Successively define $r=B(p)e_{i}$ and $s=B(p)e_{i}$ for all $(e_{i})$ , an orthonormal basis, and obtain

\begin{gathered}\begin{aligned} (Xr,r)=r^{T}Xr&=\sum_{i}(Be_{i})^{T}X(Be_{i})=\sum_{i}e_{i}^{T}B^{T}XBe_{i}=\text{trace}(B^{T}XB)=\text{trace}(A(x,p)X).\end{aligned}\end{gathered}

Therefore, we can write

\begin{gathered}\begin{aligned} \text{trace}(A(x,p)X)+\text{trace}(A(y,q)Y)&=(XB(p)e_{i},B(p)e_{i})+(YB(q)e_{i},B(q)e_{i})\\ &\leq\mu_{1}|B(p)e_{i}-B(q)e_{i}|^{2}+\mu_{2}\left(|B(p)e_{i}|^{2}+|B(q)e_{i}|^{2}\right)\\ &=\mu_{1}\text{trace}\left(\left(B(p)-B(q)\right)^{T}\left(B(p)-B(q)\right)\right)+\mu_{2}\left(G(x)+G(y)\right).\\ \end{aligned}\end{gathered}

We now focus on reformulating the first term, we start by decomposing $A(x,p)$ as follows

A(x,p)=\begin{bmatrix}\frac{p_{1}}{|p|}&-\frac{p_{2}}{|p|}\\ \frac{p_{2}}{|p|}&\frac{p_{1}}{|p|}\\ \end{bmatrix}\begin{bmatrix}0&0\\ 0&G(x)\\ \end{bmatrix}\begin{bmatrix}\frac{p_{1}}{|p|}&\frac{p_{2}}{|p|}\\ -\frac{p_{2}}{|p|}&\frac{p_{1}}{|p|}\\ \end{bmatrix}=\begin{bmatrix}\frac{p_{1}}{|p|}&-\frac{p_{2}}{|p|}\\ \frac{p_{2}}{|p|}&\frac{p_{1}}{|p|}\\ \end{bmatrix}\begin{bmatrix}0&0\\ 0&\sqrt{G(x)}\\ \end{bmatrix}\begin{bmatrix}0&0\\ 0&\sqrt{G(x)}\\ \end{bmatrix}\begin{bmatrix}\frac{p_{1}}{|p|}&\frac{p_{2}}{|p|}\\ -\frac{p_{2}}{|p|}&\frac{p_{1}}{|p|}\\ \end{bmatrix}

so we have $A=BB^{T}$ where

B(p)=\begin{bmatrix}0&-\frac{p_{2}}{|p|}\sqrt{G(x)}\\ 0&\frac{p_{1}}{|p|}\sqrt{G(x)}\\ \end{bmatrix}.

Using this we compute

\begin{gathered}\begin{aligned} \text{trace}\left(\left(B(p)-B(q)\right)^{T}\left(B(p)-B(q)\right)\right)=\left|\frac{p}{|p|}\sqrt{G(x)}-\frac{q}{|q|}\sqrt{G(y)}\right|^{2}.\end{aligned}\end{gathered}

Substituting this in the overall trace sum we have

\begin{gathered}\begin{aligned} \text{trace}(A(x,p)X)+\text{trace}(A(y,q)Y)\leq\mu_{1}\left|\frac{p}{|p|}\sqrt{G(x)}-\frac{q}{|q|}\sqrt{G(y)}\right|^{2}+2\mu_{2}\theta.\\ \end{aligned}\end{gathered}

as $G(x)<\theta$ ( $G$ is bounded) for all $x\in\Omega$ . Focussing on the first term in this expression we compute

\begin{gathered}\begin{aligned} \left|\frac{p}{|p|}\sqrt{G(x)}-\frac{q}{|q|}\sqrt{G(y)}\right|^{2}&=\left|\frac{p}{|p|}\sqrt{G(x)}-\frac{p}{|p|}\sqrt{G(y)}+\frac{p}{|p|}\sqrt{G(y)}-\frac{q}{|q|}\sqrt{G(y)}\right|^{2}\\ &\leq 2\left(\sqrt{G(x)}-\sqrt{G(y)}\right)^{2}+2G(y)\left|\frac{p}{|p|}-\frac{q}{|q|}\right|^{2}\\ &\leq 2\left(\sqrt{G(x)}-\sqrt{G(y)}\right)^{2}+8\theta\rho(p,q)^{2}\\ \end{aligned}\end{gathered}

where $\rho=\min\left(\frac{|p-q|}{\min(|p|,|q|)},1\right)$ . This uses inequality $\left|\frac{p}{|p|}-\frac{q}{|q|}\right|^{2}\leq 2\rho(p,q)$ (see [15, 16, 17, 18, 24, 35]). We now note that $g(s)=\frac{1}{1+s^{2}}$ is Lipschitz continuous with Lipschitz constant $\frac{3\sqrt{3}}{8}$ .

Note. In the Geodesic Model we fix $G(x)=g(|\nabla z|)$ . Therefore, assuming $G(x)$ and $\sqrt{G(x)}$ as Lipschitz requires us to assume that the underlying $z$ is a smooth function [16]. Thankfully, $z$ is typically provided as a smoothed image after some filtering (e.g. Gaussian smoothing) and we can assume regularity of $z$ .

Remark 10.

It is less clear that $\sqrt{G(x)}$ is Lipschitz, we now prove it explicitly. Firstly, it is relatively easy to prove that

\sqrt{G(x)}-\sqrt{G(y)}\leq\frac{2}{3\sqrt{3}}\bigg{|}\left|\nabla z(x)\right|-\left|\nabla z(y)\right|\bigg{|}

by letting $K(s)=\sqrt{g(s)}$ and we find $\sup\limits_{s}|K^{\prime}(s)|=\frac{2}{3\sqrt{3}}$ . We now need to prove that $|\nabla z(x)|$ is Lipschitz also. Take $h(x)=|\nabla z(x)|$ , then by a remark in [16], we can conclude $\exists\,\zeta<\infty$ such that

\bigg{|}\left|\nabla z(x)\right|-\left|\nabla z(y)\right|\bigg{|}\leq\zeta|x-y|

and so $\sqrt{G(x)}$ is Lipschitz with constant $\frac{2}{3\sqrt{3}}\zeta$ .

After some computations we obtain

\left|\frac{p}{|p|}\sqrt{G(x)}-\frac{q}{|q|}\sqrt{G(y)}\right|^{2}\leq 2\left(\frac{2}{3\sqrt{3}}\zeta\right)^{2}\left|x-y\right|^{2}+8\theta\rho(p,q)^{2}=\frac{8}{27}\zeta^{2}\left|x-y\right|^{2}+8\theta\rho(p,q)^{2}.

Following the results in [15, 16, 17, 18, 24, 35] we have

\left|\nabla G(x)-\nabla G(y)\right|\left|p\right|<\kappa|p||x-y|\leq\kappa\max(|p|,|q|)|x-y|.

so overall

\langle\nabla G(x),p\rangle-\langle\nabla G(y),q\rangle\leq\kappa\max(|p|,|q|)|x-y|+\eta|p-q|\\

where $|\nabla G(y)|<\eta<\infty$ . Finally, we note that $-\left(|p|-|q|\right)=|q|-|p|\leq\Big{|}|q|-|p|\Big{|}\leq|p-q|$ . If we now write

\begin{gathered}\begin{aligned} -\left(F(t,x,u,p,X)-F(t,y,u,q,-Y)\right)=&\mu\left(\text{trace}(A(x,p)X)+\text{trace}(A(y,q)Y)\right)\\ &+\mu\left(\langle\nabla G(x),p\rangle-\langle\nabla G(y),q\rangle\right)\\ &-\left(|p|-|q|\right)k(u)-|p|f(x)+|q|f(y)\\ \leq&~\mu\mu_{1}\left(\frac{8}{27}\zeta^{2}|x-y|^{2}+8\theta\rho(p,q)^{2}\right)+2\mu\mu_{2}\theta\\ &+\mu\kappa\max(|p|,|q|)|x-y|+\mu\eta|p-q|\\ &-\left(|p|-|q|\right)\left(k(u)+2\max\limits_{x\in\Omega}f(x)\right)\\ \leq&~\mu\mu_{1}\left(\frac{8}{27}\zeta^{2}|x-y|^{2}+8\theta\rho(p,q)^{2}\right)+2\mu\mu_{2}\theta\\ &+\mu\kappa\left(\max\left(|p|,|q|\right)+1\right)|x-y|+\mu\eta|p-q|+C_{1}|p-q|.\end{aligned}\end{gathered}

where $C_{1}=\max\limits_{x\in\Omega}\left(k(u)+2\max\limits_{x\in\Omega}f(x)\right)$ (we must assume $k(u),f(x)$ are bounded). Hence we have

\begin{gathered}\begin{aligned} F(t,x,u,p,X)&-F(t,y,u,q,-Y)\geq\\ &-\max\left\{\frac{8}{27}\zeta^{2}\mu,8\mu\theta,2\mu\theta,\mu\eta+C_{1},\mu\kappa\right\}\left[\mu_{1}\left(|x-y|^{2}+\rho(p,q)^{2}\right)+\mu_{2}\right.\\ &\hskip 158.99377pt+|p-q|+|x-y|\left(\max(|p|,|q|)+1\right)\Big{]}\\ \end{aligned}\end{gathered}

and setting $\omega_{R}=\max\left\{\frac{8}{27}\zeta^{2}\mu,8\mu\theta,2\theta,\eta+C_{1},\mu\kappa\right\}R$ , this is a non-decreasing continuous function, maps $[0,\infty)\rightarrow[0,\infty)$ and $\omega_{R}(0)=0$ as required. We have proven that condition (I7) is satisfied.

References

[1] R. Adams and L. Bischof. Seeded region growing. IEEE Transactions on Pattern Analysis and Machine Intelligence, 16(6):641–647, 1994.
[2] N. Badshah and Ke Chen. Multigrid Method for the Chan-Vese Model in Variational Segmentation. Communications in Computational Physics, 4(2):294–316, 2008.
[3] N. Badshah and Ke Chen. On two multigrid algorithms for modeling variational multiphase image segmentation. IEEE Transactions on Image Processing, 18(5):1097–1106, 2009.
[4] N. Badshah and Ke Chen. Image selective segmentation under geometrical constraints using an active contour approach. Communications in Computational Physics, 7(4):759–778, 2010.
[5] N. Badshah, Ke Chen, H. Ali, and G. Murtaza. A coefficient of variation based image selective segmentation model using active contours. East Asian J. Appl. Math., 2:150–169, 2012.
[6] Xue Bai and Guillermo Sapiro. Geodesic matting: A framework for fast interactive image and video segmentation and matting. International journal of computer vision, 82(2):113–132, 2009.
[7] Dhruv Batra, Adarsh Kowdle, Devi Parikh, Jiebo Luo, and Tsuhan Chen. Interactive co-segmentation of objects in image collections. Springer Science & Business Media, 2011.
[8] Vicent Caselles, Ron Kimmel, and Guillermo Sapiro. Geodesic Active Contours. International Journal of Computer Vision, 22(1):61–79, 1997.
[9] Francine Catté, Pierre-Louis Lions, Jean-Michel Morel, and Tomeu Coll. Image selective smoothing and edge detection by nonlinear diffusion. SIAM Journal on Numerical analysis, 29(1):182–193, 1992.
[10] Tony F. Chan, Selim Esedoglu, and Mila Nikolova. Algorithms for Finding Global Minimizers of Image Segmentation and Denoising Models. SIAM Journal on Applied Mathematics, 66(5):1632–1648, 2006.
[11] Tony F. Chan and Luminita A. Vese. Active contours without edges. IEEE Transactions on Image Processing, 10(2):266–277, 2001.
[12] M. G. Crandall, H. Ishii, and P.-L. Lions. User’s guide to viscosity solutions of second order partial differential equations. ArXiv Mathematics e-prints, June 1992.
[13] J. Duan, W. Lu, Z. Pan, and L. Bai. New second order mumford-shah model based on ${\Gamma}$ -convergence approximation for image processing. Infrared Physics and Technology, 76:641–647, 2016.
[14] D. G. Gordeziani and G.V. Meladze. Simulation of the third boundary value problem for multidimensional parabolic equations in an arbitrary domain by one-dimensional equations. USSR Computational Mathematics and Mathematical Physics, 14:249–253, 1974.
[15] Christian Gout and Carole Le Guyader. Segmentation of complex geophysical structures with well data. Computational Geosciences, 10(4):361–372, 2006.
[16] Christian Gout, Carole Le Guyader, and Luminita Vese. Image segmentation under interpolation conditions. Preprint CAM-IPAM, page 44 pages, 2003.
[17] Christian Gout, Carole Le Guyader, and Luminita A. Vese. Segmentation under geometrical conditions using geodesic active contours and interpolation using level set methods. Numerical Algorithms, 39(1-3):155–173, 2005.
[18] Laurence Guillot and Maïtine Bergounioux. Existence and uniqueness results for the gradient vector flow and geodesic active contours mixed model. arXiv preprint math/0702255, 2007.
[19] Jia He, Chang-Su Kim, and C-C Jay Kuo. Interactive Segmentation Techniques: Algorithms and Performance Evaluation. Springer Science & Business Media, 2013.
[20] Hitoshi Ishii and Moto-Hiko Sato. Nonlinear oblique derivative problems for singular degenerate parabolic equations on a general domain. Nonlinear Analysis: Theory, Methods and Applications, 57(7):1077 – 1098, 2004.
[21] Paul Jaccard. The distribution of the flora in the alpine zone.1. New Phytologist, 11(2):37–50, 1912.
[22] Michael Kass, Andrew Witkin, and Demetri Terzopoulos. Snakes: Active contour models. International Journal of Computer Vision, 1(4):321–331, 1988.
[23] M. Klodt, F. Steinbrücker, and Daniel Cremers. Moment Constraints in Convex Optimization for Segmentation and Tracking. Advanced Topics in Computer Vision, pages 1–29, 2013.
[24] Carole Le Guyader. Imagerie Mathématique: segmentation sous contraintes géométriques Théorie et Applications. PhD thesis, INSA de Rouen, 2004.
[25] Carole Le Guyader and Christian Gout. Geodesic active contour under geometrical conditions: Theory and 3D applications. Numerical Algorithms, 48(1-3):105–133, 2008.
[26] Chunxiao Liu, Michael Kwok-Po Ng, and Tieyong Zeng. Weighted variational model for selective image segmentation with application to medical images. Pattern Recognition, pages –, 2017.
[27] T. Lu, P. Neittaanmäki, and X-C. Tai. A parallel splitting up method and its application to navier-stokes equations. Applied Mathematics Letters, 4(2):25 – 29, 1991.
[28] Calvin R Maurer, Rensheng Qi, and Vijay Raghavan. A linear time algorithm for computing exact euclidean distance transforms of binary images in arbitrary dimensions. IEEE Transactions on Pattern Analysis and Machine Intelligence, 25(2):265–270, 2003.
[29] D. Mumford and J. Shah. Optimal approximation of piecewise smooth functions ans associated variational problems. Commu. Pure and Applied Mathematics, 42:577–685, 1989.
[30] Thi Nhat Anh Nguyen, Jianfei Cai, Juyong Zhang, and Jianmin Zheng. Robust interactive image segmentation using convex active contours. IEEE transactions on image processing : a publication of the IEEE Signal Processing Society, 21(8):3734–43, 2012.
[31] S. Osher and J. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79(1):12–49, 1988.
[32] S. M. Patel and J. N. Dharwa. Accurate detection of abnormal tissues of brain MR images using hybrid clustering and deep learning based classifier methods of image segmentation. International Journal of Innovative Research in Science, Engineering and Technology, 6, 2017.
[33] Pietro Perona and Jitendra Malik. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on pattern analysis and machine intelligence, 12(7):629–639, 1990.
[34] Alexis Protiere and Guillermo Sapiro. Interactive image segmentation via adaptive weighted distances. IEEE Transactions on image Processing, 16(4):1046–1057, 2007.
[35] Lavdie Rada. Variational models and numerical algorithms for selective image segmentation. PhD thesis, University of Liverpool, 2013.
[36] Lavdie Rada and Ke Chen. A new variational model with dual level set functions for selective segmentation. Communications in Computational Physics, 12(1):261–283, 2012.
[37] Lavdie Rada and Ke Chen. Improved Selective Segmentation Model Using One Level-Set. Journal of Algorithms and Computational Technology, 7(4):509–540, 2013.
[38] James A Sethian. A fast marching level set method for monotonically advancing fronts. Proceedings of the National Academy of Sciences, 93(4):1591–1595, 1996.
[39] Jack Spencer and Ke Chen. A Convex and Selective Variational Model for Image Segmentation. Communications in Mathematical Sciences, 13(6):1453–1472, 2015.
[40] L. Vincent and P. Soille. Watersheds in digital spaces: an efficient algorithm based on immersion simulations. IEEE Trans. Pattern Analysis Machine Intell., 13(6):583–598, 1991.
[41] G. T. Wang, W. Q. Li, S. Ourselin, and T. Vercauteren. Automatic brain tumor segmentation using cascaded anisotropic convolutional neural networks. https://arxiv.org/pdf/1709.00382.pdf, preprint, 2017.
[42] G. T. Wang, W. Q. Li, M. A. Zuluaga, R. Pratt, P. A. Patel, M. Aertsen, T. Doel, A. L. David, J. Deprest, S. Ourselin, and T. Vercauteren. Interactive medical image segmentation using deep learning with image-specific fine-tuning. CoRR, abs/1710.04043, 2017.
[43] J. Weickert, B. M. Ter Haar Romeny, and M. A. Viergever. Efficient and reliable schemes for nonlinear diffusion filtering. IEEE Transactions on Image Processing, 7(3):398–410, 1998.
[44] J. P. Zhang, Ke Chen, and B. Yu. A 3D multi-grid algorithm for the CV model of variational image segmentation. International Journal of Computer Mathematics, 89(2):160–189, 2012.
[45] Feng Zhao and Xianghua Xie. An overview of interactive medical image segmentation. Annals of the BMVA, 2013(7):1–22, 2013.
[46] W. Zhu, X.-C. Tai, and T. F. Chan. Image segmentation using euler’s elastica as the regularization. Journal of Scientific Computing, 57:414–438, 2013.

A Convex Geodesic Selective Model for Image Segmentation

Abstract

1. Introduction

2. Review of Variational Segmentation Models

2.1. Global Segmentation

2.2. Selective Segmentation

3. Proposed Convex Geodesic Selective Model

3.1. Computing the Geodesic Distance Term 𝒟M​(x,y)\mathcal{D}_{M}(x,y)

3.2. Comparing Euclidean and Geodesic Distance Terms

Remark 1 (Computational Complexity.).

3.3. Improvements to Geodesic Distance Term

3.4. The new model and its Euler-Lagrange equation

4. An additive operator splitting algorithm

5. Existence and Uniqueness of the Viscosity Solution

Theorem 2 (Theorem 3.1 [20]).

5.1. Existence and uniqueness for the Geodesic Model

Remark 3.

Theorem 4 (Theory for the Geodesic Model).

5.2. Generalisation to other related models

Corollary 5.

Remark 6.

Remark 7.

6. Numerical Results

Remark 8 (A note on M5 and M6).

Remark 9.

7. Conclusions

Acknowledgements

Appendix A — Proof that Condition (I7) Holds in Theorem 4

Remark 10.

References

3.1. Computing the Geodesic Distance Term $\mathcal{D}_{M}(x,y)$