Eikonal Regionbased Active Contours for Image Segmentation
Abstract
The minimal path model based on the Eikonal partial differential equation (PDE) has served as a fundamental tool for the applications of image segmentation and boundary detection in the passed three decades. However, the existing minimal pathsbased image segmentation approaches commonly rely on the image boundary features, potentially limiting their performance in some situations. In this paper, we introduce a new variational image segmentation model based on the minimal path framework and the Eikonal PDE, where the regionbased functional that defines the homogeneity criteria can be taken into account for estimating the associated geodesic paths. This is done by establishing a geodesic curve interpretation to the regionbased active contour evolution problem. The image segmentation processing is carried out in an iterative manner in our approach. A crucial ingredient in each iteration is to construct an asymmetric Randers geodesic metric using a sufficiently small vector field, such that a set of geodesic paths can be tracked from the geodesic distance map which is the solution to an Eikonal PDE. The object boundary can be delineated by the concatenation of the final geodesic paths. We invoke the Finsler variant of the fast marching method to estimate the geodesic distance map, yielding an efficient implementation of the proposed Eikonal regionbased active contour model. Experimental results on both of the synthetic and real images exhibit that our model indeed achieves encouraging segmentation performance.
Keywords:
Regionbased active contours minimal path Randers metric image segmentation Finsler variant of the fast marching method Eikonal partial differential equation∎
1 Introduction
Image segmentation is a crucial and fundamental task in the applications of image analysis, computer vision and medical imaging. Among existing approaches for image segmentation, active contour models established upon a solid theoretical background provide an effective way to exploit the image edgebased features and a broad variety of regionbased homogeneity measures, hence are capable of obtaining desired segmentation results in various situations. The active contour models provide a natural way to exploit continuous curves to delineate image structures of interest. These curves can be derived through the minimization processing of a problemdependent energy functional.
The original active contour model proposed by kass1988snakes allows to combine the image data and the curve regularization for boundary detection and image segmentation. Since then, a broad variety of edgebased active contour models have been considered to address, for example, the initialization and convergence problems arisen by the original snakes model. The classical geodesic active contour model (caselles1997geodesic; yezzi1997geometric) is a geometric variant of the original snakes model. The energy functional for these models are derived by simplifying the energy of (kass1988snakes). An important merit for the new energy invoked by (caselles1997geodesic; yezzi1997geometric) is that it is independent to the curve parameterization yielding a geometry quantity. kimmel2003regularized proposed a new geometric active contour model based on the alignment measures to take the advantages from the curve directionality. melonakos2008finsler established a general geodesic active contour functional using Finsler metrics. In order to enlarge the capture range of the active contours, many approaches take into account an external force defined through a dense vector field to build the respective curve evolution equation. Significant examples for these external forces include the approaches derived from the Chamfer distance (cohen1993finite), the dense extension of the image gradients (xu1998snakes; li2007active), and the physical theorems (jalba2004cpm; xie2008mac). paragios2004gradient provides a way to combine an external force based on the gradient vector flow (xu1998snakes) and the geodesic active contours flow, which is able to address the convergence problem in some extent. Instead of solving the curve evolution problem in a continuous setting, boykov2003computing addressed the problem of finding geodesic paths in a discrete domain based on the graph cut theorem and the CauchyCrofton formula. mishra2011decoupled proposed an efficient active contour model posed on prescribed graphs, where the optimal contours are searched by the Viterbi algorithm.
On the other hand, the regionbased active contour models in common differ them to the edgebased approaches by an additional term of the homogeneity measure within each region. A pioneering approach is the original MumfordShah functional proposed by (mumford1989optimal), which involves a similarity measure estimated from a piecewise smooth data fitting function and a curve length regularization term. Finding the solutions to the original MumfordShah functional has motivated a series of approaches such as (tsai2001curve; vese2002multiphase; bresson2007fast). Basically, the MumfordShah functional assumes that the image data can be approximated by a piecewise smooth function. Such an assumption is not always reasonable especially when the gray level distributions of the target are complicated. Many efforts have been dedicated to address this problem, such as the Bayes generalizations of the MumfordShah functional which exploit adequate probability density functions to establish the regionbased homogeneity criteria (zhu1996region; brox2009local; paragios2002geodesic) and the local cluster method (lialevel2011) which takes into account an external bias field to cope with the gray level inhomogeneities. The parametric statistical active contour models such as the region competition model (zhu1996region) assume that the image gray levels admit a prescribed parameterdependent probability distribution in each region, where the parametric Gaussian distributions are often utilized to define the regionbased similarity measure. The class of the nonparametric active contour models (ni2009local; michailovich2007image) using the probability distribution functions of the image gray levels is able to remove the dependency on these statistical parameters. Both parametric and nonparametric models come with advantages and disadvantages and we refer to (cremers2007review) for a more detailed review of the statistical active contour models.
In essence, the piecewise constant active contour models (chan2001active; cohen1997avoiding) are built upon a simplified and practical variant of the full MumfordShah functional. These variant models impose a homogeneous prior on the image data within each region. In general, the similarity measure used in (chan2001active; cohen1997avoiding) is constructed via the type distance between the image gray levels and the fitting data. In contrast, the type distance estimation was taken into account in (kimmel2003fast; bresson2007fast; jung2017piecewise) to define the regionbased homogeneity penalization terms. The active contour approaches based on the pairwise similarity or dissimilarity criteria (sumengen2006graph; bertelli2008variational; jung2012nonlocal) are capable of handling more complicated segmentation situations but requiring high computation complexity to estimate the associated gradient descent flows. sundaramoorthi2007sobolev; sundaramoorthi2009new defined the associated gradient descent flows for curve evolution through the Sobolevtype inner product in a infinite dimensional Riemannian manifold, leading to a more robust curve evolution scheme.
The active contour models mentioned above are not exhaustive, other interesting approaches include the hybrid active contour models simultaneously integrating the regionbased homogeneity penalization and the edgebased features (paragios2002geodesic; kimmel2003regularized; zach2009globally), the shape priorbased active contour models (cremers2002diffusion; bresson2006variational), and the multiregion/multilabel segmentation models (dubrovina2015multi; chambolle2012convex; bae2011global; brox2006level).
1.1 Minimal Pathbased Image Segmentation Models
In the basic formulation of the geodesic active contour approach (caselles1997geodesic; yezzi1997geometric), the minimization of the energy, which is referred to as the weighted curve length, is implemented by the curve evolution method relying on the gradient descent flow. Unfortunately, there is no guarantee for this scheme to find the globally minimizing curve. In order to address this problem, cohen1997global introduced a new minimal path model based on the Eikonal PDE, where a minimal geodesic path between two prescribed points can be recovered from the respective geodesic distance map.
Most of the existing minimal pathbased image segmentation models, also including the proposed one, aims to exploit geodesic paths derived from a geodesic distance map to depict the boundaries of the target. For instances, the saddle points detected from the geodesic distance map are often served as landmark points to track geodesic paths, the concatenation of which can be used to depict the object contour (cohen1997global; cohen2001multiple; mille2015combination). appleton2005globally proposed a method which can generate a closed curve that globally minimizes the weighted curve length, providing that this optimal curve surrounds a prescribed point. A keypoint detection method was introduced in (benmansour2009fast), which can be efficiently applied to the tasks of boundary detection and closed contour extraction. It only requires a single source point for initialization and is able to iteratively add new source points. In each iteration, a new keypoint can be detected and a geodesic path linking one of the existing keypoints to the latest one can be tracked. The authors in (chen2017global) proposed to connect a set of prescribed points of unknown orders through the curvaturepenalized geodesic paths. These paths are tracked in an orientationlifted space and are derived using the image boundary features extracted by a steerable filter.
Actually, these minimal pathsbased image segmentation approaches commonly define the datadriven geodesic metrics using boundary information such as image gradients. As a consequence, the corresponding minimal paths, either isotropic or anisotropic, essentially fall into the image edgedriven limitation, despite their relevant applications for image segmentation in conjunction with regionbased homogeneity criteria, see (mille2015combination; appia2011active) for such examples. However, the minimal paths derived from the image boundary features (e.g. image gradients) sometimes may fail to depict the boundaries of the target, especially for the cases that the target object involves complex structures or the image gradients are not reliable. In this paper, we propose a new minimal path model for regionbased active contours evolution and image segmentation. The corresponding geodesic metric has a Randers form constructed by the regionbased homogeneity term and/or the anisotropic boundary features.
In Fig. 1, we take the combination of paths model proposed by mille2015combination as an example to illustrate the limitations of using edgedriven minimal paths for segmentation. Basically, the paths combination model is an interactive segmentation method for which a set of vertices on the target boundary should be provided. In Fig. 1a, we show a natural image with four vertices indicated by red dots. One of the crucial points for the paths combination model lies at the construction of a potential used for estimating the geodesic distance maps. In Fig. 1b, we show the segmentation contour (visualized by blue lines) from the paths combination model. One can point out that some portions of the desired boundary are not accurately followed, since in these regions the boundaries are not well characterized by edge appearance features. Finally, the paths tracking using the proposed model, which blend the benefits from both region and edgebased information, are able to follow the desired boundary, as depicted in Fig. 1c.
1.2 Overview of the Proposed Framework
In summary, the main contribution in this paper is threefold:

Firstly, we establish the relationship between the general regionbased active contours model and the minimal geodesic path framework in conjunction with the Eikonal PDE, utilizing an anisotropic and asymmetric Randers metric. The proposed approach, regarded as the Eikonal regionbased active contour model, allows to blend the benefits from both of the regionbased homogeneity penalty term and the image gradient features including both the magnitude and the directions.

As a second contribution, we present the analysis for the existence of the aforementioned Randers metric, relying on the existence of a sufficiently small vector field derived from the regionbased homogeneity penalization term. A crucial point underlying this analysis is the search of the solution to a minimization problem that is constrained by a divergence equation. Moreover, for numerical considerations, we also propose two methods to efficiently estimate the approximate solution to that PDEconstrained minimization problem.

Finally, we propose two methods to exploit a set of piecewise geodesic paths associated to the proposed Randers metric for shape deformation and image segmentation. Both methods rely on a set of ordered vertices and each pair of successive vertices defines the source and end points for each geodesic path. The set of these ordered vertices can be either specified by the user or can be generated through a point sampling scheme.
The preliminary short versions of this work were first presented in the conferences (chen2016finsler; chen2017anisotropic), upon which more contributions such as the theoretical analysis and experimental results have been added.
1.3 Paper Outline
In Section 2, we briefly introduce the minimal path models and the corresponding Eikonal PDEs and gradient descent ordinary equations, including the Riemannian cases and the asymmetric Randers case. The preliminary on the hybrid active contour energy functionals and on the iterative minimization scheme for these functionals are presented in Section 3. In Section 4, we transform the hybrid active contour energy functional to a weighted curve length associated to a Randers metric. This transformation relies on the solution to a divergence equation. In Section 5, we introduce two methods to exploit the constructed Randers metrics and the associated minimal paths for the application of image segmentation. The numerical implementation details and the experimental results are presented in Sections 6 and 7, respectively. The conclusion can be found in Section 8.
2 Background on Minimal Path Models
2.1 Riemannian Minimal Paths
The original minimal path model introduced by cohen1997global provides an efficient way to determine the globally minimizing curves over an open and bounded domain . In this model, the weighted length of a Lipschitz continuous curve can be measured independently to the parameterization by
(1) 
where is the firstorder derivative of the curve with respect to the parameter , is the Euclidean norm of , and the function is a potential which is supposed to take low values around image features of interest.
Let be the set of all Lipschitz continuous curves (of finite length) . A minimal path (or a geodesic path) linking the given source point to a target point , is defined as a curve that globally minimizes to the weighted curve length . In the context of image segmentation, the potential is often defined as a deceasing function of the image gradient magnitude (caselles1997geodesic; cohen1997global). Nevertheless, a geodesic path favours to pass through the positions with high values of gradient magnitude, which usually correspond to the object boundaries. The weighted curve length defined in Eq. (1) in essence is measured through an isotropic Riemannian metric, yielding that the length is independent to the directions . In contrast, the minimal path models based on the general Riemannian metrics (bougleux2008anisotropic; jbabdi2008accurate) are capable of taking the advantages from the datadriven anisotropy features (which can be derived from image gradients) when tracking the geodesic paths.
Denote by the set of all the positive symmetric definite matrices of size . A general Riemannian metric can be constructed via a tensor field as follows:
(2) 
where stands for the Euclidean scalar product on , i.e. .
The weighted length of any curve associated to a Riemannian metric can be measured by
(3) 
The length is a simplification of the general case . If the tensor field is proportional to the identity matrix of size , the Riemannian metric gets to be the isotropic Riemannian case:
As a consequence, the weighted curve length has the same form to the length measured through .
The minimization of the weighted curve length (2.1) with respect to a fixed source point is regraded as a way of searching for a geodesic distance map , where for each point the value of is equal to the minimal curve length between and
(4) 
where is a geodesic path connecting the source point to the target point .
The geodesic distance map admits the anisotropic Eikonal PDE, which is referred to as the static HamiltonianJacobBellman PDE
(5) 
where is the standard Euclidean gradient of .
The geodesic path can be tracked by means of the solution to a gradient descent ordinary differential equation (ODE) on the geodesic distance map . Let be a geodesic path parameterized by its arclength. One can solve the following ODE
(6) 
until the source point is reached, where is the arclength parameter of the path . Finally, the desired geodesic path can be generated by reparameterizing and reversing the geodesic path .
2.2 Randers Minimal Paths
The minimal path model, which is not restricted to the Riemannian scope, can be generalized to the asymmetric Finsler case, where a Finsler metric can be naturally characterized through asymmetric norm (bao2012introduction; melonakos2008finsler). The utilization of Finsler metrics allows to take the asymmetry enhancement into consideration (melonakos2008finsler; zach2009globally) for the computation of minimal paths, where the asymmetry information may lead to promising segmentation performance in some tasks, as discussed in (zach2009globally).
In the scenario of image segmentation and boundary detection, the Randers metric (randers1941asymmetrical) is a widely considered case of the Finsler metrics. Specifically, a Randers metric can be constructed through a tensor field (see Eq. (2)) and a sufficiently small vector field as follows:
(7) 
Basically, the Randers metric should follow the positive definiteness condition, i.e. the inequality holds for any point and any vector . The following Theorem 2.1 yields a simple characterization of the positive definiteness condition of the Randers metric .
Theorem 2.1 (Positive Definiteness)
A Randers metric defined in Eq. (7) is positively defined if and only if
(8) 
In particular, if the tensor field is set to , where is an identity matrix, then the inequality (8) gets to be
(9) 
Proof
For any point and any vector , if , we clearly have Ley us choose , which yields that for all points
which immediately implies .
Conversely, one has by CauchySchwartz inequality
(10) 
Hence for any vector , one has
(11) 
yielding that
∎
The weighted curve length measured through a Randers metric can be defined by
(12) 
Given a source point , the geodesic distance map with respect to a Randers metric satisfies the following Eikonal PDE
(13) 
The equation (13) can be also regarded as the general Eikonal PDE with respect to any Finsler metric . In particular, for a Randers metric of the form (7), the Eikonal PDE (13) gets to be the following equivalent form (mirebeau2019riemannian; chen2018fast)
(14) 
One can point out that when the vector field , the Eikonal PDE (14) becomes the Riemannian Eikonal equation formulated in Eq. (5).
The gradient descent ODE for tracking a geodesic path associated to a Randers metric is defined by
(15) 
where is a vectorvalued function defined as follows
(16) 
Again, the target geodesic path can be generated by reparameterizing and reversing the solution to the gradient descent ODE (15).
Minimal Path Examples.
In Fig. 2, we illustrate three minimal path computation examples respectively associated to the isotropic Riemannian, anisotropic Riemannian and Randers metrics. The isotropic metric in Fig. 2b is endowed by a potential set as a constant function of , such that the geodesic path (white line) between the source (red dot) and end points (black dot) is a straight segment. The construction for the anisotropic Riemannian metric and the Randers metric rely on a unit vector field . Each vector is tangent to the circle
In this case, the corresponding geodesic path favours to follow the directions as shown in Fig. 2c. The Randers metric used in Fig. 2d is constructed as follows:
One can observe that the obtained Randers geodesic path tries to travel along the directions opposite to these arrows as much as possible.
2.3 Overview on the Combination of Paths Model
The paths combination model proposed in (mille2015combination) applies edgebased minimal paths for interactive image segmentation. The initialization for this model relies on a set of fixed points for . These points are distributed along the target boundary in a clockwise order. Basically, the goal of the combination of paths model is to search for a simple and closed curve by minimizing an energy functional
(17) 
where is a simple curve and are scalarvalued parameters which control the related importance between the three energy terms. Specifically, the simplicity term measures the amount of the selftangency and selfintersection of the curve . The edgebased term is measured using an edge indicator along the curve , and the regionbased one is set as the Bhattacharyya coefficient (michailovich2007image) of the image gray level or color probability distribution functions inside and outside . We refer to (mille2015combination) for more details on the energy .
One of the crucial steps for the combination of paths model lies at construction of a set of open admissible curves, such that each curve travels from a point and to the successive one . All of the admissible curves between each pair of successive points are tracked from a single geodesic distance map with source points and . This is done by performing the partial fast marching scheme (deschamps2001fast; chen2019minimal), where the front propagation will expand from both source points simultaneously. The scheme leads to a meeting interface of the two propagating fronts associated to the respective source points. The interface is a set made up of all equidistant points to and . One can detect several saddle points from the meeting interface. Each saddle point can generate two geodesic paths from the combined action map and an admissible curve is built using those two geodesic paths by a concatenation operation. In other words, each saddle point between two successive points is supported to be passed by an admissible curve.
The considered geodesic distance map, referred to as a combined action map (mille2015combination), is the viscosity solution to the Riemannian Eikonal PDE (5) by setting , where is a potential which can be computed from the magnitude of the image gradients. In the scenario of image segmentation and boundary detection, the minimal paths from the geodesic distance map which depend only on the image edge information, favour to travel along the region of strong edge appearance features. This sometimes leads to unexpected image segmentation results, especially when the values of image gradients magnitude close to the targets are insufficiently salient. In order to illustrate the drawbacks of the combination of paths model, we show in Fig. 3 the admissible curves between two pairs of successive points. In this figure, the prescribed points are identical to the ones used in Fig. 1a. As depicted in Fig. 3a, the potential takes low values around the image features of interest and high values otherwise. From Figs. 3b and 3c, one can point out that the admissible curves fail to travel along the desired boundaries but along their neighbouring ways, due to the lack of the reliability of images gradients as the indicator of target boundaries. This gives arise to the motivation to induce a new minimal path model whose metrics relying not only the image edges but also the regionbased homogeneity criteria to find suitable image segmentation.
3 Hybrid Active Contour Formalism
In this paper, a shape is defined as an open and bounded subset with a rectifiable boundary. Moreover, we suppose that each shape is measurable and simply connected (without holes). The characteristic function of a shape can be defined by
(18) 
We denote by the set of all functions defined almost everywhere on the domain and taking only the values and , i.e. the characteristic functions of all the subsets of the domain .
3.1 Riemannian Metricmeasured Curve Length for Regularization
In general, a typical regionbased active contour energy functional should be comprised of a regionbased homogeneity penalization term and a regularization term such that
(19) 
where is an arbitrary shape and is a constant which is a parameter determining the importance of . The term is the total variation norm of the function which equals to the perimeter of the boundary . It is also the 1dimensional Haussdoff measure of . bresson2007fast proposed to apply a weighted total variation norm as the regularization term, which in essence corresponds to the weighted curve length (see Eq. (1)) associated with an isotropic metric. In this case, the active contour energy functional involving a weighted total variation norm is extended to the hybrid version (kimmel2003regularized; paragios2002geodesic; sagiv2006integrated; zach2009globally).
The hybrid active contour models are capable of taking advantages of both regionbased and boundarybased image features to enhance the segmentation. In this paper, we take the weighted curve length associated to a Riemannian metric defined in Eq. (2.1) as the regularization term, yielding an active contour energy functional
(20) 
where is a Lipchitz continuous map that parameterizes the boundary .
The tensor field for the Riemannian metricmeasured curve length encompasses both the magnitude and the directions of the image gradient features
3.2 Iterative Scheme for the Minimization of the General Regionbased Active Contour Model
In regionbased active contour models, the associated homogeneity measure is commonly defined through a functional , embedding with the regional similarity or dissimilarity information. For a a given shape , the functional is said to be differentiable at , if and only if there exists a scalarvalued function such that for all of the admissible functions
(21) 
Instead of being defined over the space , the functional often has a natural extension to the linear space , see Appendix A for such examples. Then the function is the gradient of the energy functional at and can be derived from the Gâteaux derivative of along the direction such that
(22) 
Let be a shape that is close to the fixed shape . By choosing , the regionbased functional can be reformulated by its firstorder expansion
(23) 
Indeed the norm stands for the area of the symmetric difference between the shapes and .
For the sake of simplification, we define a new regionbased functional
(24) 
By integrating Eqs. (20), (23) and (24), the active contour energy functional can be reformulated as
(25) 
Iterative Curve Evolution Solution. When working with the shape evolution scheme, the minimization of the active contour energy functional (25) can be carried out by iteratively solving the following two subproblems:

Finding a shape to minimize the energy with respect to a given shape ;

Updating the shape gradient by choosing , i.e. setting that .
Note that the minimization of the active contour energy by repeatedly solving the aforementioned two subproblems is expected to generate a shape to delineate the target object. This shape deformation or shape evolution algorithm will terminate whenever the Hausdorff distance between the two successive shapes is small enough.
Variational image segmentation methods can be backed by two celebrated existence results, for the full piecewise smooth MumfordShah model (mumford1989optimal) and its piecewise constant approximation models (zhu1996region; cohen1997avoiding; chan2001active), see (morel2012variational) for their proofs. These works also encompass the more difficult problem of multiregion segmentation, and describe their triple point junctions. On the other hand, these results indeed limit their scope to only two specific types of regionbased homogeneity penalty terms , related to either piecewise smooth approximation, or piecewise constant approximation with a gradient squared penalization. The authors are not aware of a theory sufficiently general to accommodate the various nonlinear and nonlocal region terms typically found in applications of (cremers2007review; jung2012nonlocal) based on e.g. their regularity properties. Nevertheless, the existence of the minimizers for these variants is generally regarded as a consequence of the two fundamental cases treated. Finally, let us also mention the existence results in the context of fuzzy (li2010variational), or convex (chambolle2012convex), relaxations of the segmentation problem.
4 Regionbased Eikonal Active Contour Model
In this section, we introduce a new regionbased Ã¥Eikonal active contour model for image segmentation. The main objective of the proposed model is to address the hybrid active contours problem through the Randers Eikonal PDE framework as well as the associated minimal geodesic paths.
Given a shape , we start the description of the proposed Eikonal active contour model by considering a new minimization problem
(26) 
where is a set of admissible curves determined through the given shape . The construction for will be introduced in Section 4.1.
4.1 Construction of the Set
The construction of the set relies on the tubular neighbourhood region of the shape boundary and the tool of the winding number of closed curves.
Tubular Neighbourhood
In essence, the tubular neighbourhood is an open and bounded narrowband domain surrounding the boundary , which can be defined by
(27) 
where is the radius of the tubular neighbourhood and is the Euclidean distance between the points and . The tubular neighbourhood can be regraded as the union of all the disk structures of radius , each of which is centred at the curve . In Fig. 4b, we show an example for the tubular neighbourhood with respect to a shape that is illustrated in Fig. 4a.
In practice, the construction of the tubular neighbourhood can be implemented by means of searching for the level set of a Euclidean distance map which satisfies the following isotropic Eikonal PDE
(28) 
Nevertheless, the construction of the tubular neighbourhood can be efficiently implemented through the isotropic fast marching algorithm (sethian1999fast). The use of a constant function as the potential for the Eikonal PDE (28) generates a symmetric neighbourhood with respect to the shape boundary . Moreover, one can obtain an asymmetric tubular neighbourhood by slightly modifying (28), as introduced in Section 6.1.2.
Winding Number about a Point
Suppose that is a simple and closed curve. Let be a point such that does not to pass through the point . The winding number, or the index of the closed curve associated to the point counts the number of times that rotates around in a given order, that is clockwise order or counterclockwise order. More specifically, the winding number is a nonpositive integer if the curve is parameterized in a clockwise order, and otherwise a nonnegative integer (krantz2012handbook).
For a point which is located outside the simple curve , the winding number is equal to as illustrated in Fig. 5a. In Figs. 5b and 5c, the point is enclosed by the simple and closed curve . This generates that in case is parameterized in a counterclockwise order as shown in Fig. 5b and otherwise, as illustrated in Fig. 5c.
With these definitions in hands, the set associated to a given shape can be expressed as follows
(29) 
where is a set that encompasses all of the simple and closed curves subject to .
By the definition of the set and for each shape such that , we can obtain
(30) 
which yields the following relation
(31) 
Basically, the set defines the searching space for the solution to the minimization problem formulated in Eq. (4).
4.2 Reformulation of the Hybrid Active Contour Energy Functional via a Randers Metric
In this section, we introduce the core of the proposed Eikonal active contours model, that is to transform the functional in Eq. (24) to a weighted curve length associated to an asymmetric Randers metric, providing that the shape is fixed. This transformation requires the construction of a sufficiently small vector field derived from a divergence equation.
For a fixed shape , the functional can be naturally reformulated as
(32) 
where is a scalar value which can be expressed by
(33) 
One can see that the scalar depends only on the characteristic function .
We then focus on the second term of in Eq. (32) which satisfies
Taking the relation (31) into consideration, we obtain
(34) 
where the function is the characteristic function of the tubular neighbourhood . Similar to the use of (see Eq. (33)), the second term of (34) is also independent of the characteristic function . Thus we define a new quantity by
Let be a vector field which solves the divergence equation over the space
(35) 
For a shape with boundary , the combination of Eqs. (32) and (34) yields
(36)  
(37) 
where in Eq. (36) is the uint outward normal to the curve and is the clockwise rotation matrix with angle . Note that equation (36) is obtained from the divergence theorem and equation (37) is derived due to the fact that
(38) 
We denote by a new vector field that is generated by rotating using the matrix in a clockwise order
(39) 
Using the equation (37), one has
(40) 
where is a function with a Randers form as expressed in Eq. (7). It can be formulated for any point and for vector by
(41) 
The function is said to be a Randers metric iff it obeys the positive definiteness condition (8). More specifically, the tensor field and the vector field should be comparable such that
(42) 
Given a shape , we have shown the possibility of transforming the functional to the weighted curve length through the Randers metric . A crucial issue to be addressed for this transformation lies at the existence of a sufficiently small solution to the divergence equation (35), i.e., the existence of such a Randers metric with compatible components and . In next Section, we will present the existence analysis for the metric in case is fixed.
4.3 Existence Analysis
We show in Proposition 1 that there exists a vector field which satisfies the positive definiteness condition (8).
Proposition 1 (Existence)
Let and with be two constants. Let be an open and bounded domain, and let be a vector field which solves the PDEconstrained minimization problem
(43) 
where is a scalarvalued function. Then one has the following inequality
(44) 
where the scalar values and are respectively set as
and is the Lebesgue measure of the set .
Proof
The solution to the PDEconstrained minimization problem stated in Eq. (43) is , where solves the following Poisson equation on :
Expressing the solution in terms of the Green kernel over the domain , we thus obtain for all
(45) 
Therefore, the vector field can be expressed for all as
(46) 
Let be the radius defined by , so that the disk with origin has the same area as the domain . The Hölder’s inequality and a rearrangement inequality yield
(47) 
Evaluating the right hand side of Eq. (47), we have for all such that
which concludes the proof. ∎
4.4 Variational Formulationbased Solution to the PDEconstrained Minimization Problem
In practice, it is difficult to find the solution to the minimization problem (43) over . Instead, we solve that problem over the neighbourhood : finding a vector field which solves
(48) 
where is a scalarvalued function.
Proposition 2
The vector field solving the minimization problem in Eq. (48) admits the variational formulation: finding such that for all in the same spaces, one has
(49)  
(50) 
Proof
For this optimization problem, we first define the augmented functional
(51) 
where is a Lagrange multiplier (function). For the second term in Eq. (51), using integration by parts states that
(52) 
By applying Green’s theorem to the first term of Eq. (4.4) and using the assumption , one has
(53) 
The augmented functional (51) can be reformulated as
(54) 
The firstorder directional derivatives of the augmented functional along and can be respectively expressed as
and
∎
In this paper, the numerical solutions to the variational formulas in Eqs. (49) and (50) are solved by using a finite differences discretization on the regular grid with being the grid size. We also store the values of the functions and on a staggered grid in order to improve the accuracy of the gradient operator.
In Fig. 6b, we show the vector field on a synthetic image, where is estimated by solving the PDEconstrained problems (48). In Fig 6a, the boundary of the given shape is indicated by the red line. In this experiment, the gradient is derived from the piecewise constantstype homogeneity term
4.5 Vectorvalued Convolutionbased Solution to the PDEconstrained Minimization Problem
In this section we make use of convolution method to search for the approximate solution to the minimization problem (43) based on a vectorvalued kernel.
Let us firstly rewrite the closed form (46) by a convolution operator
(55)  
(56) 
where is a vectorvalued kernel defined as follows
(57) 
For numerical implementation, we need a truncated version of the kernel in Eq. (57) by
(58) 
where the scalar parameter is the radius of the truncated domain for the kernel . In Fig. 7, we show an example for a vectorvalued kernel with maximal radius . Without loss of generality, the kernel can be decomposed to two scalarvalued components and such that for any point , the vector can be denoted by . Then the computation of the solution using Eq. (59) can be implemented by
(59) 
which involves two discrete linear convolution operations. Nevertheless, the vector field can be efficiently solved by applying the fast Fourier transform and the inverse fast Fourier transform.
We note that the authors in the literature (li2007active) also proposed an external force computation method by convolving an edge appearance feature map with a vectorvalued kernel, leading to a vector field including edgebased information. While in this section the proposed method is able to estimate a vector field involving regionbased homogeneity information.
4.6 Summary
Let be a given shape with being the sufficiently small tubular neighbourhood of . Let be a vector field which solves the PDEconstrained minimization problem (43). For a given tensor field , the proposition 1 implies that there exists a Randers metric of the form (41) that satisfies the positive definiteness condition (42). As a consequence, the minimization problem (4) is equivalent to the following minimal path computation problem
(60) 
Obviously the minimizer is a closed geodesic path with respect to the Randers metric .
Finding the image segmentation through curve evolution method amounts to iteratively addressing the minimal path computation problem (60) and the subproblem (II). More specifically, the solution to the former problem is a closed curve tracked from the corresponding geodesic distance map. cohen1997global investigated the detection method for simple and closed curves associated to isotropic Riemannian metrics, for which the crucial step is to search for a saddle point from the geodesic distance map. Then one can track two distinct geodesic paths between the source point and the detected saddle point. Unfortunately, this method is difficult to be straightforwardly applied to the asymmetric Randers case. Alternatively, we consider to solve a simplified piecewise geodesic path tracking problem, instead of the original one (60). The basic idea is to reconstruct each evolving curve as the concatenation of a set of disjoint geodesic paths derived using the constructed Randers metric involving regionbased homogeneity information. Hence, the proposed Eikonal active contour model is actually carried out through a disjoint geodesic paths extraction scheme which will be detailed in next section.
5 Disjoint Randers Geodesic Paths Extraction Scheme for Image Segmentation
In this section, we provide two methods for the extraction of disjoint geodesic paths to construct the evolving curve during the curve evolution. Since the computation of a geodesic path requires a source point and a target point, we exploit a set of successive vertices such that each pair of successive points will be taken as a source point and an end point for geodesic path computation.
5.1 Prescribed Pointsbased Scheme: Initialization
The first disjoint geodesic paths extraction scheme depends on a set of userprescribed points for tracking geodesic paths, as in the combination of paths method (mille2015combination). These points are supported to be distributed along the target boundary in a clockwise order and are not allowed to move in the course of curve evolution. Thus the proposed method can be regarded as an interactive image segmentation model.
Let be the set of () points indexed by . In the initialization stage, the goal is to generate a simple and closed curve made up of a series of disjoint open curves, each of which connects a pair of successive prescribed points. We illustrate in Fig. 8a an example of a set of successive points and the corresponding disjoint curves linking them.
Simply connecting each pair of two successive points by a straight segment may potentially yield a curve with selfintersection or selftangency structures in some cases, as depicted in Fig. 9a. This issue can be solved manually by replacing some segments by polyline to remove the unexpected structures, leading to high requirement on user interaction. In order to segment an object in a minimally interactive manner, we exploit the variants of the combination of paths model (mille2015combination) to configure the initialization for the prescribed pointsbased scheme. The combination of paths model integrates the regionbased homogeneity penalty, the edgebased term, and the contour simplicity measure to get a closed curve that minimizes the target energy , as formulated in Eq. (17). The resulting closed curve is expected to be simple, agreeing with the requirement on the initialization of the prescribed pointsbased scheme. By modifying the target energy , we obtain a polygon construction method and an edge simplicity method for the initialization of the prescribed pointsbased scheme, which will be detailed in the following.
Polygon Construction Method. The first method is to construct a polygon as the initial curve passing through all the ordered points . For this purpose, we replace the image datadriven terms and by a Euclidean curve length term , leading to a new functional