A Convex Geodesic Selective Model for Image Segmentation

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.


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 M = {(x i , y i ) ∈ Ω, 1 ≤ i ≤ k} where Ω ⊂ R 2 is the image domain [4,5,17].Nguyen et al. [30] considered marker sets M and 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 M. This was enforced with a distance function which is zero at 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 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 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 M as a distance penalty term, however we propose replacing this with the edge-weighted geodesic distance from 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.

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 Ω ⊂ R 2 .

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 where the foreground Ω 1 is the subdomain to be segmented, the background Ω 2 = Ω\Ω 1 and µ, λ 1 , λ 2 are fixed non-negative parameters.The values c 1 and c 2 are the average intensities of z(x, y) inside Ω 1 and Ω 2 respectively.We use a level set function an idea popularised by Osher and Sethian [31]) and reformulate (1) as with H ε (φ) a smoothed Heaviside function such as H ε (φ) = 1  2 + 1 π arctan( φ ε ) for some ε, we set ε = 1 throughout.We solve this in two stages, first with φ fixed we minimise F CV with respect to c 1 and c 2 , obtaining and secondly, with c 1 and c 2 fixed we minimise (2) with respect to φ.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.

Selective Segmentation
Selective segmentation models make use of user input, i.e. a marker set M of points near the object or objects to be segmented.Let M = {(x i , y i ) ∈ Ω, 1 ≤ i ≤ k} be such a marker set.The aim of selective segmentation is to design an energy functional where the segmentation contour Γ is close to the points of 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 Γ g(|∇z(x, y)|)dΓ.
The term g(|∇z(x, y)|) is an edge detector, one example is g(s) = 1/(1 + βs 2 ) with β a tuning parameter.It is common to smooth the image with a Gaussian filter G σ where σ is the kernel size, i.e. use g(|∇ (G σ * z(x, y)) |) 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 D(x, y) into this integral, i.e. the integrand is D(x, y)g(|∇z|).The distance term is a penalty on the distance from M, this model encourages the contour to be near to the set 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 where µ, λ 1 , λ 2 , γ are fixed non-negative parameters.There is freedom in choosing the distance term D(x, y), see [37] for some examples.A 1 is the area of the polygon formed from the points of M and A 2 = |Ω| − 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 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 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, α ∈ [0, 1] and minimisation is constrained, requiring φ ∈ [0, 1], so F NG (φ) 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 where µ, λ 1 , λ 2 , θ are fixed non-negative parameters.Note that the regulariser of this model differs from the Rada-Chen model (4) as the distance function D(x, y) has been separated from the edge detector term and is now a standalone penalty term D E (x, y).The authors use normalised Euclidean distance D E (x, y) from the marker set 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 (φ, c 1 , c 2 ) have the same stationary solutions as for with the minimisation constrained to u ∈ [0, 1].This is a constrained convex minimisation which can be reformulated to an unconstrained minimisation using an exact penalty term ν(u) := max{0, 2|u − 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 ν ε (u) to ν(u) given by and perform the unconstrained minimisation of 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.[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

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 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 where 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 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.

Computing the Geodesic Distance Term D M (x, y)
The geodesic distance from the marker set M is given by D M (x, y) = 0 for (x, y) ∈ M and  where f (x, y) is defined later on with respect to the image contents.
If f (x, y) ≡ 1 (i.e.|∇D 0 M (x, y)| = 1) then the distance penalty D M (x, y) is simply the normalised Euclidean distance 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 this gives us the desired property that in areas where |∇z(x, y)| ≈ 0, the distance function increases by some small ε D ; here image z(x, y) is scaled to [0, 1].At edges, |∇z(x, y)| is large and the geodesic distance increases here.We set value of β G = 1000 and ε 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.

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 θ, 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 intensitythis penalty can also be quite high if our marker set is small and doesn't cover the whole object.

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 D M (x, y) is O(N log(N)) for an image with N pixels.This is is only marginally more complex than computing the Euclidean distance which has O(N) complexity [28].

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 |∇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 M with low penalty.As the geodesic distance only uses marker set M for its initial condition (see (11)), this can result in objects far from 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 where G σ is a Gaussian convolution with standard deviation σ.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: where μ, ι are non-negative parameters (we fix throughout μ = 10 −3 , ι = 5 × 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 where We update all pixels once per iteration and solve the PDE in (11) with f (x, y) replaced by 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.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 D E (x, y) incorporated somehow.
Our solution is to modify the term f 1 (x, y) from ( 16) to In Figure 5 the effect of this is clear, as ϑ increases, the distance function resembles the Euclidean distance more.We use value ϑ = 10 −1 in all experiments as it adds a reasonable penalty to pixels far from the marker set.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 AM.We propose to use a geodesic distance map from  the set AM denoted by D AM (x, y) which penalises pixels near to the set AM and doesn't add any penalty to those far away.We could naïvely choose D AM (x, y) = 1 − DGAM (x, y) where DGAM (x, y) is the normalised geodesic distance from AM. However this puts a large penalty on those pixels inside the object we actually want to segment (as DGAM (x, y) to those pixels is small).To avoid this problem, we propose the following anti-marker distance term where α is a tuning parameter.We choose α = 200 throughout.This distance term ensures rapid decay of the penalty away from the set 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 D AM (x, y) decays rapidly from AM, we do require that the anti-marker set be close to the blurred edge and away from the object we desire to segment.

The new model and its Euler-Lagrange equation
The Proposed Geodesic Model.Putting the above 3 ingredients together, we propose the model min where D G (x, y) = (D M (x, y) + D AM (x, y)) /2 and D M (x, y) is the geodesic distance from the marker set M. We compute 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 and the minimisation with respect to u (with c 1 and c 2 fixed) gives the PDE in Ω, where we replace |∇u| with |∇u| ε 2 = u 2 x + u 2 y + ε 2 to avoid zero denominator; we choose ε 2 = 10 −6 throughout.We also have Neumann boundary conditions ∂u ∂n = 0 on ∂Ω where 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.

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 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 and after discretisation, we can rewrite this as [43] For notational convenience we write G = G(u).The matrix A 1 (u) can be obtained as follows and similarly to [37,39], for the half points in G we take the average of the surrounding pixels, e.g.
. 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 ν ε (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 ν ε (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 ν ε (u) around u = 0 and u = 1 and group the terms into the constant and linear components in u, we can respectively write ).We actually find that b 0 (ε) = b 1 (ε) and denote the linear term as b from now on.Therefore, for a change in u of δu around u = 0 and u = 1, we can approximate the change in ν ε (u) by b • δu.To focus on the jumps, define the interval in which ν ε (u) jumps as and refine the linear function by bk i,j = b, u k i,j ∈ I ζ 0, else.
Using these we can now offset the change in ν ε (u k ) by changing the formulation (21) to or in AOS form u k+1 = u k + τµ∇ • (G(u k )∇u k+1 ) − τα bk u k+1 + τα bk u k − f k which, following the derivation in [39], can be reformulated as where Bk = diag(τα bk ).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 where the matrix Q 2 does have unit row sum, however the matrix is not always symmetrical.We can guarantee convergence for ζ = 0.5 (in which case Q 2 must be symmetrical) but we desire to use a small ζ to give a small interval I ζ .We find experimentally that convergence is achieved for any small value of ζ, 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 ζ and I + Bk is a matrix with all diagonal entries 1 + τα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 ζ = 0.01.

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 Ω is a bounded domain with C 1 boundary.
Update u k to u k+1 using the AOS scheme (22).end for u * ← u k .
From the work of [12,20], we have the following Theorem for analysing the solution of a partial differential equation of the form F(x, u, Du, D 2 u) = 0 where is the set of n × 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 := x for the vector of a general point in R n .Theorem 2 (Theorem 3.1 [20]).Assume that the following conditions (C1)-( C2) and (I1)-(I7) hold.Then for each u 0 ∈ C(Ω) there is a unique viscosity solution u ∈ C([0, T) × Ω) of ( 23) and ( 24) satisfying (25).

Conditions (I1)-(I7).
Assume Ω is a bounded domain in R n with C 1 boundary.
Here F * and F * denote, respectively, the upper and lower semi-continuous envelopes of F, which are defined on [0 , where C 1,1 is the Hölder functional space.
(I5) For each x ∈ R n the function p → B(x, p) is positively homogeneous of degree one in p, i.e.B(x, λp) = λB(x, p) for all λ ≥ 0 and p ∈ R n \{0}.
Here n(x) denotes the unit outward normal vector of Ω at x ∈ ∂Ω.

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 |∇u| ε 2 , obtaining the nonlinear PDE where G(x, ∇z) = g(|∇z(x, y)|).We can rewrite this as where Theorem 4 (Theory for the Geodesic Model).The parabolic PDE ∂u ∂t + F(t, x, u, Du, D 2 u) = 0 with u 0 = u(0, x) ∈ C(Ω), F as defined in (28) and Neumann boundary conditions has a unique solution u = u(t, x) in C([0, T) × Ω).
• Spencer-Chen (Convex) [39] Then if we define a PDE of the general form we have a unique solution u ∈ C([0, T) × Ω) 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 ν ε (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. 2 Remark 6. Theorem 4 and Corollary 5 also generalise to cases where G(x) = 1 1+β|∇z| 2 and to G(x) = D(x)g(|∇z|) where 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ν Ω H ε (u) dΩ − A 1 .

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]; 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 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 for µ, λ 1 , λ 2 non-negative fixed parameters, α and ν ε (u) as defined in (7) and ω 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.
Quantitative Comparisons.To measure the quality of a segmentation, we use the Tanimoto Coefficient (TC) (or Jaccard Coefficient [21]) defined by where GT is the 'ground truth' segmentation and ũ 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.
Following [10] we 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 ω enforces these terms at all edges in the image and near 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 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.(i) (ii) (iii) (iv) 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 µ, λ 1 , λ 2 , θ and ε 2 .In all tests we set µ = 1, which is simply a rescaling of the other parameters, and we set λ = λ 1 = λ 2 .In the first example, in Figure 12, we compare the TC value for various λ and θ values for segmentation of a bone in a knee scan.We see that the segmentation is very good for a larger range of θ and λ 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 ε 2 values with fixed parameters λ = 2 and θ = 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 ε 2 tested, and in the case of the knee segmentation the results are almost identical for any ε 2 < 10 −6 , above which the quality slowly deteriorates.
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 15, 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.

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 ε 2 Knee Segmentation (Figure 12) Circle Segmentation (Figure 13) 10  (i)  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.
discussions at the early stages of this work.The second author is grateful to the UK EPSRC for the grant EP/N014499/1.Therefore, we can write trace(A(x, p)X) + trace(A(y, q)Y) = (XB(p)e i , B(p)e i ) + (YB(q)e i , B(q)e i )

Note.
In the Geodesic Model we fix G(x) = g(|∇z|).Therefore, assuming G(x) and 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. ζ.

Figure 1 :
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.

Figure 2 :Figure 3 :
Figure 2: Examples of images showing the problems discussed and the resulting geodesic distance maps.Column 1shows the lack of robustness to noise, column 2 shows that outside the patient we have unreasonably low distance penalty, column 3 shows how the blurred edge under the aorta leads to the distance term being very low throughout the heart.

Figure 5 : 3 . 3 :
Figure 5: Displayed is D M (x, y) using f 2 (x, y) for various ϑ values.The marker set is the same as that used in Figure 3.

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.

Figure 9 :
Figure 9: Test 1 setting: (i) Image 1; (ii) Image 1 with marker and anti-marker set shown in green and pink respectively; (iii) Test Image 2; (iv) Image 2 with marker set shown.

Figure 10 :
Figure 10: Visual comparison of M1 -M7 results for Test Image 1. M1 segmented part of the object, M2 -M4 failed to segment the object, M5 gave a reasonable result (though not accurate) and, M6 and M7 correctly segmented the object.

Figure 11 :
Figure 11: Visual comparison of M1 -M7 results for Test Image 2. M1 segmented part of the object, M2 -M4 failed to segment the object, M5, M6 and M7 correctly segmented the object.

( a )
Original Image.(b) Ground Truth Segmentation.(c) M3 TC values for various λ and θ values.(d) M6 TC values for various λ and θ values.(e) M7 TC values for various λ and θ values.

Figure 14 :
Figure 14: Three further test results obtained using our Geodesic Model M7, all with parameters θ = 5, λ = 5.The first row shows the original image with the marker set (plus anti-marker set), the second row the final segmentation result and the final row shows the residual history.

Figure 16 :
Figure 16: Test 4 on M7's marker set (θ = 5, λ = 3).Row 1 shows the original image with 3 marker points, the normalised geodesic distance map and the final segmentation result.Row 2 shows the original image with 1 marker point, the normalised geodesic distance map and the final segmentation result.Clearly the second and third columns are the same for different marker points.Thus M7 is robust.

= µ 1
A -Proof that Condition (I7) Holds inTheorem 4    Using the assumption in(26), we write(Xr, r) + (Ys, s) = r T Xr + s T Ys ≤ µ 1 r T s T I |r − s| 2 + µ 2 |r| 2 + |s| 2 .Note that matrix A from (29) is a real symmetric matrix and decomposes asA = 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(Xr, r) = r T Xr = ∑ i (Be i ) T X(Be i ) = ∑ i e T i BT XBe i = trace(B T XB) = trace(A(x, p)X).

Table 1 :
The Tanimoto Coeffcient for various ε 2 values, segmenting the images in Figures12 and 13.