Trajectory Grouping with Curvature Regularization for Tubular Structure Tracking

Let $\Omega\subset\mathbb{R}^{2}$ be an open and bounded image domain and let $\operatorname{Lip}([0,1],\Omega)$ denote the set of curves $\gamma:[0,1]\to\Omega$ with Lipschitz continuity. An energy functional, or the weighted curve length, of a curve $\gamma\in\operatorname{Lip}([0,1],\Omega)$ encoding curvature penalization can be formulated as

where $\kappa:[0,1]\to\mathbb{R}$ stands for the curvature of $\gamma$, and $\gamma^{\prime}=d\gamma/du$ is the first-order derivative of $\gamma$. One can see that the energy (1) involves two cost functions $\mathfrak{C}_{0}$ and $\psi$. Specifically, the function $\mathfrak{C}_{0}:\Omega\times\mathbb{R}^{2}\to\mathbb{R}^{+}$ can be derived from the image data using a steerable filter. Basically, $\mathfrak{C}_{0}(x,\mathbf{v})$ is expected to have low values if the point $x$ is close to a tubular centerline and the direction $\mathbf{v}$ is proportional to the orientation that a tubular structure should have at $x$. Furthermore, for the FE and FSR minimal path models, the cost function $\psi$ can be respectively formulated as $\psi:=\psi_{\rm FE}$ and $\psi:=\psi_{\rm FSR}$ such that

where $\beta\in\mathbb{R}^{+}$ is a scalar parameter that controls the importance of the curvature $\kappa$.

Let $\tilde{\Omega}=\Omega\times\mathbb{S}^{1}$ be an orientation-lifted space of higher dimension, where $\mathbb{S}^{1}=[0,2\pi)$ is an interval with periodic boundary condition, and denote by $\mathbf{n}_{\theta}=(\cos\theta,\sin\theta)^{\top}$ a unit vector of an angle $\theta\in\mathbb{S}^{1}$. As discussed in [27, 4], a key idea for minimizing the energy (1) is to lift a planar curve $\gamma\in\operatorname{Lip}([0,1],\Omega)$ to the space $\tilde{\Omega}$ in conjunction with a parametric function $\tau:[0,1]\to\mathbb{S}^{1}$ defined being such that

yielding that

The orientation lifting of $\gamma$ yields a new curve $\tilde{\gamma}=(\gamma,\tau)$ subject to Eq. (2) and its first-oder derivative obeys $\tilde{\gamma}^{\prime}(u)=(\gamma^{\prime}(u),\tau^{\prime}(u))$ for any $u\in[0,1]$. Note that each point $(x,\theta)\in\tilde{\Omega}$ lying in an orientation-lifted curve $\tilde{\gamma}=(\gamma,\tau)$ has three coordinates, where the first two coordinate $x\in\Omega$ indicates the physical positions and the third coordinate $\theta\in\mathbb{S}^{1}$ characterizes the tangent vector $\gamma^{\prime}(u)$. With these definitions in hands, the minimization of the energy (1) can be efficiently addressed by seeking the minimal curve length of orientation-lifted curves measured through an orientation-lifted Finsler metric $\mathcal{F}_{\epsilon}:\tilde{\Omega}\times\mathbb{R}^{3}\to\mathbb{R}^{+}$, which implicitly encodes the curvature terms. Specifically, for any point $\tilde{x}=(x,\theta)\in\tilde{\Omega}$ and any vector $\tilde{\mathbf{u}}=(\mathbf{u},\nu)\in\mathbb{R}^{3}$, a general form of $\mathcal{F}_{\epsilon}$ can be expressed as

where $\mathfrak{C}:\tilde{\Omega}\to\mathbb{R}^{+}$ is an orientation-dependent cost function derived from image data, subject to $\mathfrak{C}(x,\theta)=\mathfrak{C}_{0}(x,\mathbf{n}_{\theta})$. The computation for the cost function $\mathfrak{C}$ is presented in Section II-B. In addition, the metric $\mathfrak{F}_{\epsilon}$ only involves the curvature penalty, thus independent to image data.

Specifically, the metric $\mathfrak{F}_{\epsilon}$ for the FE minimal path model [27] reads

While for the FSR minimal path model [4], one has

For an orientation-lifted curve $\tilde{\gamma}$, the energy of $\mathcal{L}_{\epsilon}(\tilde{\gamma})$ measured by the metric $\mathcal{F}_{\epsilon}$ can be expressed as follows:

As discussed in [27, 4], the minimum of the energy $\mathcal{L}_{0}(\gamma)$ in Eq. (1) can be well approximated by the length of a geodesic path associated to the metric $\mathcal{F}_{\epsilon}$ as $\epsilon\to 0$, providing that $\mathfrak{C}_{0}(x,\mathbf{n}_{\theta})=\mathfrak{C}(x,\theta)$. In this way, the minimization problem of the energy $\mathcal{L}_{0}(\gamma)$ is transferred to minimizing the new energy $\mathcal{L}_{\epsilon}(\tilde{\gamma})$, where the later problem can be efficiently addressed by the eikonal PDE framework.

The geodesic distance map very often lends itself for the minimization of the length $\mathcal{L}_{\epsilon}$. For a fixed source point $\tilde{a}\in\tilde{\Omega}$, the geodesic distance map defines a minimal curve length for each point $\tilde{x}\in\tilde{\Omega}$

It is well known that the geodesic distance map $\mathcal{U}_{\tilde{a}}$ satisfies the Eikonal equation such that $\mathcal{U}_{\tilde{a}}(\tilde{a})=0$ and for any orientation-lifted point $\tilde{x}\in\tilde{\Omega}\backslash\{\tilde{a}\}$ we have

The eikonal equation (9) can be solved by using the state-of-the-art Hamiltonian fast marching method [5], in terms of a Hamiltonian reformulation of Eq. (9) . A geodesic path $\mathcal{C}_{\tilde{a},\tilde{x}}$ linking from $\tilde{a}$ to $\tilde{x}$ can be derived by re-parametering the solution $\mathcal{C}$ (which is also a geodesic path) to a gradient descent ordinary differential equation (ODE) on the geodesic distance map $\mathcal{U}_{\tilde{a}}$

For two given points with tangents, the minimal paths associated to the data-driven curvature-penalized metric $\mathcal{F}_{\epsilon}$ attempt to keep smooth, since $\mathcal{F}_{\epsilon}$ implicitly encodes curvature penalization as regularization. In next section, we present the construction details for the curvature-penalized metric $\mathcal{F}_{\epsilon}$.

The function $\mathfrak{C}$ can be estimated from the image data via a steerable filter [27]. In the context of tubular structure tracking, the optimally oriented flux (OOF) filter [33] is an effective tool for extracting geometry features from images, which will be adopted in this paper. For clarity, we assume that the tubular structures are supported to have locally lower intensities than background.

Let $G_{\sigma}$ be a Gaussian kernel with standard deviation $\sigma$ and let $\{\partial_{x_{i}x_{j}}G_{\sigma}\}_{i,j}$ be the Hessian matrix of the kernel $G_{\sigma}$. The response of the OOF filter on a gray level image $I:\Omega\to\mathbb{R}$ at a point $x$ and a scale $r\in[R_{\rm min},R_{\rm max}]$ is a symmetric matrix of size $2\times 2$

where “$\ast$” stands for the convolution operator. The function $\chi_{r}:\Omega\to\{0,1\}$ is the indicator for a disk of radius $r$

The eigenvalues of $\Psi(\cdot)$, denoted by $\lambda_{1}(\cdot)$ and $\lambda_{2}(\cdot)$ with assumption $\lambda_{1}(\cdot)\leq\lambda_{2}(\cdot)$, can be used to characterize the appearance of tubular structure centerlines. As in [33], one can compute a scalar-valued function $\zeta:\Omega\to\mathbb{R}^{+}_{0}$, referred to as vessel score map, such that $\zeta(x)$ is derived from $\lambda_{1}(x,\cdot)$ at the optimal scale

By the eigenvalues $\lambda_{1}(\cdot)$, one can also derive an optimal scale map $\rho:\Omega\to[R_{\rm min},R_{\rm max}]$ as follows:

The score $\zeta(x)$ indicates the likelihood of a point $x$ belonging to a tubular structure. As in [33, 3], the direction of a tubular structure at a point $x$ can be characterized by $\mathbf{v}_{2}(x)$ which is the eigenvector of $\Psi(x,\rho(x))$ corresponding to the eigenvalue $\lambda_{2}(x,\rho(x))$. The tool of the orientation scores extends the confidence map $\Psi$ to a multi-orientation space. Denoting by $\psi_{\rm os}:\Omega\times\mathbb{S}^{1}\to\mathbb{R}^{+}$ the orientation scores associated to an image, we have

where $\mathbf{n}_{\theta}^{\perp}=(-\sin\theta,\cos\theta)$ is a unit vector that is perpendicular to $\mathbf{n}_{\theta}$.

Then the image data-driven function $\mathfrak{C}$ used in Eq. (4) can be computed as

where $\alpha\in\mathbb{R}^{+}$ is a weighting parameter on the image data.

In many applications, tubular structures involved in the images may consist of multiple tree structures. In image analysis applications, a fundamental task is to trace a path to delineate a curvilinear structure between two given points from the entire tree structures or from a complicated background. The classical tubular minimal path models [3, 24, 22] likely tend to travel along the regions with strong appearance features or directly pass through the background regions, yielding the shortcuts or short branches combination problems. Alternatively, exploiting minimal paths with curvature regularization for tracing tubular structures may reduce the risk of the above problems in some extent. However, modeling a tubular centerline as a globally optimal curvature-penalized geodesic path is not always suitable in practical applications. Moreover, the computation complexity for the curvature-penalized minimal path models are too high for real-time applications, since we expect the tracking time to be comparable to the user interaction time. Examples for illustration of the shortcuts and short branches combination problems can be seen in Figs. 1 and 2. In Fig. 1, the goal is to extract an artery vessel between two points from a retinal image patch. The results show that both the anisotropic tubular model and the FSR model suffer from the shortcuts and short branches combination problems, as depicted in Figs. 1b and 1c.

Moreover, we make a test in Fig. 2 with a goal to delineate a spiral from a synthetic image interrupted by noise. In this experiment, we evaluate the anisotropic tubular model [3], the progressive model [26], the FSR model and the proposed one. The results from the existing models suffer from the shortcuts problem due to the effects from the noise. In order to overcome these problems mentioned above, we propose a new minimal path model for minimally interactive tubular structure tracing in conjunction with curvature-penalized minimal paths and trajectories grouping. The resulting path from the proposed model can avoid these problems as much as possible, as depicted in Fig. 1d and in column $5$ of Fig. 2.

Basically, tubularity segmentation is regarded as a way of classifying all points into either tubular structures or background. The segmentation is usually represented by a binary mapping, which can be generated by many tubular structure segmentation approaches as reviewed in literature [1, 2]. In this paper, we implement a simple method to segment the tubular structures by directly thresholding a vessel score map $\zeta$, see Eq. (13). An example for the visualization of the score map $\zeta$ that is derived from a retinal image patch can be seen in Fig. 3b.

Following that we apply the morphological filters on tubularity segmentation to get the skeleton structure of entire tubular structure network. In this paper, we suppose that each trajectory is a connected set involving skeleton points. In order to obtain disjoint trajectories, we remove all the branch points from the computed tubular skeleton structures, where a branch point is a skeleton point connecting at least three trajectories. We illustrate an example for the disjoint trajectories in Fig. 3c, where each trajectory is tagged by different colors. Moreover, in order to reduce the computation complexity, we remove the trajectories for which the length (in grid points) is lower than a given thresholding value.

The proposed model for tracing tubular structure centerlines is established on a graph $\mathcal{G}=(V,E)$, where $V$ represents the set of nodes and $E$ stands for the set of edges between connected nodes. We denote by $e_{ij}\in E$ an edge joining two adjacent nodes $\vartheta_{i}$, $\vartheta_{j}\in V$. Each edge $e_{ij}$ is assigned a weight $\omega_{ij}\in\mathbb{R}^{+}_{0}$. For the sake of simplicity, we leverage the convention that the weight $\omega_{ij}=+\infty$ implies the node $\vartheta_{i}$ is disconnected to $\vartheta_{j}$. In tubular structure tracking, one may expect to get the same tracking result by choosing either of the two prescribed points as the source point and taking the remaining one as end point. For this purpose, we focus on the indirect graph, which means that for each pair of edges $e_{ij}$ and $e_{ji}$, one has $\omega_{ij}=\omega_{ji}$.

As in [30], the graph $\mathcal{G}$ is constructed by associating a node $\vartheta_{i}\in V$ to a trajectory $\mathcal{T}_{i}$. The edge set $E$ characterizes the connection between two trajectories $\mathcal{T}_{i}$ and $\mathcal{T}_{j}$ for any $i\neq j$. For a fixed trajectory $\mathcal{T}_{i}$, the joint trajectories $\mathcal{T}_{j}$ can be detected by a tubular neighbourhood $M_{i}\subset\Omega$ surrounding $\mathcal{T}_{i}$. A possible way is to perform front propagation associated to a suitable metric emanating from $\mathcal{T}_{i}$ such that $M_{i}$ can be identified by thresholding the propagated geodesic distance. However, this method may increase the computation complexity of the entire algorithm. Alternatively, we make use of the method proposed in [30] which uses Euclidean distance. Specifically, the trajectory $\mathcal{T}_{i}$ is first prolonged from its two endpoints along the respective tangents, in order to generate an prolonged trajectory $T_{i}$. Then we build a regular tubular neighbourhood $M_{i}$ for the extended trajectory by

where $\tau$ is a given thresholding value. A trajectory $\mathcal{T}_{j}$ is said to be adjacent to the fixed $\mathcal{T}_{i}$ for any $j\neq i$, if the prolongated trajectories $T_{j}$ satisfies

Given two vertices, optimizing the graph $\mathcal{G}$ by Dijkstra’s algorithm [31] yields a shortest path that consists of a series of ordered trajectories. Once the sets $V$ and $E$ for the graph $\mathcal{G}$ are built, one can point out that the weight $\omega_{ij}$ dominates the tracing results. The weight $\omega_{ij}$ measures the cost of bridging the gap between a pair of adjacent trajectories $\mathcal{T}_{i}$ and $\mathcal{T}_{j}$. In contrast to [30] which exploits straight segments, we complete the gap between a pair of adjacent trajectories by curvature-penalized geodesic paths.

The generation of a set of orientation-lifted trajectories constitutes the first step in the proposed model. Basically, a trajectory $\mathcal{T}_{i}\subset\Omega$ likely passes through the centerline of a tubular structure. Thus, a point $x\in\mathcal{T}_{i}$ can be assigned to a pair of orientations $\theta_{x}\in[0,\pi]$ and $\theta_{x}+\pi$ which characterize the directions of the centerline at $x$. Such an orientation $\theta_{x}$ can be computed via the orientation scores $\psi_{\rm os}$ as follows

Accordingly, a trajectory $\mathcal{T}_{i}$ can be lifted to the space $\tilde{\Omega}$, leading to a new subset $\widetilde{\mathcal{T}}_{i}\subset\tilde{\Omega}$ as a union of two sets

We say that two orientation-lifted trajectories $\widetilde{\mathcal{T}}_{i}$ and $\widetilde{\mathcal{T}}_{j}$ for any $i\neq j$ are adjacent if their physical projections $\mathcal{T}_{i}$ and $\mathcal{T}_{j}$ are adjacent.

In the proposed model, we expect to choose a set of trajectories $\mathcal{T}_{i}$ to constitute a shortest path which delineates the target tubular structure, providing that two vertices are given. In principle, the weights $\omega_{ij}$ should be as low as possible if two adjacent trajectories $\mathcal{T}_{i}$ and $\mathcal{T}_{j}$ lie at the same tubular structure. For a fixed trajectory $\mathcal{T}_{i}$ lying in the target structure, the construction for the weights must allow to differentiate between a trajectory that belongs to the same structure with $\mathcal{T}_{i}$ and one belongs to another structure. In many scenarios, tubular structures appear to be locally smooth. Thus it is reasonable to estimate the edge weights $\omega_{ij}$ using curvature-penalized geodesic distance. The smoothness property of the curvature-penalized minimal paths agrees with the requirement for connecting the gaps between two adjacent vessel segments.

For a fixed trajectory $\widetilde{\mathcal{T}}_{i}$, we consider a geodesic distance map $\mathcal{U}_{i}$ with respect to $\widetilde{\mathcal{T}}_{i}$ such that

where $\mathcal{L}_{\epsilon}(\tilde{\gamma})$ is the weighted curve length of $\tilde{\gamma}$ associated to the metric $\mathcal{F}_{\epsilon}$, see Eq. (7). As discussed in Section II, the geodesic distance map $\mathcal{U}_{i}$ admits the eikonal PDE (9) with a boundary condition $\mathcal{U}_{i}(\tilde{x})=0$ for any point $\tilde{x}\in\widetilde{\mathcal{T}}_{i}$.

By the definition (III-C), we can estimate a distance value $\mathcal{D}_{i,j}$ that is the minimal weighted curve length between the fixed trajectory $\widetilde{\mathcal{T}}_{i}$ and a neighbouring one $\widetilde{\mathcal{T}}_{j}$

Once the distance map $\mathcal{U}_{i}$ is computed, a geodesic path $\mathcal{C}_{i,j}\in\operatorname{Lip}([0,1],\tilde{\Omega})$ such that $\mathcal{C}_{i,j}(0)\in\widetilde{\mathcal{T}}_{i}$ and $\mathcal{C}_{i,j}(1)\in\widetilde{\mathcal{T}}_{j}$ can be tracked by performing the gradient descent ODE on the distance map $\mathcal{U}_{i}$, see Eq. (10).

The minimal length $\mathcal{D}_{i,j}$ associated to the metric $\mathcal{F}_{\epsilon}$ encodes the image data-driven function $\mathfrak{C}$ defined in Eq. (4), due to the use of the metric $\mathcal{F}_{\epsilon}$. Accordingly, the length $\mathcal{D}_{i,j}$ encodes the tubular appearance features. This may introduce bias to the cost of connecting the trajectory $\mathcal{T}_{i}$ to any neighbour $\mathcal{T}_{j}$, especially when the target structures weaker than its neighbouring ones. In order to remove the effect from the image data, we take into account a new weighted length $\mathfrak{D}_{i,j}$ measured along the geodesic path $\mathcal{C}_{i,j}$ derived from the distance map $\mathcal{U}_{i}$. Let us denote $\mathcal{C}_{i,j}:=(\mathscr{C}_{i,j},\tau_{i,j})$ where $\mathscr{C}_{i,j}(u)$ indicates the spatial positions of the orientation-lifted geodesic path $\mathcal{C}_{i,j}$, and the parametric function $\tau_{i,j}(u)$ determines the tangent $\mathscr{C}^{\prime}_{i,j}$, see Eq. (2). We expect that the quantity $\mathfrak{D}_{i,j}$ is explicitly dependent to the curvature of $\mathscr{C}_{i,j}$, but independent to the image data $\mathfrak{C}$. Thus the quantity $\mathfrak{D}_{i,j}$ can be estimated by

where $\kappa_{i,j}$ denotes the curvature of the physical projection $\mathscr{C}_{i,j}$. The parameter $\beta\in\mathbb{R}^{+}$ is a constant controlling the importance of the curvature. Note that the estimation of $\mathfrak{D}_{i,j}$ by Eq. (23) allows to give more importance to the curvature penalization by setting a larger value to parameter $\beta$ than the one used for defining $\mathcal{D}_{i,j}$. Such a quantity $\mathfrak{D}_{i,j}$ can be computed directly through the path $\mathcal{C}_{i,j}$. In this paper, we also consider an alternative way for estimating the quantity $\mathfrak{D}_{i,j}$ relies on a new map $\mathcal{E}_{i}:\tilde{\Omega}\to\mathbb{R}^{+}_{0}$. For each point $\tilde{x}\in\tilde{\Omega}$, we can obtain a geodesic path $\mathcal{C}=(\mathscr{C},\tau)$ such that $\mathcal{C}(0)\in\widetilde{\mathcal{T}}_{i}$ and $\mathcal{C}(1)=\tilde{x}$. Similar to Eq. (23), we set

where $\kappa$ represents the curvature of the physical projection $\mathscr{C}$. We denote by $\tilde{x}_{i}^{*}$ an optimal point such that $\mathcal{D}_{i,j}=\mathcal{U}_{i}(\tilde{x}_{i}^{*})$, yielding that $\mathfrak{D}_{i,j}=\mathcal{E}_{i}(\tilde{x}_{i}^{*})$. Note that both of the maps $\mathcal{U}_{i}$ and $\mathcal{E}_{i}$ can be computed simultaneously, without tracking the geodesic paths $\mathcal{C}$, as discussed in [23, 34]. In the following experiments, we apply the map $\mathcal{E}_{i}$ for the estimation of the value $\mathfrak{D}_{i,j}$.

Similar to the computation of the geodesic path $\mathcal{C}_{i,j}$, we can generate a different geodesic path $\mathcal{C}_{j,i}=(\mathscr{C}_{j,i},\tau_{j,i})$ traveling from the trajectory $\mathcal{T}_{j}$ to $\mathcal{T}_{i}$ subject to $\mathcal{C}_{j,i}(0)\in\widetilde{\mathcal{T}}_{j}$ and $\mathcal{C}_{j,i}(1)\in\widetilde{\mathcal{T}}_{i}$. By Eq. (23), we can obtain a weighted curve length $\mathfrak{D}_{j,i}$ associated to $\mathcal{C}_{j,i}$.

As a consequence, we define the weight $\omega_{ij}$ for the edge $e_{ij}$ as follows

The definition in Eq. (25) imposes $\omega_{ij}=\omega_{ji}$, where $\omega_{ij}$, $\omega_{ji}$ are the respective weights for the edges of $e_{ij}$ and $e_{ji}$. Following that, the path completing the gap between $\mathcal{T}_{i}$ and $\mathcal{T}_{j}$ can be chosen as the minimal one between $\mathcal{C}_{i,j}$ and $\mathcal{C}_{j,i}$ in the sense of the weighted curve length formulated in Eq. (23). In the following, the geodesic paths for connecting the gaps are referred to as bridge paths.

We conduct the numerical experiments with both qualitative and quantitative comparison to the progressive minimal path model with bending constraint (Progressive) [26], the anisotropic tubular minimal path model (Aniso) [3], the curvature-penalized minimal path model with a FSR metric [4], the grouping method with Euclidean distance and angles (Group-Angle) [30] on both synthetic and real images. We refer to the proposed model as Group-FE (resp. Group-FSR) if the edge weights $\omega_{ij}$ of the graph $\mathcal{G}$ is estimated by the FE metric (resp. FSR metric).

The Progressive minimal path model invokes a dynamic isotropic metric by taking into account nonlocal path features, which introduces the bending constraint into the fast marching fronts propagation scheme, in order to avoid the shortcuts and short branches combination problems. The path feature is computed using a backtracked truncated geodesic path whose length is fixed as $10$ grid points. The thresholding value that defines the maximally admissible bending degree is set to $0.9$, see [26] for more details. In the following experiments, we use the same model as [3] to construct the anisotropic Riemannian metric for the Aniso model. The anisotropic ratio for this metric is fixed to $10$ for all the tests.

The FSR and FE minimal path models descried in Section. II are also considered in the experiments. For the FSR and FE metrics, we fix the parameter $\epsilon=0.1$ to achieve good balance between the accuracy and the computation complexity. The parameters $\alpha$ and $\beta$ are weighted parameters which control the importance of the image data and the curvature penalization, respectively. We set $\alpha=5$ and $\beta=20$ for all the tests. In addition, both metrics are also used to estimate the edge weights $\omega_{ij}$ for the proposed Group-FE and Group-FSR models. For the Group-Angle model, the same edge set $V$ as the proposed Group-FE model is adopted to construct the graph. We exploit the identical method as in  [30] for the estimation of the edge weights $\omega_{ij}$.

In the OOF filter, the Gaussian kernel $G_{\sigma}$ with the standard deviation $\sigma$ is used to slightly smooth of the input image. The interval $[R_{\rm min},R_{\rm max}]$ is the range of the possible radii. We set $\sigma=1.5$, $R_{\rm min}=1$, and $R_{\rm max}=8$. For numerical consideration, we normalize $\psi_{\rm os}$ to the range $[0,1]$.

Finally, all the experiments are performed on a standard 6-core Intel Core i7 of 3.2GHz architecture with 64Gb RAM.

In this work, we compare the proposed model with the state-of-the-art minimal path-based tracing algorithms.

In Fig. 5, we show the tubular centerline tracing results on our synthetic images. The average size of the four images is $474\times 321$ grid points. The gray levels of these synthetic images are normalized to the range $[0,1]$ interrupted by additive Gaussian noise of normalized standard derivation $\sigma_{n}=0.15$. In each image, there exist two tubular structures of low gray levels, along which the directions vary smoothly. Our goal is to extract the one between two given points. In rows $1$ and $2$, one can see that the source and end points are close to one another and the target structure has a long Euclidean length. In rows $3$ and $4$, the target structures have strong tortuosity appearance. In Fig. 5, columns $2$ to $5$ illustrate the results from the Aniso, Progressive, FSR and Group-Angle models, where we have observed the occurrence of the shortcuts and short branches combination problems. While the paths generated by the proposed Group-FSR and Group-FE models can correctly trace the objective vessels, as depicted in columns $6$ and $7$columns $5$ and $6$, due to the use of smooth geodesic paths for the connection of gaps between trajectories.

In Fig. 6, we qualitatively compare the proposed Group-FSR model¹¹1We observed that tracing results from the proposed Group-FE model are almost identical to those from Group-FSR model. For better visualization, we only show the results from the Group-FSR model. to the Aniso, Progressive, FSR and Group-Angle models on from the DRIVE [35] and IOSTAR [36] datasets. The original retinal images are RGB color images, where we only use the green channel of each test image [36] for the computation of the cost function $\mathfrak{E}$, see Eq. (16). The goal is to trace an artery vessel between the given points indicated by dots. The corresponding ground truth is depicted in column $1$. We note that each model is performed in the entire retinal image and we only show the image patch in gray level containing the target vessel for well visualization. In retinal images, the main challenge is that artery vessels appear weaker and are often close to or even cross over veins of strong visibility. In this experiment, some portions of the target vessels in general appear to be strongly tortuous. This may introduce difficulties to the Progressive and FSR models. Moreover, the strong neighbours of the target vessels will highly affect the results of the Aniso model. Finally, using straight segments to link adjacent trajectories for vessels of high length may accumulate the approximation errors. In columns $2$ to $5$ of Fig. 6, we can see that the artery vessel tracing results from the existing models again suffer from the shortcuts and short branches combination problems. In contrast, the proposed Group-FSR model passed by the correct way, see column $6$.

In Figs. 7 and 8, we present the tubular structure centerline tracing results from the Aniso, Progressive, FSR, Group-Angle and the proposed Group-FSR models on road and river images. In these experiments, the color images are first converted into gray levels in a preprocessing step, such that the cost function $\mathfrak{E}$ is extracted from these gray images. However, we still draw the tubular structure tracking results on the original color images for better visualization. In Fig. 7, we show the qualitative results on aerial images of road networks [20]. The main difficulty usually lies at the complicated background such as buildings which may also produce strong tubular appearance features. From columns $2$ to $5$ of Fig. 7, one can observe the shortcuts occur when implementing the existing models. While for the results from the proposed model, the objective structures can be correctly traced thanks to the proposed edge weights construction way. In Fig. 8, we demonstrate the extracted paths indicated by red lines on road and river patches from satellite images. The target tubular structures surrounded by complicated background also appear very long and have strong tortuosity. The tracing results derived from the existing results are shown in columns $2$ to $5$, from which we observe paths with shortcuts. In contrast, the proposed Group-FSR model can follow the road and river structures successfully.

The qualitative comparisons on synthetic images, retinal vessel images and natural images presented in Section IV-B prove that the proposed Group-FSR and Group-FE models indeed obtain promising results. In order to evaluate the proposed models in a rigorous and convincing manner, here we run the quantitative comparisons on the synthetic images and retinal images, respectively. The accuracy that measures the performance of the tested models is carried out by a score $\mathcal{J}$ defined as follows:

where $S$ is the set of grid points passed through by the evaluated paths, $GT$ denotes the region of the ground truth, and $\#|S|$ stands for the elements involved in the set $S$. The accuracy score $\mathcal{J}$ is ranged within the interval [0,1], where higher values of $\mathcal{J}$ means better performance on tracing tubular structures.

We first present the performance of the compared models on the synthetic images shown in Fig. 5. The ground truth $GT$ is set as a binary segmentation of all pixels corresponding to the target structures. The accuracy scores $\mathcal{J}$ computed from the evaluated paths are shown in Table I. One can point out that proposed Group-FSR and Group-FE models achieve the best performance among all the considered models.

In Table II, we show the quantitative evaluation results on the retinal images from the datasets of DRIVE [35] and IOSTAR [36]. Tracing artery vessels from retinal images providing that only few user-provided points are given is a challenging task, thus can well measure the performance of the considered models. For a quantitative evaluation, we compare our proposed Group-FSR and Group-FE models to the Aniso, Progressive, FSR and Group-Angle models, in the sense of the accuracy score $\mathcal{J}$ defined in Eq. (26). We have extracted the ground truth regions for each tested artery vessel from the artery-vein ground truth images. Note that each artery-vein ground truth image classifies each grid point belonging to either artery, vein or background. We made use of $394$ artery vessels sampled from two datasets, where most of the major artery vessels in each retinal image are taken into account for evaluation. Among all the tests, we provide only $2$ points for $321$ vessels and the remaining cases require $3$ points. Again, we note that each individual artery vessel is tracked using the entire image. In each test, all the evaluated models are under the same user-supplied points. In Table. II, we illustrate the average accuracy scores for different models. In this table, we observe that the proposed Group-FSR and Group-FE models demonstrate large gaps to the other tested models, proving the advantages of the proposed models. In addition, we can see that in each dataset the average values of $\mathcal{J}$ for the proposed Group-FSR and Group-FE models are similar to each other and achieve around $12\%$ higher than the Group-Angle model. Note that we utilized the same node and edge sets to construct the graphs for the Group-Angle model and the proposed models, except for the computation of the edge weights.

Finally, we show more statistical results on the accuracy scores $\mathcal{J}$ for all the tested models through the tool of box plots, as depicted in Fig. 9. The left and right columns in Fig. 9 correspond to the results on the DRIVE and IOSTAR datasets, respectively. The results shown in this figure again prove that the proposed models outperform the compared state-of-the-art models in both robustness and accuracy.