# An adaptive geodesic voting method for curvilinear tree structure extraction

Hongda Liang, Da Chen, Tao Chen, Li Liu, Jiong Zhang, Laurent D. Cohen

► To cite this version:

Hongda Liang, Da Chen, Tao Chen, Li Liu, Jiong Zhang, et al. An adaptive geodesic voting method for curvilinear tree structure extraction. Sixteenth International Conference on Graphics and Image Processing (ICGIP 2024), Nov 2024, Nanjing, China. 10.1117/12.3057740. hal-05344786

HAL Id: hal-05344786

https://hal.science/hal-05344786v1

Submitted on 3 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# An Adaptive Geodesic Voting Method for Curvilinear Tree Structure Extraction

Hongda Liang¹, Da Chen¹, Tao Chen², Li Liu³, ⁴*, Jiong Zhang², Laurent D. Cohen⁵

¹Shandong Artificial Intelligence Institute, Qilu University of Technology (Shandong Academy of Sciences), Jinan 250014, China

²Laboratory of Advanced Theranostic Materials and Technology, Ningbo Institute of Materials Technology and Engineering, Chinese Academy of Sciences, Ningbo, China

³Department of Rehabilitation Medicine, Tong Ren Hospital, Shanghai Jiao Tong University School of Medicine, 200336, Shanghai, China.

⁴Institute of Rehabilitation, Shanghai Jiao Tong University School of Medicine, 200025, Shanghai, China.

⁵University Paris Dauphine, PSL Research University, CNRS, UMR 7534, CEREMADE, 75016 Paris, France.

*corresponding author: liuli19910101@163.com

## ABSTRACT

Geodesic voting method is known as a powerful tool for extracting curvilinear structures, which is able to find a tree structure from a single point. However, this method may fail to generate accurate results in complex scenarios such as complex network-like structures, due to the limitation of single source point. In order to solve this problem, we propose an adaptive curvature-penalized geodesic voting method where multiple source points with geometric voting constraint can be used for constructing the voting score map. In addition, we exploit the introduced adaptive geodesic voting method for the task of retinal vessel tracking, in conjunction with a deep learning-based junction points detection procedure. Experimental results on both synthetic images and retinal images prove the efficiency of the introduced adaptive geodesic voting method.

**Keywords:** Curvilinear Tree Structure Tracking, Geodesic voting, Junction detection

## 1. INTRODUCTION

The extraction of curvilinear structures is a fundamental task in the fields of image analysis. Significant applications include blood vessel tracking from medical images, roads detection from aerial images and find optimal trajectories in motion planning and robotic navigation. In the past decades, models on addressing the curvilinear structure extraction have been extensively studied. Despite great advances around this issue, accurately extracting complex curvilinear tree structures remain a challenging task, due to their elongated shapes [16].

Approaches based on the curvilinear detectors are a type of simple yet powerful methods [8], which are capable of transforming the images to process to an enhanced map. These models can take advantages of the geometry of curvilinear structures, thus have been widely exploited the curvilinear feature descriptors. Typical examples on curvilinear detectors include [6,12,18,25]. Models [3,24] driven based deep learning techniques refer to the curvilinear structure segmentation as a binary classification problem. Thus the strong ability of the neural network in features representation allows these models to extract satisfactory results even in some difficult tasks. However, it is difficult for both types of approaches mentioned above to benefit from the connectivity prior of curvilinear structures, thus very often leading to segmentation featuring gaps.
The optimal paths are taken as an efficient curvilinear descriptor, which can overcome that problem, since a curvilinear structure can be naturally modeled as a continuous curve. Significant examples include the active contour models [17,26] and the minimal paths models based on the Hamiltonian-Jacobi-Bellman (HJB) equation [20], relying on an energy minimization framework. Cohen and Kimmel [4] connect the computation of globally optimal paths to the HJB equation, giving various applications of minimal paths in image analysis. From then on, many efforts have been devoted to exploiting minimal paths to address curvilinear structures from various images [1,7,9]. However, most of these minimal paths-based models can be categorized as an interactive segmentation method, where the endpoints of a curvilinear structure should be prescribed manually. The geodesic voting method [22] uses a set of sampled points as the endpoint candidates, from which a large number of minimal paths can be obtained. The basic idea is that if most of these sampled points are inside or nearby the curvilinear structures, then the density of the respective minimal paths can be used to measure the appearance of curvilinear structures. The geodesic voting method can extract a complete curvilinear tree structure from a single source point. Unfortunately, it is difficult for this model to get satisfactory results in complicated scenarios, due to the limitations of single source point.

In this work, we propose an adaptive geodesic voting method which can use multiple source points to enhance the curvilinear structure extraction results. We introduce a geometric voting constraint in order to remove the voting scores from unexpected minimal paths. In addition, in contrast to the original model which only considers isotropic curvilinear features, the introduced model utilizes a second-oder model where the path curvature is taken into account. Fig.1.1 illustrates examples on the geodesic voting method on synthetic images. In columns 1 and 2, the tests are performed in the same image but with different numbers of source points. We can observe that the result in column 1 suffers from gaps while the result of the introduced adaptive voting method with multiple source points is satisfactory. In column 3, we perform the adaptive geodesic voting method in a noisy image, where one can also observe accurate extraction result.

The remaining of this paper is organized as follows. Section 2 presents the background on computing elastica geodesic paths and section 3 presents the core of this work: the adaptive geodesic voting method. The numerical experimental results and the conclusion are given in Sections 4 and 5.

<figure><img src="image_2.png" /><figcaption>Figure 1.1. Results of the geodesic voting method. The top row shows synthetic images with green dots indicating the planar points of the source points. The second row illustrates the extraction results, with single source point (column 1), multiple source points (column 2) and multiple source points in noisy image (column 3).</figcaption></figure>
## 2. BACKGROUND ON THE ELASTICA GEODESIC MODEL

The second-order elastica geodesic model [2,13,15] is designed to search for the global minimization of the classical Euler-Mumford elastica problem [19]. Let $\Omega \subset \mathbb{R}^2$ be an open and bounded space which we referred to as the image domain and let $\mathbb{S}^1 := \mathbb{R} \setminus 2\pi\mathbb{Z}$ denote the circle. The elastica geodesic model considers a curvature-regularized energy of a smooth curve $\gamma: [0,1] \to \Omega$, reading as

$$ Length(\gamma) := \int_{0}^{1} \alpha(\gamma(t), \eta(t)) (1 + (\xi\kappa(t))^2) dt, \qquad (2.1) $$

where $\alpha: \mathbb{M} := \Omega \times \mathbb{S}^1 \to \mathbb{R}^+$ is a data-driven cost, $\eta: [0,1] \to \mathbb{S}^1$ is the turning angles of $\gamma$ subject to $\forall t \in [0,1]$, $\dot{n}(\eta(t)) = (\cos(\eta(u)), \sin(\eta(u)))\gamma'(t) \parallel \gamma'(t) \parallel$, $\xi > 0$ is a weighting constant controlling the relative importance of the curvature $\kappa: [0,1] \to \mathbb{R}$ of $\gamma$. Using the turning angles $\eta$, the curvature $\kappa$ can be reformulated as $\kappa(t) = \dot{\theta}(t)/\|\gamma'(t)\|$.

The minimization of the length in Eq. (2.1) is carried out in the HamiltonianJacobi-Bellman (HJB) equation framework, via a degenerated metric $\mathcal{F}: \mathbb{M} \times \mathbb{R}^3 \to [0, \infty]$ for any orientation - lifted point $\mathbf{x} = (x, \theta) \in \mathbb{M}$ and any vector $\dot{\mathbf{x}} = (\dot{x}, \dot{\theta}) \in \mathbb{R}^3$

$$ \mathcal{F}(\mathbf{x}, \dot{\mathbf{x}}) = \begin{cases} \| \dot{\mathbf{x}} \| + \frac{(\xi \dot{\theta})^2}{\| \dot{\mathbf{x}} \|}, & \text{if } \dot{\mathbf{x}} = \dot{n}(\theta) \| \dot{\theta} \|, \\ \infty, & \text{otherwise,} \end{cases} $$

In the elastica geodesic model, given a fixed source point $\mathbf{s}$ and any target point $\mathbf{x}$, the minimization of Length($\gamma$) yields a minimal action map $\mathcal{U}$ such that $\mathcal{U}(\mathbf{x})$ denotes the minimal curve length from $\mathbf{s}$ to $\mathbf{x}$, i.e.

$$ \mathcal{U}(\mathbf{x}) = \inf\{\text{Length}(\gamma) \mid \gamma \in \text{Lip}([0,1], \mathbb{M}), \gamma(0) = \mathbf{s}, \gamma(1) = \mathbf{x}\}, $$

where $\text{Lip}([0,1] \to \mathbb{M})$ is the set of all Lipschitz continuous curves $\gamma: [0,1] \to \mathbb{M}$. Fortunately, the estimation of $\mathcal{U}$ can be efficiently obtained by solving the HJB equation such that $\mathcal{U}(\mathbf{s}) = 0$, and for any $\mathbf{x} \in \mathbb{M} \setminus \{\mathbf{s}\}$

$$ \mathcal{H}(\mathbf{x}, \nabla \mathcal{U}(\mathbf{x})) = \frac{1}{2} \alpha(\mathbf{x})^2, \qquad (2.2) $$

where $\nabla \mathcal{U}$ is the euclidean gradient of the map $\mathcal{U}$, and where $\mathcal{H}: \mathbb{M} \times \mathbb{R}^3 \to [0, \infty]$ is the Hamiltonian of the metric $\mathcal{F}$ with a closed form [13]

$$ \mathcal{H}(\mathbf{x}, \dot{\mathbf{x}}) = \frac{1}{8} \left( \langle \dot{\mathbf{x}}, (\dot{n}(\theta), 0) \rangle + \sqrt{\langle \dot{\mathbf{x}}, (\dot{n}(\theta), 0) \rangle^2 + \xi^{-2} \hat{\theta}^2} \right)^2, $$

for any point $\mathbf{x} = (x, \theta) \in \mathbb{M}$ and for co-vector $\dot{\mathbf{x}} = (\dot{x}, \dot{\theta}) \in \mathbb{R}^3$ where we recall that $\dot{n}(\theta) = (\cos \theta, \sin \theta)$ Then a minimal path $\mathcal{G}: [0, T] \to \mathbb{M}$, or a geodesic, where $T := \mathcal{U}(\mathbf{x})$, is generated as the solution to the following gradient descent ordinary differential equation (ODE) on the minimal action map $\mathcal{U}$

$$ \mathcal{G}'(t) = \partial_2 \mathcal{H}(\mathcal{G}(t), \nabla \mathcal{U}(\mathcal{G}(t))), \quad \text{s.t.} \quad \mathcal{G}(0) = \mathbf{x}, \mathcal{G}(T) = \mathbf{s}. \qquad (2.3) $$

The minimal path $\mathcal{G}$ can be easily re-parameterized so as to get a new minimal path $\mathcal{G}_{s,x} \in \text{Lip}([0,1], \mathbb{M})$.

Numerically, the HJB equation (2.2) can be efficiently addressed by the state-of-the-art Hamiltonian fast marching method (HFM) [13]. The HFM adopts a monotone-advanced front propagation scheme in a single-pass way, allowing to terminate the algorithm when all the target points are reached. In addition, we adopt the robust scheme proposed in [14] for numerically solving the ODE (2.3).
### 3. ADAPTIVE GEODESIC VOTING METHOD WITH MULTIPLE SOURCE POINTS

In this section, we introduce the core of this work, named the adaptive curvature-regularized geodesic voting method with multiple source points.

#### 3.1 Constructing the Sets of Source and Target Points

In this section, we introduce an adaptive geodesic voting method allowing multiple source points to set up the model. We start by defining the set of orientation-lifted points which are taken as the source points of the geodesic voting method. We denote by $\mathcal{S} := \{\tilde{s}_1, \cdots, \tilde{s}_K\} \cup \{\tilde{s}_1, \cdots, \tilde{s}_K\}$ the set of the given source points for the introduced adaptive geodesic voting method. Each orientation-lifted source point $s_i$ (resp. $\tilde{s}_i$) consists of a planar point $s_i \in \Omega$ and an angle $\vartheta_i \in \mathbb{S}^1$ (resp. $\vartheta_i + \pi$), i.e. $s_i = (s_i, \vartheta_i)$ (resp. $\tilde{s}_i = (s_i, \vartheta_i + \pi)$). The planar points $s_i$ can be prescribed manually or by an application-specified algorithm such as vessel junction detection method. The angles $\vartheta_i$ associated to each planar point $s_i$ points to the tangent directions of the curvilinear structure at $s_i$. Note that the model allows that a planar point $s_i$ corresponds to more than two angles, due to the potential junction structures. Basically, the angles associated to $s_i$ are regarded as the local minimizers of $\alpha(s_i, \cdot)$. We say that $(s_i, \vartheta_k) \in \mathcal{S}$ and $(s_i, \vartheta_k + \pi) \in \mathcal{S}$ if $\vartheta_k$ satisfies $\alpha(s_i, \vartheta_k - \epsilon) < \alpha(s_i, \vartheta_k) \ge \alpha(s_i, \vartheta_k + \epsilon)$, i.e $\alpha(s_i, \vartheta_k)$ is a peak in a small local neighbourhood $[\vartheta_k - \epsilon, \vartheta_k + \epsilon]$.

Furthermore, let $\wp := \{\mathfrak{p}_1, \mathfrak{p}_2, \cdots, \mathfrak{p}_N\} \cup \{\tilde{\mathfrak{p}}_1, \tilde{\mathfrak{p}}_2, \cdots, \tilde{\mathfrak{p}}_N\}$ be a set of sampled points, which provide the target points for the introduced adaptive geodesic voting method. Again, the point $\mathfrak{p}_i$ has the same planar point $p_i$ with the point $\tilde{\mathfrak{p}}_i$, but opposite directions $\dot{n}(\theta_i)$ and $-\dot{n}(\theta_i)$. The angles $\theta_i$ corresponding to the planar points $p_i$ are estimated via $\theta_i := \arg \min_{\theta \in [0, 2\pi]} \alpha(p_i, \theta)$. The method for calculating all of the planar points $\{p_i\}_{1 \le i \le N}$ is as follows: First, the Hessian matrix is computed for the 2D image, defined as:

$$H = \begin{pmatrix} I_{xx} & I_{xy} \\ I_{xy} & I_{yy} \end{pmatrix}$$

where $I_{xx}, I_{xy}$, and $I_{yy}$ are the second derivatives of the image in the x and y directions. The Hessian matrix is then decomposed to obtain eigenvalues $\lambda_1$ and $\lambda_2$ (assuming $|\lambda_1| \le |\lambda_2|$), and the vesselness measure is calculated based on these eigenvalues:

$$\text{Vesselness} = \begin{cases} 0, & \text{if } \lambda_2 > 0, \\ \exp\left(-\frac{R_B^2}{2\beta^2}\right)\left(1 - \exp\left(-\frac{S^2}{2c^2}\right)\right), & \text{otherwise.} \end{cases}$$

where $R_B = \frac{|\lambda_1|}{|\lambda_2|}$ and $S = \sqrt{\lambda_1^2 + \lambda_2^2}$. By setting different scale ranges and ratios, multi-scale analysis is performed on the image, and the maximum response is extracted to form a single scale vesselness image. To determine the centerlines of vessel structures, the vesselness image is sorted, and a threshold is determined based on a preset threshold ratio. The vesselness image is then binarized. Next, the skeleton extraction algorithm [21][28] is applied to extract the centerlines of the vessels. Finally, the points on the skeleton are taken as the planar points $\{p_i\}_{1 \le i \le N}$.

#### 3.2 Adaptive Geodesic Voting Method with Multiple Source Points

The classical geodesic voting method [22] investigates only one single source point for computing the voting paths, very often leading to shortcut problems. We solve this issue by introducing a new curvature-penalized geodesic voting method featuring a multiple source points mechanism with a geometric voting constraint. Algorithm 1 presents the main steps of the introduced adaptive geodesic voting method.

The geodesic voting score map is regarded as a core ingredient in the geodesic voting method [22]. In the basic definition of the set $\mathcal{S}$, each planar point $s_i$ may correspond to multiple orientation-lifted points included in the set $\mathcal{S}$. For this purpose, we split $\mathcal{S}$ into a family of subsets in terms of planar points

$$\mathcal{S}_i := \{s = (s, \vartheta_k) \in \mathcal{S}, k \ge 2\}, \quad (3.1)$$
where each set $\mathfrak{F}_i$ contains all the orientation-lifted points $s$ such that they have the identical planar point $s$ but different angles. The final geodesic voting map, denoted by $\Psi: \Omega \to \mathbb{Z}_0^+$, can be computed as the summation of all $\psi_i$ for $1 \le i \le K$:

$$ \Psi(x) = \sum_i \psi_i(x), \quad \forall x \in \Omega. \qquad (3.2) $$

In our work, we suppose that each set $\mathfrak{F}_i$ leads to a voting score map $\psi_i: \Omega \to \mathbb{Z}^+ \cup \{0\}$. At each planar point $p_j \in \Omega$ which is the planar point involved in $\wp$, we expect to track only one minimal path $G_{i,j} = (\gamma_{i,j}, \eta_{i,j})$ with $G_{i,j}(0) \in \mathfrak{F}_i$ and $G_{i,j}(0) \in \{p_j, \tilde{p}_j\}$ to contribute to the voting score map $\psi_j$. This is done by firstly addressing the following minimal curve length problem

$$ \inf\{\text{Length}_{\mathcal{F}}(\Gamma) \mid \Gamma \in \text{Lip}([0,1], \mathbb{M}), \Gamma(0) \in \mathfrak{F}_i, \Gamma(1) \in \{p_j, \tilde{p}_j\}\} \qquad (3.3) $$

where $\text{Length}_{\mathcal{F}}(\Gamma)$ is the weighted curve length of $\Gamma$ associated with the metric $\mathcal{F}$ and cost function $\alpha$:

$$ \text{Length}_{\mathcal{F}}(\Gamma) = \int_{0}^{1} \alpha(\Gamma(t)) \mathcal{F}(\Gamma(t), \Gamma'(t)) dt. \qquad (3.4) $$

Using the HJB equation framework as discussed in Section 2, the minimization of (3.3) can be carried out by the minimal action map $\mathcal{U}_i$ which is the viscosity solution to the HJB equation

$$ \begin{cases} \mathcal{H}(\mathbf{x}, \nabla \mathcal{U}_i(\mathbf{x})) = \frac{1}{2} \alpha(\mathbf{x})^2, & \forall \mathbf{x} \in \mathfrak{F}_i \\ \mathcal{U}_i(\mathbf{x}) = 0, & \forall \mathbf{x} \in \mathfrak{F}_i. \end{cases} \qquad (3.5) $$

Specifically, once the minimal action map $\mathcal{U}_i$ is estimated, from the target points $\{p_j, \tilde{p}_j\}$ we can obtain two minimal paths by reparameterizing the solution to the gradient descent ODE (2.3) on the minimal action map $\mathcal{U}_i$. Among them, the target minimal path $G_{i,j} = (\gamma_{i,j}, \eta_{i,j})$ is the one with minimal value of $\mathcal{U}_i$. In this case, by the minimal action map $\mathcal{U}_i$, one can compute $K$ minimal paths minimal path $\{G_{i,j}\}_{1 \le j \le K}$, where $K$ is the half number of the points in $\wp$.

In principle, if the planar point $p_j$ is inside a curvilinear structure, the path $G_{i,j}$ delineates the curvilinear structure whose endpoints are $\gamma_{i,j}(0)$ and $\gamma_{i,j}(1)$ associated with tangent directions $\dot{n}(\eta_{i,j}(0))$ and $\dot{n}(\eta_{i,j}(1))$, respectively. However, for the case that the curvilinear structures are very complex and the two points $p_j$ and $s_i$ are too far to each other, the planar curve $\gamma_{i,j}$ very often suffers from the shortcut problem. In order to reduce the possibility of such an issue, we introduce a geometric voting constraint. Specifically, we suppose that only a minimal path $G_{i,j}$ whose target point $G_{i,j}(1) = p_j$ will contribute to the voting score map $\psi_j$ if the distance between $p_j$ and $s_i$ is sufficiently low. In this work, this distance is measured using the curvature-weighted length of $G_{i,j}$, that is defined as

$$ \mathcal{L}(G_{i,j}) = \int_{0}^{1} \sqrt{\gamma'_{i,j}(t)^2 + \beta \eta'_{i,j}(t)^2} dt = \int_{0}^{1} \sqrt{1 + \beta \kappa_{i,j}(t)^2} \| \gamma'_{i,j}(t) \| dt \qquad (3.6) $$

where $\kappa_{i,j}: [0,1] \to \mathbb{R}$ is the curvature of the curve $\gamma_{i,j}$ and $\beta > 0$ is a constant.

Eventually, the voting score map $\psi_i$ can be estimated by

$$ \psi_i(x) = \sum_{j=1}^{N} \delta_x(\gamma_{i,j}) \chi_{\mathcal{L}(G_{i,j})>\lambda} \qquad (3.7) $$

where $N$ is the total number of points included in the set $\wp$ of target points, $\delta_x(\gamma_{i,j})$ is the voting detection function [22] such that $\delta_x(\gamma_{i,j}) = 1$ if $\exists t \in [0,1], \gamma_{i,j}(t) = x$ and $\delta_x(\gamma_{i,j}) = 0$ otherwise, and $\chi_{\mathcal{L}}(G_{i,j}) > \lambda$ is defined as

$$ \chi_{\mathcal{L}(G_{i,j})>\lambda} := \begin{cases} 1, & \text{if } \mathcal{L}(G_{i,j}) > \lambda, \\ 0, & \text{otherwise,} \end{cases} $$
with $\lambda > 0$ being a thresholding value. By Eq. (3.2), the target voting score map $\Psi$ can be computed, where a high value $\Psi(x)$ indicates a high possibility that the planar point $x \in \Omega$ is located in the centerline of a curvilinear structure. Therefore, after a thresholding procedure to the voting score map $\Psi$, we can obtain the tracking result of a complex curvilinear three structure.

<caption>**Algorithm 1:** Adaptive Geodesic Voting with Multiple Source Points</caption>

**Input:** The metric $\mathcal{F}$, the sets $\mathcal{S}, \wp$ and a constant $\lambda$.
**Output:** The voting score map $\Psi$
**Initialization:**
  Split the set $\mathcal{S}$ into $M$ subsets $\mathfrak{I}_i$ with $1 \le i \le M$.
  Set $i \leftarrow 1$ and $\Psi(x) \leftarrow 0, \forall x \in \Omega$.
**for** $i \le M$ **do**
  Compute the minimal action map $U_i$ using $\mathfrak{I}_i$ as the source point set.
  Set $j \leftarrow 1$ and $\psi_i(x) \leftarrow 0, \forall x \in \Omega$.
**for** $j \le K$ **do**
  Compute the minimal path $G_{i,j}$.
**if** $\mathcal{L}(G_{i,j}) > \lambda$ **then**
  Compute the voting score map $\psi_i$ using Eq. (3.7).
**end**
Set $j \leftarrow j + 1$.
**end**
Set $\Psi(x) \leftarrow \Psi(x) + \psi_i(x), \forall x \in \Omega$ and set $i \leftarrow i + 1$.
**end**

## 3.3 Numerical Implementation

### 3.3.1 The Computation of the Minimal Action Maps

The numerical implementation of the introduced adaptive geodesic voting method are established over a regular grid $M_h = M \cap (h\mathbb{Z}^2 \times h\mathbb{Z}/2\pi\mathbb{Z})$ with grid scale $h = 2\pi/N_\theta$, where $N_\theta$ is the discrete angles along the orientation dimension and where $h\mathbb{Z}/2\pi\mathbb{Z} := \{0, h, 2h, \dots, (N_\theta - 1)h\}$. In addition, the discretized image domain $\Omega_h$ is denoted as $\Omega_h = \Omega \cap h\mathbb{Z}^2$.

The computation of the minimal action map $U_i$ is the core in the proposed model, implemented using the HFM method [13]. The HFM method invokes the following discretized system of the HJB equation (3.5)

$$ \sum_{1 \le l \le L} \omega_l(\mathbf{x}) \sum_{1 \le j \le J} \rho_{jl}(\mathbf{x}) \left( \frac{u_i(\mathbf{x}) - u_i(\mathbf{x} - h\dot{\mathbf{e}}_{jl}(\mathbf{x}))}{h} \right)_+^2 = 1, \quad (3.8) $$

where $\omega_l(\mathbf{x}) > 0, \rho_{jl}(\mathbf{x}) > 0$ are positive weights, $\dot{\mathbf{e}}_j l(\mathbf{x}) \in \mathbb{Z}^3$ are offsets with integer coordinates, $(a)_+ = \max\{0, a\}$ denotes the positive part of a scalar value $a$. The weights $\omega_l(\mathbf{x})$ are computed using the Fejer rule. The weights $\rho_{jl}(\mathbf{x})$ and the offsets $\dot{\mathbf{e}}_j l(\mathbf{x})$ are estimated using Selling's decomposition of positive quadratic forms. For more details we refer to [13]. The offsets $\dot{\mathbf{e}}_j l(\mathbf{x})$ forms a neighbourhood system $\mathcal{S}$ adaptive to the metric $\mathcal{F}$. Specifically, one has $\mathcal{S}(\mathbf{x}) = \{\dot{\mathbf{e}}_{jl}(\mathbf{x}) \in \mathbb{Z}^3, 1 \le j \le J, 1 \le l \le L\}$ such that $(\mathbf{x} - h\dot{\mathbf{e}}_j l) \in M_h$ is a neighbour point of $\mathbf{x}$. The course of the HFM can be seen in Algorithm 2.

The computation of the curvature-weighted length (3.6) can be estimated once the paths $G_{i,j}$ are tracked. Also, it is can be obtained in an accumulation way during the computation of the map $U_i$, using the upwind gradient of $U_i$, see [10].

### 3.3.2 The Computation of the Cost Function $\alpha$

The cost function $\alpha$ characterizes the appearance features of curvilinear structures, such that $\alpha(x, \theta)$ should have a high value if the planar point $x$ is passed by a curvilinear structure and $\dot{n}(\theta)$ is collinear to the its tangent direction. In
our work, the cost $\alpha$ is computed using Frangi filter [5]. Let $G_r$ be a Gaussian kernel of variance $r^2$ and let $H_r(x)$ be the Hessian matrix of the Gaussian kernel $G_r$. The output of the Frangi filter can be expressed as $(H_r * I)(x)$, where $I: \Omega \to \mathbb{R}$ is the image to process, provided that the gray levels of curvilinear structures are locally darker than background. At a planar point $x$ belonging to a curvilinear structure, the matrix $Q(x, r) = (H_r * I)(x)$ have two eigenvalues $\lambda_1(x, r)$, $\lambda_2(x, r)$ and two corresponding eigenvectors $v_1(x, r)$, $v_2(x, r)$, subject to $\lambda_1(x, r) \le \lambda_2(x, r)$. The analysis $Q(x, r)$ yields a multi-scale vesselness map $f(x, r) = g(\lambda_1(x, r), \lambda_2(x, r))$ where the function $g$ is defined as [5]

$$g(a, b) = \begin{cases} 0, & \text{if } b < 0 \\ \exp(-A(a, b))(1 - \exp(-B(a, b))), & \text{otherwise.} \end{cases}$$

where $A(a, b) = a^2/(2u^2b^2)$ and $B(a, b) = (a^2 + b^2)/(2v^2)$ with $u, v$ being two constants. Let $\tau(x) = \arg\max_r f(x, r)$ be an optimal scale map. In this case, $v_1(x, \tau(x))$ indicates direction that a curvilinear structure should have at $x$. Then an orientation scores used in this work is defined as $\mathcal{P}(x, \theta) = f(x, \tau(x))|\langle \dot{n}(v_1(x, \tau(x))), \dot{n}(\theta)\rangle|$. As a result, the cost function $\alpha$ is expressed as

$$\alpha(x, \theta) = \exp(-\zeta \mathcal{P}(x, \theta)),$$

where $\zeta > 0$ is a constant.

<table><caption>**Algorithm 2:** Computing the Map $\mathcal{U}_i$ via the HFM Method</caption><tr><td>Input:</td><td>a metric $\mathcal{F}$ and a source point $s \in \Omega_h$.</td></tr><tr><td>Output:</td><td>Minimal action map $\mathcal{U}_i$.</td></tr><tr><td>Initialization:</td><td>Set $\mathcal{U}_i(\mathbf{x}) \leftarrow 0, \forall \mathbf{x} \in \mathfrak{F}_i$ and $\mathcal{U}_i(\mathbf{x}) \leftarrow \infty$, else.<br>Tag each grid point $\mathbf{x} \in \mathcal{M}_h$ as TRIAL.<br>Construct the set $\mathcal{T}$ involving all TRIAL points.<br>Construct the stencils $\mathcal{S}$ associated to the metric $\mathcal{F}$.</td></tr><tr><td>while $\mathcal{T} \neq \emptyset$ do</td><td>Set $\mathbf{x}_* \leftarrow \arg\min_{\mathbf{x} \in \mathcal{T}} \mathcal{U}_i(\mathbf{x})$.<br>Tag $\mathbf{x}_*$ as Accepted and set $\mathcal{T} \leftarrow \mathcal{T} \setminus \{\mathbf{x}_*\}$.<br>for each TRIAL point $y$ such that $\mathbf{x}_* \in \mathcal{S}(y)$ do<br>Set $\mathcal{U}_i(y)$ as the solution to the discretized HJB equation (3.8).</td></tr><tr><td>end</td><td></td></tr><tr><td>end</td><td></td></tr></table>

# 4. NUMERICAL EXPERIMENTS

In this section, we conducted experiments using the adaptive geodesic voting method with multiple source points on both synthetic image, highway image, satellite image and fundus image so as to verify the availability of the introduced method.

## 4.1 Application to Retinal Vessel Network Extraction

The retinal vessels appear as a network structure consisting of complex vessel branches such as bifurcations and crossings. In order to generate accurate extraction results, we first exploit a procedure to detect the junction points which connect more than one vessel segments. These junction points are then taken as the planar points to identify the set $\mathfrak{S}$ of source points.

In the past years, significant advances have been made in junction detection tasks [23,27] based on the deep learning techniques. Among them, we choose the Attention O-Net model [27] as the detector for junction points of retinal vessels, since this model is free to image segmentation. Basically, it involves two branches, where the first one is called junction detection branch (JDB) that utilizes heatmaps as labels, and the second one is a local enhancement branch (LEB) that
considers adaptive radius labels. For the JDB, the corresponding heatmap is created using a Gaussian kernel. The LEB applies a distance transform [11] to the ground truth vessel maps to ascertain the approximate radius of each junction point. Subsequently, the label is formulated based on the positions of the junctions and their respective radii. In this section, we train the junction detection model using the training set of the DRIVE dataset [27]. We process the images to be analyzed through the network, obtaining a set of planar positions {$s_i$}_i that locate the junctions. This set is then taken as the planar points to generate the set $\mathcal{S}$.

## 4.2 Application to Satellite Road Extraction

To validate the generality of our method, we captured several patches from Google Earth, including highways synthesized by Google Earth and actual satellite images of roads. We manually provide some junctions through an interactive program, obtaining a set of planar positions {$s_i$}_i. This set is then taken as the planar points to generate the set $\mathcal{S}$.

## 4.3 Results on Curvilinear Structure Extraction

In Fig. 4.1 a, the green dot is the planar position of the orientation-lifted source point. The red arrows characterize the angles associated to the planar position, using the method presented in Section 3.1. Fig. 4.1 b illustrates the points by red dots sampled using the farthest point sampling scheme. They are then used for constructing the target point set $\rho$.

<figure><img src="image_3.png" /><figcaption>Figure 4.1. a: The red arrows indicate the directions associated to the planar position (green dot). b: The red dots are the target points.</figcaption></figure>

Fig. 4.2 illustrates the extraction results on a retinal image patch, respectively: with single source point (column 2), multiple source points without geometric voting constraint (column 3) and multiple source points with geometric voting constraint (column 4).

Fig. 4.3 illustrates the extraction results on a highway image patch, respectively: with single source point (column 2), multiple source points without geometric voting constraint (column 3) and multiple source points with geometric voting constraint (column 4).

Fig. 4.4 illustrates the extraction results on a satellite image patch, respectively: with single source point (column 2), multiple source points without geometric voting constraint (column 3) and multiple source points with geometric voting constraint (column 4).

By testing our method in different scenarios, one can observe that the adaptive geodesic voting method indeed obtain promising results.
<figure><img src="image_4.png" /><figcaption>Figure 4.2 Vessel tree structure extraction from retinal image patches. Row1, column1: Retinal vessel image patch. Row2, column1: The green dots indicates the detected vessel junction points. Other columns in the top row illustrates the voting score maps Ψ, with single source point (column 2), multiple source points without geometric voting constraint (column 3) and multiple source points with geometric voting constraint (column 4). Other columns in the second row illustrates the extraction results, with single source point (column 2), multiple source points without geometric voting constraint (column 3) and multiple source points with geometric voting constraint (column 4).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4.3 Curvilinear Structure Extraction from highway image patch. Row1, column1: Highway image patch. Row2, column1: The green dots indicates the junction points. Other columns in the top row illustrates the voting score maps Ψ, with single source point (column 2), multiple source points without geometric voting constraint (column 3) and multiple source points with geometric voting constraint (column 4). Other columns in the second row illustrates the extraction results, with single source point (column 2), multiple source points without geometric voting constraint (column 3) and multiple source points with geometric voting constraint (column 4).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4.4 Curvilinear Structure Extraction from satellite image patch. Row1, column1: Satellite image patch. Row2, column1: The green dots indicate the junction points. Other columns in the top row illustrate the voting score maps Ψ, with single source point (column 2), multiple source points without geometric voting constraint (column 3) and multiple source points with geometric voting constraint (column 4). Other columns in the second row illustrates the extraction results, with single source point (column 2), multiple source points without geometric voting constraint (column 3) and multiple source points with geometric voting constraint (column 4).</figcaption></figure>

## 5. CONCLUSION

In this work, we introduce an adaptive geodesic voting model, encoding a geometric voting constraint, and allowing to incorporate multiple source points into the computation of voting score map that indicates the possibility of each point belonging to a curvilinear structure. In addition, we show that the integration of the adaptive geodesic voting method and a junction points detector indeed can extract promising results from images with complex curvilinear structures. Future work will be devoted to designing more types of geometric voting constraints in conjunction with deep learning techniques.

## REFERENCES

[1] F. Benmansour and L. D. Cohen. Tubular structure segmentation based on minimal path method and anisotropic enhancement. Int. J. Comput. Vis., 92(2):192–210, (2011)

[2] D. Chen, J.-M. Mirebeau, and L. D. Cohen. Global minimum for a Finsler elastica minimal path approach. Int. J. Comput. Vis., 122(3):458–483, (2017)

[3] Mingfei Cheng, Kaili Zhao, Xuhong Guo, Yajing Xu, and Jun Guo. Joint topology-preserving and feature-refinement network for curvilinear structure segmentation. In Proc. ICCV, pages 7147–7156, (2021)

[4] L. D. Cohen and R. Kimmel. Global minimum for active contour models: A minimal path approach. Int. J. Comput. Vis., 24(1):57–78, (1997)

[5] A. F. Frangi, W. J. Niessen, K. L. Vincken, and M. A. Viergever. Multiscale vessel enhancement filtering. In Proc. MICCAI, pages 130–137, (1998)

[6] J. Hannink, R. Duits, and E. Bekkers. Crossing-preserving multi-scale vesselness. In Proc. MICCAI, pages 603–610, (2014)

[7] V. Kaul, A. Yezzi, and Y. Tsai. Detecting curves with unknown endpoints and arbitrary topology using minimal paths. IEEE Trans. Pattern Anal. Mach. Intell., 34(10):1952–1965, (2012)

[8] Jonas Lamy, Odyssée Merveille, Bertrand Kerautret, and Nicolas Passat. A benchmark framework for multiregion analysis of vesselness filters. IEEE Trans. Med. Imag., 41(12):3649–3662, (2022)

[9] H. Li and A. Yezzi. Vessels as 4-D curves: Global minimal 4-D paths to extract 3-D tubular surfaces and centerlines. IEEE Trans. Med. Imag., 26(9):1213–1223, (2007)
[10] Li Liu, Da Chen, Minglei Shu, Baosheng Li, Huazhong Shu, Michel Paques, and Laurent D Cohen. Trajectory grouping with curvature regularization for tubular structure tracking. IEEE Trans. Image Processing, 31:405–418, (2021)

[11] Min Liu, Weixun Chen, Chao Wang, and Hanchuan Peng. A multiscale ray-shooting model for termination detection of tree-like structures in biomedical images. IEEE Trans. Med. Imag., 38(8):1923–1934, (2019)

[12] O. Merveille, H. Talbot, L. Najman, and N. Passat. Curvilinear structure analysis by ranking the orientation responses of path operators. IEEE Trans. Pattern Anal. Mach. Intell., 40(2):304–317, (2018)

[13] J.-M. Mirebeau. Fast-marching methods for curvature penalized shortest paths. J. Math. Imag. Vis., 60(6):784–815, (2018)

[14] J.-M. Mirebeau. Riemannian fast-marching on Cartesian grids, using Voronoi's first reduction of quadratic forms. SIAM J. Numer. Anal., 57(6):2608–2655, (2019)

[15] Jean-Marie Mirebeau and Jorg Portegies. Hamiltonian fast marching: A numerical solver for anisotropic and non-holonomic eikonal pdes. Image Processing On Line, 9:47–93, (2019)

[16] S. Moccia, E. De Momi, S. El Hadji, and L. S. Mattos. Blood vessel segmentation algorithms—Review of methods, datasets and evaluation metrics. Comput. Methods Programs Biomed., 158:71–91, (2018)

[17] V. Mohan, G. Sundaramoorthi, and A. Tannenbaum. Tubular surface segmentation for extracting anatomical structures from medical imagery. IEEE Trans. Med. Imag., 29(12):1945–1958, (2010)

[18] Stefano Moriconi, Maria A Zuluaga, H Rolf Jäger, Parashkev Nachev, Sébastien Ourselin, and M Jorge Cardoso. Inference of cerebrovascular topology with geodesic minimum spanning trees. IEEE Trans. Med. Imag., 38(1):225–239, (2018)

[19] David Mumford. Elastic and computer vision. in Algebraic Geometry and Its Applications (C. L. Bajaj, Ed.), Springer-Verlag, New York, pages 491–506, (1994)

[20] G. Peyré, M. Péchaud, R. Keriven, and L. D. Cohen. Geodesic methods in computer vision and graphics. Foundations and Trends in Computer Graphics and Vision, 5(3–4):197–397, (2010)

[21] Ta-Chih Lee, Rangasami L. Kashyap and Chong-Nam Chu "Building skeleton models via 3-D medial surface/axis thinning algorithms." Computer Vision, Graphics, and Image Processing, 56(6):462–478, (1994)

[22] Y. Rouchdy and L. D. Cohen. Geodesic voting for the automatic extraction of tree structures. Methods and applications. Comput. Vis. Image Underst., 117(10):1453–1467, (2013)

[23] Yinghui Tan, Min Liu, Weixun Chen, Xueping Wang, Hanchuan Peng, and Yaonan Wang. DeepBranch: Deep neural networks for branch point detection in biomedical images. IEEE Trans. Med. Imag., 39(4):1195–1205, (2019)

[24] Giles Tetteh, Velizar Efremov, Nils D Forkert, Matthias Schneider, Jan Kirschke, Bruno Weber, Claus Zimmer, Marie Piraud, and Bjørn H Menze. Deepvesselnet: Vessel segmentation, centerline prediction, and bifurcation detection in 3-d angiographic volumes. Frontiers in Neuroscience, 14:592352, (2020)

[25] Jierong Wang and Albert CS Chung. Higher-order flux with spherical harmonics transform for vascular analysis. In Proc. MICCAI, pages 55–65. Springer, (2020)

[26] Y. Wang, A. Narayanaswamy, and B. Roysam. Novel 4-D open-curve active contour and curve completion approach for automated tree structure extraction. In Proc. CVPR, pages 1105–1112, (2011)

[27] Yuqiang Zhang, Min Liu, Fuhao Yu, Tieyong Zeng, and Yaonan Wang. An o-shape neural network with attention modules to detect junctions in biomedical images without segmentation. IEEE J. Biomed. Health Inform., 26(2):774–785, (2021)

[28] Kerschnitzki, Kollmannsberger et al., "Architecture of the osteocyte network correlates with bone material quality." Journal of Bone and Mineral Research, 28(8):1837–1845, (2013)

## Authors' background

*This form helps us to understand your paper better, the form itself will not be published. You can delete it when you prepare the final paper if it's accepted.*

<table><thead><tr><td>Your Name</td><td>Title*</td><td>Research Field</td><td>Personal website</td></tr></thead></table>
<table><tr><td>Hongda Liang</td><td>Master student</td><td>Medical Image Processing</td><td></td></tr><tr><td>Da Chen</td><td>associate professor</td><td>Medical Image Processing</td><td></td></tr><tr><td>Tao Chen</td><td>Master student</td><td>Medical Image Processing</td><td></td></tr><tr><td>Li Liu</td><td>lecture</td><td>Medical Image Processing</td><td></td></tr><tr><td>Jiong Zhang</td><td>associate professor</td><td>Medical Image Processing</td><td></td></tr><tr><td>Laurent D. Cohen</td><td>full professor</td><td>Image Processing</td><td></td></tr></table>

*Title can be chosen from: master student, PhD candidate, assistant professor, lecture, senior lecture, associate professor, full professor