Keywords: Circular geodesic, fast marching method, image segmentation, tangent-constrained prior, curvature penalization

Highlights↩︎

• Unified tangent-constrained curvature-penalized geodesic models: We introduce a family of tangent-constrained Finslerian metrics that extend classical curvature-penalized geodesic models (including Reeds–Shepp, Dubins, Euler–Mumford elastica and curvature prior elastica models). The proposed formulation enforces hard tangent constraints while remaining fully compatible with curvature regularization.

• HJB PDE formulation with efficient numerical solvers: The resulting Hamilton–Jacobi–Bellman equations can be efficiently solved using variants of Hamiltonian fast marching methods, preserving single-pass computational complexity and enabling practical application to large-scale image segmentation tasks.

• Practical segmentation pipeline with deep-learning-based initialization: We present an automatic segmentation framework that leverages neural network pre-segmentation to extract intrinsic shape representatives, which are then used to construct tangent priors and initialize circular geodesic models. Extensive experiments on synthetic, natural, and medical images demonstrate that the proposed approach significantly improves robustness to weak boundaries and complex object geometries, outperforming classical geodesic models as well as representative graph-based and learning-based baselines.

1 Introduction↩︎

Image segmentation remains a fundamental challenge in computer vision and medical imaging [1]–[6]. Despite significant advances, existing segmentation approaches are still struggling with weak boundaries, noise, complex object geometries and image intensity distributions. Both global and local shape constraints have proven essential for achieving suitable segmentations across various scenarios [7]. In practice, shape constraints serve as effective geometry priors in deep learning-based models, energy minimization-based geometric models and hybrid frameworks combining both paradigms.

Energy minimization is a flexible framework for addressing image segmentation, capable of encoding various types of image features, geometric regularization, and shape priors into a unified objective functional. Existing energy minimization models can be broadly divided into two categories: graph-based models (e.g., [8]–[12]) and variational models (e.g., [13]–[16]). Graph-based segmentation models cast the problem as a discrete optimization problem, which can be solved efficiently using algorithms such as graph cuts, random walks or power watershed [17]–[19]. However, they typically approximate geometric quantities using finite-order clique potentials defined over a discrete lattice, often leading to metrication artifacts and limiting the incorporation of high-order geometric priors [20]–[22]. In contrast, variational models address segmentation within a PDE framework, providing a continuous, geometry-preserving representation that naturally accommodates curvature, other high-order geometric regularization terms and shape priors, without being restricted to discrete structures. The proposed segmentation method falls into the latter category, leveraging the PDE-based formulation to integrate curvature penalization and local shape constraints in a principled manner.

The minimal path model [23] was first introduced to obtain the global minimization of a path energy associated with a suitable isotropic Riemannian metric, establishing the connection between the PDE and the differential geometry. It was originally investigated to track image edges and tubular structures. These applications have inspired a series of works to address the limitations of this original model, such as the integration of anisotropic and/or asymmetric metrics [24]–[29]. Among those geodesic approaches, a significant research line is the curvature-penalized models, which incorporate second-order curvature terms for regularization, and have demonstrated promising results in several practical applications [30]–[34]. These models must address a fundamental issue of incorporating a second-order curvature regularization within an inherently first-order HJB PDE framework. The orientation-lifting strategy provides an elegant resolution by introducing an auxiliary angular dimension, thereby expressing path curvature as a ratio of first-order derivatives. Existing second-order geodesic models either employ curvature as a regularization term or leverage it for tracking circular minimal paths with a convexity shape prior. However, they cannot impose mandatory directional constraints on the path tangents to accommodate more flexible, spatially varying shape constraints. Fig. 1 illustrates a fundamental limitation inherent to classical curvature-penalized models when applied to complicated segmentation scenarios. We consider examples of segmenting a spiral using three representative models: the Reeds-Shepp forward (RSF) model [30], the Dubins model [31] and the Euler-Mumford elastica model [34]. As shown in Figs. [fig_intro-rsf] to [fig_intro-em], all three models produce geodesics that take undesirable shortcuts. This shortcutting behavior arises because these models lack any mechanism to enforce directional consistency with the underlying shape geometry. In contrast, Figs. [fig_intro-tc-rsf] to [fig_intro-tc-em] present the segmentation results obtained by integrating tangent-constrained priors into the same curvature-penalized geodesic frameworks. The resulting models⁷ successfully delineate the complete spiral boundary without shortcuts, demonstrating that the proposed tangent constraints effectively resolve this limitation while preserving the benefits of curvature regularization.

In this work, we propose a unified geodesic framework that integrates tangent-constrained priors and curvature penalization within a single PDE formulation. The tangent constraints are derived from the ISR and formulated as orientation-dependent costs in a high-dimensional space, enabling the use of efficient single-pass Hamiltonian Fast-Marching methods [31], [35] or the GPU-accelerated solvers for the HJB PDE [36]. Compared to the convexity shape prior, the proposed tangent-constrained priors are less restrictive and thus applicable to a broader range of segmentation scenarios. To the best of our knowledge, the most relevant works to our model are the star convexity-embedded approaches [37]–[41]. Among these, the star shape constraint is either incorporated into graph-based models, or integrated with neural networks to guide the segmentation output, both of which have achieved notable success. However, the graph-based models are susceptible to metrication errors and anisotropic bias, while neural network-based approaches rely on soft relaxations of geometric constraints, offering no theoretical guarantee that the output satisfies the star convexity property. Moreover, neither class of models incorporates curvature penalization, which plays a crucial role in the presence of strong noise, weak edges and complex image content. This provides the motivation for this work, where the main contributions are two-fold:

• We propose a unified HJB PDE framework for computing tangent-constrained and curvature-penalized circular optimal paths of minimum arrival time, where the tangent priors are derived from the ISR that serves as local shape constraints. We present an efficient method that transforms the vector-valued tangent priors as a scalar-valued field over the orientation-lifted space, allowing to invoke the well-established numerical solvers (e.g. Hamiltonian fast marching method [31] or fast iterative method [36]) for computing satisfactory minimal paths.

• A practical pipeline that leverages deep learning-based pre-segmentation to extract the ISR (e.g., skeletons or interior landmarks) and to provide automatic initialization for image segmentation based on the proposed geodesic models. This enables automatic segmentation while ensuring that the resulting circular geodesics encode both image features and tangent constraints, thereby effectively preventing unexpected shortcuts in challenging scenarios.

The manuscript is organized as follows. Section 2 reviews the curvature-penalized minimal path models based on the static HJB PDE framework. Section 3 presents the core of the proposed model: the unified tangent-constrained and curvature-penalized geodesic framework. Section 4 describes the applications of the introduced geodesic models to image segmentation. Section 5 reports experimental results and Section 6 concludes the paper.

2 Curvature-penalized Minimal Paths↩︎

In this section, we revisit the curvature-penalized minimal path models upon the HJB PDE framework, which are posed in an orientation-lifted space [30], [31], [33], [34].

2.1 Orientation Lifting for Path Tangents and Curvature↩︎

To represent curvature within the energy minimization framework, the open bounded and connected domain \(\Omega \subset \mathbb{R}^{2}\) is lifted to an orientation-lifted space \(\mathbb{M}:= \Omega \times \mathbb{S}^{1}\), where \(\mathbb{S}^1 := \mathbb{R}/2\pi\mathbb{Z}\) is an angle space with period boundary condition. Any point \(\mathbf{x}\in \mathbb{M}\) is a pair \((x, \theta)\) consisting of a spatial coordinate \(x \in \Omega\) and an angular coordinate \(\theta \in \mathbb{S}^{1}\). The tangent space at any point \(\mathbf{x}\) is given by \(\mathbb{E}:= \mathbb{R}^{2} \times \mathbb{R}\), with tangent vectors denoted by \({\dot{\mathbf{x}}}= (\dot{x}, {\dot{\theta}})\).

A smooth planar curve \(\gamma:[0,1] \to \Omega\) with a non-vanishing velocity (i.e., \(\|\gamma^\prime(\cdot)\|\neq 0\)) can be lifted to a curve \(\Gamma:[0,1] \to \mathbb{M}\), defined as \(\Gamma = (\gamma, \eta)\). The angular component \(\eta(t) \in \mathbb{S}^{1}\) represents the orientation of the curve \(\gamma\) for any \(t\in[0,1]\). Consequently, the tangent vector of the planar curve is constrained to be aligned with its orientation: \[\label{eq:first95order95gamma} \gamma^\prime(t) = \|\gamma^\prime(t)\| \mathop{\mathrm{\mathbf{n}}}(\eta(t)),\tag{1}\] where \(\mathop{\mathrm{\mathbf{n}}}(\theta) = (\cos\theta, \sin\theta)\) denotes the unit vector associated with the angle \(\theta\). The first-order derivative of the orientation-lifted curve is then given by \(\Gamma^\prime(t)= \big(\gamma^\prime(t), \eta^\prime(t)\big),\,\forall t\in[0,1]\).

This orientation-lifted representation allows for a direct formulation of the curvature \(\kappa(t)\) of \(\gamma\). It is defined as the ratio of two first-order terms \[\kappa(t) := \frac{\eta^\prime(t)}{\|\gamma^\prime(t)\|},\] allowing to incorporate the curvature terms into the first-order static HJB PDE framework.

2.2 Curvature-penalized Geodesic Models↩︎

The bending energy of a regular curve \(\gamma : [0,1] \to \Omega\) is formulated to include a penalty on its curvature, which generalizes the Euclidean curve length as follows \[\label{eq:energy} \int_{0}^{1} \psi(\gamma(t), \eta(t))\; \mathcal{C}(\xi\kappa(t)) \|\gamma^\prime(t)\|\, dt,\tag{2}\] where \(\psi:\mathbb{M}\to \mathbb{R}^{+}\) is an image data-driven cost function, \(\mathcal{C}:\mathbb{R}\to [0,\infty]\) is a cost function of curvature dependent on the specific geodesic models⁸, and \(\xi\in\mathbb{R}^+\) is a constant modulating the strength of curvature penalization. The term \(\eta(t)\) denotes the turning angle of \(\gamma(t)\) satisfying the assumption 1 .

In the orientation-lifted space, this energy corresponds to the length of the lifted path \(\Gamma(t) = (\gamma(t), \eta(t))\): \[\mathcal{L}(\Gamma) := \int_{0}^{1} \psi(\Gamma(t))\mathcal{F}(\Gamma(t),\Gamma^\prime(t)) \, dt,\] where \(\mathcal{F}:\mathbb{M}\times \mathbb{E}\to [0,\infty)\) is an orientation-lifted Finsler metric. For a point \(\mathbf{x}=(x,\theta) \in \mathbb{M}\) and a tangent vector \(\dot{\mathbf{x}}=(\dot{x},\dot{\theta}) \in \mathbb{E}\), the metric is defined as: \[\label{eq:metric} \mathcal{F}(\mathbf{x}, \dot{\mathbf{x}}) = \begin{cases} \mathcal{C}\bigl(\xi {\dot{\theta}}/\|\dot{x}\|\bigr)\|\dot{x}\|, & \text{if~} \dot{x}= \mathop{\mathrm{\mathbf{n}}}(\theta)\|\dot{x}\|, \\ \infty, & \text{otherwise}. \end{cases}\tag{3}\] The infinite penalty enforces the constraint that the spatial velocity \(\dot{x}\) must be aligned with the orientation vector \(\mathop{\mathrm{\mathbf{n}}}(\theta)\).

The geodesic path between a source point \(\mathbf{p}\in \mathbb{M}\) and a target point \(\mathbf{x}\in\mathbb{M}\) is found by computing the minimal action, also known as the geodesic distance map \(\mathcal{U}_{\mathbf{p}}\): \[\label{eq:distance95map} \mathcal{U}_{\mathbf{p}}(\mathbf{x}) = \inf_{\Gamma \in \mathrm{Lip}([0,1], \mathbb{M})} \left\{ \mathcal{L}(\Gamma) \,;\, \Gamma(0) = \mathbf{p},\;\Gamma(1) = \mathbf{x} \right\}.\tag{4}\] This distance map \(\mathcal{U}_{\mathbf{p}}\) is the unique viscosity solution of a static Hamilton-Jacobi equation. The equation is defined in terms of the Hamiltonian \(\mathcal{H}\) associated with the metric \(\mathcal{F}\): \[\label{eq95HJBPDE95Curvature} \mathcal{H}(\mathbf{x},d\mathcal{U}_{\mathbf{p}}(\mathbf{x}))= \frac{1}{2}\,\psi(\mathbf{x})^{2},\quad\forall \mathbf{x}\in \mathbb{M}\backslash \{\mathbf{p}\},\tag{5}\] subject to \(\mathcal{U}_{\mathbf{p}}(\mathbf{p})=0\) and an outflow boundary condition on \(\partial \mathbb{M}\). Here, \(d\mathcal{U}_{\mathbf{p}}\) represents the differential of the geodesic distance map \(\mathcal{U}_\mathbf{p}\). The Hamiltonian \(\mathcal{H}\) is obtained through the Legendre-Fenchel duality of the metric \(\mathcal{F}\): \[\label{eq:hamiltonian} \mathcal{H}(\mathbf{x},{\hat{\mathbf{x}}}):= \sup_{{\dot{\mathbf{x}}}\in \mathbb{E}} \left\{ \left\langle \hat{\mathbf{x}}, \dot{\mathbf{x}} \right\rangle - \frac{1}{2}\,\mathcal{F}(\mathbf{x}, \dot{\mathbf{x}})^{2} \right\},\tag{6}\] for any \(\mathbf{x}=(x,\theta) \in \mathbb{M}\) and any \({\hat{\mathbf{x}}}=({\hat{x}},{\hat{\theta}}) \in \mathbb{R}^{2} \times \mathbb{R}\).

Once the geodesic distance map \(\mathcal{U}_{\mathbf{p}}\) is computed, the geodesic path \(\mathcal{G}_{\mathbf{p},\mathbf{x}}\) from the source \(\mathbf{p}\) to any target point \(\mathbf{x}\in \mathbb{M}\) can be recovered by a backtracking scheme. This is implemented by solving a gradient descent ordinary differential equation (ODE) backward in time from the target, reading for any \(t \in (0, T]\) \[\mathcal{G}'(t) = \mathbf{V}\left(\mathcal{G}(t)\right),\quad \mathbf{V}(\mathbf{x})= \partial_2\mathcal{H}(\mathbf{x},d\mathcal{U}_{\mathbf{p}}(\mathbf{x}))\] till the source point \(\mathbf{p}\) is reached, i.e. \(\mathcal{G}(T)=\mathbf{p}\) is achieved, where \(T:=\mathcal{U}_{\mathbf{p}}(\mathbf{x})\) is the arrival time and where \(\mathbf{V}\) is the geodesic flow. Finally, the target geodesic path \(\mathcal{G}_{\mathbf{p},\mathbf{x}}\) can be obtained by reparameterizing \(\mathcal{G}\) such that \(\mathcal{G}_{\mathbf{p},\mathbf{x}}(0)=\mathbf{p}\) and \(\mathcal{G}_{\mathbf{p},\mathbf{x}}(1)=\mathbf{x}\). Note that \(\partial_2\mathcal{H}(\mathbf{x},{\hat{\mathbf{x}}}):=\partial\mathcal{H}(\mathbf{x},{\hat{\mathbf{x}}})/\partial{\hat{\mathbf{x}}}\) denotes the partial differential of the Hamiltonian \(\mathcal{H}\) with respect to its second argument.

3 A HJB PDE Framework for Computing Curvature-penalized Geodesics with Tangent-Constrained Priors↩︎

In this section we present a unified minimal path framework that incorporates both path tangent constraints and curvature penalization. We begin by introducing the formulation of tangent constraints interpreted in the orientation-lifted space \(\mathbb{M}\), followed by the development of three distinct geodesic models that integrate these constraints with curvature penalization.

3.1 Constructing the Tangent-Constrained Priors↩︎

In image segmentation, a crucial challenge lies at the efficient representation of object shapes. Among various representations, the shape skeleton or interior landmark points are known to offer compact yet expressive descriptions. Hereafter, the skeleton or the interior landmark points are referred to as ISR. It encodes the internal geometric organization of a region, thus can be taken as an important cue to derive shape appearance features. In order to extend the effective scope to the whole domain \(\Omega\), a common way is to construct a distance map emanating from the skeleton structure [37]–[40]. In the remainder of this subsection, we discuss the computation for the geometric constraints using ISR.

3.1.1 Distances from the ISR↩︎

Let \(\Im\subset \Omega\) denote the ISR, from which one can compute a weighted distance map \(\mathcal{D}_{\Im}\) using the Eikonal PDE [23], [42] \[\label{eq95PriorHJB} \begin{cases} \tilde{\mathcal{H}}(\mathbf{x}, d\mathcal{D}_{\Im}(\mathbf{x}))=\frac{1}{2},&\forall \mathbf{x}\in\Omega\backslash\Im, \\ \mathcal{D}_{\Im}(\mathbf{x})=0,&\forall\mathbf{x}\in\Im \end{cases}\tag{7}\] where \(\tilde{\mathcal{H}}\) is a Hamiltonian. In particular for the Riemannian case, the Hamiltonian \(\tilde{\mathcal{H}}\) can be expressed as \[\tilde{\mathcal{H}}(\mathbf{x},{\hat{\mathbf{x}}})=\left\langle{\hat{\mathbf{x}}},\mathfrak M(\mathbf{x}){\hat{\mathbf{x}}}\right\rangle,\] where \(\mathfrak M\) is a tensor field such that \(\mathfrak M(\mathbf{x})\) is a symmetric definite positive matrix of size \(2\times2\), and where \(\langle\cdot,\cdot\rangle\) is the standard Euclidean scalar product in \(\mathbb{R}^2\). The tensor field \(\mathfrak M\) can be derived from the image data, for instance using image gradients [37], or from the shape of prescribed segmentation mask. In particular, the geodesic distance map \(\mathcal{D}_\Im\) can be reduced to a Euclidean distance map by setting \(\mathfrak M\equiv\mathop{\mathrm{Id}}\), where \(\mathop{\mathrm{Id}}\) is the identity of size \(2\times2\). The Eikonal PDE 7 can be efficiently solved using the fast marching method [42]–[44].

Let \(\tilde{\mathfrak} q:\Omega\to\mathbb{R}^2\) be a unit vector field subject to \(\|\tilde{\mathfrak} q\|\equiv1\), which can be naturally defined using the gradient of the distance map \(\mathcal{D}_{\Im}\) \[\tilde{\mathfrak} q(\mathbf{x}):= \begin{cases} \displaystyle \frac{\partial_2\tilde{\mathcal{H}}(\mathbf{x},d\mathcal{D}_{\Im}(\mathbf{x}))}{\|\partial_2\tilde{\mathcal{H}}(\mathbf{x},d\mathcal{D}_{\Im}(\mathbf{x}))\|},&\forall\mathbf{x}\in\Omega\backslash\Im,\\ \mathbf{0},&\text{otherwise} \end{cases}\] which leads to a rotated unit vector field \(\mathfrak q\) such that \[\label{eq95tangentPriors} \mathfrak q(\mathbf{x})=\mathbf{M}(\pi/2)\;\tilde{\mathfrak} q(\mathbf{x}),\tag{8}\] where \(\mathbf{M}(\pi/2)\) is a counter-clockwise rotation matrix with respect to an angle \(\pi/2\). Hereafter, we refer to the vector field \(\mathfrak q\) as the tangent priors.

3.1.2 Tangent-constrained priors in the orientation-lifted space \(\mathbb{M}\)↩︎

We introduce a first-order constraint to limit the path tangents. For a regular curve \(\gamma:[0,1]\to\Omega\), we assume that its tangent \(\gamma^\prime(t)\) satisfies the following acute alignment inequality \[\label{eq95TangentConstraint} \bigl\langle\gamma^\prime(t),\mathfrak q(\gamma(t)) \bigr\rangle \geq \rho(\gamma(t))\|\gamma^\prime(t)\|,\quad \forall t\in[0,1],\tag{9}\] where \(\rho:\Omega\to[0,1]\) is a scalar-valued function. The inequality in 9 imposes a mandatory constraint to the path tangent, i.e., the angle between the vectors \(\gamma^\prime(t)\) and \(\mathfrak q(\gamma(t))\) for any \(t\in[0,1]\) is required to be acute, since the values of \(\rho\) are limited to the range \([0,1]\).

The orientation lifting of a regular curve \(\gamma\) states that the tangent \(\gamma^\prime(t)\) can be characterized by its turning angles \(\eta\) via  1 . In this case, the inequality 9 can be rewritten as \[\label{eq95OLAcuteConstraint} \langle\mathop{\mathrm{\mathbf{n}}}(\eta(t)),\mathfrak q(\gamma(t)) \rangle\geq \rho(\gamma(t)).\tag{10}\] It allows to define a binary-valued function \(\mathfrak C:\mathbb{M}\to\{0,\infty\}\) encoding the inequality 10 \[\label{eq95AcuteConstraint} \mathfrak C_{\rho,\mathfrak q}(\mathbf{x})= \begin{cases} 0,&\langle\mathop{\mathrm{\mathbf{n}}}(\theta),\mathfrak q(x) \rangle\geq \rho(x)\\ \infty,&\text{otherwise} \end{cases}\tag{11}\] for any point \(\mathbf{x}=(x,\theta)\in\mathbb{M}\), which can be used to build a HJB PDE framework for computing minimal paths with both curvature penalization and tangent constraint.

3.2 Tangent-constrained Curvature-penalized Geodesic Models↩︎

In this section, we introduce the core of this work: a HJB PDE framework for computing minimal paths with both curvature penalization and tangent constraint, based on the acute alignment constraint defined in 9 10 . In conjunction with \(\mathfrak C_{\rho,\mathfrak q}\) formulated in 11 , we introduce a new type of geodesic metric \(\mathfrak F_{\rho,\mathfrak q}:\mathbb{M}\times\mathbb{R}^3\to[0,\infty]\), encoding the tangent constraint and reading for any point \(\mathbf{x}=(x,\theta)\) and for any vector \({\dot{\mathbf{x}}}=(\dot{x},{\dot{\theta}})\) as \[\begin{align} \tag{12} \mathfrak F_{\rho,\mathfrak q}(\mathbf{x},{\dot{\mathbf{x}}})&=\begin{cases} \left(\mathcal{C}\bigl(\xi {\dot{\theta}}/\|\dot{x}\|\bigr)+\mathfrak C_{\rho,\mathfrak q}(\mathbf{x})\right)\|\dot{x}\|, & \text{if~} \dot{x}= \mathop{\mathrm{\mathbf{n}}}(\theta)\|\dot{x}\|,\\ \infty,& \text{otherwise} \end{cases}\\ \tag{13} &=\begin{cases} \mathcal{F}(\mathbf{x},{\dot{\mathbf{x}}}),&\text{if~}\langle\mathop{\mathrm{\mathbf{n}}}(\theta),\mathfrak q(x) \rangle\geq \rho(x)\\ \infty,&\text{otherwise}, \end{cases} \end{align}\] where \(\mathcal{F}\) is a curvature-penalized geodesic metric dependent on the models, as described in Appendix 7. For any point \(\mathbf{x}=(x,\theta)\), the geodesic metric \(\mathfrak F_{\rho,\mathfrak q}(\mathbf{x},{\dot{\mathbf{x}}})\) is assigned the value \(\infty\) if the prescribed acute alignment constraint is not satisfied.

The computation for minimal paths under the introduced metric \(\mathfrak F_{\rho,\mathfrak q}\) can be performed within the HJB PDE framework. The corresponding Hamiltonian of the metric \(\mathfrak H_{\rho,\mathfrak q}:\mathbb{M}\times\mathbb{R}^3\to[0,\infty]\) can be constructed using \(\mathcal{H}\) of \(\mathcal{F}\), see 6 . This is to say \[\mathfrak H_{\rho,\mathfrak q}(\mathbf{x},{\hat{\mathbf{x}}}):=\mathcal{H}(\mathbf{x},{\hat{\mathbf{x}}}),\quad \text{if~}\langle\mathop{\mathrm{\mathbf{n}}}(\theta),\mathfrak q(x) \rangle\geq \rho(x),\] for any point \(\mathbf{x}=(x,\theta)\in\mathbb{M}\) and any co-vector \({\hat{\mathbf{x}}}=({\hat{x}},{\hat{\theta}})\).

Given a source point \(\mathbf{p}\) and a target point \(\mathbf{q}\) located in \(\mathbb{M}\), the minimal path \(\mathcal{G}_{\mathbf{p},\mathbf{q}}\) between them can be generated using a geodesic distance map \(\mathfrak U_{\mathbf{p}}\) by solving a HJB PDE associated with the Hamiltonian \(\mathfrak H_{\rho,\mathfrak q}\), i.e. the distance map \(\mathfrak U_{\mathbf{p}}\) satisfies that \(\mathfrak U_\mathbf{p}(\mathbf{p})=0\) and otherwise \[\mathfrak H_{\rho,\mathfrak q}(\mathbf{x},d\mathfrak U_\mathbf{p}(\mathbf{x}))=\frac{1}{2}\psi(\mathbf{x})^2,\quad \forall\mathbf{x}\in\mathbb{M}\backslash\{\mathbf{p}\},\] whose solution can be obtained by the Hamiltonian fast marching method [31].

Alternatively, we introduce here a practical way for computing the geodesic distance map \(\mathfrak U_\mathbf{p}\) by reconstructing its supporting domain using the inequality 10 : \[\mathbb{M}_{\mathfrak q,\rho}:=\bigl\{\mathbf{x}=(x,\theta)\in\mathbb{M}~|~\langle\mathop{\mathrm{\mathbf{n}}}(\theta),\mathfrak q(x) \rangle\geq \rho(x)\bigr\}.\] Then we estimate the viscosity solution to the following HJB PDE using the domain \(\mathbb{M}_{\mathfrak q,\rho}\) \[\label{eq95HJB95ours} \begin{cases} \mathcal{H}(\mathbf{x},d\mathfrak U_\mathbf{p}(\mathbf{x}))=\frac{1}{2}\psi(\mathbf{x})^2,&\quad \forall\mathbf{x}\in\mathbb{M}_{\mathfrak q,\rho}\backslash\{\mathbf{p}\},\\ \mathfrak U_\mathbf{p}(\mathbf{p})=0.& \end{cases}\tag{14}\]

In this work, we use the method introduced in [45] for computing the data-driven cost \(\psi\), which encodes both of the image gradients and the image region-based appearance features.

4 Applications to Robust Image Segmentation↩︎

In this section, we present a practical method that applies the proposed geodesic framework with tangent-constrained priors and curvature penalization to image segmentation. While the framework in 3 supports interactive segmentation where users manually specify the ISR, we develop here a fully automatic approach by leveraging deep learning-based pre-segmentation to extract the ISR.

The proposed segmentation method mainly consists of three subsequent stages: (1) generating the ISR from neural network pre-segmentation (4.1), (2) constructing the image data-driven cost \(\psi\) (4.2), and (3) computing the final segmentation contours via the adaptive circular geodesic models [45], [47], [48] (4.3). Fig. 2 illustrates the complete segmentation method.

4.1 Generation of the ISR↩︎

In our work, the ISR are first constructed to initialize the tangent-constrained geodesic model. These representatives are produced through a pre-segmentation procedure using a pre-trained neural network. While various architectures such as PolarMask [49] and other prevalent segmentation models [2] are suitable, we employ the DeepLabv3+ [50] model. Its encoder-decoder architecture equipped with atrous spatial pyramid pooling effectively captures multi-scale context and preserves boundary details, thereby ensuring the generation of topologically consistent initial masks for the subsequent extraction of ISR.

From the probability map output by the network, we obtain a binary segmentation mask by thresholding. We then extract the largest connected component to represent the object of interest. The centroid of this region, denoted as \(x_C\), is computed to serve as the reference center. Additionally, we extract the morphological skeleton \(\mathcal{S}\) of the binary mask, which provides a topological representation of the object’s shape and is used to define the tangent constraint field.

4.2 Data-driven Cost from Image Features↩︎

As discussed in 3, the HJB PDE 14 of the proposed model involves two components, the tangent-constrained priors and the image data-driven cost \(\psi\). Following the literature [45], the data-driven cost \(\psi\) in our work integrates both the image gradients-based and region-based appearance features. In particular, the region-based appearance features are dependent on the evolving shape [45], [51]. We denote by \(U\subset\Omega\) a tubular neighbourhood that surrounds the boundary of the evolving shape. For each point \(\mathbf{x}= (x, \theta) \in\mathbb{M}\), the data-driven cost function, defined in an orientation-scores manner, is formulated as \[\psi(x,\theta)= \begin{cases} \exp\left(\alpha\,\tilde{\psi}(x,\theta)\right),&\forall x\in U\\ \infty,&\text{otherwise}, \end{cases}\] where \(\alpha>0\) is a weighting parameter. The function \(\tilde{\psi}\) combines both types of feature as follows \[\label{eq95SubCosts} \tilde{\psi}(x,\theta)=\psi_{\rm grad}(x,\theta)+\frac{\mu\langle\omega(x),\mathop{\mathrm{\mathbf{n}}}(\theta)\rangle}{\max_{y\in U}\|\langle\omega(y),\mathop{\mathrm{\mathbf{n}}}(\theta) \rangle\|},\tag{15}\] where \(\psi_{\rm grad}\) is the orientation score from the image gradients and \(\omega\) is a vector field of sufficiently small magnitude, respectively formulated in Appendices 8 and 9. The first term in Eq. 15 captures the edge probability map through the structure tensor \(\mathcal{W}\), while the second term encodes regional preference via the vector field \(\omega\), with \(\mu>0\) balancing their relative contributions. Note that the vector field \(\omega\) and the tubular neighbourhood will be updated according to the evolving shapes.

4.3 Final Segmentation with Circular Geodesic Model↩︎

The circular geodesic model (CGM) is a practical method to extract closed planar geodesic curves within the image domain \(\Omega\), which provides a feasible geodesic-based solution to the task of image segmentation. In contrast to the original interactive context, where a point \(\mathbf{p}\) on the target boundary and its tangent \(\theta_p\) are user-defined, our implementation automates this initialization using the pre-segmentation results.

The seed point \(z\) required by the CGM is set to the centroid \(x_C\) provided by the pre-segmentation stage. The source point \(\mathbf{p}\) on the target boundary is automatically selected as the point on the contour of the pre-segmentation mask that is farthest from the centroid \(z\). This selection strategy typically places \(\mathbf{p}\) at an extremity of the object, facilitating the robust propagation of the geodesic front. The corresponding initial orientation \(\theta_p\) is estimated from the tangent vector of the contour at \(\mathbf{p}\).

A ray line connecting the centroid \(z\) and the source point \(\mathbf{p}\) is introduced as a cut to locally disconnect the domain, transforming the problem into finding a path from \(\mathbf{p}\) to itself in the cut domain. The skeleton \(\mathcal{S}\) is also utilized as a set of internal obstacles to ensure the generated curve properly encloses the object’s core structure. Minimal paths are then computed using the Hamiltonian fast marching method over the orientation-lifted domain.

5 Experimental Results↩︎

In this section, we evaluate the performance of the proposed geodesic models that unify tangent-constrained priors with curvature penalization. We first conduct a study to validate the effectiveness of the tangent-constrained priors when integrated with different curvature penalization metrics, including the RSF, Dubins, and Elastica models. Next, we demonstrate how these priors improve segmentation accuracy when using neural network pre-segmentation results as initialization. Finally, we compare the geodesic models against two state-of-the-art methods that also benefit from pre-segmentation: a graph-based model with convexity constraints (Graph-Convexity) and the Segment Anything Model (SAM).

5.1 Study on Tangent-Constrained Priors↩︎

In Fig. 3, we evaluate the influence of tangent constraint strictness on the minimal paths generated by the RSF, Dubins, and Elastica models. The strictness is governed by the scalar parameter \(\rho\) in the acute alignment condition (Eq. 9 ), where increasing \(\rho\) from \(0\) to \(1\) narrows the admissible angular sector centered on the tangent prior \(\mathfrak q\). Fig. 3 compares the baseline models (without tangent constraints, column (a)) against their tangent-constrained counterparts with maximum angular deviations of \(\pi/3\) and \(\pi/6\), as illustrated in columns (b) and (c). In the absence of constraints or when constraints are loose (column (b)), all models suffer from severe shortcutting, failing to track the hairpin contour and instead cutting directly across the gap between the start (green) and end (red) points. As the constraints become stricter, the paths are forced to align better with the boundary; notably, with a \(\pi/6\) constraint, all models successfully recover the true contour.

We further evaluate the impact of tangent-constrained priors on image segmentation using the proposed framework. Fig. 4 provides a qualitative comparison between standard curvature-penalized models (RSF, Dubins, Elastica) and their tangent-constrained counterparts (TC-RSF, TC-Dubins, TC-Elastica) on images from the MSRA-B dataset [52]. In the first two rows, the objects exhibit weak internal image gradients, resulting in low values of the data-driven cost \(\psi\) within the object interior. Consequently, the unconstrained geodesic models (columns (b), (d), (f)) tend to take shortcuts through the object’s body to minimize the total energy. In contrast, by incorporating the tangent constraint into the metric \(\mathfrak F_{\rho,\mathfrak q}\) (Eq. 12 ), the proposed models restrict the admissible path tangents to align with \(\mathfrak q\). This effectively prevents the geodesics from traversing the object interior, forcing them to follow the boundary even where the image edge evidence is weak (columns (c), (e), (g)). The third row illustrates a scenario where the object contains salient internal structures that act as local traps. Without tangent constraints, the geodesics are distracted by these internal boundaries, failing to encompass the entire object. However, the tangent priors successfully guide the optimization, ensuring the recovered paths respect the global topology of the object and avoid getting trapped in local minima.

5pt

5.2 Evaluation on Interactive Segmentation↩︎

We evaluate the proposed tangent-constrained geodesic models on the interactive image segmentation task. In this setting, the priors required for segmentation are provided through user interaction. Specifically, for our method, the user provides the skeleton (serving as the ISR to generate the tangent priors \(\mathfrak q\) via Eq. 8 ) along with the source point \(\mathbf{p}\) and its tangent direction \(\theta_p\) for the circular geodesic model, as described in Sec. 4.3.

For comparison, we select the graph-based Convexity Shape Prior (Graph-Convexity) method [53], a representative shape-aware interactive segmentation approach that enforces global convexity constraints through graph cuts optimization. This method requires user-provided scribbles indicating both foreground and background regions. To ensure a fair comparison, we set the foreground scribble identical to the skeleton used by our method, thereby sharing the same interior shape representation.

Fig. 5 presents a qualitative comparison between Graph-Convexity and the proposed tangent-constrained geodesic models (TC-RSF, TC-Dubins, TC-Elastica). The results reveal that Graph-Convexity, constrained by limited user input, frequently fails to capture the true object contour—either under-segmenting parts of the object, including extraneous background regions, or producing imprecise boundary delineation. In contrast, our geodesic models leverage the tangent-constrained metric \(\mathfrak F_{\rho,\mathfrak q}\) (Eq. 12 ) to restrict admissible path tangents, enabling accurate segmentation even with minimal user guidance.

We further conduct quantitative evaluation using the Dice similarity coefficient (DSC): \[\label{eq:dice} \mathrm{Dice}(S, GT) = \frac{2 \, |S \cap GT|}{|S| + |GT|},\tag{16}\] where \(S\) and \(GT\) denote the segmentation result and the ground truth, respectively. Experiments are performed on 20 images from the MSRA-B dataset [52], with results summarized in Tab. 1. The proposed tangent-constrained models achieve substantially higher mean Dice scores (TC-RSF: 0.9817, TC-Dubins: 0.9829, TC-Elastica: 0.9799) compared to both Graph-Convexity (0.7359) and the unconstrained geodesic models (RSF: 0.7790, Dubins: 0.6893, Elastica: 0.7157). Notably, the standard deviations of our models are an order of magnitude smaller, indicating consistently reliable performance across diverse object shapes. The unconstrained geodesic models exhibit high variance due to their susceptibility to shortcutting through low-gradient object interiors, as discussed in the previous subsection. These results demonstrate that incorporating the tangent-constrained priors effectively prevents such degenerate solutions, yielding robust segmentation with minimal user input.

5.3 Evaluation on Automated Segmentation↩︎

In this subsection, we evaluate the proposed framework in the automated segmentation setting, where all priors are derived from neural network pre-segmentation without user interaction. As described in Sec. 4.1, a DeepLabv3+ network provides the initial segmentation from which we extract the skeleton \(\mathcal{S}\) (serving as the ISR) to construct the tangent priors \(\mathfrak q\) and initialize the CGM.

We compare against two representative methods: the graph-based Graph-Convexity [53] and the Segment Anything Model (SAM) [54]. To ensure fair comparison, we adapt these originally interactive methods to utilize the same pre-segmentation. Specifically, for Graph-Convexity, the foreground scribbles are derived from the skeleton \(\mathcal{S}\). For SAM, positive prompt points are sampled from this skeleton. Since both baseline methods require background information (unlike our approach, which only requires the ISR), we generate background scribbles (for Graph-Convexity) and negative prompt points (for SAM) by randomly sampling from the region external to the pre-segmentation mask. All geodesic-based models share identical initialization parameters, including the source point \(\mathbf{p}\) and the cut configuration for the CGM.

Fig. 6 presents qualitative results on natural images from the MSRA-B dataset. The convexity constraint in Graph-Convexity proves overly restrictive, failing to capture non-convex boundary details when object shapes are complex (columns (a)–(b)). SAM, while capable of handling complex topologies, often struggles with precise boundary delineation, leading to under-segmentation or loss of fine structural details (columns (c)–(d)). In contrast, our TC-RSF model (columns (e)–(f)) leverages the tangent-constrained metric \(\mathfrak F_{\rho,\mathfrak q}\) (Eq. 12 ) combined with curvature regularization, enabling precise contour tracking even for objects with intricate shapes.

We further evaluate on medical images using the Left Atrium dataset from the Medical Segmentation Decathlon (MSD) [55]. As shown in Fig. 7, the low contrast and noise inherent in medical imaging pose significant challenges. Graph-Convexity struggles to localize boundaries accurately despite the target’s relatively convex shape, likely due to weak edge evidence in the data-driven cost \(\psi\). SAM tends toward over-segmentation, erroneously including adjacent tissues due to ambiguous intensity boundaries. Our TC-RSF model demonstrates superior robustness by leveraging the tangent priors \(\mathfrak q\) to guide the geodesic through low-contrast regions where image gradients alone are insufficient.

Table 2 presents the quantitative results on both datasets, computed using 100 selected representative images. The proposed tangent-constrained models consistently achieve the highest Dice scores: TC-Dubins leads on natural images (0.9739) while TC-RSF leads on medical images (0.9307). Compared to their unconstrained counterparts, the tangent-constrained models yield improvements of 5–9% (e.g., TC-RSF vs.RSF: 5.9% on MSRA-B, 6.2% on MSD LA), confirming that the tangent priors \(\mathfrak q\) effectively prevent shortcut paths. The standard deviations are also reduced by an order of magnitude (e.g., TC-Dubins: 0.0192 vs.Dubins: 0.2064), reflecting the stability provided by the hard constraint \(\mathfrak C_{\rho,\mathfrak q}\) (Eq. 11 ). Notably, all three tangent-constrained models outperform the DeepLabv3+ network that provides their initialization, demonstrating that the geometric constraints in the HJB PDE formulation (Eq. 14 ) effectively refine neural network predictions. In contrast, Graph-Convexity and SAM perform poorly, particularly on medical images (0.4258 and 0.7926, respectively), due to overly restrictive convexity assumptions and lack of domain-specific training.

6 Conclusion↩︎

In this work, we proposed a unified framework that incorporates tangent constraints as spatially varying shape constraints into curvature-regularized geodesic models within the HJB PDE formulation. By enforcing tangent admissibility in the orientation-lifted domain, our model restricts geodesic trajectories to prescribed angular sectors while remaining fully compatible with curvature penalization. The resulting tangent-constrained metric can be efficiently optimized using an adapted Hamiltonian Fast Marching algorithm, enabling the extraction of circular geodesics and closed contours for image segmentation. Experimental results on natural and medical images confirm that the proposed approach improves robustness under weak or noisy boundaries, mitigates topological shortcuts, and achieves competitive or superior segmentation performance compared to deep learning baselines and classical geodesic models without priors.

Future work will be devoted to extending the framework to 3D volumetric segmentation in conjunction with the framework [56], and developing end-to-end learning pipelines for joint extraction of tangent priors and the circular geodesic computation.

7 Curvature-penalized Geodesic Models↩︎

Existing second-order geodesic models [30]–[34] employ a curvature-dependent regularization term to construct the corresponding bending energy functional, see 2 . Typically, those models differ from each other mainly in the curvature cost \(\mathcal{C}\) used. Specifically, the cost \(\mathcal{C}\) for the Reeds-Shepp forward (RSF) model [30], the Dubins (D) model [31] and the Euler-Mumford elastica (EM) model [34] can be respectively formulated as \[\begin{align} &\mathcal{C}_{\rm RSF}(a)=\sqrt{1+(a)^2},\quad \mathcal{C}_{\rm D}(a)= \begin{cases} 1,&\text{if~}|a|\leq 1\\ \infty,&\text{otherwise} \end{cases},\\ &\mathcal{C}_{\rm EM}(a)=1+(a)^2, \end{align}\] for any scalar value \(a\in\mathbb{R}\).

In addition, the curvature prior elastica (CP) model [33] invokes the same curvature cost as the EM elastica model, except that the CP model utilizes a shifted curvature term \(\tilde{\kappa}\) defined such that \[\tilde{\kappa}(t)=\kappa(t)-\varsigma(\gamma(t),\eta(t)),\quad \forall t\in[0,1],\] where \(\kappa\) is the curvature of a regular curve \(\gamma\), and \(\varsigma:\mathbb{M}\to\mathbb{R}\) is a data-driven function which is also referred to as curvature priors. For a point \(x\) located in a centerline of tubular structure, \(\varsigma(x,\theta)\) is the estimated curvature of a fitting curve (to that centerline) whose path tangent is \(\mathop{\mathrm{\mathbf{n}}}(\theta)\). We refer the reader to [33] for the computation of curvature priors \(\varsigma\).

8 Orientation Scores from Image Gradients↩︎

Image edge and anisotropy features are extracted from image gradients. To enhance robustness to noise, the gradients are computed from a smoothed version of the image. First, the image \(I\) is smoothed by convolving it with a Gaussian kernel \(G_\sigma\) with standard deviation \(\sigma\), resulting in \(G_\sigma * I\). The core of this method is the Jacobian matrix, \(\mathcal{J}\), which contains the partial derivatives of the smoothed image. For a color image \(I = (I_1, I_2, I_3)\), the Jacobian at a point \(x\in\Omega\) is a matrix of size \(2 \times 3\): \[\mathcal{J}(x) = \begin{pmatrix} \partial_a G_\sigma * I_1 & \partial_a G_\sigma * I_2 &\partial_a G_\sigma * I_3 \\ \partial_b G_\sigma * I_1 & \partial_b G_\sigma * I_2 & \partial_b G_\sigma * I_3 \end{pmatrix},\] where \(\partial_a G_\sigma\) (resp. \(\partial_b G_\sigma\)) is the differential of the kernel \(G_\sigma\) along the axis \(a\) (resp. axis \(b\)).

For a grayscale image, the Jacobian simplifies to a \(2 \times 1\) vector of the smoothed image’s partial derivatives: \[\mathcal{J}(x) = (\partial_a G_\sigma * I, \partial_b G_\sigma * I)^\top(x).\] The anisotropic features, describing the orientation and shape of edges, are encapsulated in a tensor field \(\mathcal{W}\). At each point \(\mathbf{x}\), \(\mathcal{W}(\mathbf{x})\) is a \(2 \times 2\) symmetric positive-definite matrix \[\mathcal{W}(x) = \mathcal{J}(x)\mathcal{J}(x)^\top + \mathop{\mathrm{Id}}\] where \(\mathop{\mathrm{Id}}\) is the identity of size \(2\times2\), ensuring that the matrix \(\mathcal{W}(x)\) is positive-definite for any point \(x\in\Omega\).

The image-driven cost function \(\psi_{\rm grad}\) is formulated to guide the minimal path towards salient image edges. It incorporates the anisotropic feature from \(\mathcal{W}\): \[\label{eq95gradient95cost} \psi_{\rm grad}(\mathbf{x}) =\frac{\left\langle \mathop{\mathrm{\mathbf{n}}}(\theta)^{\perp}, \mathcal{W}(x) \mathop{\mathrm{\mathbf{n}}}(\theta)^{\perp} \right\rangle^{1/2}}{\max_{(y,\vartheta)\in\mathbb{M}}\left\langle \mathop{\mathrm{\mathbf{n}}}(\vartheta)^{\perp}, \mathcal{W}(x) \mathop{\mathrm{\mathbf{n}}}(\vartheta)^{\perp} \right\rangle^{1/2}}\tag{17}\] for any point \(\mathbf{x}=(x,\theta) \in \mathbb{M}\), where \(\mathop{\mathrm{\mathbf{n}}}(\theta) = (\cos\theta, \sin\theta)\) is a unit vector associated with angle \(\theta\), \(\mathop{\mathrm{\mathbf{n}}}(\theta)^{\perp} = (-\sin\theta, \cos\theta)\) is its perpendicular vector, and \(\langle \cdot, \cdot \rangle\) denotes the standard Euclidean scalar product. The term \(\left\langle \mathop{\mathrm{\mathbf{n}}}_{\theta}^{\perp}, \mathcal{W}(x) \mathop{\mathrm{\mathbf{n}}}_{\theta}^{\perp} \right\rangle\) serves as an orientation score. This score is maximized when the vector \(\mathop{\mathrm{\mathbf{n}}}_{\theta}^{\perp}\) is aligned with the image gradient direction which is the principal direction of \(\mathcal{W}(x)\). This, in turn, encourages the path’s tangent vector, \(\mathop{\mathrm{\mathbf{n}}}_{\theta}\), to align with the orientation of the edge, where the cost \(\psi_{\rm grad}(\mathbf{x})\) is minimal. This design is advantageous for segmentation, as it generates a valley of low potential along image edges, thereby allowing to attract the minimal paths to pass through the desired edge points.

9 Orientation Scores from Region-based Appearance Models↩︎

Beyond image gradient features, region-based image appearance features provide complementary cues for guiding geodesic paths for delineating the target boundary enclosing regions with complicated image content. We adopt the region competition model with Gaussian mixture model (GMM) as an example for quantifying the likelihood of image intensities belonging to foreground or background regions. Note that other region-based appearance models such as image color histograms-based models can be considered, and we refer the reader to [45], [51] for more detail.

Consider a simple closed curve \(\gamma:[0,1]\to\Omega\) that divides the domain \(\Omega\) into two disjoint regions: \(\mathcal{R}_1\) (interior) and \(\mathcal{R}_2\) (exterior) with \(\Omega=\mathcal{R}_1\cup \mathcal{R}_2\). For each region, we define a data fidelity function \(\rho_i:\Omega\to\mathbb{R}\) (for \(i=1,2\)) based on a Gaussian mixture model (GMM), where \(\rho_i(x)=-\log\,P_i(I(x)~|~\Theta_i)\) represents the negative log-likelihood of the image intensity \(I(x)\) under the probability density function \(P_i\) with the parameters \(\Theta_i\).

To incorporate these regional statistics into the geodesic framework, we utilize the method introduced in [51] to encode the difference \(\rho_1-\rho_2\) into a vector field \(\omega:\mathbb{R}^2\to\mathbb{R}^2\) of sufficiently small magnitude. This is implemented by addressing the following PDE-constrained minimization problem \[\min~\int_{\mathbb{R}^2}\|\omega(x)\|^2\,dx,~s.t.~\mathop{\mathrm{curl}}\omega=(\rho_1-\rho_2)\chi_U,\] where \(\chi_U\subset\Omega\) denotes the indicator function of a bounded tubular neighbourhood \(U \subset \Omega\). In other words, among all feasible solutions to the PDE \(\mathop{\mathrm{curl}}\omega=(\rho_1-\rho_2)\chi_U\), we choose the one minimizing the energy \(\int_{\mathbb{R}^2}\|\omega(x)\|^2dx\), which admits a closed-form solution via convolution.

References↩︎

---

1. cdi@qlu.edu.cn↩︎

2. liuli184791@sjtu.edu.cn↩︎

3. jinglin.zhang@sdu.edu.cn↩︎

4. zhenjiangli1987@163.com↩︎

5. chenda@ceremade.dauphine.fr↩︎

6. cohen@ceremade.dauphine.fr↩︎

7. Corresponding to their classical counterparts, the proposed models are referred to as the TC-RSF, TC-Dubins and TC-Elastica, respectively.↩︎

8. In Appendix 7, we show different types of curvature-penalized geodesic models and their corresponding curvature cost functions \(\mathcal{C}\).↩︎