Vessel tree extraction using radius-lifted keypoints searching scheme and anisotropic fast marching method

Abstract

Geodesic methods have been widely applied to image analysis. They are particularly efficient to extract a tubular structure, such as a blood vessel, given its two endpoints in a 2D or 3D medical image. We address here a more difficult problem: the extraction of a full vessel tree structure given a single initial root point, by growing a collection of keypoints or new initial source points, connected by minimal geodesic paths. In this article, those keypoints are iteratively added, using a new detection criteria, which utilize the weighted geodesic distances with respect to a radius-lifted Riemannian metric, the standard Euclidean curve length and a path score. Two main weaknesses of classical keypoints searching approach are that the weighted geodesic distance and the Euclidean path length do not take into account the orientation of the tubular structure or object boundaries, due to the use of an isotropic geodesic Riemannian metric, and suffer from a leakage problem. In contrast, we use an anisotropic geodesic Riemannian metric, and develop new criteria for selecting keypoints based on the path score and automatically stopping the tree growth. Experimental results demonstrate that our method can obtain the expected results, which can extract vessel structures at a finer scale, with increased accuracy.

Introduction

Minimal path models have been widely applied to the field of image processing, ¹ like tubular structure and object boundary extraction, since the pioneer work proposed by Cohen and Kimmel (Cohen–Kimmel model), ² where the tubular structures or the object boundaries are extracted under the form of minimal paths with respect to a isotropic Riemannian metric $\mathcal{P}$ over the image domain. The geodesic Riemannian metric can be isotropic ² if the metric $\mathcal{P}$ depends only on the positions, or anisotropic in the sense that the metric $\mathcal{P}$ also depends on the path orientation.^3,4 It can be defined on the physical space (image domain), or involve additional abstract variables such as the extracted vessel radius.^4,5 Once this Riemannian metric is properly defined, the single-pass fast marching methods^6–8 are the favoured methods to estimate geodesic distances, ² because of their low computational complexity. Finally, the minimal path can be extracted from the geodesic distance map, also called the action map, by simply solving an ordinary differential equation (ODE).

Modern geodesic methods for image processing involve two domains: the physical image space $\overline{\Omega }\subset ℝn,n=2$ or 3, and a higher dimensional and model dependent geodesic domain $\Omega \subset ℝm,m\ge n$ . The classical Cohen–Kimmel approach ² simply sets $\Omega :=\overline{\Omega }$ . However, for tubular structure extraction tasks, especially for extracting the centrelines and boundaries at the same time, the Cohen–Kimmel model does not work well. More recent and accurate feature extraction methods, pioneered by Li and Yezzi ⁵ and including this work, thus define Ω as the product of physical space $\overline{\Omega }$ with a parameter space $[R\text{min},R\text{max}]$ representing vessel radius, i.e. $\Omega =\overline{\Omega }\times [R\text{min},R\text{max}]$ , so that $m=n+1$ ; vessel location and size are thus extracted and estimated in a single unified step. In this article, we focus us on the 2D vessel tree extraction such that we fix the image domain $\overline{\Omega }\subset ℝ2$ and the radius lifted domain $\Omega =\overline{\Omega }\times [R\text{min},R\text{max}]\subset ℝ3$ .

The curve length $\ell (\gamma )$ , weighted by Riemannian metric $\mathcal{P}(·,·)$ , is defined by:

$\ell (\gamma ):={\int }_0^1\mathcal{P}(\gamma (t),\stackrel{·}{\gamma }(t))\text{d}t$

(1)

where $\stackrel{·}{\gamma }=\frac{d}{\mathrm{dt}}\gamma (t)$ . Fix the Riemannian metric $\mathcal{P}$ , and a source point s. The geodesic distance $\mathcal{U}(\mathbf{x})$ , or minimal action map, is the minimal weighted length of any path joining a point $\mathbf{x}\in \Omega$ from the source s:

$\mathcal{U}(\mathbf{x}):=inf\{\ell (\gamma );\gamma \in I,\gamma (1)=\mathbf{x},\gamma (0)=\mathbf{s}\}$

(2)

where ℑ denotes the collection of Lipschitz paths $\gamma :[0,1]\to \Omega$ . Under standard assumptions, there always exists at least one minimal geodesic path $\mathcal{C}\mathbf{s},\mathbf{x}$ , and only a single one for most x.

The classical Cohen–Kimmel model ² relies on the isotropic metrics of the form defined on $\overline{\Omega }$ :

$\mathcal{P}n(x,\mathbf{u}s)=(Pn(x)+\omega )\Vert \mathbf{u}s\Vert ,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\forall x\in \overline{\Omega },\forall \mathbf{u}s\in ℝ2$

(3)

where $Pn:\overline{\Omega }\to ℝ+$ is a potential function taking lower values around the interesting features of the images, and the positive constant w controls the smoothness of the minimal curve $\mathcal{C}$ .

Li and Yezzi ⁵ used metrics of the same form of $\mathcal{P}n$ (3), but with the above-mentioned additional radius dimension. Let $\mathbf{x}=(x,r)$ , where $x\in \overline{\Omega }$ and $r\in [R\text{min},R\text{max}]$ . Then, the radius-lifted metrics of Li-Yezzi model ⁵ can be formulated as: where $Pl:\Omega \to ℝ+$ is the radius-lifted potential function, which takes lower values along the centrelines of the tubular structures. Parameter ε is a positive constant to balance the speeds along the physical space and radius space. Both the models, Cohen and Kimmel ² and Li and Yezzi ⁵ are special instances of the more general Riemannian metric framework.

$\mathcal{P}(\mathbf{x},\mathbf{u}):=\sqrt{〈\mathbf{u},\mathcal{M}(\mathbf{x})\mathbf{u}〉}$

(4)

where $\mathcal{M}$ is a symmetric positive definite tensor field. The Cohen–Kimmel metric has the form

$\mathcal{M}(x)=(Pn(x)+\omega )2\mathbf{I}d,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\forall x\in \overline{\Omega }$

and the isotropic metric of Li and Yezzi model ⁵ is

$\mathcal{M}(\mathbf{x})=(Pl(\mathbf{x})+\omega )2\mathbf{D}\varepsilon ,\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\forall \mathbf{x}\in \Omega$

where $\mathbf{D}\varepsilon$ is a $(d+1)$ -dimensional diagonal matrix:

$\mathbf{D}\varepsilon =\text{diag}\{1,\dots ,1,\varepsilon \}$

Metrics like these which are proportional to the identity tensor $\mathbf{I}d$ , or to a diagonal matrix $\mathbf{D}\varepsilon$ , are referred to as isotropic.

However, minimal path models based on the isotropic metrics^2,5 often suffer from the shortcuts problem. The metric does not properly steer the paths along the vessels, and the minimal paths tend to take shortcuts outside the structure of interest. To solve this problem, Benmansour and Cohen ⁴ enriched these Riemannian metric constructions through the use of an anisotropic tensor field $\mathcal{M}r$ as formulated in (8), given by general symmetric positive definite matrices, which eigenvectors typically align with the orientation of the tubular structure.

Recently, a curvature penalty minimal path model was proposed by Chen et al., ⁹ where the extracted minimal paths compromise between path length and path smoothness. This work relies on an approximation of Euler’s elastica bending energy, by a Finsler metric. To automatically extract the retinal segments, Chen et al. proposed a novel region constrained minimal path model, ¹¹ which restricted the extract geodesic inside the given constrained region to overcome the overlapping extraction problem. Liao et al. ¹² proposed a curvature thresholding minimal path model, which compute the isotropic metric of point x by taking into account the curvature of a back tracked local geodesic leading by x. Chen et al. proposed a dynamic anisotropic Riemannian metric based geodesic model, ¹³ which penalized the local intensity consistency while preserving the tubular structure anisotropy. This model was applied to interactive retinal vessel extraction, with only two endpoints given, to overcome the short branches combination problem.

All the minimal path models mentioned above require user input endpoints as the prior knowledge to track the minimal paths. To reduce the user input, Benmansour and Cohen ¹⁴ proposed a new approach: a keypoints searching method to detect recursively new startpoints (or keypoints) along the expected features named, see ‘Path score based keypoints searching method’ section. Kaul et al. ¹⁵ proposed a new stopping criterion for both open and closed curve detection and a new method to compute the Euclidean curve length different to the keypoints methods proposed by Benmansour and Cohen ¹⁴ and Deschamps and Cohen. ¹⁶ Li et al. ¹⁷ proposed to detect the tubular structure using the extra radii model ⁵ and keypoints searching method. ¹⁴

In this article, we proposed a new vessel tree extraction method based on the automatic keypoints detection and state-of-the-art anisotropic fast marching method, ⁸ see Mirebeau ¹⁸ for an open source library. The pairwise distances between keypoints are fixed by a curve length threshold; classical keypoints searching method ¹⁴ tends to overlook parts of the vessel tree when the curve length threshold is small, and to leak or take shortcuts between distinct branches of the vessel tree for large curve length threshold (for details, see ‘Advantages of using a small curve length threshold and a path score judgement’ section and Fig. 3). We substantially improve the performance for small curve length threshold, by incorporating in the keypoint detection criterion, a path score computed from the optimally oriented flux vesselness map for each keypoint candidate detected by the classical definition.

In summary, our contribution is as follows:

We redefine the keypoints to make them reasonable even for small curve length threshold.

We present the extension of the curve length calculation in the anisotropic case, in itself a different contribution.

We give a stopping criterion to automatically stop the keypoints searching scheme.

The rest of this article is organized as follows: in ‘Background’ section, we give the scientific background including the optimally oriented flux filter, ¹⁹ the radius-lifted minimal path model based on the anisotropic Riemannian metric ⁴ and the numerical solver: anisotropic fast marching method. ⁸ In ‘Path score based keypoints searching method’ section, we introduce the proposed keypoint searching scheme. Finally, in ‘Experiments’ section, we give some experimental results.

Convention: if $\mathbf{x}\in \Omega$ denotes a point with extra radii dimension, then we use $x\in \overline{\Omega }$ to denote a point with the same space location but without extra radii dimension, i.e. $\mathbf{x}=(x,r)$ .

Background

Optimally oriented flux and vesselness map

The oriented flux ¹⁹ of an image I, of dimension 2, is defined by the amount of the image gradient projected along the orientation $\mathbf{n}$ flowing out from a 2D circle at point x with radius r:

$f(x;r,\mathbf{p})=\int \partial Sr(\nabla (G\sigma *I)(x+r\mathbf{n})·\mathbf{p})(\mathbf{p}·\mathbf{n})\text{d}a$

(5)

where $G\sigma$ is a Gaussian function with variance σ and n is the outward unit normal vector along $\partial Sr$ . da is the infinitesimal length on $\partial Sr$ . The oriented flux is a quadratic function: one has

$f(x;r,\mathbf{p})=\mathbf{p}\text{T}\mathcal{Q}(x,r)\mathbf{p}$

for some symmetric matrix $\mathcal{Q}(x,r)$ , which eigenvalues and eigenvectors we denote by $\lambda i(x,r)$ and $\mathbf{v}i(x,r)$ . Recalling that $\mathbf{x}=(x,r)$ , one then has:

$\mathcal{Q}(\mathbf{x})=\sum _i^2\lambda i(\mathbf{x})\mathbf{v}i{\mathbf{v}}_i^{\text{T}}(\mathbf{x})$

(6)

Law et al. ¹⁹ used the normalized sum of the non-zero eigenvalues to compute the vesselness map, which takes its largest values for points x in vessel centrelines. Indeed, the eigenvector of $\mathcal{Q}(\mathbf{x})$ oriented tangentially to the vessel is associated to a small or zero eigenvalue, whereas another eigenvector oriented transversally to the vessel is associated to large eigenvalue (in magnitude). Here, we suppose that pixels inside a vessel have stronger intensity value than the background. Without loss of generality, we assume that $\lambda 1\le \lambda 2\approx 0$ . Thus, we can define the vesselness map as,

$ℎ(x)=max\{\underset{r}{max}\{\frac{-1}{r}\lambda 1(x,r)\},0\}$

(7)

The vesselness map $ℎ(x)$ in equation (7) will have a large value if x is located inside the vessels, which means that we can take into account this map to guide the searching scheme of the keypoints as all the keypoints should inside the vessels.

Radius-lifted minimal path model with anisotropic Riemannian metric

In this article, we use the same metric tensor construction as in Benmansour and Cohen ⁴ : the metric tensor $\mathcal{M}r$ , defined in the radius lifted domain Ω with additional radii dimension for the vessel extraction, is as follows:

$\mathcal{M}r(\mathbf{x})=(\mathcal{M}s(\mathbf{x})00Pr(\mathbf{x}))$

(8)

where $\mathcal{M}s$ is a symmetric positive definite tensor with size 2 × 2:

$\mathcal{M}s(\mathbf{x})=\exp (\alpha \lambda 2(\mathbf{x}))\mathbf{v}1(\mathbf{x}){\mathbf{v}}_1^{\text{T}}(\mathbf{x})+\exp (\alpha \lambda 1(\mathbf{x}))\mathbf{v}2(\mathbf{x}){\mathbf{v}}_2^{\text{T}}(\mathbf{x})$

and P_r is a scalar function with implicit dependence on x:

$Pr(\mathbf{x})=\beta \exp (\alpha \frac{\lambda 1(\mathbf{x})+\lambda 2(\mathbf{x})}{2})$

The constant $\beta >0$ controls the fast marching propagation speed along the radius direction, and α controls the metric anisotropic ratio μ:

$\mu =\underset{\mathbf{x}\in \Omega }{max}\sqrt{\frac{\exp (\alpha \lambda 2(\mathbf{x}))}{\exp (\alpha \lambda 1(\mathbf{x}))}}$

The weighted curve length ℓ of a Lipschitz curve, $\gamma :[0,1]\to \Omega$ , is defined by

$\ell (\gamma )={\int }_0^1\sqrt{〈\stackrel{·}{\gamma }(t),\mathcal{M}r(\gamma (t))\stackrel{·}{\gamma }(t)〉}\text{d}t$

(9)

The minimal action map or geodesic distance with respect to ℓ (9) can be defined as

$\mathcal{U}(\mathbf{x})=inf\{\ell (\gamma );\gamma \in I,\gamma (1)=\mathbf{x},\gamma (0)=\mathbf{s}\}$

(10)

where ℑ denotes the collection of Lipschitz paths $\gamma :[0,1]\to \Omega$ and s is a fixed initial source point. It satisfies the anisotropic Eikonal equation:

$\{\Vert \nabla \mathcal{U}(\mathbf{x})\Vert {\mathcal{M}}_r^{-1}(\mathbf{x})=1,\text{forall}\mathbf{x}\in \Omega \setminus \mathbf{s};\mathcal{U}(\mathbf{s})=0$

(11)

Numerically, this Eikonal equation can be solved by the start-of-the-art lattice basis reduction based anisotropic fast marching method, ⁸ which will be introduced in the following. The minimal path $\mathcal{C}$ can be extracted by solving the following ODE:

$\stackrel{·}{\mathcal{C}}(t)\propto -{\mathcal{M}}_r^{-1}(\mathcal{C}(t))\nabla \mathcal{U}(\mathcal{C}(t))$

(12)

$\text{with}\mathcal{C}(\mathbf{s})=0$

(13)

using the standard Runge–Kutta numerical method or the robust method proposed in Mirebeau. ⁸

Anisotropic fast marching method using lattice basis reduction

The popular numerical solver for the minimal action map $\mathcal{U}(\mathbf{x})$ described in (2) is the fast marching method. It introduces a discretization grid Z of Ω, and for each $\mathbf{x}\in Z\cap \Omega$ , a stencil or a mesh $S(\mathbf{x})$ of a neighbourhood of x with vertices in the grid Z (with the adequate modification if x is at or near the boundary). An approximation of the minimal action map $\mathcal{U}$ is given by the solution of the following fixed point system: find $\mathcal{U}:Z\to ℝ$ such that (a) $\mathcal{U}(\mathbf{s})=0$ and (b) for all $\mathbf{x}\in Z\setminus \mathbf{s}$

$\mathcal{U}(\mathbf{x})=\underset{\mathbf{y}\in \partial S(\mathbf{x})}{min}\{\mathcal{P}(\mathbf{x},\mathbf{y}-\mathbf{x})+IS(\mathbf{x})\mathcal{U}(\mathbf{y})\}$

(14)

where I_S denotes piecewise linear interpolation on the boundary of stencil S. $\mathcal{P}$ is the anisotropic Riemannian metric defined by the radius lifted tensor $\mathcal{M}r$ (8) as follows:

$\mathcal{P}(\mathbf{x},\mathbf{y}-\mathbf{x})=\sqrt{〈(\mathbf{y}-\mathbf{x}),\mathcal{M}r(\mathbf{x})(\mathbf{y}-\mathbf{x})〉}$

(15)

The expression (14) is based on the fact that the minimal geodesic $\mathcal{C}\mathbf{s},\mathbf{x}$ , joining x to the initial source point s, needs to cross the stencil boundary $\partial S(\mathbf{x})$ at some point y; hence, it is the concatenation of a small path $\mathcal{C}\mathbf{y},\mathbf{x}$ joining x to y, of approximate length $\mathcal{P}(\mathbf{x},\mathbf{y}-\mathbf{x})$ , and of $\mathcal{C}\mathbf{s},\mathbf{y}$ , whose energy is approximated by interpolation.

For concreteness, we give a second description of (14), closer to implementation and specialized to three-dimensional domains as in the two-dimensional domains with radius lifting case. Opting for offset-based notations, we introduce the translated stencil:

$\mathcal{T}(\mathbf{x}):=\mathcal{S}(\mathbf{x})-\mathbf{x}$

which is a collection of tetrahedra whose union is a neighbourhood of the origin. A generic boundary stencil point $\mathbf{y}\in \partial S(\mathbf{x})$ as in optimization (14) can then be expressed as:

$\mathbf{y}=\mathbf{x}+\sum _{i=1}^3\eta i\mathbf{u}i$

(16)

$\mathrm{st}\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\eta i>0,\text{and}\sum _{i=1}^3\eta i=1$

(17)

where $\mathbf{u}1,\mathbf{u}2$ and $\mathbf{u}3$ are non-zero vertices of a common tetrahedron $T\in \mathcal{T}(\mathbf{x})$ . Consider the compact and convex set:

$\Theta :=\{\eta =(\eta 1,\eta 2,\eta 3);\eta i>0\text{and}\sum _{i=1}^3\eta i=1,i=1,2,3\}$

For each grid point, $\mathbf{x}\in Z$ and simplex $T\in \mathcal{T}(\mathbf{x})$ , we consider the function:

$\mathcal{J}\mathbf{x},T(\eta )=\Vert \sum _{i=1}^3\eta i\mathbf{u}i\Vert \mathcal{M}r(\mathbf{x})+\sum _{i=1}^3\eta i\mathcal{U}(\mathbf{x}+\mathbf{u}i)$

(18)

One can see that $\mathcal{J}\mathbf{x},T$ is a convex function on the set of Θ. Minimizing $\mathcal{J}\mathbf{x},T$ can be found in literature.^8,20 Then, (14) is equivalent to

$\mathcal{U}(\mathbf{x})=\underset{T\in \mathcal{T}}{min}\{\underset{\eta \in \Theta }{min}\{\mathcal{J}\mathbf{x},T(\eta )\}\}$

(19)

We denote the minimizer of $\mathcal{U}$ in equation (19) as $\eta *=({\eta }_1^{*},{\eta }_2^{*},{\eta }_3^{*})$ . The above one describes a discrete approximation of the minimal action map $\mathcal{U}$ , given by the solution of an apparently complex system of N-coupled non-linear equations of the form (14). This N-dimensional fixed point system, with $N=\#(Z)$ , can be solved using the single pass fast marching algorithm,^7,8 provided the stencil $S(\mathbf{x})$ satisfies some geometric properties depending on the local Riemannian metric $\mathcal{P}(\mathbf{x},·)$ . An adaptive construction of such stencils was introduced in Mirebeau, ⁸ which led to breakthrough improvements in terms of computation time and accuracy for strongly anisotropic geodesic energy potentials, as in our application. It invokes Lattice Basis Reduction, a tool from discrete geometry, which combines in an optimal way the geometric structure given by the Riemannian metric, and the arithmetic structure of the Cartesian discretization grid. We emphasize that these rather sophisticated techniques are hidden to the user, and point researchers interested in anisotropic fast marching to the open-source code joined to Mirebeau. ¹⁸

Path score based keypoints searching method

Euclidean curve length calculation

The minimal geodesic $\mathcal{C}\mathbf{s},\mathbf{x}$ , joining point x to the source s and minimizing the action $\mathcal{U}(\mathbf{x})$ (2), is generically unique. Its Euclidean curve length:

$\mathcal{L}(\mathbf{x}):=\mathcal{L}(\mathcal{C}\mathbf{s},\mathbf{x})$

can be defined as

$\mathcal{L}(\mathbf{x})=\int \mathcal{C}\mathbf{s},\mathbf{x}\Vert \stackrel{·}{\mathcal{C}}\mathbf{s},\mathbf{x}(t)\Vert \text{d}t$

(20)

which is a crucial ingredient of our algorithm. A natural approach to computing the Euclidean length map $\mathcal{L}$ , is to extract for each x a minimal geodesic $\mathcal{C}\mathbf{s},\mathbf{x}$ by solving the ODE (12), and then calculating its Euclidean curve length. This first approach turns out unfortunately to be too expensive in terms of computation time. In literature,^14,16 the authors proposed a fast method for approximating the curve length during the fast marching front propagation. However, their methods are based only on isotropic metric over the physical domain. In this article, we extend this method to the radius-lifted anisotropic Riemannian metric case.

Let $\mathcal{K}$ denote the initial source point collection. An approximation of $\mathcal{L}$ is given by the solution of the following fixed point problem: find $\mathcal{L}:Z\to ℝ$ such that (a) for all $\mathbf{x}=(x,r)\in \mathcal{K},\mathcal{L}(\mathbf{x})=0$ and (b) for all $\mathbf{x}\in Z\setminus \mathcal{K}$ , denoting by $\mathbf{y}*$ the point at which the minimum (14) is attained

$\mathcal{L}(\mathbf{x})=\Vert \mathbf{y}*-\mathbf{x}\Vert M2+\mathcal{L}(\mathbf{y}*)$

(21)

where $\mathcal{L}(\mathbf{y}*)$ can be obtained by linear interpolation. The norm $\Vert \mathbf{u}\Vert M2=\sqrt{〈\mathbf{u},M2\mathbf{u}〉}$ and M₂ with size 3 × 3 is defined as

$M2=(1\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}0\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}00\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}1\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}00\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}0\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}0)$

(22)

In offset-based coordinates, the minimizer of the fixed point problem formulated in (19) is denoted by $\eta *=({\eta }_1^{*},{\eta }_2^{*},{\eta }_3^{*})$ . Then, equation (21) can be expressed with respect to $\eta *$ :

$\mathcal{L}(\mathbf{x})=\Vert \sum _{i=1}^3{\eta }_i^{*}\mathbf{u}i\Vert M2+\sum _{i=1}^3{\eta }_i^{*}\mathcal{L}(\mathbf{x}+\mathbf{u}i)$

(23)

Equation (23) means that we only compute the approximated curve length of path Γ, where $\gamma =(\Gamma ,r)$ . A single pass solve is again possible: whenever the fast marching algorithm updates $\mathcal{U}(\mathbf{x})$ , simultaneously update $\mathcal{L}(\mathbf{x})$ , using the (just computed) minimizer $\eta *$ from (19).

Path score calculation

The score of a path $\gamma =(\Gamma ,r)$ is obtained by averaging the vesselness map ℎ computed from the optimally oriented flux filter, ¹⁹ see (7).

$\mathcal{P}\mathcal{S}(\Gamma ,\text{Th})=\frac{\int \Gamma ℎ(\Gamma (s))\delta (\text{Th},ℎ(\Gamma (s)))\text{d}s}{\int \Gamma \delta (\text{Th},ℎ(\Gamma (s)))\text{d}s}$

(24)

The threshold parameter $\text{Th}>0$ and selector $\delta (·,·)$ are used to eliminate irrelevant parts of the path

$\delta (\text{Th},d)=\{1,\text{if}d<\text{Th};0,\text{otherwises}$

(25)

For each $\mathbf{x}=(x,r)$ , we denote by $\mathcal{P}\mathcal{S}(x,\text{Th})=\mathcal{P}\mathcal{S}(\Gamma ,\text{Th})$ , the path score associated to the geodesic $\gamma =(\Gamma ,r)$ joining x to the initial source point s.

Keypoint definition with a path score

Algorithm 1 Vessel Tree Extraction

Input: Metric $\mathcal{M}r$ , initial source point s, curve length threshold λ, path score thresholds $(\text{Th}i)1\le i\le k$ , vesselness map ℎ, stencil S.

Output: Minimal path γ, keypoint set $\mathcal{K}$ .

Initialization:

For each point x, set $\mathcal{U}(\mathbf{x})=\mathcal{L}(\mathbf{x})=+\infty$ ; Tag x as Far.

Set $\mathcal{U}(\mathbf{s})=\mathcal{L}(\mathbf{s})$  = 0; Tag s as Trial.

Set $\mathcal{K}\to \{\mathbf{s}\}$ .

Set $i\to 1$ .

Set If Stop = False.

While(If Stop == False)

1: Find $\mathbf{x}\text{min}=(x\text{min},r\text{min})$ , the Trial point which minimizes $\mathcal{U}$ .

2: Tag $\mathbf{x}\text{min}$ as Accepted.

3: for All y such that $\mathbf{x}\text{min}\in S(\mathbf{y})$ and $\mathcal{U}(\mathbf{y})>\mathcal{U}(\mathbf{x}\text{min})$ do

4: Compute $\mathcal{U}\text{new}(\mathbf{y})$ using (19).

5: Compute $\mathcal{L}\text{new}(\mathbf{y})$ using (23).

6: Tag y as Trial.

7: if $\mathcal{U}\text{new}(\mathbf{y})<\mathcal{U}(\mathbf{y})$ , then.

8:  Set $\mathcal{U}(\mathbf{y})\to \mathcal{U}\text{new}(\mathbf{y})$ .

9:  Set $\mathcal{L}(\mathbf{y})\to \mathcal{L}\text{new}(\mathbf{y})$ .

10: end if

11: end for

12: if $\mathcal{L}(\mathbf{x}\text{min})\ge 3\lambda$ , then ⊳ ▹ Path score threshold reduction.

13: if i = k, then

14:  If Stop → True.  ▹ Stopping criterion.

15: else

16:   $i\to i+1$ .

17: end if

18: for all points $\mathbf{z}\in Z$ passed by the minimal paths γ do

19:  Set $\mathcal{U}(\mathbf{z})\to 0$ and $\mathcal{L}(\mathbf{z})\to 0$ .

20:  Tag z as Trial.

21: end for

22: else if $\mathcal{L}(\mathbf{x}\text{min})\ge \lambda$ , then

23: Track the minimal path γ_x from $\mathbf{x}\text{min}$ .

24: if $min\{\mathcal{P}\mathcal{S}(x\text{min},\text{Th}i-1),ℎ(x\text{min})\}>\text{Th}i$ , then

25:   $\mathcal{K}\to \mathcal{K}\cup \{\mathbf{x}\text{min}\}$ .

26:   $\gamma \to \gamma \cup \{\gamma x\}$ .

27:  Set $\mathcal{U}(\mathbf{x}\text{min})\to 0$ and $\mathcal{L}(\mathbf{x}\text{min})\to 0$ .

28:  Tag $\mathbf{x}\text{min}$ as Trial.

29: end if

30: end if

Our algorithm combines three main ingredients:

• The classical keypoint searching scheme. ¹⁴

• State-of-the-art anisotropic fast marching ⁸ with integrated geodesic curve length computation.

• Vessel tree extraction based on the original and new keypoints detection criteria as well as the stopping criteria.

In practice, our keypoints detection on method is embedded within the inner loop of the fast marching algorithm, which it augments with several robust criteria for keypoint detection, adaptation of a path score threshold and termination.

Classical keypoints searching method

The approximated geodesic distance map $\mathcal{U}$ to the currently extracted tree structure is estimated using the fast marching algorithm update scheme: initialization and steps 1–11 of the loop in Algorithm 1. Following the dynamic programming principle, image pixels are tagged as either Trial or Accepted. The Trial point $\mathbf{x}\text{min}$ currently minimizing $\mathcal{U}$ is tagged Accepted (i.e. frozen), and the value $\mathcal{U}(\mathbf{y})$ at neighbouring points y is suitably updated. In addition, line 9, we estimate the geodesic curve length $\mathcal{U}(\mathbf{y})$ and Euclidean curve length $\mathcal{L}(\mathbf{y})$ , and potentially tag y as Trial. The classical keypoints searching method (KPSM) ¹⁴ adds the latest Accepted point $\mathbf{x}\text{min}$ to the set $\mathcal{K}$ of keypoints as soon as $\mathcal{L}(\mathbf{x}\text{min})\ge \lambda$ , where λ is the user chosen curve length threshold. The reason of choosing such a point as keypoint is that among all the fast marching front points with the same minimal action map value, a point q globally maximizing the Euclidean curve length map $\mathcal{L}$ will be located in the centreline of the tubular structure.^14,15,17,21

Keypoints searching scheme based on the path score

A point $\mathbf{x}=(x,r)$ is marked as a new keypoint if its Euclidean curve length satisfies $\mathcal{L}(\mathbf{x})\in [\lambda ,3\lambda )$ , and if additionally it obeys

$min\{\mathcal{P}\mathcal{S}(x,\text{Th}i-1),ℎ(x)\}>\text{Th}i$

(26)

where ℎ is the vesselness map (7).

The key idea behind this detection criterion is that, among all points for which the action map $\mathcal{U}$ is within a given bound, the point which maximizes the Euclidean curve length map $\mathcal{L}$ should be inside the tubular structure. Indeed, the large ratio $\mathcal{L}(\mathbf{x})/\mathcal{U}(\mathbf{x})$ reflects the fact that the geodesic γ, joining x to a previous keypoint, has a small action $\mathcal{P}(\gamma ,\stackrel{·}{\gamma })$ in average. Thus, by construction of the metric, this geodesic must lie on the vessel centerlines and be aligned with the vessel orientations.

Here, we consider the minimal action map $\mathcal{U}$ for a tubular structure tree with the initial point at the root of this tree. For all the points y in the level set $C=\{\mathbf{y};\mathcal{U}(\mathbf{y})=c\}$ , where c > 0 is a constant, there may be many local maximums of $\mathcal{L}$ , some of which are the intersections of the level set C and the tubular structure branch centrelines. Thus, we select the maximal local maximum, which is inside a tubular structure judging by equation (26). In our keypoints searching method, one or several path score thresholds $\text{Th}1\ge \text{Th}2\ge \dots \ge \text{Th}k$ are given as input; generally $k\in \{1,2\}$ , and by convention $\text{Th}0=+\infty$ . If the currently active point (the latest Accepted point) $\mathbf{x}\text{min}$ has an excessive estimated curve length $\mathcal{L}(\mathbf{x}\text{min})\ge 3\lambda$ , see line 12, then the next threshold $\text{Th}i+1$ is activated to discover finer, less visible vessel structures; unless i = k in which case the method ends. If the currently active point $\mathbf{x}\text{min}$ has an appropriate estimated curve length $\lambda \le \mathcal{L}(\mathbf{x}\text{min})<3\lambda$ , then it is considered as a candidate new keypoint (i.e. a potential node of the vessel tree structure), see line 22. We extract the geodesic γ linking $\mathbf{x}\text{min}$ to the already extracted tree structure, and evaluate its relevance using a path score.

The selection test line 24 requires the path score to exceed the current selection threshold $\text{Th}i$ . The selector $\delta (T,·)$ appearing in (24) is applied with $T=\text{Th}i-1$ , so as to push the keypoint detection towards finer and less visible structures, when i > 1, and leave the vicinity of the tree extracted with the previous threshold $\text{Th}i-1$ . Using a hierarchy of successive thresholds, an anisotropic metric, allows in the end to reliably handle a much smaller curve length threshold λ than the classical KPSM, ¹⁴ without leaking outside of the vessel or tubular structure, but staying right in its centreline. In Figure 1, we show the intermediate keypoint searching results using the proposed method. In this figure, the green point is the user-input initial source point, and the red points are the searched keypoints with two path score threshold values $\text{Th}1$ and $\text{Th}2$ . Cyan contours are the boundaries of the tubular structures. In Figure 1(a)–(c), we show the detected one, two and seven keypoints with the corresponding centrelines and contours. In Figure 1(d)–(f), we show the optimal minimal action maps $\mathcal{U}\text{opt}$ of $\mathcal{U}$ with respect to (a)–(c), where $\mathcal{U}\text{opt}$ is defined as

$\mathcal{U}\text{opt}(x):=\underset{r\in [R\text{min},R\text{max}]}{min}\{\mathcal{U}(x,r)\}$

(27)

Figure 1.

Steps of keypoints searching scheme. (a) The first keypoint is found; (b) two keypoints are found; (c) seven keypoints are found. Cyan contours are the boundaries of the tubular structure and red lines are the centrelines. (d)–(f) are the optimal distance maps $\mathcal{U}\text{opt}$ for corresponding to the images of (a)–(c).

Stopping criterion

The keypoint detection process must be stopped after all vessel branches have been explored, but before spurious artefacts appear in the reconstructed vessel tree, which requires an adequate stopping criterion. The proposed algorithm terminates when all provided path score thresholds $(\text{Th}i)i=1k$ have been used, and the value of Euclidean curve length is $\mathcal{L}(\mathbf{x}\text{min})>3\lambda$ , where x is the latest accepted point in the fast marching propagation. In Figure 2, we show the complete keypoints searching result. Green point is the initial source point. Red points are the keypoints with respect to path score threshold $\text{Th}1$ . Blue points are the keypoints associated to the path score threshold $\text{Th}2$ . From Figure 2, we can see that no leakage occurs, even on weak branches, which are reliably extracted by our algorithm.

Figure 2.

Keypoints searching result with two path score threshold values $\text{Th}1$ and $\text{Th}2$ . The green dot indicates the given initial source point. Red dots are the keypoints searched using the path score threshold value $\text{Th}1$ . Blue dots denote the keypoints by the path score threshold value $\text{Th}2$ . Red curves are the centrelines and boundaries.

Figure 3.

Comparison between our algorithm and the classical KPSM. (a) and (b) are the results from the classical KPSM with curve length threshold λ = 26 and λ = 60, respectively. (c) and (d) are the results from our algorithm with λ = 26 and λ = 60, respectively.

Semi-automatic parameter setting

The proposed algorithm requires a few parameters: a curve length threshold λ, and a collection of thresholds $(\text{Th}i)i=1k$ , one initial source point. The curve length threshold λ should be slightly more than twice the largest radius of the vessels to be detected. The path score thresholds collection $(\text{Th}i)i=1k$ is chosen depending on the test case, which are automatically selected as quantiles of the vesselness map distribution on image pixels. In the following experiments, we find that k = 1 or k = 2 are enough to extract the expected vessel tree. Finally, the initial source point can be user provided or automatically selected as the point which maximizes the value of the vesselness map ℎ.

Experiments

Advantages of using a small curve length threshold and a path score judgement

Small curve length thresholds are a priori desirable when extracting vessel trees, since they favour the discovery of small structures and avoid the extraction of inadequate shortcuts linking different tree branches. Unfortunately, the KPSM ¹⁴ suffers from a leaking problem with small curve length threshold: before the main vessel branches are extracted, multiple irrelevant keypoints are detected outside the vessel structure of interest. This problem, which is mainly caused by noise and intensity inhomogeneities, is avoided with our new keypoint detection criterion involving a path score, as illustrated in Figure 3. In Figure 3, with a large curve length threshold λ = 60, the two methods produce similarly inaccurate results: some small branches are missed, and some undesirable shortcuts between different branches are extracted. In Figure 3(a), with a small curve length threshold λ = 26, the classical KPSM ¹⁴ keypoints leak and begin to accumulate outside the structure; in contrast, our path score based keypoints searched method accurately detects the vessel tree and then automatically stops the searching procedure as shown in Figure 3, which can illustrate the advantages of using both path score constraint and small curve length threshold value. Note that in this experiment, we only use one path score threshold value for the keypoints searching scheme.

We illustrate in Figure 4 the results of classical KPSM and our method, with the curve length threshold λ = 12, which empirically is the best for these two methods on this image. For the classical KPSM, we specify a certain number of keypoints to stop the keypoints searching scheme and, for our method, it is stopped automatically. For the proposed method as shown in Figure 4(b), we compute the path score threshold $\text{Th}1$ such that 10% pixels having the highest vesselness values among all pixels are chosen, in other words $\{x\in \overline{\Omega };ℎ(x)\ge \text{Th}1\}$ collects 10% of the image pixels. For $\text{Th}2$ , we use the 12% quantile. On this retinal image, one can see that the classical KPSM suffers from the leakage problem in at least three places due to the gray level inhomogeneities as shown in Figure 4(a). However, for our method, no leakage happens, see Figure 4(b). Note that in the following experiments with retinal image, we only show the centrelines and keypoints for better visualization. In fact, our algorithm also extracts the vessel radii, hence the vessel walls, as illustrated in Figure 2.)

Figure 4.

Comparison of the classical KPSM and the proposed algorithm in real retinal image. (a) Result from the classical KPSM; (b) result of our algorithm. The leakage problem occur at least in three places in (a) from the classical KPSM.

As discussed in ‘Advantages of using a small curve length threshold and a path score judgement’ section, many vessels are missed with a large curve length threshold value, resulting in the extraction of numerous irrelevant shortcuts, as shown in Figure 5. In this experiment, we use the same retinal image and path score thresholds as in Figure 4(b), except that in Figure 5, we use a large curve length threshold λ = 40. One can see that many finer vessels are missed and some short cuts occurs when compared with Figure 4(b).

Figure 5.

Keypoints searching result from our algorithm with curve length threshold λ = 40.

In Figure 6, we demonstrate the impact of the path score threshold to our algorithm in the same retinal image as used in Figure 4(b). In Figure 6, we used only one large path score threshold $\text{Th}1$ by specifying the 6% highest vesselness map value among all the pixels. Less keypoints are found due to the larger path score threshold comparing to Figure 4(b).

Figure 6.

Keypoints searching result from our algorithm with a large path score threshold.

The main drawback of the proposed path score based keypoints searching method is that when handling the tubular structure tree with loops, as shown in Figure 4(b), some small tubular segments will be missed due to the existence of the loops and the fast marching scheme.

Conclusions

In this article, we propose a new keypoints based tubular structure tree extraction method using anisotropic fast marching, and introducing of a path score selection procedure in the keypoints detection criterion. We also show the possibility that the keypoints searching scheme can be automatically stopped by only providing a set of path score thresholds and a Euclidean curve length threshold. These ingredients allow our method to search keypoints separated by small Euclidean curve lengths, leading to better extraction results compared with the classical KPSM. ¹⁴ Numerical experiments illustrate these improvements on two MRA images and one retinal image. The future work is to extend our approach for 3D vessel extraction and to validate it on a large data set.