Eikonal Region-based Active Contours for Image Segmentation

Eikonal Region-based 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 paths-based 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 region-based 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 region-based 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 region-based active contour model. Experimental results on both of the synthetic and real images exhibit that our model indeed achieves encouraging segmentation performance.

Keywords:
Region-based active contours minimal path Randers metric image segmentation Finsler variant of the fast marching method Eikonal partial differential equation
5

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 edge-based features and a broad variety of region-based 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 problem-dependent 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 edge-based 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 Cauchy-Crofton 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 region-based active contour models in common differ them to the edge-based approaches by an additional term of the homogeneity measure within each region. A pioneering approach is the original Mumford-Shah 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 Mumford-Shah functional has motivated a series of approaches such as (tsai2001curve; vese2002multiphase; bresson2007fast). Basically, the Mumford-Shah 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 Mumford-Shah functional which exploit adequate probability density functions to establish the region-based 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 parameter-dependent probability distribution in each region, where the parametric Gaussian distributions are often utilized to define the region-based 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 Mumford-Shah 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 region-based 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 Sobolev-type inner product in a infinite dimensional Riemannian manifold, leading to a more robust curve evolution scheme.

Figure 1: An example for interactive image segmentation using the combination of paths model (mille2015combination) and the proposed model. a The original image with user-provided vertices indicated by red dots. The image is obtained form the Grabcut dataset (rother2004grabcut). b The result from the combination of paths model. c The result obtained using the proposed model

The active contour models mentioned above are not exhaustive, other interesting approaches include the hybrid active contour models simultaneously integrating the region-based homogeneity penalization and the edge-based features (paragios2002geodesic; kimmel2003regularized; zach2009globally), the shape prior-based active contour models (cremers2002diffusion; bresson2006variational), and the multi-region/multi-label segmentation models (dubrovina2015multi; chambolle2012convex; bae2011global; brox2006level).

1.1 Minimal Path-based 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 path-based 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 curvature-penalized geodesic paths. These paths are tracked in an orientation-lifted space and are derived using the image boundary features extracted by a steerable filter.

Actually, these minimal paths-based image segmentation approaches commonly define the data-driven 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 edge-driven limitation, despite their relevant applications for image segmentation in conjunction with region-based 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 region-based active contours evolution and image segmentation. The corresponding geodesic metric has a Randers form constructed by the region-based 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 edge-driven 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 edge-based 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 region-based 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 region-based active contour model, allows to blend the benefits from both of the region-based 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 region-based 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 PDE-constrained 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 first-order 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 data-driven 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 Hamiltonian-Jacob-Bellman 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 arc-length. One can solve the following ODE

(6)

until the source point is reached, where is the arc-length parameter of the path . Finally, the desired geodesic path can be generated by re-parameterizing and reversing the geodesic path .

Figure 2: Minimal geodesic paths derived through different geodesic metrics. (a) The vector field used to construct the anisotropic Riemannian metric and the asymmetric Randers metric. (b) to (d) Geodesic paths obtained respectively using the isotropic metric, the anisotropic Riemannian metric and the asymmetric Randers metric. The red dots are the source points and the black dots are the end points

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 Cauchy-Schwartz 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 vector-valued function defined as follows

(16)

Again, the target geodesic path can be generated by re-parameterizing 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 circle6 passing through , as depicted in Fig. 2a. In Fig. 2c, we exploit the vector field and the identity matrix to build a Riemannian tensor field

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.

Figure 3: An example for admissible curves which are denoted by lines of different colors. a The potential derived from the image gradients. The potential takes low values around the image boundaries. b-c The admissible curves corresponding to two pairs of successive points

2.3 Overview on the Combination of Paths Model

The paths combination model proposed in (mille2015combination) applies edge-based 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 scalar-valued parameters which control the related importance between the three energy terms. Specifically, the simplicity term measures the amount of the self-tangency and self-intersection of the curve . The edge-based term is measured using an edge indicator along the curve , and the region-based 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 region-based 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 Metric-measured Curve Length for Regularization

In general, a typical region-based active contour energy functional should be comprised of a region-based 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 1-dimensional 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 region-based and boundary-based 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 metric-measured curve length encompasses both the magnitude and the directions of the image gradient features7. As a consequence, the use of the regularization term allows to benefit from the image anisotropy features to enhance the image segmentation (kimmel2003regularized). Note that for images where the gradient information at edge positions is unreliable, the Euclidean curve length of should be taken as the regularization term. In fact, the length denotes the Euclidean curve length of if we set (recall that is the identity matrix).

3.2 Iterative Scheme for the Minimization of the General Region-based Active Contour Model

In region-based 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 scalar-valued 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 region-based functional can be reformulated by its first-order 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 region-based 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:

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

  2. 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 Mumford-Shah 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 multi-region segmentation, and describe their triple point junctions. On the other hand, these results indeed limit their scope to only two specific types of region-based 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 non-linear and non-local 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.

Figure 4: (a) Ah shape and its boundary are demonstrated. (b) The tubular neighbourhood for
Figure 5: Examples for the winding number about a point , noted by black dots, for a simple and closed curve indicated by black lines. (a) The point is outside of the curve and the winding number . (b) The curve is parameterized in a counter-clockwise order and . (c) The curve is parameterized in a clockwise order and

4 Region-based Eikonal Active Contour Model

In this section, we introduce a new region-based å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 non-positive integer if the curve is parameterized in a clockwise order, and otherwise a non-negative 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 counter-clockwise 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 PDE-constrained minimization problem

(43)

where is a scalar-valued 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 PDE-constrained 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. ∎

By Proposition 1, the values of can be bounded by the area of the tubular neighbourhood . This means that given a sufficiently small tubular neighbourhood , one can obtain a vector field satisfying the condition (42).

Figure 6: Visualization for the vector field estimated by the variational formulation-based method and the variational formulation-based method. a The original image and the given shape enclosed by the red line. b The norms are superimposed over the original image. c Visualization for the vector field by a color coding way. d 3D plot for the norms

4.4 Variational Formulation-based Solution to the PDE-constrained 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 scalar-valued 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 first-order 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 PDE-constrained 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 constants-type homogeneity term8 (chan2000active). In Fig. 6d, the norms are superimposed on the original image and the 3D plot of the norms are illustrated in the first and third columns. The vector field is presented in Fig. 6c through the color coding method. As depicted Figs. 6b and 6d, the magnitude of can be regarded as the altitude of a mountainous region over the tubular neighbourhood . By evaluating the closed form (46) associated to a piecewise constant-type gradient , one can point out that the mountain ridges likely lie close to the zero-crossing regions of or close to the boundaries of the tubular neighbourhood. In addition, as the curve evolution tends to stabilize, the evolving contour usually passes through the zero-crossing regions of the gradient.

Figure 7: An example for the truncated vector-valued convolution kernel . The length of the arrow at is positively proportional to the norm of the corresponding vector . The black dot indicates the zero vector . The maximal radius for this truncated kernel is set as for the purpose of visualization

4.5 Vector-valued Convolution-based Solution to the PDE-constrained 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 vector-valued kernel.

Let us firstly rewrite the closed form (46) by a convolution operator

(55)
(56)

where is a vector-valued 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 vector-valued kernel with maximal radius . Without loss of generality, the kernel can be decomposed to two scalar-valued 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 vector-valued kernel, leading to a vector field including edge-based information. While in this section the proposed method is able to estimate a vector field involving region-based 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 PDE-constrained 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 region-based 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

Figure 8: Fixed points-based disjoint geodesic paths extraction scheme. a The red dots represent the prescribed points and the black solid line indicates the input curve . b The tubular neighbourhood . c The decomposed subdomains which are tagged by different colours. d-f The black solid lines indicate the geodesic paths generated within each subdomain , respectively

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 Points-based Scheme: Initialization

The first disjoint geodesic paths extraction scheme depends on a set of user-prescribed 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.

Figure 9: Initialization for the proposed scheme. The prescribed points are indicated by red dots. a Initial curve with self-intersection issue. b Initial curve comprised of four disjoint lines obtained through the polygon construction method
Figure 10: Edge simplicity method. a The image gradients-based potential. b Initial curve obtained using the edge simplicity method
Figure 11: The course of the curve evolution based on the prescribed points-based disjoint geodesic paths extraction scheme. a Initialization. The prescribed points are indicated by red dots. b-c The intermediate steps during the curve evolution. d The final simple and closed curve comprised of a series of disjoint geodesic paths

Simply connecting each pair of two successive points by a straight segment may potentially yield a curve with self-intersection or self-tangency 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 points-based scheme. The combination of paths model integrates the region-based homogeneity penalty, the edge-based 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 points-based scheme. By modifying the target energy , we obtain a polygon construction method and an edge simplicity method for the initialization of the prescribed points-based 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 data-driven terms and by a Euclidean curve length term , leading to a new functional