Anisotropic tubular minimal path model with fast marching front freezing scheme

In this work, we introduce an anisotropic minimal path model based on a new Riemannian tensor integrating the crossing-adaptive anisotropic radius-lifted tensor ﬁeld and the front freezing indicator by appearance and path features. The non-local path feature only can be obtained during the geodesic distance computation process by the fast marching method. The predeﬁned criterion derived from path feature is able to steer the front evolution by freezing the point causing high bending of the geodesic to solve the shortcut problem. We performed qualitative and quantitative experiments on synthetic and real images (including retinal vessels, rivers and roads) and compare with the minimal path models with classical anisotropic Riemannian metric and dynamic isotropic metric, which demonstrated the proposed method can detect desired targets from complex tubular tree structures


Introduction
Tubular structure extraction is a crucial task in many computer vision applications. Various tubular structure extraction algorithms and techniques have been studied, which can be roughly divided into two categories: supervised and unsupervised methods [1,2]. Supervised approaches usually achieve better per- 5 formance due to the prior information from the training components to decide whether a pixel should belong to a vessel or not. These models construct feature vectors in order to train the respective classifiers [3,4] or apply the neural network [5,6] with the help of manual label. However, it is very time-consuming to find discrimination features for constructing powerful models by training. Be- 10 sides, the obtained models always meet for particular requirement. Those unsupervised methods, without manual labeling information, usually have lower computation complexity and more simple procedure [2], such as the pattern recognition techniques [7,8], model-based approaches [9, 10], tracking-based algorithms [11,12] and so on. 15 The minimal path technique [13] has been adopted to detect tubular structures thanks to its efficiency and global optimality [14,15,16,17]. It is often taken as an interactive tool such that users are allowed to provide manual intervention such as source points, which constrains the minimal path to delineate the tubular structures. In essence, a geodesic path is a curve linking two points derived by globally minimizing a weighted curve length [13]. It is defined by integrating a geodesic metric P : Ω × R d → R + , where Ω ⊂ R d stands for the open and bounded domain and d is the dimension of Ω, along a regular path γ: where γ : [0, 1] → Ω ⊂ R d is a regular curve of Lipschitz continuity. The minimal path models have been extensively exploited for tubular structure segmentation, where both the centerline positions and thickness of the tubular structures are required [15,16,17]. The isotropic tubular minimal path model [15] only makes use of the curve position information, which may increase the risk of shortcuts. 20 More general anisotropic Riemannian metric [16,18] is proposed to segment tubular structures in conjunction of the pre-detected vessel anisotropy features.
However, this extended tubular minimal path model still has difficulties when dealing with complicated situation such as extract a weak vessel from a tree structure. Moreover, some researchers try to improve the classical tubular min- 25 imal path models by reducing the user intervention [19,20,21]. For these models, only one (or several) initial root point(s) is (are) necessary to extract a tubular tree structure. Recently, the higher order properties of the geodesic curves are exploited for shortest path estimation as well as tubular structure segmentation. Specifically, Ulen et al. [22] proposed a curvature and torsion 30 regularized shortest path model by computing the curvature and the torsion.
Based on the Eikonal equation framework, the curvature-penalized minimal path models [23,24,25] are established over an orientation-lifted space such that the orientation dimension can be used to represent the curve tangents.
Controlling the course of the fast marching front propagation is a solution for 35 the above problem. Criterions are defined to make the front propagation stop in undesired directions or promote in specific directions [26,27]. In contrast to the minimal path models using static geodesic metrics which are fixed during fast marching front propagation, the dynamic minimal path models allow to update the respective geodesic metrics in the course of geodesic distance estimation.

40
In [28], the authors proposed a new fast marching front propagation scheme in order to address the shortcut and short branches combination problems. The basic idea is to freeze the front points of which the respective local path features violate the pre-defined criteria. In this case the fast marching fronts will not pass through the points corresponding to undesired geodesic paths. Chen et al. [29] 45 proposed an adaptively Riemannian metric using a dynamic metric construction way in order to penalize the appearance feature coherence property. In this case the obtained geodesic paths favour to pass by the region of slow-varying tubular features. Geoffrey et al. [30] created a connectivity metric to determine the reasonable pathways of connection. Windheuser et al. [31] estimated the 50 curvature information by computing the angles between each pair of adjacent edges. Krueger et al. [32] incorporated the curvature estimation into a greedy core algorithm to improve the efficiency. However, these models track geodesic paths only in the spatial space where the vessel thickness information cannot be extracted.

55
In this paper, we focus on the dynamic minimal path-based tubular structure segmentation, for which the geodesic metrics are established over a radius-lifted space. The existing tubular minimal path models [15,16,17] often suffer from the shortcut and short branches combination problems, especially for the case that the target is weakly-defined and crosses or parallels with strong ones, see In this figure, the target vessel has low gray level contrast, while its neighbour vessels are strongly-defined (i.e., high gray level contrast). In Fig. 1(c), we show the minimal path extraction result obtained by the anisotropic tubular minimal path model [16]. One can observe short branches combination occur.
In other words, the obtained minimal path prefers to pass through strong seg-65 ments belonging to different vessels. In order to simultaneously obtain the centerlines and thickness of tubular structures as well as to overcome the shortcut and short branches combination problems, we proposed an anisotropic minimal path model for tubular structure segmentation via a dynamic radius-lifted anisotropic metric and state-of-the-art anisotropic fast marching method [33].

70
The motivation of this work is that the anisotropic tubular minimal path model 4 [16] cannot take into account path features for distance estimation. Using the path features detected during the front propagation for minimal path computation is able to obtain desired results even in complicated situation. The main contribution in this paper lies at the introduction of a dynamic anisotropic Rie-75 mannian metric in the radius-lifted space using back-tracked geodesic paths. In The paper outline is shown as follows. In Section 2, we briefly introduce the background for the anisotropic radius-lifted Riemannian metric. Section 3 80 describes the main contribution of this work including the proposed dynamic Riemannian metric in a radius-lifted space and its numerical implementation in conjunction with the fast marching method. Section 4 shows the experimental results and the final section 5 gives the conclusion.

85
In this section, we give the background on anisotropic radius-lifted minimal path model [16] for tubular structure segmentation.

Tubular Geometry Features Computation
In this section, we first consider the extraction of the local tubular geometry features from the image I : Ω → R through the optimally oriented flux (OOF) filter [34]. Note that the local geometry features can also be extracted using other steerable filters such as [34,35,36]. The output of the OOF filter is estab- is the radius space 1 . The output of the OOF filter at the position x and scale r is a symmetric matrix of size d × d as follows: 1 In the spaceΩ, a pointx ∈Ω is an ordered pair (x, r) where x ∈ Ω ⊂ R 2 and r ∈ [r min , rmax].

5
where (∂ i,j G σ ) is the Hessian matrix of the Gaussian kernel G σ with variance σ. χ Sr is a characteristic function of a disk with radius r.

90
The eigenvalues λ i (i = 1, 2, · · · , d) extracted using the OOF filter are the values of the oriented flux along the corresponding eigenvectors v i as following: Without loss of generality, we assume that λ 1 (·) ≤ λ 2 (·) ≤ · · · λ d (·). In this paper, we only consider the 2D tubular segmentation problem, which means that d = 2. For the case that the tubular structures have lower gray levels than background, the vector v 1 (x, r * ) indicates the orientation of the tubular structure at the point x inside the vessel structures, where r * is the optimal 95 vessel scale of point x.

Radius-lifted Minimal Paths for Tubular Segmentation
We denote by S + d the set of symmetric positive definite matrices of size d × d (d = 2, 3) and let Lip([0, 1],Ω) be the set of Lipchitz continuous curves γ : [0, 1] →Ω. The radius-lifted minimal path model is proposed in [15] for 100 tubularity segmentation. Benmannsour and Cohen [16] generalize this isotropic model [15] to the anisotropic case. The basic idea for both models is to lift a spatial geodesic to the radius-lifted spaceΩ by adding one abstract dimension to represent the vessel radii.
The path length L for such a radius-lifted curve γ ∈ Lip([0, 1],Ω) can be measured through a Riemannian tensor field M scale :Ω → S + 3 , where ·, · represents the Euclidean scalar product. The tensor field M scale is constructed [16] by the eigenvalues λ i and eigenvectors v i , where M aniso :Ω → S + 2 is a tensor field associated to the spatial anisotropy with size 2×2 and P scale :Ω → R + is a scalar function. Both can be constructed by combining the eigenvalues and eigenvectors of the OOF response F in Eq. (2) as follows: where α ∈ R and β ∈ R + are two parameters which control the regularization 105 of the spatial and radius dimensions, respectively. Note that the value of α is negative for dark-vessel and bright-background case and α > 0, otherwise.
In the Riemannian minimal path model, once the metric tensor is given, one can define the relevant geodesic distance Uŝ(x) as the minimal curve length betweenx and the fixed source pointŝ as follows where the values of Uŝ can be regarded as the arrival times of a front propagating with oriented velocity related to the metric tensor M −1 scale . The geodesic distance map Uŝ from the source pointŝ satisfies the anisotropic Riemannian Eikonal equation with boundary condition Uŝ(ŝ) = 0, where u M = √ u T M u for any matrix M ∈ S + 3 . A geodesic pathĈx ,ŝ linkingx to the source pointŝ can be tracked from the pointx by solving the following ordinary differential equation (ODE) till the source pointŝ is reached withCx ,ŝ (0) =x. The geodesicCx ,ŝ is parameterized by its arc-length. The  Overview. The main goal in this section is to establish a new anisotropic geodesic metric in a radius-lifted space in order to avoid the shortcut and short branches combination problems as shown in Fig. 1. In contrast to [16] for which only the local geometry information are considered to establish the geodesic metrics, the front freezing-based minimal path model [28] constructs the geodesic metrics by taking into account the nonlocal information such as the local path.
During the geodesic distance front propagation, the front points for which the features do not satisfy the given criterion will be frozen, where the features for each front point are extracted through two extra points. This scheme proposed in [28] is able to avoid the shortcut problem in some extent. However, it is restricted to the case of centerline detection and cannot take advantages of the path orientation due to the isotropic nature of the used metrics in [28]. In this section, we extend this front frozen scheme to a radius-lifted space through an anisotropic Riemannian metric. This is done by constructing a new Riemannian tensor field M dyn :Ω → S + 3 during the geodesic distance computation which is actually carried out in a front advancing procedure. The invoked tensor field M dyn consists of two ingredients: the crossing-adaptive anisotropic radius-lifted tensor field M adap :Ω → S + 3 and the front freezing indicator δ :Ω → {1, ∞}. It can be expressed by We will explain the construction details for the tensor field M adap in Section 3.1 115 and the scalar indicator δ in Section 3.2, respectively.

Computation of the Crossing-Adaptive Tensor Field
By applying the OOF filter to an image involving vessel tree structures, one can obtain an optimal direction v 1 (x, r) at each point (x, r) belonging to the tree structures through eigen-decomposition to the response of the filter OOF We suppose that inside the tubular structures the intensities are lower than background, thus the corresponding eigenvalue λ 2 can be used to characterize the appearance feature of the vessels, by which we define a vesselness map The value of ζ(x) is derived from λ 2 at the optimal scale η(x) ∈ [r min , r max ], where the η : Ω → [r min , r max ] is a map defined by Based on the optimal scale map η and the eigenvectors v 1 , we can define the tubular feature vector p : Ω → R 2 that indicates the orientation for a tubular structure The tubular direction p(x) well describes the vessel orientation at non-crossing point x. However, for crossing points, denoted by y, there will be at least two vessels across one another. This means that the highly anisotropic tensors problem is to adaptively remove the anisotropy from the tensors M aniso (y, ·), which can be done by invoking the orientation scores [2] or by structure tensors [37]. In this section, we construct our crossing-adaptive radius-lifted tensor field through the tool of structure tensor.
The shortcut and short branches combination problems usually occur where the vessel crosses its neighbours, since the eigenvector of a crossing point usually indicates the orientation of the vessel with strong appearance feature. So the speed computed from the anisotropic geodesic metric is slower along the weak vessel direction than along the strong one. To solve the shortcut and short branches combination problems, we make use of a crossing-adaptive tensor field M adap :Ω → S + 3 established in the radius-lifted spaceΩ, the anisotropy of which is kept at single-vessel region and removed or reduced at crossing points.

It can be constructed by
where α is a positive constant which has been used in Eqs. (6) and (7). The tensor field T smooth :Ω → S + 2 is computed via a Gaussian kernel G p with standard derivative p and the identity matrix I d of size 2 × 2 where ∈ R + is a sufficiently small constant to generate a regular matrix M adap (·). Note that for any radius values r 1 , r 2 ∈ [r min , r max ] where r 1 = r 2 , one always has T smooth (x, r 1 ) = T smooth (x, r 2 ). The scalar-valued function : Ω → R + 0 is a weighted function over the image domain Ω. In this paper, we propose three ways to compute the function , where the first way is to invoke the vesselness map ζ (see Eq. (13)).
Alternatively, the weighted function can be computed as a binary-valued function, which depends on a threshold value Th ∈ R + and also on the vesselness map ζ The last method utilizes the skeleton of the vessel mask. The function is used 130 in order to reduce the influence from the regions outside the vessel structures.
In other words, the information inside the vessel region dominates the tensor field T smooth .
In Eq. (18) the crossing point, the tensor field T smooth is impacted greater by the vessel with bigger radius. Note that the main difference between T smooth and the structure tensor field used in [37] lies at the existence of the function . As discussed above, it can reduce or avoid the effects derived from the non-vessel regions.
This means that we are able to make use of a Gaussian filter with more flexible 150 variance values p.
We express the tensor filed T −1 smooth by One can see that the Gaussian filtering is operated on the feature vectors p(y)

Computation of the Front-Freezing Indicator
In this section, we construct the proposed front-freezing indicator δ with non-local path features derived from the local truncated geodesic paths to solve the shortcut and short branches combination problems. The non-local feature is  10)) on the obtained geodesic distance map Uŝ, which has already been computed within the region where the front visited [28]. During the front propagation, one can easily track a radius-lifted geodesic pathCx ,ŝ parameterized by its arc-length by solving the corresponding gradient descent ODE (10). Sincẽ Cx ,ŝ lies in the radius-lifted space, we haveCx ,ŝ = (G x,s , η) whereG x,s (·) ∈ Ω is the projected path and η(·) ∈ [r min , r max ]. We first detect a pointq = (q, r q ) fromCx ,ŝ such that its projected path |G x,q | = Γ, where |G x,q | stands for the Euclidean length of G x,q and Γ is a positive constant. Note thatq is the first extra point and the second one, referred to asm = (m, r m ), is defined as the middle point of the part of the pathCx ,q . The process to search the two extra points are shown in Fig. 4(a). The bending measure K is estimated by the angle between two vectors related to the two extra points q, m ∈ Ω as follows [28] K where ·, · denotes the scalar product and ||·|| represents the norm of the vector.

170
The values of K are ranged at [−1, 1]. The estimation of K is done during the fast marching front propagation, which will be described in Section 3.3. We mention that a small value of K implies that the corresponding geodesic path likely appears with a sharp turning. Findx min = (x min , r min ), the T rial point which minimizes Uŝ;

11:
for allŷ ∈ N * (x min ) and L(ŷ) = Accepted do where K 0 is a given threshold. The bending measure K(x) > K 0 which yields δ = 175 1 means the fast marching front is propagated as usual. By using the indicator δ, the shortcut and short branches combination problems can be avoided as shown in Fig. 1. Our method with an appropriated threshold is able to detect the desired vessel although the neighbouring vessel has stronger appearance feature.  Uŝ(x, r) The obtained geodesic map H s (x) is shown on the original vessel images as in

Discussion
The proposed method in essence constructs the geodesic metric in a dynamic way and we discuss the differences between our method and two state-of-the-art 235 dynamic minimal path approaches [28,29].
• In work [28], instead of freezing the front points. Sometimes one needs to carefully control the importance between the tubular confidence information and the coherence penalization. Secondly, the dynamic metric in [29] is still built in the spatial domain 3 while our method can detect both centerlines and thickness in one step. Finally, we analyze the implementation difference 260 on the adaptive establishment. In [29], the authors build the dynamic metric using a reference point form a truncated geodesic path vector to distinguish the correct vessel direction in crossing points. In contrast, the basic idea in the the proposed work for the considered crossing-adaptive anisotropy is to reduce or remove the anisotropy property at the crossing 265 points and to keep the anisotropy property at the single vessel points, done by the tool of structure tensor.

Experimental Results and Discussion
In this section, we conduct the numerical experiments on qualitative and quantitive comparison between state-of-the-art minimal path models and the  compare with the ArR metric [16], the variant of the isotropic dynamic minimal path model (dIsoM) [28] and the appearance feature coherence-penalized Riemannian metric with adaptively anisotropy (AFC) [29]. Note that for the 20 sake of fair comparison, the dIsoM metric-driven minimal paths are extracted by submitting the scalar-valued function P scale defined in Eq. (7) to the geodesic 295 distance propagation scheme with front freezing procedure as depicted in Algorithm 1.
Evaluation Score. In order to evaluate the proposed dynamic minimal path models quantitatively, we consider an accuracy measurement to validate the segmentation results. Let S represent the segmentation region from the considered models and G denotes the region from ground truth data. In addition, # · be the number of grid points within the set. Thus the measurement is defined as One can point out that the measurement R is ranged within the interval [0, 1], where R = 1 means that the vessel segmentation is exactly identical to the ground truth. When the measurement R is used in the vessel centerlines com-300 parison, the ground truth region G is the dilated centerline from the ground truth by a disk of radius 2 (in grid points) and S represents the set of grid points located in the detected centerlines.

Experiments on Synthetic Images
In this section, we compare the proposed dArR metric to the ArR met-  as described in Table. 1, one can point out that the proposed dArR metric-driven minimal paths are indeed in sensitive to the effects from noise. The minimal paths from the dArR metric can benefit from the use of the front freezing scheme thus can improve the accuracy and robustness of the anisotropic fast marching 320 propagation, leading to favourable centerline and boundary extraction results.
In Fig. 8, we show the minimal paths driven by the dArR metric as well as the AFC metric. In this figure, there are two crossing structures with almost identical gray levels, where the target is longer than another one in the sense of Euclidean curve length [29]. The geodesic path derived using the AFC metric 325 fails to detect the target as shown in Fig. 8  from the image data. Thus in case the target vessels are smooth (without sharp 335 turning), the dArR metric can obtain desired minimal paths.

Experiments on Retinal Images
In this section, we validate the proposed minimal path model on 40 retinal patches obtained from the DRIVE dataset [38] and 40 retinal patches from the IOSTAR dataset [2]. Each patch includes a retinal artery vessel which is near In Fig. 9, we illustrate the effect of the caArR metric for which the anisotropy is reduced at the crossing points. In column (a), the yellow lines indicate the anisotropy feature vectors of the tubular structures detected by the OOF filter. 350 We can observe that the anisotropy feature vectors located at the cross section are approximately collinear to the directions of the vessel of strong appearance feature. This may lead to incorrect geodesic tangents on these crossing points In column (b), we can see that the short branches combination problem occurs.

360
The minimal paths prefer to combine vessel segments with strong appearance features. While in the column (c), the paths obtained from the caArR metric can get the target vessels correctly. Finally, the segmented vessels by the caArR metric are illustrated in column (d).
In Fig. 11, we show the comparison results of the proposed dArR to the results from the long vessel crossing or near with strong ones. The third row shows the detect results from the image with long and tortuosity vessels. The last two rows denote detected geodesics of vessel regions with strong tortuosity.
We observe that the classical ArR metric suffered from short combination prob- racy results for the vessel regions with respect to different metrics are shown in Table. 2. The detection accuracy of the vessel centerlines is shown in Table. 3.
From these quantitive results, we can see that the proposed metrics indeed overcomes the shortcuts and short branches combination problems and outperform 380 all the other compared metrics.
The dIsoM 4 and dArR metrics (also including the dArR-M and dArR-S metrics) are able to take into account the curve bending measure to freeze front points so as to avoid shortcuts and short branches combination problems. From Tables. 2 and 3, we can see that these metrics achieve better vessel segmenta-385 4 As mentioned above, dIsoM metric is a variant of the one proposed in [28].  We analyze the computation time of each component of the proposed method.

405
There are mainly three steps in the proposed method including the vesselness computation, the crossing-adaptive anisotropic radius-lifted tensor field construction and the fast marching front propagation. The front-freezing indicator is computed during the fast marching front propagation. The results are shown in Table. 4.

Experiments on Satellite Images
In this section, we conduct the experiments on satellite images with rivers, bridges and streets from Google Earth to evaluate the proposed metric.  In Fig. 14, we illustrate the qualitative comparison for the ArR metric, the caArR metric and the dArR metric, which are shown from the second to the 415 fourth columns. The first row shows the segmentation results from the image with river crossing several bridges whose intensities are similar to these in the background region. The ArR and caArR metrics-driven minimal paths miss the river direction from the middle located bridge, where the bending of the geodesic path varies strongly. The bending constraint is able to effectively handle this 420 problem. The second row illustrates the segmentation results of a satellite image with two neighbouring rivers and complex background. We aim to detect the smaller one by given two points. We observe that the correct river region cannot be recognized by the ArR metric. The dArR-based minimal paths are capable of tracking the favourable river correctly without being effected by the stronger

Conclusion
In this paper, there are mainly three contributions of the proposed method to overcome the shortcut and short branches combination problems in tubular structure detection. Firstly, we propose a crossing-adaptive anisotropic tensor 430 field in the radius-lifted space, which is used to reduce the in favourable effects from the incorrect vessel anisotropy features at the crossing points. This is done by reducing the anisotropy in the crossing points adaptively and simultaneously keeping the high anisotropy in the single vessel points. Secondly, we propose a new way to reduce the influence from the region outside the vessel structure.

435
Finally, we show a way for exploiting front freezing scheme for the fast marching algorithm in a radius-lifted space. We show the practical application of the proposed dynamic minimal path model on synthetic images, retinal images and satellite images. The results demonstrate that our method indeed outperforms classical minimal path methods with static geodesic metrics and the dynamic 440 isotropic metric. In the future work, the proposed method will be extended to 3D tubular structure detection application and try to reduce the manually intervention.
valuable comments, which greatly improved this manuscript. This work was