From Active Contours to Minimal Geodesic Paths: New Solutions to Active Contours Problems by Eikonal Equations

TL;DR

This paper reviews a minimal path framework based on Eikonal PDEs that unifies active contours and elastica models for image analysis tasks like segmentation and boundary detection, offering efficient numerical solutions.

Contribution

It introduces new minimal path models using Riemannian and Randers metrics that solve a wide range of active contour problems and elastica curves within a unified framework.

Findings

Effective boundary detection and segmentation results.

Broad applicability to various image analysis tasks.

Efficient computation using Eikonal solvers like fast marching.

Abstract

In this chapter, we give an overview of part of our previous work based on the minimal path framework and the Eikonal partial differential equation (PDE). We show that by designing adequate Riemannian and Randers geodesic metrics the minimal paths can be utilized to search for solutions to almost all of the active contour problems and to the Euler-Mumford elastica problem, which allows to blend the advantages from minimal geodesic paths and those original approaches, i.e. the active contours and elastica curves. The proposed minimal path-based models can be applied to deal with a broad variety of image analysis tasks such as boundary detection, image segmentation and tubular structure extraction. The numerical implementations for the computation of minimal paths are known to be quite efficient thanks to the Eikonal solvers such as the Finsler variant of the fast marching method.

Equations

$$\nabla I_{\sigma }={\left({\partial }_x{G}_{\sigma }\ast I,{\partial }_y{G}_{\sigma }\ast I\right)}^T,$$

$$g=\Vert \nabla I_{\sigma }\Vert .$$

$${\mathbf{W}}_{\sigma }(\mathbf{x})=\left(\begin{array}{ccc}{\partial }_x{G}_{\sigma }\ast I_1,& {\partial }_x{G}_{\sigma }\ast I_2,& {\partial }_x{G}_{\sigma }\ast I_3\\ {\partial }_y{G}_{\sigma }\ast I_1,& {\partial }_y{G}_{\sigma }\ast I_2,& {\partial }_y{G}_{\sigma }\ast I_3\end{array}\right)(\mathbf{x}).$$

$$g=\sum _{i=1}^3\sqrt{({\partial }_x{G}_{\sigma }\ast I_i{)}^2+({\partial }_y{G}_{\sigma }\ast I_i{)}^2}.$$

$$E_{\mathrm{snake}}(\mathcal{C})={\int }_0^1\left({\alpha }_1\Vert {\mathcal{C}}^{\prime }(u){\Vert }^2+{\alpha }_2\Vert {\mathcal{C}}^{\prime \prime }(u){\Vert }^2+P(\mathcal{C}(u))\right)du,$$

$${\mathcal{C}}^{\prime }(u):=\frac{d}{du}\mathcal{C}(u),{\mathcal{C}}^{\prime \prime }(u):=\frac{d}{du}{\mathcal{C}}^{\prime }(u),$$

$$\frac{\partial {\mathcal{C}}_{\tau }}{\partial \tau }=\underset{\text{Internal Forces}}{\underbrace{{\alpha }_1{\mathcal{C}}_{\tau }^{\prime \prime }-{\alpha }_2{\mathcal{C}}_{\tau }^{\prime \prime \prime \prime }}}-\underset{\text{External Force}}{\underbrace{\nabla P({\mathcal{C}}_{\tau })}},$$

$${\mathbf{F}}_{\mathrm{balloon}}={\mathcal{N}}_{\mathcal{C}},$$

$$\frac{\partial {\mathcal{C}}_{\tau }}{\partial \tau }={\alpha }_1{\mathcal{C}}_{\tau }^{\prime \prime }-{\alpha }_2{\mathcal{C}}_{\tau }^{\prime \prime \prime \prime }+{\alpha }_3{\mathbf{F}}_{\mathrm{balloon}}({\mathcal{C}}_{\tau })-\frac{\nabla P({\mathcal{C}}_{\tau })}{\Vert \nabla P({\mathcal{C}}_{\tau })\Vert },$$

$${\mathcal{L}}_{\mathrm{geo}}(\mathcal{C})={\int }_0^{1}P(\mathcal{C}(u))\Vert {\mathcal{C}}^{\prime }(u)\Vert du.$$

$$\frac{\partial {\mathcal{C}}_{\tau }}{\partial \tau }=P({\mathcal{C}}_{\tau }){\kappa }_{\tau }{\mathcal{N}}_{\tau }-⟨\nabla P({\mathcal{C}}_{\tau }),{\mathcal{N}}_{\tau }⟩{\mathcal{N}}_{\tau },$$

$$\frac{\partial {\mathcal{C}}_{\tau }}{\partial \tau }\propto {\kappa }_{\tau }{\mathcal{N}}_{\tau },$$

$${\mathcal{L}}_{\varrho }(\mathcal{C})={\int }_0^1\varrho (⟨V_{\mathrm{align}}(\mathcal{C}),\mathcal{N}⟩)\Vert {\mathcal{C}}^{\prime }\Vert du,$$

$$\mathscr{L}_{\rm align}(\mathcal{C})=\int_{0}^{1}\big{\langle}\nabla I(\mathcal{C}),\mathcal{N}\big{\rangle}\,\|\mathcal{C}^{\prime}\|du.$$

$$\frac{\partial {\mathcal{C}}_{\tau }}{\partial \tau }=\Delta I_{\sigma }({\mathcal{C}}_{\tau }){\mathcal{N}}_{\tau },$$

$$\varrho (⟨\nabla I(\mathcal{C}),\mathcal{N}⟩)=|⟨\nabla I(\mathcal{C}),\mathcal{N}⟩|,$$

$$\mathscr{L}^{+}_{\rm align}(\mathcal{C})=\int_{0}^{1}\left|\big{\langle}\mathcal{N},\nabla I_{\sigma}(\mathcal{C})\big{\rangle}\right|\,\|\mathcal{C}^{\prime}\|du.$$

$$\frac{\partial {\mathcal{C}}_{\tau }}{\partial \tau }=\mathrm{sign}(⟨\nabla I_{\sigma },{\mathcal{N}}_{\tau }⟩)\Delta I_{\sigma }({\mathcal{C}}_{\tau }){\mathcal{N}}_{\tau }.$$

$$E_{\mathrm{MS}}(f,K)=\alpha {\int }_{\Omega }(I-f{)}^{2}d\mathbf{x}+\hat{\alpha }{\int }_{\Omega \backslash K}\Vert \nabla f{\Vert }^{2}d\mathbf{x}+\mathrm{Per}(K),$$

$$E_{\mathrm{ACWE}}(\mathcal{C},{\mu }_{\mathrm{in}},{\mu }_{\mathrm{out}})=\alpha {\int }_{A_{\mathcal{C}}}(I-{\mu }_{\mathrm{in}}{)}^{2}d\mathbf{x}+\alpha {\int }_{\Omega \backslash A_{\mathcal{C}}}(I-{\mu }_{\mathrm{out}}{)}^{2}d\mathbf{x}+{\mathcal{L}}_{\mathrm{Euclid}}(\mathcal{C}),$$

$$H(a)=\{\begin{array}{ll}1,& a\ge 0,\\ 0,& a<0.\end{array}$$

$${\mathcal{E}}_{\mathrm{ACWE}}(\varphi ,{\mu }_{\mathrm{in}},{\mu }_{\mathrm{out}})=$$

$$+\alpha {\int }_{\Omega }(I-{\mu }_{\mathrm{out}}{)}^2(1-H(\varphi ))d\mathbf{x}$$

$$+{\int }_{\Omega }\delta (\varphi )\Vert \nabla \varphi \Vert d\mathbf{x},$$

$${\mu }_{\mathrm{in}}^{\ast }=\frac{{\int }_{\Omega }IH(\varphi )d\mathbf{x}}{{\int }_{\Omega }H(\varphi )d\mathbf{x}},{\mu }_{\mathrm{out}}^{\ast }=\frac{{\int }_{\Omega }I(1-H(\varphi ))d\mathbf{x}}{{\int }_{\Omega }(1-H(\varphi ))d\mathbf{x}}.$$

$$\underset{\varphi }{\min }{\mathcal{E}}_{\mathrm{ACWE}}(\varphi ,{\mu }_{\mathrm{in}}^{\ast },{\mu }_{\mathrm{out}}^{\ast }),$$

$$\frac{\partial \varphi }{\partial \tau }={\delta }_{ϵ}(\varphi )\left[\alpha (I-{\mu }_{\mathrm{out}}^{\ast }{)}^2-\alpha (I-{\mu }_{\mathrm{in}}^{\ast }{)}^2-\mathrm{div}\left(\frac{\nabla \varphi }{\Vert \nabla \varphi \Vert }\right)\right]$$

$${\mathcal{L}}_{\mathrm{IR}}(\mathcal{C})={\int }_0^1(c_1+P(\mathcal{C}(u))\Vert {\mathcal{C}}^{\prime }(u)\Vert du$$

$${\mathcal{G}}_{\mathbf{s},\mathbf{x}}=\underset{\mathcal{C}\in \mathrm{Lip}([0,1],\Omega )}{\mathrm{argmin}}{\mathcal{L}}_{\mathrm{IR}}(\mathcal{C}),s.t.\{\begin{array}{l}\mathcal{C}(0)=\mathbf{s},\\ \mathcal{C}(1)=\mathbf{x}.\end{array}$$

$${\mathcal{D}}_{\mathbf{s}}(\mathbf{x})$$

$$={\mathcal{L}}_{\mathrm{IR}}({\mathcal{G}}_{\mathbf{s},\mathbf{x}}).$$

$$\{\begin{array}{ll}\Vert \nabla {\mathcal{D}}_{\mathbf{s}}(\mathbf{x})\Vert =c_1+P(\mathbf{x}),& \forall \mathbf{x}\in \Omega \backslash \{\mathbf{s}\},\\ {\mathcal{D}}_{\mathbf{s}}(\mathbf{s})=0,& \end{array}$$

Full text

1 Introduction

The original active contour model proposed by Kass et al., (1988) is a variational approach for the applications of boundary detection, image segmentation and object tracking. It is capable of locally minimizing an energy functional for the purpose of detecting an optimal and continuous curve, either open or closed, to depict the image features of interest such as image edges. The energy functional of this model can be established using the image gray levels or gradients. A broad variety of region-based active contour approaches have been intensively studied in order to develop the original active contour model (Kass et al.,, 1988). They take into account a type of more robust regional homogeneity penalization way for image segmentation. As a consequence, these region-based variant approaches are able to search for suitable solutions in more complicated situations. In the other hand, the minimal geodesic path model was first introduced by Cohen and Kimmel, (1997) in order to search for the global minimum of a simplified version of the functional used in the original active contour model (Kass et al.,, 1988). In essence, this simplified active contour energy functional can be regarded as a weighted curve length and it can be solved through the Eikonal partial differential equation (PDE). The objective of this chapter is to show how the terms that are commonly used in the active contour models can be handled using different kinds of metrics for geodesic paths, in particular Randers metrics.

In its basic formulation, the minimal geodesic path model extracts image features by a globally optimal curve between two prescribed points, usually provided by the user. A minimal path can be tracked through the solution to the Eikonal PDE. A large number of well-established Eikonal solvers have been exploited such as the fast marching methods (Sethian,, 1999; Mirebeau, 2014a, ; Mirebeau, 2014b, ; Mirebeau,, 2018). The global optimality and the efficient solutions lead to a series of successful minimal path-based applications as reviewed in (Peyré et al.,, 2010). However, the original minimal geodesic path model (Cohen and Kimmel,, 1997) depends only on the first-order derivative term of a curve. Thus this model cannot impose curvature to define the regularization term which corresponds to the curve rigidity property. Finally, in the minimal path-based image segmentation scenario, almost all of the existing geodesic metrics (Mille et al.,, 2015; Benmansour and Cohen,, 2009; Appia and Yezzi,, 2011) are derived from the image boundary features, which may limit their performance on complicated image segmentation tasks. In the remaining of this chapter, we proposed different geodesic metrics to overcome the aforementioned drawbacks suffered by the existing minimal path models.

In (Kimmel,, 2003; Kimmel and Bruckstein,, 2003), the authors proposed a listing of the active contour energy terms for image segmentation. In our chapter, we propose to revisit Kimmel’s chapter (Kimmel,, 2003), named fast edge integration, with minimal path methods to address the active contour problems.

1.1 Outline

In this chapter, we illustrate different types of geodesic metrics and their respective applications in image analysis based on the minimal path framework and the Eikonal PDE. The structure of this chapter is organized as follows:

• •

In Section 2, we present the preliminary background on the active contour models involving the approaches respectively driven by edge- and region-based terms.

• •

In Section 3, we review the existing minimal path models associated with both of the Riemannian and Randers metrics, the corresponding Eikonal PDEs for the computation of the geodesic distance as well as the gradient descent ordinary differential equations (ODEs) for tracking geodesic curves.

• •

In Section 4, we present two minimal path models for the minimization of the geometric active contour models with alignment terms (Kimmel,, 2003; Kimmel and Bruckstein,, 2003). We respectively induce a Randers metric and an anisotropic Riemannian metric from the variants of the asymmetric and symmetric alignment terms. In this case, the active contour problems with different alignment terms are transferred to the corresponding minimal path problems.

• •

In Section 5, we present a Finsler elastica minimal path model (Chen et al.,, 2017) to solve the Mumford-Euler elastica curve problem (Mumford,, 1994). This is done by using a strongly anisotropic and asymmetric Randers geodesic metric that is established over an orientation-lifted space. The orientation dimension is used to represent the curve directions. The introduced Finsler elastica minimal paths can be naturally applied for the applications of retinal imaging and boundary detection.

• •

In Section 6, we present a new method that applies the Randers minimal path model to address the region-based active contour problems (Chen et al.,, 2016, 2019). The used Randers metric is derived from the region-based homogeneity term as well as the curve length-based regularization term, relying on the solution to divergence equation-constrained minimization problem. In this case, we can make use of the Randers Eikonal PDE to solve a broad variety of region-based active contour problems.

2 Active Contour Models

The active contour models originated from the work proposed in (Kass et al.,, 1988) have obtained successful applications in the fields of boundary detection, image segmentation, object tracking and shape modelling. Basically, the goal of an active contour model can be summarized as seeking a regular parametric curve for the delineation of the image feature of interest, such as the object boundaries, by minimizing an energy functional. In general, the minimization of an active contour energy functional can be carried out by the curve evolution scheme. The crucial point for this scheme is the estimation of the gradient descent flow with adequate boundary conditions. Basically, this scheme amounts to searching for a parametric curve $\mathcal{C}_{\tau}:(u,\tau)\in[0,1]\times\mathbb{R}^{+}_{0}\mapsto\mathcal{C}_{\tau}(u)\in\Omega$, where $u\in[0,1]$ is the curve parameter and $\tau\in\mathbb{R}^{+}_{0}$ represents the (artificial) time and $\Omega\subset\mathbb{R}^{2}$ denotes the image domain. The curve $\mathcal{C}_{\tau}$ is supposed to converge to the target boundary when $\tau\to\infty$.

2.1 Edge-based Active Contour Model

2.1.1 The Computation of the Image Gradient Features

Let $I:\Omega\to\mathbb{R}$ be a gray level image. The standard Euclidean gradient of a (smoothed) image $I$ can be defined in conjunction with a Gaussian kernel $G_{\sigma}$ with variance $\sigma$:

[TABLE]

where $\ast$ denotes the convolution operator, and $\partial_{x}G_{\sigma}$, $\partial_{y}G_{\sigma}$ are the first-order partially differentials of the Gaussian kernel $G_{\sigma}$ respectively along the axes $x$ and $y$. The magnitude $g:\Omega\to\mathbb{R}^{+}_{0}$ of the image gradient $\nabla I_{\sigma}$ is a scalar-valued function that can be formulated by

[TABLE]

In the context of boundary detection and image segmentation, the magnitude $g$ often serves as the edge appearance feature map.

For a vector-valued image $\mathbf{I}=(I_{1},I_{2},I_{3})$ in a RGB color space, we first consider a Jacobi matrix $\mathbf{W}_{\sigma}(\mathbf{x})$ of size $2\times 3$ to construct the associated image boundary features for a (smoothed) color image $\mathbf{I}$, i.e.

[TABLE]

The edge appearance feature map $g$ for a color image can be estimated using the Frobenius norm of the matrix $\mathbf{W}_{\sigma}(\cdot)$. In other words, one can build the magnitude $g$ as follows:

[TABLE]

2.1.2 The original active contour model

The original edge-based active contour model or the snakes model (Kass et al.,, 1988) aims to minimize the following energy functional

[TABLE]

where $\mathcal{C}\in H^{2}([0,1],\Omega)$ is a regular curve with non-vanishing velocity and the parameters $\alpha_{1}$ and $\alpha_{2}$ are real positive constants controlling the importance for each component. The functions $\mathcal{C}^{\prime}$ and $\mathcal{C}^{\prime\prime}$ of the forms

[TABLE]

are respectively the first- and second-order derivatives of the curve $\mathcal{C}$. Note that the squared norm $\|\mathcal{C}^{\prime\prime}(\cdot)\|^{2}$ is relevant to the squared curvature of $\mathcal{C}$. The function $P:\Omega\to\mathbb{R}^{+}_{0}$, which is referred to as potential, is an image data-driven function. It is defined to have low values around the image structures of interest and high values otherwise. For instances, the potential $P$ can be set as the decreasing function of the image appearance feature map $g$ or can be derived from the image gray levels (Kass et al.,, 1988).

The curve evolution equation can be derived from the Euler-Lagrange equation of the energy functional $E_{\rm snake}$, which reads

[TABLE]

where $\mathcal{C}_{\tau}^{\prime\prime\prime\prime}(u)$ is the fourth-order derivative of the curve $\mathcal{C}_{\tau}$ with respect to $u$. The first two terms in the evolution equation (7) obtained from the curve itself is referred to as the internal forces which impose the regularity of the curve, while the last term is the external force used to attract the curve to the features (e.g. object edges) of the target.

2.1.3 Balloon Force

Cohen, (1991) introduced an external force, named the balloon force, to relax the demanding initialization requirement of the original active contours model (Kass et al.,, 1988). According to (Cohen,, 1991), the balloon force $\mathbf{F}_{\rm balloon}$ can be expressed by

[TABLE]

where $\mathscr{N}_{\mathcal{C}}$ is the outward unit normal to the curve $\mathcal{C}$. The balloon force $\mathbf{F}_{\rm balloon}$ is able to drive the curve expanding from the interior region of $\mathcal{C}$ to the desired features such as the object boundaries. Thus the initial guess should be located inside the target and allowed to be far from the boundaries.

In conjunction with the balloon force $\mathbf{F}_{\rm balloon}$, the corresponding curve evolution equation eventually turns out to be

[TABLE]

where $\alpha_{3}\in\mathbb{R}^{+}$ is a weighting parameter that controls the importance of the balloon force during the curve evolution.

2.1.4 Geodesic active contour model

Let $\operatorname{Lip}([0,1],\Omega)$ be the set of all the Lipschitz continuous curves $\mathcal{C}:[0,1]\to\Omega$. The geodesic active contour model proposed in (Caselles et al.,, 1997; Yezzi et al.,, 1997) takes into account a potential function $P$ to build a weighted curve length

[TABLE]

The weighted curve length $\mathscr{L}_{\rm geo}$ is intrinsic since it is independent of the parameterization of the curve $\mathcal{C}$. Comparing to the energy functional (5) considered in the original snakes model (Kass et al.,, 1988), the term which relies on the second-order derivative of the curve has been removed.

An extremum curve that minimizes the weighted curve length $\mathscr{L}_{\rm geo}$ is a geodesic path. In (Caselles et al.,, 1997; Yezzi et al.,, 1997), a locally minimizing curve can be generated by the curve evolution scheme. Denoting by $\kappa_{\tau}:[0,1]\to\mathbb{R}$ the curvature of $\mathcal{C}_{\tau}$, the curve evolution equation obtained through the Euler-Lagrange equation of the weighted curve length $\mathscr{L}_{\rm geo}$ reads

[TABLE]

where the operator $\langle\cdot,\cdot\rangle$ stands for the Euclidean scalar product, $\mathcal{N}_{\tau}$ is the inward unit normal to the curve $\mathcal{C}_{\tau}$, and $\nabla P$ is the Euclidean gradient of the function $P$ with respect to positions. In the homogeneous regions of an image, the gradients of the potential $P$ tend to vanish, leading to the following Euclidean curve shortening flow

[TABLE]

where $\propto$ is the positively proportional operator. When the curve $\mathcal{C}_{\tau}$ is near an edge at some time $\tau$, the values of $P(\mathcal{C}_{\tau}(\cdot))$ are low and the curve evolution will be dominated by the second term of equation (11). Therefore, it is a good choice to initialize the curve $\mathcal{C}_{\tau}|_{\tau=0}$ surrounding the target. The drawback of the solution to the original active contour energy (5) is the high possibility of being stuck in a local minimum. The level set method introduced by Osher and Sethian, (1988) is a very popular way for solving the curve evolution problem (Malladi et al.,, 1995; Caselles et al.,, 1997; Yezzi et al.,, 1997), but it still can be stopped by unexpected local minima. This is why the authors in (Cohen and Kimmel,, 1997) introduced the use of the minimal geodesic path to find the global minimum of a simplified active contour energy functional, i.e. the weighted curve length (10). The global minimum is attained by solving a nonlinear partial differential equation, called the Eikonal PDE, which will be briefly introduced in Section 3.1.

2.1.5 Geometric active contour models with alignment terms

The aforementioned geodesic active contour model (Caselles et al.,, 1997; Yezzi et al.,, 1997) as well as the fast implementation approach (Goldenberg et al.,, 2001) utilize a direction-independent potential function $P$ to measure the length of a curve. Kimmel and Bruckstein, (2003) introduced an anisotropic variant of the geodesic active contour model based on an image data-driven vector field $V_{\rm align}:\Omega\to\mathbb{R}^{2}$, leading to a new weighted curve length formula

[TABLE]

where $\varrho$ is a scalar-value function and $\langle V_{\rm align}(\mathcal{C}),\mathcal{N}\rangle$ is an alignment measure term. For the applications of boundary detection and image segmentation, the vector field $V_{\rm align}$ can be naturally built by the image gradient vector field $\nabla I_{\sigma}$.

Kimmel and Bruckstein, (2003) suggested two effective types of $\varrho$ for the alignment measure term. The first one is to choose $\varrho(\langle\nabla I_{\sigma}(\mathcal{C}),\mathcal{N}\rangle)=\langle\nabla I_{\sigma}(\mathcal{C}),\mathcal{N}\rangle$ which yields an asymmetric functional

[TABLE]

The curve evolution equation for the evolving curve $\mathcal{C}_{\tau}$ with respect to the functional (13) can be expressed by

[TABLE]

where $\Delta$ is the Laplacian operator such that $\Delta I_{\sigma}=\operatorname{div}(\nabla I_{\sigma})$. The extremum curves $\mathcal{C}_{\tau}$ as $\tau\to\infty$ corresponds to a maximizer of the functional $\mathscr{L}_{\rm align}$.

The second choice introduced by Kimmel and Bruckstein, (2003) for the function $\varrho$ is to set

[TABLE]

which yields a robust and symmetric variant of the asymmetric case (13) as follows:

[TABLE]

The curve evolution equation in this case has a form

[TABLE]

The maximization of the alignment-based functionals $\mathscr{L}_{\rm align}$ and $\mathscr{L}^{+}_{\rm align}$ can be respectively solved through the curve evolution equations (14) and (16) in conjunction with the level set framework (Osher and Sethian,, 1988). An optimal curve generated by means of the evolution equation (14) (resp. the evolution equation (16)) corresponds to the local maximum of the alignment-based functional $\mathscr{L}_{\rm align}$ (resp. the functional $\mathscr{L}^{+}_{\rm align}$). In Sections. 4.1 and 4.2, we will introduce a Randers metric and an anisotropic Riemannian metric respectively for the global minimization of the slight variants of the alignment-based functionals $\mathscr{L}_{\rm align}$ and $\mathscr{L}^{+}_{\rm align}$.

2.2 The Piecewise Smooth Mumford-Shah Model and the Piecewise Constant Reduction Model

The edge-based active contour approaches mentioned above, which take the image gradients as the edge appearance features for image boundary extraction, are quite efficient and effective models. However, the image gradients are very sensitive to the image noise and spurious edges, which also likely yield strong edge appearance features. In order to address this issue, the region-based active contour approaches exploit more complicated homogeneity criteria to establish the energy functionals. Significant examples of these region-based approaches involve the piecewise smooth Mumford-Shah model (Mumford and Shah,, 1989) as well as the associated variant approaches such as (Zhu and Yuille,, 1996; Cohen,, 1997; Chan and Vese, 2001a, ; Brox and Cremers,, 2009).

2.2.1 The Piecewise Smooth Mumford-Shah Model

Let $K$ be a closed subset of the image domain $\Omega\subset\mathbb{R}^{2}$, which is made up of a finite set of smooth curves. The goal for the original piecewise smooth Mumford-Shah model is to seek the set $K$ and a piecewise smooth function $f:\Omega\to\mathbb{R}$ by minimizing the energy functional

[TABLE]

where $\alpha$ and $\hat{\alpha}$ are two positive constants which control the relative importance between different terms, $\operatorname{Per}(K)$ is the perimeter of the set $K$, and $f$ is referred to as an image data fitting function. The first integration term defines a similarity measure by estimating the approximation errors between the image gray levels $I$ and the data fitting function $f$. The second integration is a smoothness penalty term that imposes a prior to the function $f$.

Since the original work in (Mumford and Shah,, 1989), many efforts have been devoted to seeking practical solutions to the minimization of the piecewise smooth Mumford-Shah problem (17). Among them, the explicit curve evolution scheme is a widely considered tool for minimizing the energy functional (17), which can be implemented either by the level set method (Osher and Sethian,, 1988; Chan and Vese, 2001b, ; Vese and Chan,, 2002; Tsai et al.,, 2001) or by the convex relaxation framework (Chan et al.,, 2006; Bresson et al.,, 2007; Pock et al.,, 2009).

2.2.2 Piecewise Constant Reduction of the Full Mumford-Shah Model

The piecewise constant active contour model (Chan and Vese, 2001a, ; Chan et al.,, 2000; Cohen,, 1997) is a practical variant of the full Mumford-Shah model. Instead of using a piecewise smooth function $f$ to approximate the image data, the piecewise constant model or the active contours without edges model (ACWE) assumes that the image gray levels within each region can be approximated by the mean intensity value estimated in the corresponding region.

Let $\mathcal{C}$ be a simple and closed curve. We denote by $A_{\mathcal{C}}\subset\Omega$ the open and bounded subset enclosed by the curve $\mathcal{C}$. In the case of binary-valued image segmentation scenario, the piecewise constant active contour model (Chan and Vese, 2001a, ; Cohen,, 1997) considers an energy functional

[TABLE]

where $\mu_{\rm in}$ and $\mu_{\rm out}$ are respectively the mean gray levels inside and outside the closed curve $\mathcal{C}$ and the term $\mathscr{L}_{\rm Euclid}(\mathcal{C})$ stands for the Euclidean curve length of $\mathcal{C}$ which serves as a regularization term for the active contour model. Note that in the original Chan-Vese piecewise constant model (Chan and Vese, 2001a, ), there is a term related to the area of the region $A_{\mathcal{C}}$, which has been ignored here.

2.2.3 Minimization of the Piecewise Constant Active Contour Model

Chan and Vese, 2001a proposed a solution for the minimization of the active contour energy functional $E_{\rm ACWE}$ based on the level set method (Osher and Sethian,, 1988; Zhao et al.,, 1996). Basically, the level set formulation relies on a Lipschitz function $\phi:\Omega\to\mathbb{R}$ such that the open and bounded subset $A_{\mathcal{C}}$ can be represented by $A_{\mathcal{C}}:=\{\mathbf{x}\in\Omega;\phi(\mathbf{x})>0\}$ and its boundary $\mathcal{C}$ can be identified from the zero level set of $\phi$, i.e. $\mathcal{C}=\{\mathbf{x}\in\Omega;\phi(\mathbf{x})=0\}$. In numerical implementation, the level set function $\phi$ is often set to be a signed Euclidean distance function.

Let $H$ be the Heaviside function which can be formulated by

[TABLE]

The characteristic function of the subset $A_{\mathcal{C}}$ can be defined by the Heaviside function of the corresponding level set function. Nevertheless, one can reformulate the functional (18) using the Heaviside function $H$ and the level set function $\phi$ as follows:

[TABLE]

where $\delta=H^{\prime}$ is the one-dimensional Dirac measure (Chan and Vese, 2001a, ).

As introduced in (Chan and Vese, 2001a, ), the minimization of the functional $\mathcal{E}_{\rm ACWE}$ formulated in Eq. (2.2.3) can be addressed by iteratively evaluating the following two steps:

(1). Keeping the level set function $\phi$ fixed and minimizing the functional $E_{\rm ACWE}$ with respect to the scalar values $\mu_{\rm in}$ and $\mu_{\rm out}$ will yield

[TABLE]

(2). Based on the scalar values $\mu_{\rm in}^{*},\,\mu_{\rm out}^{*}$ derived from Eq. (20), one solves

[TABLE]

with respect to the level set function $\phi$.

According to (Chan and Vese, 2001a, ), finding the numerical solution to the minimization problem (21) requires a regularized function $H_{\epsilon}$ which approximates the Heaviside function $H$ as $\epsilon\to 0$. By computing the Euler-Lagrange equation for the level set function $\phi$, one can obtain the following level set evolution equation

[TABLE]

where $\delta_{\epsilon}:=H_{\epsilon}^{\prime}$ is the regularized version of $\delta$.

3 Minimal Paths for Edge-based Active Contours Problems

3.1 Cohen-Kimmel Minimal Path Model

The minimal path model introduced by Cohen and Kimmel, (1997) is a powerful tool in the fields of image analysis and medical imaging (Cohen,, 2006; Peyré et al.,, 2010). It is designed to search for the global minimum of the weighted curve length energy (Caselles et al.,, 1997; Yezzi et al.,, 1997), measured along a curve $\mathcal{C}\in\operatorname{Lip}([0,1],\Omega)$ as follows:

[TABLE]

where $c_{1}$ is a positive constant used for regularization and the function $P:\Omega\to\mathbb{R}^{+}$ is a scalar-valued potential. In the context of boundary detection, the potential $P$ is usually set as a decreasing function of the image gradient magnitude (Cohen and Kimmel,, 1997).

Given a source point $\mathbf{s}$ and a target point $\mathbf{x}$ in the domain $\Omega$, a geodesic path or a minimal path $\mathcal{G}_{\mathbf{s},\mathbf{x}}\in\operatorname{Lip}([0,1],\Omega)$ with $\mathcal{G}_{\mathbf{s},\mathbf{x}}(0)=\mathbf{s}$ and $\mathcal{G}_{\mathbf{s},\mathbf{x}}(1)=\mathbf{x}$, is a curve that globally minimizes the length $\mathscr{L}_{\rm IR}$ among all the paths in the set $\operatorname{Lip}([0,1],\Omega)$. In other words, a geodesic path $\mathcal{G}_{\mathbf{s},\mathbf{x}}$ can be determined by

[TABLE]

Tracking such a geodesic path from the source point $\mathbf{s}$ to a target point $\mathbf{x}$ relies on a geodesic distance map $\mathscr{D}_{\mathbf{s}}:\Omega\to\mathbb{R}^{+}_{0}$. At each point $\mathbf{x}$, the value of the geodesic distance map $\mathscr{D}_{\mathbf{s}}(\mathbf{x})$ is equal to the minimal curve length from $\mathbf{s}$ to $\mathbf{x}$

[TABLE]

It is known that the geodesic distance map $\mathscr{D}_{\mathbf{s}}$ is the unique viscosity solution to the Hamiltonian-Jacobi-Bellman equation or the Eikonal equation

[TABLE]

where $\nabla\mathscr{D}_{\mathbf{s}}$ is the standard Euclidean gradient of the geodesic distance map $\mathscr{D}_{\mathbf{s}}$.

Solving a gradient descent ODE on the distance map $\mathscr{D}_{\mathbf{s}}$ can generate a geodesic path, or a minimal path $\hat{\mathcal{G}}_{\mathbf{x},\mathbf{s}}$, which connects from $\mathbf{x}$ to $\mathbf{s}$. It corresponds to a back-propagation processing from the target point $\mathbf{x}$ to the source point $\mathbf{s}$. The gradient descent ODE can be expressed by

[TABLE]

with boundary condition $\hat{\mathcal{G}}_{\mathbf{x},\mathbf{s}}(0)=\mathbf{x}$. Note that the geodesic curve $\hat{\mathcal{G}}_{\mathbf{x},\mathbf{s}}$ is parameterized by its arc-length. The gradient descent ODE (27) can be numerically solved by Heun’s or Runge-Kutta’s methods.

In Fig. 1, we show an example for the computation of a geodesic path between a given pair of prescribed points on an image from the BSDS500 dataset (Arbelaez et al.,, 2011). In Fig. 1a, the red and yellow dots indicate the source point and the end point, respectively. The corresponding geodesic distance map and the target geodesic path are shown in Figs. 1b and 1c, respectively. The potential $P$ is computed by in terms of image gradient features

[TABLE]

where $\beta$ is a positive constant and $g$ is the magnitude of the color image gradients, see Eq. (4). We apply the isotropic fast marching method (Sethian,, 1999) for the estimation of the geodesic distance map.

In Fig. 2 we show an example of applying a set of minimal geodesic paths to extract a curvilinear tree structure. In this example, in order to find the whole tree structure, we give three end points which are located at the end of each bifurcation. These end points are indicated by yellow dots as illustrated in Fig. 2a. The geodesic distance map and the associated minimal geodesic paths (represented by blue lines) are shown in Figs. 2b and 2c, repsectively. In this experiment, we build the potential $P$ using the image gray levels as follows:

[TABLE]

where the scalar value $c_{2}\in\mathbb{R}^{+}$ is a regularization factor.

3.2 Finsler and Randers Minimal Paths

A Finsler metric is a continuous map $\mathfrak{F}_{\rm F}:\Omega\times\mathbb{R}^{2}\to[0,+\infty]$ over the space $\Omega\times\mathbb{R}^{2}$. For each point $\mathbf{x}\in\Omega$, the Finsler metric $\mathfrak{F}_{\rm F}(\mathbf{x},\vec{\mathbf{u}})$ can be defined through an asymmetric norm, which is a convex and 1-homogeneous function on its second argument. Basically, the asymmetry property of a Finsler metric allows to take advantages of path directions during the computation of the geodesic distance and of the geodesic paths, as explored in (Melonakos et al.,, 2008; Chen and Cohen,, 2018).

The weighted length of a curve $\mathcal{C}\in\operatorname{Lip}([0,1],\Omega)$ can be measured with respect to a Finsler metric $\mathfrak{F}_{\rm F}$ (see Eq. (32)) by

[TABLE]

Similar to the isotropic case defined in Eq. (29), the geodesic distance map $\mathscr{D}_{\mathbf{s}}:\Omega\to\mathbb{R}^{+}_{0}$ associated to the weighted curve length $\mathscr{L}_{\rm F}$ can be formulated by

[TABLE]

In order to estimate the geodesic distance map $\mathscr{D}_{\mathbf{s}}$, one can solve the following Finsler Eikonal equation:

[TABLE]

The corresponding gradient descent ODE on the geodesic distance map $\mathscr{D}_{\mathbf{s}}$ with respect to the Finsler metric $\mathfrak{F}_{\rm F}$ can be written as

[TABLE]

The target geodesic path $\mathcal{G}_{\mathbf{s},\mathbf{x}}$ can be computed by re-parameterizing the geodesic path $\hat{\mathcal{G}}_{\mathbf{x},\mathbf{s}}$.

In the following, we focus on the Randers metric (denoted by $\mathfrak{F}_{\rm Randers}$) and the Riemannian metric (denoted by $\mathfrak{R}$), which are two particular cases of Finsler metric. The Eikonal PDEs with respect to these two types of geodesic metrics $\mathfrak{F}_{\rm Randers}$ and $\mathfrak{R}$ are presented in Sections 3.2.1 and 3.2.2, respectively.

3.2.1 Randers Minimal Paths

Let $\mathbb{S}^{+}_{2}$ be the set of all the positive definite symmetric matrices of size $2\times 2$. A Randers metric (Randers,, 1941) is comprised of a positive definite symmetric tensor field $\mathcal{M}:\Omega\to\mathbb{S}^{+}_{2}$ and a vector field $\omega:\Omega\to\mathbb{R}^{2}$

[TABLE]

where $\langle\cdot,\cdot\rangle$ stands for the standard Euclidean scalar product on $\mathbb{R}^{2}$. A Randers metric $\mathfrak{F}_{\rm Randers}$ of the form (32) should obey the positive definiteness condition over the domain $\Omega\times\mathbb{R}^{2}$. In other words, one has $\mathfrak{F}_{\rm Randers}(\mathbf{x},\vec{\mathbf{u}})>0$ for any point $\mathbf{x}\in\Omega$ and any vector $\vec{\mathbf{u}}$, if and only if

[TABLE]

For the applications of image analysis, the tensor field $\mathcal{M}$ and the vector field $\omega$ can be constructed dependently on the tasks, as formulated in Sections. 5 and 6.

The Eikonal equation (30) with respect to a Randers metric $\mathfrak{F}_{\rm Randers}$ is equivalent to the following nonlinear PDE (Mirebeau,, 2017; Chen and Cohen,, 2018)

[TABLE]

with boundary condition $\mathscr{D}_{\mathbf{s}}(\mathbf{s})=0$.

3.2.2 Riemannian Minimal Paths

The Randers metric $\mathfrak{F}_{\rm Randers}$ of the form (32) gets to be a Riemannian metric $\mathfrak{R}$ of a general form, if the vector field $\omega$ obeys that $\omega\equiv\mathbf{0}$, i.e.,

[TABLE]

Then the Eikonal equation associated to the Riemannian metric $\mathfrak{R}$ defined in Eq. (35) can be formulated by

[TABLE]

In this case, the gradient descent ODE defined in (31) gets to

[TABLE]

This ODE can be solved by Mirebeau’s ODE solver (Mirebeau, 2014a, ).

4 Minimal Paths for Alignment Active Contours

Inspired by the geometric active contour models with edge-based alignment terms (Kimmel and Bruckstein,, 2003; Kimmel,, 2003), we introduce a Randers metric and an anisotropic Riemannian metric based on the image gradients for the application of boundary detection. These metrics are induced using the variant forms of the asymmetric and symmetric edge-based alignment terms.

4.1 Randers Alignment Minimal Paths

We introduce a method, named the Randers alignment minimal path model, for boundary detection. The basic idea is to integrate the asymmetric alignment term $\mathscr{L}_{\rm align}$ formulated in Eq. (13), and a Euclidean curve length-based regularization term to establish a new edge-based energy functional

[TABLE]

where $\beta\in\mathbb{R}^{+}$ is a parameter that controls the relative importance between the regularization term $\mathscr{L}_{\rm Euclid}$ and the edge-based asymmetric align term $\mathscr{L}_{\rm align}$. Note that $\nabla I_{\sigma}$ represents the gradient vector field of the smoothed image $I$, see Eq. (1).

The objective in this section is to globally minimize $\mathcal{L}_{\rm align}$ by searching for a geodesic path between two given points. Let $\mathbf{M}$ be a rotation matrix with angle $\pi/2$ and let $\xi$ be the rotated image gradient vector field

[TABLE]

In this case, a vector $\xi(\mathbf{x})$ indicates the direction of a boundary should have at an edge point $\mathbf{x}$.

We reformulate the energy functional $\mathcal{L}_{\rm align}$ by

[TABLE]

where $\mathfrak{F}_{\rm align}$ can be expressed by

[TABLE]

Now we can see that $\mathcal{L}_{\rm align}$ has been formulated as a weighted curve length associated to a Randers metric. The function $\mathfrak{F}_{\rm align}$ is required to obey the positive definiteness condition (33) in order to globally minimize $\mathcal{L}_{\rm align}$ through the Eikonal PDE framework. One simple solution is to limit the range of the parameter $\beta$ as follows

[TABLE]

However, a small value of $\beta$ may reduce the importance of the image data, thus increasing the risk of shortcuts.

As a second choice, we make use of a new vector field $\tilde{\xi}:\Omega\to\mathbb{R}^{2}$ to replace the term $\beta\xi$ used in the Randers metric $\mathfrak{F}_{\rm align}$. The considered vector field $\tilde{\xi}$ can be expressed by

[TABLE]

where $\varphi:\mathbb{R}^{+}\to\mathbb{R}^{+}$ is a scalar function

[TABLE]

By the vector field $\tilde{\xi}$, we obtain a new Randers metric $\mathcal{F}_{\rm align}$

[TABLE]

One can see that the positive definiteness condition characterized by (33) will always hold for the Randers metric $\mathcal{F}_{\rm align}$.

Notice that through the vector field $\tilde{\xi}$, we actually minimize a new weighted curve length

[TABLE]

where the second term in the second line of Eq. (44) can be treated as a variant of the asymmetric alignment term $\beta\mathscr{L}_{\rm align}$.

The asymmetric edge feature at each point $\mathbf{x}$ still can be characterized by the vector $\tilde{\xi}(\mathbf{x})$. Specifically, a large value of $\|\tilde{\xi}(\mathbf{x})\|$ indicates a high possibility that the point $\mathbf{x}$ belongs to an image edge. Moreover, when $\mathbf{x}$ is located at an edge, the direction $\tilde{\xi}(\mathbf{x})/\|\tilde{\xi}(\mathbf{x})\|=\xi(\mathbf{x})/\|\xi(\mathbf{x})\|$ is proportional to the tangent of the image edge at point $\mathbf{x}$.

In Fig. 3b, we visualize the vector field $\tilde{\xi}$ through the color coding scheme (Baker et al.,, 2011). In this figure, we only illustrate the vectors $\tilde{\xi}(\cdot)$ within a rectangle as shown in Fig. 3a. In Fig. 3c, we show the color coding, where the direction of a vector is coded by hue while the norm is coded by saturation (Baker et al.,, 2011). The color image used in this experiment is obtained from the BSDS500 dataset (Arbelaez et al.,, 2011).

Notice that when dealing with a vector-valued image $\mathbf{I}:=(I_{1},I_{2},I_{3})$, we build the vector field $\xi$ through the following equation:

[TABLE]

The asymmetry property of the Randers metrics $\mathfrak{F}_{\rm align}$ and $\mathcal{F}_{\rm align}$ allows to exploit the asymmetric edge features for boundary detection and image segmentation. In Fig. 4, we show the boundary detection results using the Randers metric $\mathcal{F}_{\rm align}$ on a color image. In Fig. 4a, the red dot and the green dot respectively denote the given source and end points. In this experiment, the objective is to extract a geodesic path to depict the boundary of the object. The geodesic distance map superimposed in the original image is shown in Fig. 4c. Note that we adopt the Finsler variant of the fast marching method introduced by Mirebeau, 2014b as the Eikonal solver. The computation of the geodesic distance map is terminated once the end point, denoted by the yellow dot, is reached by the geodesic distance front (Deschamps and Cohen,, 2001).

In Fig. 5, we illustrate the asymmetry property of the Randers alignment minimal paths. In Fig. 5a, the geodesic path is obtained by taking the red dot as the source point and the blue dot as the end point. While in Fig. 5b, we exchange the source and end points used in Fig. 5a. One can obtain a different minimal path to the one shown in Fig. 5a. Then both minimal paths can be used to form a closed curve which depicts the complete object boundary, as shown in Fig. 5c.

In Fig. 6, we show the Randers alignment minimal path extraction results with respect to different values of the parameter $\beta$. Specifically, we increase the values of $\beta$ from Figs. 6a to 6d. In particular, we use $\beta=0$ in Fig. 6a which generates a straight line segment. One can see that a small value of $\beta$ which yields a weakly asymmetric Randers metric may generate a geodesic path with unexpected shortcuts, as shown in Figs. 6b and 6c. For a suitable value of $\beta$, the obtained minimal path can successfully depict the object boundary, as shown in Fig. 6d.

4.2 Riemannian Alignment Minimal Paths

In this section, we induce an anisotropic Riemannian metric from a variant form of the symmetric alignment term $\mathscr{L}_{\rm align}^{+}$ which is defined in Eq. (15). By integrating a regularization term and the symmetric alignment term, we first introduce an energy functional with a sufficiently small $\beta$

[TABLE]

The tensor field $\mathcal{M}_{\rm align}:\Omega\to\mathbb{S}^{+}_{2}$ can be formulated by

[TABLE]

providing that $\beta$ is sufficiently small, where $\xi=\mathbf{M}\,\nabla I_{\sigma}$ is the rotated image gradient vector field with angle $\pi/2$. In order to ensure the matrix $\mathcal{M}_{\rm align}(\mathbf{x})$ to be positive symmetric definite for any point $\mathbf{x}$, we replace the term $\beta^{2}\xi(\cdot)\xi(\cdot)^{T}$ in Eq. (46) by $\tilde{\xi}(\cdot)\tilde{\xi}(\cdot)^{T}$ defined in Eq. (41).

Hence we can obtain a new tensor field

[TABLE]

The matrix $\tilde{\mathcal{M}}_{\rm align}(\mathbf{x})$ will always be positive definite symmetric for each point $\mathbf{x}$ since the inequality $\|\tilde{\xi}(\mathbf{x})\|<1$ always holds. Finally, we can define a new Riemannian metric

[TABLE]

Following (Sapiro,, 1997; Di Zenzo,, 1986), we take into account the matrix $\mathbf{W}_{\sigma}(\mathbf{x})$ (defined in Eq. (3)) to estimate the vector $\xi(\mathbf{x})$ from a color image $\mathbf{I}$. More precisely, we set

[TABLE]

where the vector $\upsilon(\mathbf{x})$ is the eigenvector corresponding to the smaller eigenvalue of the matrix $\mathbf{W}_{\sigma}(\mathbf{x})\mathbf{W}_{\sigma}(\mathbf{x})^{T}$ and $g$ is the Frobenius norm of $\mathbf{W}_{\sigma}(\mathbf{x})$, see Eq. (4).

In Fig. 7b, we show the geodesic path associated to the Riemannian metric $\mathcal{R}_{\rm align}$ on a color image from the Grabcut dataset (Rother et al.,, 2004). The red and yellow dots respectively denote the given source and end points. In Fig. 7c, we show the corresponding geodesic distance map superimposed on the original image.

In Fig. 8, we illustrate the Riemannian alignment minimal paths with respect to different values of the parameter $\beta$. In the first row, from columns $1$ to $3$, the minimal paths indicated by the blue lines are obtained by setting $\beta=4/g_{\rm max},\,8/g_{\rm max}$ and $10/g_{\rm max}$, respectively, where $g_{\rm max}$ is a scalar value defined by

[TABLE]

In the first row, we illustrate the respective geodesic distance maps associated to different values of $\beta$. For a small value of $\beta$, we find a shortcut as shown in row $1$ and column $1$. When increasing the values of $\beta$, one can obtain desired minimal paths which pass through the boundaries of the target, as illustrated in columns $2$ to $3$ of Fig. 8.

5 Orientation-lifted Randers Minimal Paths for Euler-Mumford Elastica Problem

In this section, we denote by $\mathbb{S}^{1}=[0,2\pi)$ an orientation space with periodic boundary condition, based on which we can define an orientation-lifted domain $\mathbb{D}=\Omega\times\mathbb{S}^{1}$. In this case, each point $\hat{\mathbf{x}}=(\mathbf{x},\theta)\in\mathbb{D}$ is an ordered pair including a position $\mathbf{x}$ in the domain $\Omega$ and an orientation $\theta$ in the orientation space $\mathbb{S}^{1}$.

The original Euler-Mumford elastica bending energy functional (Mumford,, 1994; Nitzberg and Mumford,, 1990) assigns to each regular curve $\mathcal{C}\in H^{2}([0,1],\Omega)$ a curvature-dependent length which can be expressed by

[TABLE]

where $\alpha\in\mathbb{R}^{+}$ is a scalar-valued parameter which controls the importance of the curvature $\kappa:[0,1]\to\mathbb{R}$ in the functional $\mathscr{L}_{\rm Elastica}$.

In order to extract the image features such as object boundaries or tubular structure centerlines, we take into account the following data-driven elastica bending energy (Chen et al.,, 2017):

[TABLE]

The function $\Phi:\Omega\times\mathbb{R}^{2}\to\mathbb{R}^{+}$ is a data-driven function relying on both positions and directions. Around the image features of interest and along the proper directions, the orientation-dependent function $\Phi$ is supposed to have small values, and large values, otherwise. The elastica bending energy functional $\mathscr{L}_{\rm Elastica}$ is a geometric variant of the original active contour energy (5) which also takes the curvature-dependent term for regularization. The goal of both the Euler-Mumford elastica model (Mumford,, 1994) and the original active contour model (Kass et al.,, 1988) are relevant to one another, i.e. finding an optimal curve of low curvature integration value, which simultaneously tends to pass through the image features of interest.

For the sake of simplicity, we will set the data-driven function $\Phi\equiv 1$ in Sections 5.1 and 5.2, and present a Randers minimal path solution to the Euler-Mumford Elastica problem (48). The minimization of the bending energy functional $\mathcal{L}_{\rm Elastica}$ including the data-driven function will be introduced in Section 5.3.

5.1 Euler-Mumford Elastica Problem and its Finsler Metric Interpretation

The goal in this section is to search for an optimal curve which minimizes the elastica bending energy functional

[TABLE]

by the Eikonal equation-based minimal path framework. The basic idea is to establish a Finsler geodesic metric to interpret $\mathscr{L}_{\rm Elastica}$ as a weighted curve length. Since the Eikonal equation is a first-order partial differential equation, we should transform the squared curvature term to a first-order term via an auxiliary parametric function of orientations.

We denote by $\varsigma:[0,1]\to\mathbb{S}^{1}$ a parametric function over the interval $[0,1]$, where $\varsigma$ is defined being such that

[TABLE]

By a short computation, one can obtain the following expression

[TABLE]

where ⟂ is a perpendicular operator, i.e. $\vec{\mathbf{u}}^{\perp}$ is a vector perpendicular to $\vec{\mathbf{u}}$.

On the other hand, we have the following equations

[TABLE]

In terms of Eqs. (52) and (5.1), we can express the curvature $\kappa(u)$ by

[TABLE]

Now we have reformulated the curvature $\kappa(\cdot)$ as a ratio of two first-order terms $\varsigma^{\prime}(\cdot)$ and $\|\mathcal{C}^{\prime}(\cdot)\|$. By incorporating the formula in Eq. (54) to the elastica bending length $\mathscr{L}_{\rm Elastica}$, one obtains that

[TABLE]

following the constraint in Eq. (51).

Let $\Gamma=(\mathcal{C},\varsigma)\in\operatorname{Lip}([0,1],\mathbb{D})$ be a canonical orientation-lifted curve and one has $\Gamma^{\prime}(u)=(\mathcal{C}^{\prime}(u),\varsigma^{\prime}(u)),\,\forall u\in[0,1]$. Now the elastica bending energy functional $\mathscr{L}_{\rm elastica}$ of the form (55) can be rewritten by

[TABLE]

where $\mathfrak{F}_{\infty}:\mathbb{D}\times\mathbb{R}^{3}\to[0,+\infty]$ is an orientation-lifted Finsler metric. It is defined for any orientation-lifted point $\hat{\mathbf{x}}=(\mathbf{x},\theta)\in\mathbb{D}$ and any vector $\tilde{\mathbf{u}}=(\vec{\mathbf{u}},\nu)\in\mathbb{R}^{3}$ by

[TABLE]

where $\vec{\mathbf{p}}_{\theta}=(\cos\theta,\sin\theta)^{T}$ is a unit vector associated to an orientation $\theta\in\mathbb{S}^{1}$ and where $\vec{\mathbf{u}}\propto\vec{\mathbf{p}}_{\theta}$ means that the vector $\vec{\mathbf{u}}$ is positively proportional to the vector $\vec{\mathbf{p}}_{\theta}$, i.e., $\vec{\mathbf{u}}=\|\vec{\mathbf{u}}\|\vec{\mathbf{p}}_{\theta}$.

5.2 Finsler Elastica Geodesic Path for Approximating the Elastica Curve

The orientation-lifted Finsler metric $\mathfrak{F}_{\infty}$ defined in Eq. (57) involves the infinite value $+\infty$, which will introduce difficulties to the numerical computation of the geodesic distance maps and the associated minimal geodesic paths. In order to address this problem, here we consider a relaxed orientation-lifted Randers metric $\mathfrak{F}_{\lambda}$ with a scalar factor $\lambda\gg 1$, which can be expressed by

[TABLE]

where $\alpha\in\mathbb{R}^{+}$ is a positive weighing parameter as used in Eq. (57). In the following, we refer to the Randers metric $\mathfrak{F}_{\lambda}$ as the Finsler elastica metric.

The metric $\mathfrak{F}_{\lambda}$ can also be reformulated in the form (32) through a tensor field $\mathcal{M}_{\lambda}:\mathbb{D}\to\mathbb{S}^{+}_{3}$ and a vector field $\omega_{\lambda}:\mathbb{D}\to\mathbb{R}^{3}$, which can be respectively formulated for any point $\forall\hat{\mathbf{x}}=(\mathbf{x},\theta)$ as

[TABLE]

One can claim that the positive definiteness condition formulated in Eq. (33) can be satisfied:

[TABLE]

Now we give the analysis for the relation between the orientation-lifted Finsler metric $\mathfrak{F}_{\infty}$ and the Finsler elastica metric $\mathfrak{F}_{\lambda}$. As $\lambda\to\infty$, we can express the metric $\mathfrak{F}_{\lambda}$ as follows:

[TABLE]

The value of the term $(\lambda-1)(\|\vec{\mathbf{u}}\|-\langle\vec{\mathbf{u}},\vec{\mathbf{p}}_{\theta}\rangle)$ will get to [math] when the directions $\vec{\mathbf{u}}$ and $\vec{\mathbf{p}}_{\theta}$ are positively proportional to each other. Thus, one can see that $\mathfrak{F}_{\lambda}\to\mathfrak{F}_{\infty}$ pointwisely, as $\lambda\to\infty$.

The control set or the unit ball is a basic tool for visualizing the geometry distortion of a metric. For the metrics $\mathfrak{F}_{\infty}$ and $\mathfrak{F}_{\lambda}$, the sets of the unit balls can be respectively defined as

[TABLE]

and

[TABLE]

In order to characterize the sets $B_{\infty}$ and $B_{\lambda}$, we first define two scalar values $\mu_{1}$ and $\mu_{2}$ as follows (Chen et al.,, 2017)

[TABLE]

Hence any vector $\tilde{\mathbf{u}}=(\vec{\mathbf{u}},\nu)\in B_{\infty}(\hat{\mathbf{x}})$ can be characterized by the following inequalities

[TABLE]

Taking into account the Eq. (62), one obtains

[TABLE]

Thus the ball $B_{\infty}(\hat{\mathbf{x}})$ can be characterized by a flat 2D ellipse which is embedded in the 3D tangent space with the origin being on its boundary. In particular, if one sets $\alpha=1$, the unit ball $B_{\infty}(\hat{\mathbf{x}})$ gets to be a flat 2D disk of radius $1/2$ as shown in Fig. 9a.

When the factor $\lambda<\infty$, one can point out that any vector $\tilde{\mathbf{u}}=(\vec{\mathbf{u}},\nu)\in B_{\lambda}(\hat{\mathbf{x}})$ can be characterized by the following inequality

[TABLE]

where the scalar values $a_{\lambda},b_{\lambda},c_{\lambda}$ are all $1+\mathcal{O}(1/\lambda)$. Hence the ball $B_{\lambda}(\hat{\mathbf{x}})$ is an ellipsoid and is almost flat in the direction $\vec{\mathbf{p}}_{\theta}^{\perp}$ due to the large value of $\lambda/2$, as shown in Fig. 9b for which we set $\alpha=1$. One can see that the unit ball $B_{\lambda}(\hat{\mathbf{x}})$ will converge to $B_{\infty}(\hat{\mathbf{x}})$ in the sense of the Haussdorf distance, as $\lambda\to\infty$.

In Fig. 10b, we show the minimal paths derived from the metric $\mathfrak{F}_{\lambda}$. We denote by the red and black dots the source and end positions, respectively. In Fig. 10a, the directions at each position are indicated by arrows.

5.3 Data-driven Finsler Elastica Metric

In order to apply the Finsler elastica minimal paths for image segmentation and tubular structure centerline extraction, we should take into account the image data carried out by the orientation score $\psi:\mathbb{D}\to\mathbb{R}^{+}_{0}$ which is defined over the orientation-lifted domain $\mathbb{D}$.

The data-driven Finsler elastica metric $\mathcal{F}_{\lambda}$ can be expressed by

[TABLE]

or

[TABLE]

where $\beta$ is a positive constant. From Eqs. (65) and (66), we can see that the data-driven function $\Phi$ (see Eq. (49)) is replaced by the orientation score-based functions $\exp(-\beta\psi(\hat{\mathbf{x}})/\|\psi\|_{\infty})$ or $1+\beta\psi(\hat{\mathbf{x}})/\|\psi\|_{\infty}$. Note that in the experiments of this section, we make use of the metric defined in (65) for the data-driven Finsler elastica minimal path computation.

The computation of the orientation score $\psi$ should be task-dependent. We respectively denote by $\psi_{\rm edge}$ and $\psi_{\rm tube}$ the orientation score $\psi$ for image segmentation and for tubular structure centerline extraction, each of which can be estimated by steerable filters (Freeman and Adelson,, 1991).

Let $\mathbf{I}=(I_{1},I_{2},I_{3}):\Omega\to\mathbb{R}^{3}$ be a RGB color image with three channels. For boundary detection and image segmentation applications, we use the $m$-order canny-like steerable filter (Jacob and Unser,, 2004) with $m$ an odd integer to build the orientation score $\phi_{\rm edge}$. The $m$-order canny-like steerable kernel $\mathbf{G}^{\theta}_{m}$ can be expressed by

[TABLE]

where the function $G_{\sigma}$ is a Gaussian kernel with variance $\sigma$. The coefficients $f_{a,b}$ are independent of the positions $\mathbf{x}$ but relying on the orientation $\theta$ and we refer to (Jacob and Unser,, 2004) for the computation of $f_{a,b}$. Based on the kernels $\mathbf{G}^{\theta}_{m}$, the orientation score $\psi_{\rm edge}$ can be constructed by

[TABLE]

The orientation score $\psi_{\rm tube}$ can be computed through a multi-scale tubular structure enhancing filter such as the classical vesselness filter (Frangi et al.,, 1998) or the optimally oriented flux (OOF) filter  (Law and Chung,, 2008). In this section, we use the OOF filter to derive the orientation score $\psi_{\rm tube}$, where the kernel $\mathbf{F}_{r}$ at a scale $r$ can be expressed as

[TABLE]

where $\mathbf{1}_{r}$ is the indicator function of a disk structure with radius $r$. In other words, the function $\mathbf{1}_{r}$ can be defined as a step function relying on a radius $r$ such that

[TABLE]

Then we can define a multi-scale response function at position $\mathbf{x}$ and scale $r$ for a given image $I:\Omega\to\mathbb{R}$ through the kernel $\mathbf{F}_{r}$

[TABLE]

Note that the response $\mathbf{Q}(\mathbf{x},r)$ is a symmetric matrix of size $2\times 2$. Let us denote by $\eta_{1}(\mathbf{x},r),\,\eta_{2}(\mathbf{x},r)$ the two eigenvalues of the matrix $\mathbf{Q}(\mathbf{x},r)$ obeying that $\eta_{2}(\cdot)\geq\eta_{1}(\cdot)$. We suppose that the gray levels inside the tubular structure are lower than background. In this case, for each point $\mathbf{x}$ we can estimate an optimal scale map $R_{*}:\Omega\to\mathbb{R}^{+}$ by

[TABLE]

Then for a point $\mathbf{x}$ inside the vessel region, the eigenvalues of the matrix $\mathbf{Q}(\mathbf{x},R_{*}(\mathbf{x}))$ satisfy that $\eta_{2}(\mathbf{x},R_{*}(\mathbf{x}))\gg\eta_{1}(\mathbf{x},R_{*}(\mathbf{x}))\approx 0$.

The orientation score $\psi_{\rm tube}$ can be computed by

[TABLE]

where $\vec{\mathbf{n}}_{\theta}=(-\sin\theta,\cos\theta)^{T}$ is a unit vector associated to an orientation $\theta\in\mathbb{S}^{1}$.

The minimal geodesic paths derived from the data-driven Finsler elastica metric should be as smooth as possible and simultaneously follow the desired image features. In Fig. 11, we illustrate an example for boundary detection using the data-driven Finsler elastica metric. In Fig. 11a, the arrows indicate the directions at the given positions (denoted by red and yellow dots). In Fig. 11b, the obtained minimal geodesic path indicated by blue line can detect an image edge of high Euclidean length, which satisfies the properties of the Finsler elastica minimal paths as discussed above.

In Fig. 12, we show the geodesic paths derived from the data-driven Finsler elastica metric on a synthetic image. The arrows in Figs. 12a and 12c denote the corresponding directions assigned to the given positions indicated by the red and yellow dots. The geodesic paths shown in Figs. 12b and 12d are computed using the same source and end positions but opposite directions. One can see that these geodesic paths tend to pass through the centerline of the target tubular structure, and simultaneously try to avoid sharp turnings as much as possible.

In Fig. 13, we apply the Finsler elastica minimal path model for retinal vessel centerline tracking. In this experiment, only the positions are provided. Each position will be assigned two directions corresponding to the orientations maximizing the orientation score $\psi_{\rm tube}$ over the space $[0,2\pi)$. In this case, the source position (resp. end position) will generate two orientation-lifted source (resp. end) points. Hence we will obtain four pairs of orientation-lifted source and end points, which correspond to four geodesic paths. Among them, the geodesic path with the smallest geodesic distance value will be chosen as the desired path, see Fig. 13b. For the purpose of comparison, we show the minimal paths respectively derived from the isotropic Riemannian metric (Li and Yezzi,, 2007) and the anisotropic Riemannian metric (Benmansour and Cohen,, 2011), as shown in Figs. 13c and Fig. 13d. One can point out that these paths from the Riemannian metrics travel the segments belonging to different vessels, while the path from the Finsler elastica metric can track the correct vessel. Note that both of the Riemannian metrics used in this experiment are constructed by the OOF filter (Law and Chung,, 2008).

6 Randers Minimal Paths for Region-based Active Contours

In this section, we introduce a new minimal path model (Chen et al.,, 2016, 2019) for solving the active contour problem involving a region-based homogeneity penalization term, using the minimal path and Eikonal PDE framework.

6.1 Hybrid Active Contour Model

A hybrid active contour model (Cohen,, 1997; Kimmel and Bruckstein,, 2003; Paragios and Deriche,, 2002; Sagiv et al.,, 2006) invokes an energy functional that is made up of a region-based homogeneity penalization functional $\mathscr{H}_{\rm region}$ and an edge-based weighted curve length $\mathscr{L}_{\rm IR}$ (see Eq. (23)), i.e.

[TABLE]

where $\alpha\in\mathbb{R}^{+}$ is a parameter that controls the relative importance of the two terms $\mathscr{H}_{\rm region}$ and $\mathscr{L}$, $\mathcal{C}\in\operatorname{Lip}([0,1],\Omega)$ is a simple and closed curve and $\Omega$ represents the image domain. The function $\mathbf{1}_{A_{\mathcal{C}}}:\Omega\to\{0,1\}$ is the characteristic function of the open and bounded subset $A_{\mathcal{C}}$ enclosed by $\mathcal{C}$, which can be defined by

[TABLE]

The weighted curve length $\mathscr{L}_{\rm IR}$ plays the role for regularization, which can be measured using an edge-based potential $P:\Omega\to\mathbb{R}^{+}$. In this section, we make use of the following method to construct the potential

[TABLE]

where $\beta\in\mathbb{R}^{+}$ is a contrast parameter and the function $g$ represents the magnitude of the image gradient vector field, as defined in Eqs. (2) and (4). One can see that by Eq. (72) the potential $P$ obeys that $\inf_{\mathbf{x}\in\Omega}\{P(\mathbf{x})\}=1$ and appears to have low values around the image edges.

Let $\tilde{\mathcal{C}}\in\operatorname{Lip}([0,1],\Omega)$ be a fixed simple and closed curve that is parameterized in a clockwise order. Supposing that the region $A_{\mathcal{C}}$ is close to $A_{\tilde{\mathcal{C}}}$ and by the differentiability of the region-based functional $\mathscr{H}_{\rm region}$, we can obtain the following approximation (Chen et al.,, 2019)

[TABLE]

where $\rho_{\tilde{\mathcal{C}}}$ is the gradient of the functional $\mathscr{H}_{\rm region}$ at $\mathbf{1}_{A_{\tilde{\mathcal{C}}}}$ and $c\in\mathbb{R}$ is a scalar value that is independent to the characteristic function $\mathbf{1}_{A_{\mathcal{C}}}$.

One can point out that only the second term in the right hand side of Eq. (73) relies on the characteristic function $\mathbf{1}_{A_{\mathcal{C}}}$ of the region $A_{\mathcal{C}}$. For convenience, let us denote this term by

[TABLE]

In our model, solving the hybrid active contour problem through the curve evolution scheme amounts to iteratively searching for a family of successive simple and closed curves $\mathcal{C}_{k}$ (indexed by $k$) such that the final curve $\mathcal{C}_{\infty}$ can be used to depict the target boundary.

In the $k$-th iteration ($k\geq 1$), the input is a curve $\mathcal{C}_{k}$ obtained from the last iteration and the output is a new simple and closed curve $\mathcal{C}_{k+1}$. By setting $\tilde{\mathcal{C}}:=\mathcal{C}_{k}$, we can obtain the gradient $\rho_{\mathcal{C}_{k}}$ of the functional $\mathscr{H}_{\rm region}$ at $\mathbf{1}_{\mathcal{C}_{k}}$. Then the curve $\mathcal{C}_{k+1}$ can be generated by solving

[TABLE]

where $\Xi_{\mathcal{C}_{k}}$ is a set of simple and closed curves. The construction of $\Xi_{\mathcal{C}_{k}}$ relies on $\mathcal{C}_{k}$ and will be introduced in the next section. In order to solve the problem (75) by the minimal path model and the Eikonal PDE, we transform the region-based functional $\mathcal{J}_{\mathcal{C}_{k}}(\mathbf{1}_{A_{\mathcal{C}}})$ into a weighted curve length associated to a Randers metric using the divergence theorem, providing that the curve $\mathcal{C}_{k+1}$ is chosen from the set $\Xi_{\mathcal{C}_{k}}$.

6.2 A Randers Metric Interpretation to the Hybrid Energy

In the $k$-th iteration, let $U_{\mathcal{C}_{k}}=\{\mathbf{x}\in\Omega;d(\mathbf{x},\mathcal{C}_{k})<r\}$ be a tubular neighbourhood of the simple and closed curve $\mathcal{C}_{k}$ with radius $r$, where $d(\mathbf{x},\mathcal{C})=\min_{u}\|\mathbf{x}-\mathcal{C}(u)\|$ denotes the Euclidean distance between $\mathbf{x}$ and $\mathcal{C}$. Then the set $\Xi_{\mathcal{C}_{k}}$, which is made up of simple and closed curves, can be defined by

[TABLE]

For any simple and closed curve $\mathcal{C}\in\Xi_{\mathcal{C}_{k}}$, one has $A_{\mathcal{C}}\backslash U_{\mathcal{C}_{k}}=A_{\mathcal{C}_{k}}\backslash U_{\mathcal{C}_{k}}$. This means that the region $A_{\mathcal{C}}$ can be decomposed as

[TABLE]

In this case, we can reformulate the region-based functional $\mathcal{J}_{\mathcal{C}_{k}}$ as follows:

[TABLE]

where $\mathbf{1}_{U_{\mathcal{C}_{k}}}$ is defined as the characteristic function of the tubular neighbourhood $U_{\mathcal{C}_{k}}$, $\mathscr{N}$ is the outward normal to the curve $\mathcal{C}$, and $\vartheta_{\mathcal{C}_{k}}:\mathbb{R}^{2}\to\mathbb{R}^{2}$ is a vector field which satisfies the following divergence equation over the domain $\mathbb{R}^{2}$

[TABLE]

One can always find such a vector field $\vartheta_{\mathcal{C}_{k}}$ due to the existence of the solution to the divergence equation (77).

In Eq. (76), the second term is independent to the curve $\mathcal{C}$ and thus we only take into account the first term in the following calculation. Let $\mathbf{M}$ be the rotation matrix of angle $\pi/2$. Integrating the weighted curve length $\mathscr{L}_{\rm IR}$ and the first term in Eq. (76), we obtain that

[TABLE]

where the metric $\mathfrak{F}_{\mathcal{C}_{k}}$ has a form of

[TABLE]

The positive definiteness condition (33) for $\mathfrak{F}_{\mathcal{C}_{k}}$ requires that

[TABLE]

However, the requirement (80) is difficult to satisfy. Hence we consider a new inequality formulated by

[TABLE]

Now we have induced a Randers metric formulated in Eq. (79), which embeds the region-based homogeneity information as well as the image gradient features. This gives us the possibility of applying the minimal path framework to solve the region-based image segmentation problem. In next section, we will introduce a new vector field to approximate $\vartheta_{\mathcal{C}_{k}}$, such that the inequality (81) will always hold.

6.3 Practical Implementations

In this section, we first introduce a method for computing the vector field $\vartheta_{\mathcal{C}_{k}}$ used in the $k$-th iteration. In order to satisfy (81), the vector field $\vartheta_{\mathcal{C}_{k}}$ is expected to be as small as possible. Since the next evolving curve lies within the tubular neighbourhood $U_{\mathcal{C}_{k}}$, we consider to find a vector field $\vartheta_{\mathcal{C}_{k}}:U_{\mathcal{C}_{k}}\to\mathbb{R}^{2}$ by solving the following minimization problem on the domain $U_{\mathcal{C}_{k}}$

[TABLE]

As discussed in (Chen et al.,, 2016, 2019), the value of $\sup_{\mathbf{x}}\|\vartheta_{\mathcal{C}_{k}}(\mathbf{x})\|$ is bounded by the area of the tubular neighbourhood $U_{\mathcal{C}_{k}}$. Since the tubular neighbourhood $U_{\mathcal{C}_{k}}$ acts as the search space for the evolving curves, reducing the area of $U_{\mathcal{C}_{k}}$ may lead the proposed model being stuck in unexpected local minima. In order to address this issue, we make use of a new vector field $\varpi_{\mathcal{C}_{k}}:U_{\mathcal{C}_{k}}\to\mathbb{R}^{2}$ derived by the nonlinear mapping $\varphi$ defined in Eq. (42) such that

[TABLE]

where $\tilde{\alpha}\in\mathbb{R}^{+}$ is a constant.

Replacing the vector field $\alpha\vartheta_{\mathcal{C}_{k}}$ (see Eq. (79)) by $\varpi_{\mathcal{C}_{k}}$, one can obtain a new Randers metric:

[TABLE]

In each iteration $k$, the weighted curve length measured along a curve $\mathcal{C}\in\Xi_{\mathcal{C}_{k}}$ can be formulated by

[TABLE]

For any point $\mathbf{x}\in U_{\mathcal{C}_{k}}$ the vector $\vartheta_{\mathcal{C}_{k}}(\mathbf{x})$ is positively proportional to $\varpi_{\mathcal{C}_{k}}(\mathbf{x})$ in the sense of the magnitude and of the direction. This gives the relevance between the original region-based active contour problem and the minimal path computation associated to the geodesic metric $\mathcal{F}_{\mathcal{C}_{k}}$. More analysis for the reasonability of the use of the Randers metric $\mathcal{F}_{\mathcal{C}_{k}}$ can be found in (Chen et al.,, 2019).

6.4 Application to Image Segmentation

We make use of the piecewise constant homogeneity term (Chan and Vese, 2001a, ) to build the Randers metric $\mathcal{F}_{\mathcal{C}_{k}}$ in Eq. (84) and illustrate how to apply $\mathcal{F}_{\mathcal{C}_{k}}$ for image segmentation. The objective is to search for a family of successive $\mathcal{C}_{k}$, each of which is simple and closed, to depict the object boundary, as $k\to\infty$.

The region-based homogeneity term in the piecewise constant variant of the Mumford-Shah model can be expressed as follows

[TABLE]

Recall that in the $k$-th iteration of the curve evolution scheme, the input is the curve $\mathcal{C}_{k}$ ($k\geq 1$) and the output is $\mathcal{C}_{k+1}$. Note that when $k=1$, $\mathcal{C}_{1}$ is the initial curve. We can estimate the gradient $\rho_{\mathcal{C}_{k}}$ of the region-based functional $\mathscr{H}_{\rm region}$, which can be formulated by

[TABLE]

where $\mu_{\rm in}[\mathcal{C}_{k}]$ and $\mu_{\rm out}[\mathcal{C}_{k}]$ are respectively the mean gray levels inside and outside the curve $\mathcal{C}_{k}$.

Based on the gradient $\rho_{\mathcal{C}_{k}}$ and the tubular neighbourhood $U_{\mathcal{C}_{k}}$, one can solve the PDE-constrained minimization problem (82) to obtain the vector field $\vartheta_{\mathcal{C}_{k}}$ and also to obtain $\varpi_{\mathcal{C}_{k}}$. Finally, the Randers metric $\mathcal{F}_{\mathcal{C}_{k}}$ can be constructed using Eq. (84).

Interactive image segmentation method. By solving the gradient descent ODE, one can obtain an open geodesic path between two points within a given domain. However, for image segmentation, the objective is to seek simple and closed curves. Thus we need to solve an issue: how to derive a closed and simple curve $\mathcal{C}_{k+1}$ in each iteration $k$, providing that the Randers metric $\mathcal{F}_{\mathcal{C}_{k}}$ is given. In this section, we consider an interactive segmentation scheme to obtain the evolving curve $\mathcal{C}_{k+1}$, which depends on a set $\{\mathbf{q}_{i}\}_{1\leq i\leq m}$ of $m$ ($m\geq 3$) prescribed vertices distributed on the target boundary in a clockwise order. These vertices will be fixed in the course of curve evolution. Each evolving curve $\mathcal{C}_{k}$ is supposed to pass through each of these vertices in a clockwise order. The initial curve $\mathcal{C}_{1}$ is a polygon which is generated using the vertices $\mathbf{q}_{i}$. The simple and closed curve $\mathcal{C}_{k+1}$ can be generated by the end-to-end concatenation of a set of geodesic paths obtained using $\mathcal{F}_{\mathcal{C}_{k}}$. Each of these geodesic paths, denoted by $\mathcal{G}_{i,k+1}$ with $1\leq i\leq m$, connects a vertex $\mathbf{q}_{i}$ to its successive one.

Now we present the method for computing the geodesic paths $\mathcal{G}_{i,k+1}$ to obtain the curve $\mathcal{C}_{k+1}$. Recall that the input curve $\mathcal{C}_{k}$ passes through all of the fixed vertices $\{\mathbf{q}_{i}\}_{1\leq i\leq m}$. In this case, we are able to decompose $\mathcal{C}_{k}$ into a set $\{\mathcal{G}_{i,k}\}_{i}$ of $m$ subpaths. Each path $\mathcal{G}_{i,k}$, defined over the integral $[0,1]$, links a vertex $\mathbf{q}_{i}$ to its successive one $\mathbf{q}_{i+1}$ for $i<m$ or to $\mathbf{q}_{1}$ for $i=m$. Following that the tubular neighbourhood $U_{\mathcal{C}_{k}}$ can be decomposed into a family of disjoint subdomains $\mathscr{U}_{i}\subset U_{\mathcal{C}_{k}}$, each of which can be regarded as a narrowband region (with radius $r$) involving the open curve $\mathcal{G}_{i,k}$. In essence, these regions $\mathscr{U}_{i}$ ($i=1,2,\cdots m$) can be identified by using the Voronoi index map over the domain $U_{\mathcal{C}_{k}}$ associated to the open curves $\mathcal{G}_{i,k}$ ($1\leq i\leq m$) as introduced in (Chen et al.,, 2019). In other words, a point $\mathbf{x}\in\mathscr{U}_{i}$ implies that $d(\mathbf{x},\mathcal{G}_{i,k})<d(\mathbf{x},\mathcal{G}_{j,k}),\,\forall j\neq i$.

With these definitions in hands, we can track a geodesic path $\mathcal{G}_{i,k+1}$ in the corresponding region $\mathscr{U}_{i}\cup\{\mathbf{q}_{i},\mathbf{q}_{i+1}\}$ if $i<m$, or $\mathscr{U}_{m}\cup\{\mathbf{q}_{m},\mathbf{q}_{1}\}$ if $i=m$. The use of these regions $\mathscr{U}_{i}\cup\{\mathbf{q}_{i},\mathbf{q}_{i+1}\}$ for geodesic path computation is to avoid the curve self-crossing issue and also to reduce the computation time. The decomposition of the tubular neighbourhood $U_{\mathcal{C}_{k}}$ can be efficiently implemented by the Fast Marching method with a Voronoi index map estimation procedure. Once all of the geodesic paths $\mathcal{G}_{i,k+1}\,(1\leq i\leq m)$ are generated, we can construct a simple and closed curve $\mathcal{C}_{k+1}$ by the concatenation of the paths $\mathcal{G}_{i,k+1}$ in an end-to-end manner. Finally, the curve $\mathcal{C}_{k+1}$ is regarded as the output of the $k$-th iteration.

Note that for a pair of successive vertices $(\mathbf{q}_{i},\mathbf{q}_{i+1})$, the extracted geodesic path $\mathcal{G}_{i,k+1}$ is the globally optimal curve for the length $\operatorname{Length}_{\mathcal{F}_{\mathcal{C}_{k}}}$ (see Eq. (85)) within the domain $\mathscr{U}_{i}\cup\{\mathbf{q}_{i},\mathbf{q}_{i+1}\}$.

In Fig. 14, we show an example for this interactive image segmentation procedure providing that a set of $3$ ordered vertices are given. These vertices are illustrated by red dots, see Fig. 14a. In Fig. 14b, the blue straight segment lines connecting each pair of red dots represent the initial curve. In Figs. 14c to 14e, we show the intermediate curve evolution results, where the blue contour is generated by the concatenation of a set of the Randers minimal paths. The final curve is shown in Fig. 14f. Again, the vertices denoted by red dots are fixed during the curve evolution.

Remark 1**.**

The method presented in this section assumes that the vertices are provided by user and are fixed during the curve evolution. As a matter of fact, in each iteration $k$, these vertices can be sampled from the input curve $\mathcal{C}_{k}$. In this case, one can initialize the proposed model by providing a closed curve and the vertices will evolve in the course of the curve evolution processing (Chen et al.,, 2016).

Remark 2**.**

Note that a variant of the edge-based balloon active contour model can also be addressed by the framework proposed in this section. The balloon force $\mathbf{F}_{\rm balloon}$ in Eq. (8) can be obtained by minimizing the following term

[TABLE]

yielding a curve evolution flow

[TABLE]

Recall that $\mathscr{N}_{\tau}$ is the outward normal to $\mathcal{C}_{\tau}$.

Now we can establish a variant form of the isotropic geodesic active contour functional (Caselles et al.,, 1997; Yezzi et al.,, 1997) with the balloon force term

[TABLE]

In this case, for a given simple and closed curve $\tilde{\mathcal{C}}$, the gradient of $E_{\rm balloon}$ thus turns out to be a constant function, i.e. $\rho_{\tilde{\mathcal{C}}}\equiv 1$. This edge-based active contour model thus can be solved by the Randers metric-based model proposed in this section.

Furthermore, an alternative choice for the regularization term in Eqs. (71) and (89) can be constructed using an anisotropic Riemannian metric $\mathfrak{R}$ (see Eq. (35)), which allows to take into account the path directions for image segmentation (Chen and Cohen,, 2017; Chen et al.,, 2019).

7 Conclusion

In this chapter, we review the methods for the construction of the Riemannian and Randers metrics in terms of various active contour terms. We show that by these metrics the edge- and region-based active contour problems and the Euler-Mumford elastica problem can be efficiently solved by the minimal path framework based on the Eikonal partial differential equation. Moreover, we also exploit the relevant applications of minimal paths associated to the proposed metrics in image analysis such as image segmentation, boundary detection as well as tubular structure centerline tracking, which are able to blend the benefits from the minimal path framework and the active contour models. The well-established fast marching methods, which are the Eikonal solvers, allow efficient and practical implementations for these applications.

Acknowledgement

The authors thank Dr. Jean-Marie Mirebeau from Université Paris-Saclay for his fruitful discussion and regular collaboration. This research has been partially funded by Roche pharma (project AMD_short) and by a grant from the French Agence Nationale de la Recherche ANR-16-RHUS-0004 (RHU TRT_cSVD). Figures 9 are reprinted by permission from Springer Nature Customer Service Centre GmbH: Springer Nature, International Journal of Computer Vision, Global minimum for a Finsler elastica approach, Da Chen, Jean-Marie Mirebeau and Laurent D. Cohen (Chen et al.,, 2017).