Automatic tracking of vessel-like structures from a single starting point

TL;DR

This paper introduces an efficient method for fully tracking vascular networks from a single starting point using a spherical sampling approach and network flow optimization, applicable to synthetic and real medical images.

Contribution

It presents a novel iterative tracking algorithm based on a min-cost flow model that handles bifurcations and paths with minimal user interaction.

Findings

Achieved over 98% accuracy on synthetic blood vessel datasets.

Demonstrated robustness to parameter variations.

Successfully segmented vascular and nerve fiber networks in real images.

Abstract

The identification of vascular networks is an important topic in the medical image analysis community. While most methods focus on single vessel tracking, the few solutions that exist for tracking complete vascular networks are usually computationally intensive and require a lot of user interaction. In this paper we present a method to track full vascular networks iteratively using a single starting point. Our approach is based on a cloud of sampling points distributed over concentric spherical layers. We also proposed a vessel model and a metric of how well a sample point fits this model. Then, we implement the network tracking as a min-cost flow problem, and propose a novel optimization scheme to iteratively track the vessel structure by inherently handling bifurcations and paths. The method was tested using both synthetic and real images. On the 9 different data-sets of synthetic blood vessels, we achieved maximum accuracies of more than 98\%. We further use the synthetic data-set to analyse the sensibility of our method to parameter setting, showing the robustness of the proposed algorithm. For real images, we used coronary, carotid and pulmonary data to segment vascular structures and present the visual results. Still for real images, we present numerical and visual results for networks of nerve fibers in the olfactory system. Further visual results also show the potential of our approach for identifying vascular networks topologies. The presented method delivers good results for the several different datasets tested and have potential for segmenting vessel-like structures. Also, the topology information, inherently extracted, can be used for further analysis to computed aided diagnosis and surgical planning. Finally, the method's modular aspect holds potential for problem-oriented adjustments and improvements.

Tables

Table 1. Table 1: OPF database results. The table shows the results obtained for each exam available in the website. The metric is the one made available for the competition, therefore we have a straight comparison with other methods.

	Exam1	Exam3	Exam4	Exam5	Exam6	Exam7	Exam8	Exam9
k-MST Türetken11	–	–	0.865	–	0.898	–	0.722	–
HGD-QMIP Türetken12	–	–	0.923	–	0.911	–	0.722	–
Proposed method	0.800	0.818	0.745		0.833	0.843	0.692	0.327

Equations

$$A=Ls$$

$$M(l)=\lceil l\alpha \rceil +1$$

$$r(m,l)=ls\cos (\frac{\pi }{2}-\frac{m}{l})$$

$$\begin{array}{r}N(m,l)=\lceil \frac{2\pi r(m,l)}{s}\rceil =\lceil 2\pi l\cos (\frac{\pi }{2}-\frac{m}{l})\rceil \text{for}l\in \{\mathrm{1..}L\}\text{and}m\in \{\mathrm{1..}(M(l)-1)\}\end{array}$$

$$a(m,l)=\frac{2\pi }{ml}$$

$$p(n,m,l)=\left[\begin{array}{c}x(n,m,l)\\ y(n,m,l)\\ z(n,m,l)\end{array}\right]=\left[\begin{array}{c}\cos (n\frac{a(m,l)}{r(m,l)})\cos (ls)\\ \sin (n\frac{a(m,l)}{r(m,l)})\cos (ls)\\ \cos (ls)\end{array}\right]$$

$$\mathcal{R}(\text{D})=\left[\begin{array}{ccc}\frac{-D_y}{|V_1|}& \frac{D_x}{|V_1|}& 0\\ \frac{-D_z{D}_x}{|V_2|}& \frac{-D_z{D}_y}{|V_2|}& \frac{D_x{D}_x+D_y{D}_y}{|V_2|}\\ D_x& D_y& D_z\end{array}\right]$$

$$p^{\star }(n,m,l,D,c)=\mathcal{R}(\text{D})p(n,m,l)+c$$

$$W_m=\underset{R_{in},D,c}{\min }\frac{{\sum }_{p\in V_{in}}{G}_{out}(I_p)+{\sum }_{p\in (V_{out}-V_{in})}{G}_{in}(I_p)}{\Vert V_{out}\Vert }$$

$$\mathcal{V}\ast =\underset{\mathcal{V}}{\mathrm{argmax}}P(\mathcal{V}|\mathcal{O})=\underset{\mathcal{V}}{\mathrm{argmax}}\prod _i{P}_{\text{det}}({\mathbf{o}}_i|\mathcal{V})P(\mathcal{V}),$$

$$P(\mathcal{V})$$

$$C_{\text{t}}(i,j)$$

$$\mathcal{V}\ast$$

$$\mathcal{F}\ast$$

$$f_{\text{in}}(i)+f_{\text{det}}(i)={\sum }_j{f}_{\text{t}}(i,j){\sum }_j{f}_{\text{t}}(j,i)=f_{\text{out}}(i)+f_{\text{det}}(i)$$

$$C_{\text{toll}}(i,j)=K_{\text{toll}}\cdot \exp \left(\frac{-\Vert {(\mathbf{o}}_i-{\mathbf{o}}_j)\Vert }{D_{\text{radius}}}\right)$$

Full text

††journal: Computerized Medical Imaging and Graphics

\reserveinserts

18

1 Introduction

The identification of vascular networks is a topic of general interest in medical image analysis and in particular for diagnosis of problems related to the vascular system such as cerebrovascular accidents or thrombosis. While efficient algorithms for single vessel tracking exist, the segmentation of complete vascular networks is still a challenge, mainly due to huge search space involved. Indeed, most solutions are only locally optimal, very computationally intensive or both.

We refer the interested reader to lesagemia2009 for a detailed review of methods to tackle the vascular segmentation problem. Many of these works are based on a propagating structure emanating from a given starting point using different techniques such as level-sets lorigo2001; luboz2005 or minimal paths deschamps2001; wink2000. Other approaches are based on particle filters florin2005; lesage2008, Markov chain processes lacoste2006, statistical methods for tubular structures wang2012 and multiple hypothesis testing frimanmia2010. Recent works turetkenneuro2011; turetkencvpr2012 have shown that global optimization can be reliable as it finds all the vessels of the structure jointly. In turetkenneuro2011, authors find a set of vessel candidate points and then solve the $k$-minimum spanning tree ($k$-MST) problem to find the final structure, while in turetkencvpr2012, the vessel detection model needs to be specifically trained for each type of data. The latter assumes a well behaved distribution of vessel points and a careful training.

In this paper we present a novel Linear Programming (LP) solution for tracking vascular networks through a fast iterative method that tracks each vessel branch independently and is able to handle bifurcations. Our approach does not impose any restrictions to the form of vessels’ offshoots. It relies on two simple assumptions: vessels are nearly cylindrical, having circular or elliptic cross sections, and present a tree like structure. Contrary to most methods proposed thus far, it only requires one starting point to track a full vessel network with arbitrary form.

1.1 Motivation and Graph Modeling

The use of directed graphs to segment tree-like structures (such as vascular networks) is somewhat intuitive. Considering vessel point detections as graph nodes, a vascular network can be structured as a directed graph, which allows the segmentation of its branches through the analysis of graph paths, as depicted in Figure 1.

Not very intuitive is the problem of using directed graphs to represent vascular structures in images with voxels distributed uniformly in a given direction, such as CT scans. Vessels usually change direction continuously and hardly follow the direction of a single axis. Since directed graphs consist of nodes organized in interconnected layers following a certain direction, this behaviour is problematic. The use of a single axis to define the directed graph levels, would lead to the impossibility of segmenting vessels whose direction goes in favor of the axis at some point, but against it later. In other words, the directed graph modeling does not allow a path passing through a sequence of levels N and N+1 to return to any node on level N. Figure 2 depicts the problem. The yellow branch cannot be entirely represented since it crosses some planes more than once.

This paper proposes a way to overcome this difficulty through the use of multiple graphs, each one modelling a portion of the vascular network that conforms the implied sequential spatial arrangement of vessel detection. These local graphs are created iteratively and their nodes are associated to spatial positions expressed in a convenient coordinate system, as illustrated in figure 3. This approach provides, at least locally, a model in which the problem described in figure 2 does not occur, so that the modelling through directed graphs becomes effective. The proposed sampling model is formally defined in section 2.2.

Our contribution is three-fold:

Oriented conical sampling method using vessel point detections to allow tracking of each vessel along its axis.

• 2.

New formulation of the vascular network tracking problem as a min-cost flow problem to track local vascular structures from a single seed point while dealing with bifurcations.

• 3.

Novel optimization scheme that iteratively tracks the full vascular network and guarantee anatomical vascular properties.

Our paper is structured as the following: in section 2 we explain in detail our method, in section 3 we show and discuss the results, and in section 4 we present our conclusions.

2 Materials and Methods

2.1 General Algorithm Description

The proposed method starts at a single user-defined vessel point, called initial seed. From this seed, emerging branches of the vascular network are detected, each one deriving a new seed and possibly new branches. Subsequent iterations are executed until the full vascular network is segmented, as outlined in figure 4.

The full procedure involves the following steps:

Initial seed point definition. In this step a initial seed is defined to be used at the starting point for the whole vascular network to be found. The definition of this point can be done manually by the user or using any automatic procedure. 2. 2.

Sampling model creation. Here a cloud of sample points is created within a conical neighborhood with apex at a given seed point. 3. 3.

Vessel point detection. In this step we test all sample points within the cloud to select those very likely to belong to a vessel. 4. 4.

Vascular network identification. Here we use the vessel point candidates selected in the previous step to track a local vascular network using a proposed network flow analysis approach. 5. 5.

New seeds definition. In this step new seed points are created from the vascular network branches found in the previous step. 6. 6.

Repeat steps 2 to 5 until there is no new seed to be evaluated.

Each of the steps from step 2 on is explained in detail in the following sections.

2.2 Sampling Model Creation

The sampling method is important for the procedure outlined in previous section, as it determines the accuracy as well as the computational efficiency of our proposal. In this section we present our procedure for this matter. It is important to mention that any other sampling procedures that fulfill the requirements of generating a nearly uniform distribution of sample points over a conical neighborhood emerging from a seed point, could be used.

Figure 5 describes the sampling model. The cylinder at the bottom of Fig. 5a represents a cylindrical elementary vessel section centered at the current seed, which is taken as origin of a coordinate system having $(x,y,z)$ as its orthonormal basis, whereby the unitary vector $z$ is aligned with the vessel central axis.

A conical neighborhood with apex in the origin (see fig. 5a) will be explored for vessel points. It is defined by two parameters, namely:

the aperture angle $\alpha$, and

• 2.

the height $A$.

Instead of examining all voxels in the neighborhood, only points over spherical calottes centered at the origin with radii equal to $ls$ are considered (see Figure 6), where $s$ is a user-defined parameter specifying the separation between any adjacent calottes, and $l$ is an integer number between 1 and $L$, the number of calottes given as a further input parameter. Hence, the height $A$ is determined by equation 1:

[TABLE]

Over each calotte $l$, $M(l)$ circumferences are determined as shown in Figure 5b and defined in equation 2, so that the minimum distance between any two points lying on consecutive circumferences over the same calotte is at most equal to $s$ (actually a value close to it).

[TABLE]

Notice that the number of circumferences determined over a calotte increases with the calotte radius. The radius $r(m,l)$ of the $m^{th}$ circumference over the $l^{th}$ calotte (see Figure 6) is given by equation 3:

[TABLE]

Notice that for $m=0$, the circumference boils down to a single point lying on the cone central axis.

Sample points are uniformly distributed along each circumference, whereby the arc determined by two adjacent sample points is at most equal to $s$. Thus, the number $N(m,l)$ of sample points uniformly distributed along the $m^{th}$ circumference of the $l^{th}$ calotte is given by equation 4:

[TABLE]

Hence, the actual arc length $a(m,l)$ determined by adjacent sample points along the $m^{th}$ circumference of the $l^{th}$ calotte will be given by equation 5:

[TABLE]

Recall that the unitary vector $\mathbb{z}$ is aligned with the conic central axis, and $\mathbb{x,y,z}$ form an orthonormal basis. It can be easily demonstrated that the coordinate vector $\mathbb{p}(n,m,l)$ of the $n^{th}$ sample point over the $m^{th}$ circumference of the $l^{th}$ calotte can be given by equation 6:

[TABLE]

where $r(m,l)$ and $a(m,l)$ are given by equations 3 and 5, respectively.

Once the basic spatial distribution of sample points has been defined, we need a geometric transformation that places the point cloud properly on the neighborhood of the seed being considered. Basically, the conic points cloud must be rotated so that the cone central axis, aligns with vessel’s central axis at the seed point position, and its origin shifted to the seed point position.

This is done by first multiplying the coordinates given in equation 6 by a 3x3 rotation matrix $\mathcal{R}(\textbf{D})$ having as its third row the transpose of the unitary vector $\textbf{D}=[Dx,Dy,Dz]^{\intercal}$ aligned with the vessel central axis at the seed, and then, by adding the result to the seed coordinate vector. The first and second rows of $\mathcal{R}(\textbf{D})$ are arbitrary, provided that $\mathcal{R}(\textbf{D})$ is unitary. Hence, we define $\mathcal{R}(\textbf{D})$ as:

[TABLE]

where $V_{1}$ is an arbitrary unitary vector orthogonal to D and $V_{2}=\textbf{D}\times V_{1}$.

This procedure creates a canonical cloud of points, at each seed point, a cloud of points in this shape is created by translating and rotating this model accordingly fitting the seed point position and direction.

Finally, the coordinate vector $\mathbb{p}^{\star}$ of sample points associated to a seed located at $\mathbb{c}$, whose vessel section is oriented in the direction defined by the unitary vector D is given by:

[TABLE]

where $\mathcal{R}(\textbf{D})$ is given by equation 7 and $\mathbb{p}(n,m,l)$ is defined in equation 6.

The spatial distribution of sample points in the proposed model (figure 5) presents some interesting characteristics:

The probability that a voxel belongs to the vessel network diminishes for layers with larger radii.

• 2.

The area of a layer grows with its distance to the seed. This complies with the expected increase of vessel network spreading as one moves away from the seed.

• 3.

The cone opening angle and height are related to the potential of vessel to change direction suddenly.

Once the sample point coordinates have been defined, each of them is evaluated as a vessel point candidate, using the vesselness metric proposed in the next section.

2.3 Vessel Point Detection

This section describes the procedure to assess how well sample points qualify as vessels. It involves two sequential steps: (a) computation of a metric, hereafter called vesselness, to assess how well a sample point qualifies as part of a vessel; (b) selection of points based on local characteristics, including vesselness. These steps are described in detail in the following subsections.

2.3.1 Vesselness Computation

Recalling section 1.1, sample points are associated to vessel hypotheses according to a determinate vessel model. Different vessel models can be found in the literature, such as elliptical cross-sections models florin2005, spheres rossignac2007 and template models frimanmia2010; worz2007.

In this work, we propose a vesselness measurement consisting of two concentric cylinders $V_{in}$ and $V_{out}$, comprising respectively a vessel section and its outer neighborhood, as shown in figure 7. Let $R_{in}$ and $R_{out}$, be the radii of $V_{in}$ and $V_{out}$, respectively, $\mathbb{c}$ the common central point and $D$ the unitary vector representing the orientation of their common central axis. Let also $R_{out}$ be such that the volume of $V_{in}$ is equal to the volume of $V_{out}$ - $V_{in}$. The cylinder height is defined as two times the inner radius, and any value around this does not impact much the final outcome.

It is further assumed that the voxels’ intensities inside and outside the vessel can be appropriately represented by two distinct Gaussian distributions, denoted respectively $G_{in}(I_{p})$ and $G_{out}(I_{p})$, where $I_{p}$ stands for the intensity of the voxel at coordinate $p$. The mean and standard deviation of each Gaussian distribution are computed as the mean value and standard deviation of the voxels inside each corresponding volume $V_{in}$ and $V_{out}$ - $V_{in}$.

The vesselness is then given by:

[TABLE]

The optimization implied in the vesselness computation can be performed by different stochastic methods. In this work we used Differential Evolution diffev for this purpose.

As byproducts of vesselness computation we obtain the optimum position $\mathbb{c}$ , the radius $R_{in}$ and the direction $D$ of the hypothesized vessel section. In particular, the optimization procedure starts at a initial location for the vessel point and looks for the most accurate position in its neighborhood.

2.3.2 Vessel Point Candidate Selection

Vessel point candidates are selected from sample points in two steps:

Vessel candidates with intensity mean value differing greatly from the mean intensity value observed at the cone apex point are discarded and not even submitted to the next step, to save processing time. 2. 2.

Points selected in step 1 are submitted to vesselness computation. The optimal cylinder found by the implied optimization process is validated according to the following rules:

1. (a)

The mean intensity value of the inner cylinder must be higher than the mean intensity value of the outer volume (see figure 7 for details). This rule models the assumption that vessels appear as roughly bright cylindrical structures in CT images. This rule is applied during the cylinder fit optimization process. 2. (b)

The vesselness of a selected sample point must be higher than a fraction of the vesselness at the seed point. The variable controlling this proportion is user defined and works as a sensibility parameter, which allows for finding weaker vessels (at the cost of adding noise) or just the vascular branches with stronger vesselness information. 3. (c)

Vessel points with too small cylinder radius are discarded, assuming these cases are noise. The threshold for the radius is configurable and set as the image spacing by default (vessels with radius lower than image spacing are hardly detectable).

Having detected the vessel points among sampled points, we structure them as a directed graph as detailed in section 2.2 and analyze it for tracking the vascular network, as explained in the next section.

2.4 Vascular Network Identification

The vascular network tracking for a given collection of vessel point detections is composed of two steps: vascular network tracking and vascular network validation.

2.4.1 Vascular Network Tracking

The idea for detecting vascular networks is to build and analyze local graphs for which the nodes represent vessel detections (as found in Section 2.3), using the sample points coordinates and their respective vesselness value. These nodes are connected to neighboring nodes both from adjacent layers (connected by transition edges) and from the same layer (connected by toll edges). Each edge represents the relation between two nodes using a determinate cost, as explained further in this section.

Thereby, the matching problem can be approached as a minimum-cost network flow problem: finding the optimal set of vessels can be solved by sending flow through the graph so as to minimize its total cost.

Let $\mathcal{O}=\{{\bf o}_{i}\}$ be a set of vessel detections comprehended by the conical sampling volume associated to the current seed (recall Section 2.3) whereby ${\bf o}_{i}=(x_{i},y_{i},z_{i})$ represented by its 3D coordinates. A vessel branch that emanates from the current seed brings about a set of vessel detections $V_{k}=\{{{\bf o}_{k_{1}}},{{\bf o}_{k_{2}}},\cdots,{{\bf o}_{k_{h}}}\}$, laying on a sequence of layers $k_{1},k_{2},\cdots k_{h}$, where $k_{i}$ is a natural number and $|k_{i+1}-k_{i}|\leq 1\quad\textrm{for all}\quad i=1,...,h-1$, as seen in Figure 5. Let $\mathcal{V}=\{V_{k}\}$ be an arbitrary set of vessel branches. Then, the set of vessel branches $\mathcal{V}*=\{V_{k}\}$ composing the vascular network can be defined as the set of vessel branches that best explains the detections in $\mathcal{O}$. Thus, $\mathcal{V}*$ can be found by maximizing the posterior probability of $\mathcal{V}$ given the set of vessel detections $\mathcal{O}$.

Assuming that the detections are $i.i.d.$, the objective function is expressed as:

[TABLE]

where $P_{\textrm{det}}({\bf o}_{i}|\mathcal{V})$ is the likelihood of $k_{1},k_{2},\cdots k_{h}$.

In order to reduce the space of $\mathcal{V}$, we make the assumption that vessel branches do not overlap, i.e., a detection cannot belong to more than one vessel). Assuming that vessel branches are $i.i.d.$ as well, leads to the decomposition:

[TABLE]

for each vessel $\{V_{k}\}$, represented by an ordered chain. $P_{\textrm{in}}({\bf o}_{i})$ and $P_{\textrm{out}}({\bf o}_{i})$ is the probability that a vessel branch starts and ends at detection ${\bf o}_{i}$. $P_{\textrm{t}}({\bf o}_{i}|{\bf o}_{j})$ is the probability that ${\bf o}_{j}$ is followed by ${\bf o}_{i}$ in the vessel branch.

To solve the objective function implied in equation 10, we linearize it using the network flow paradigm, where probabilities of connected nodes are represented by means of costs and flows. According to this paradigm, maximizing the probabilities of nodes in a directed network is equivalent to finding the flows that minimize the total flow over the network, given a set of predefined costs.

Let $G=(\mathcal{N},\mathcal{E})$ be a directed network with costs associated to every edge $e_{i}\in\mathcal{E}$. An example of such a network is shown in Fig. 9. It contains two special nodes, the source $S$ and the sink $T$. In our proposition, all flow that goes through the graph starts at the $S$ node and ends at the $T$ node, and each flow represents a vessel candidate $V_{k}$. Each vessel detection ${\bf o}_{i}$ is represented by two nodes, the beginning node $b_{i}\in\mathcal{N}$ and the end node $e_{i}\in\mathcal{N}$ (see Fig. 9), and a detection edge connecting them.

More specifically, each vessel point detection is represented in the graph by a pair of nodes $b_{i}$ and $e_{i}$ connected by an edge with cost $C_{det}(i)$ associated to the vesselness measurement for the sampled point. Also, a pair of nodes ($e_{i}$, $b_{i}$) representing a vessel point is connected to the vessel points laying in the next layer, say ($e_{j}$, $b_{j}$), by so called transition edges with cost $C_{t}(i,j)$ related to the Euclidean distance between points $i$ and $j$. In this model the optimal path corresponds to the minimum cost, which combines vesselness and the distance between nodes along the path.

Below we detail four of the five different types of edges present in the graphical model. The edges associated to $C_{toll}(i,j)$ are going to be explained later.

Transition edges. Transition edges connect end nodes $e_{i}$ over a given layer to beginning nodes $b_{j}$ of the next layers (orange edges in Fig. 9), with cost $C_{\textrm{t}}(i,j)$ and flow $f_{\textrm{t}}(i,j)=1$, if ${\bf o}_{i}$ and ${\bf o}_{j}$ belong to $V_{k}$, and [math] otherwise. The cost associated to transition edges relates to the distance of adjacent vessel detections. Since nearby points are more likely to belong to the same vessel, we define the costs to be a decreasing function of the distance between neighboring vessel detections, assuming a distance $D_{\textrm{max}}$:

[TABLE]

Detection edges. These edges (plotted in green in Fig. 9) connect the beginning node $b_{i}$ and end node $e_{i}$, with flow $f_{\textrm{det}}(i)=1$ if ${\bf o}_{i}$ belongs to $V_{k}$, and [math] otherwise. It is worth noting that, if all edge costs are positive, the solution of the minimum-cost problem is the trivial null flow. The trick in our model is to represent each vessel detection ${\bf o}_{i}$ with two nodes and a detection edge in between with negative cost $C_{\textrm{det}}(i)=\log\left(1-P_{\textrm{det}}({\bf o}_{i})\right)$. The higher the likelihood of a detection $P_{\textrm{det}}({\bf o}_{i})$ the more negative the cost of the detection edge. Hence true detections are likely to be in the path of the flow so as to minimize the total cost.

Entrance and exit edges. Entry edges (purple in Fig. 9) connect the source $S$ to all end nodes $e_{i}$, with cost $C_{\textrm{in}}(i)=1$ and flow $f_{\textrm{in}}(i)$. Similarly, exit edges connect the start node $b_{i}$ with sink $T$, with cost $C_{\textrm{out}}(i)=1+C_{\textrm{det}}(i)$, to ensure that a trajectory ends at a detection with high probability, and flow $f_{\textrm{out}}(i)$. The flows are 1 if the trajectory $V_{k}$ starts/ends at ${\bf o}_{i}$. It is interesting to notice that these constants were arbitrarily set, since the flow optimization would find the same paths independently on their specific values.

Now equation 10 can be rewritten in terms of costs and flows, as usual in the network flow paradigm. First, we apply negative log-likelihoods to the equation 10 to derive equation 13:

[TABLE]

then, denoting each of the probability terms in equations 13 and 11 as a product of costs and flows, we derive equation 14, which computes the set of flows $\mathcal{F}*$ that minimizes the global network flow. Since maximizing the probabilities in the graph is equivalent to finding the nodes that minimize the total flow over the network, we find $\mathcal{V}*$ indirectly: the set of $V_{k}$ vessels of the optimal vascular network $\mathcal{V}*$ corresponds to the set of paths of minimum flow in $\mathcal{F}*$. All the related symbols are intrinsically related to the graph shown in Figure 9.

[TABLE]

subject to the following constraints:

Edge capacities: we assume that each detection belongs either to one vessel or to none, i.e., flows can assume 0 or 1 values, or $0\leq f(i)\leq 1$ in its linearly relaxed form.

• 2.

Flow conservation at the nodes: the entering flow of a node equals its exiting flow.

[TABLE]

By now we have fully defined a linear program for finding $\mathcal{F}*$ and consequently $\mathcal{V}*$. This linear problem can be solved using many of the Linear Programming solvers proposed thus far by the optimization community, such as Simplex or k-shortest paths LP. However, equation 14 is very computationally intensive. We propose a fast iterative procedure to compute $\mathcal{V}*$, which finds the global solution for each vessel sequentially and also identifies the network topology.

The first vessel $V_{1}$ is found by solving Eq. 14, allowing entering flow only at the seed point (the cone apex) and setting the maximum flow going out of node $T$ to be 1. A further condition is imposed in order to avoid multiple paths representing a single vessel. Specifically, the vessel must be distant enough from $V_{1}$, in order to be included in $\mathcal{V}*$. To represent such condition in the graph structure a new type of edge is proposed, which connects vessel detections. Edges of this type are associated to a penalty cost, hereby called “toll”. They are represented by the thick black edges in Fig. 9. This ”toll” cost is defined as:

[TABLE]

where $K_{\textrm{toll}}$ stands for a weight, and ${D_{\textrm{radius}}}$ stands for a maximum penalty distance in millimeters.

For all vessel detections ${\bf o}_{j}$ which are at a distance ${D_{\textrm{radius}}}$ or less to any vessel detection found thus far in vessel branch $V_{1}$, we compute a corresponding $C_{\textrm{toll}}$, which will be included to Eq. 14 as a positive penalty. We privilege this way the detection of branches other than $V_{1}$, since too similar paths will be penalized through the summation of toll costs. We found out empirically that W $K_{\textrm{toll}}=5$ is a proper choice.

For detecting bifurcations and their corresponding branches a simple solution is adopted. We update the flow conditions, allowing any point of the detected $V_{k}$ to be the starting point of a new vessel branch, as shown in Fig. 10. Formally this is done by allowing $f_{in}$ to range between 0,1 for such vessel points, and by setting it to zero for all the other vessel detections. This solution handles bifurcations, and new branches are found iteratively. $f_{in}$ and toll costs are updated accordingly, until the cost of finding a new vessel branch becomes positive, i.e., until the positive costs outweigh the negative costs in the minimum flow corresponding to the new path.

2.4.2 Vascular Network Validation

A final vascular network validation procedure is performed to ensure that the final outcome (figure 11) conforms to vascular anatomy, specifically:

Vessel branches cannot be too close to each other along their length, in which case, the detected branches are very likely to represent the same vessel. This condition is implemented by imposing a minimum distance between branches found at each iteration. If this value is lower than a threshold, one of them is discarded. 2. 2.

In order to avoid loops, vessel branches are not permitted to reconnect to the so far segmented vascular network. This is implemented by computing the minimum distance between the branches found and the already segmented network. If this value is lower than a threshold, the branch forming a loop is discarded.

It is important to mention that those threshold values are relative to the estimated radius. In this way, the absolute threshold values change dynamically during the segmentation process, so as to avoid misconnections, while still keeping the ability to find small vessels.

2.5 Next seeds definition

The final step consists in defining new seeds from the detected branches to be used in the iterative method. A simple procedure is defined to find the best seed of each branch: we seek for the vessel point detection with the best vesselness value in the ending part (30% last points) of each branch found, and define it as a new seed, as depicted in figure 12.

Each new seed derives a new vascular network tracking process at its location, until no more seeds are found. An overview of the proposed method is shown in Algorithm 1.

3 Results and Discussion

3.1 Results

Experiments have been conducted to evaluate the proposed method. The scarcity of public data sets with reference information has prevented a thorough and objective experimental comparison among (semi) automatic methods for vascular segmentation proposed thus far. In fact, to our best knowledge, for this reason none of these works managed to perform an objective performance comparison with alternative approaches for vessel network segmentation.

Some authors tried to circumvent this difficulty either by relying on synthetic data (e.g. frimanmia2010; worz2007; WongC07; schap07) or by limiting themselves to visual evaluations (e.g. worz2007; WongC07; Manniesing07; schap07). The analysis reported hereafter follows one or the other strategy depending on the characteristics of the test data.

In some experiments we used a data set created for conference challenges, which contain references. In those cases the following accuracy metrics were used (detailed in cls2009):

The overlap Dice similarity index for 3D volumes.

• 2.

Root mean squared (RMS) distance between reference and segmented 3D surfaces.

• 3.

Hausdorff distance between reference and segmented 3D surfaces.

Concerning the processing time, it is important to mention that the iterative nature of the method can derive very variable amount of time to run depending on the complexity of the segmented network. Just to give an idea of processing time, an iteration with 5000 points can take between 5 and 10 minutes approximately. Depending on the complexity of a network, the number of cones needed in the iterative process to segment the whole network vary a lot.

It is worth mentioning that the optimization procedure implied in the vesselness computation was done using the differential evolution diffev method.

In the following we report the experiments conducted to evaluate the proposed method.

3.1.1 Synthetic data 1

The first experiment was based on a data set available in Macedo10. It consists of a simple planar vascular like structure modeled as sinusoidal shapes with bifurcations. The data set is affected by Gaussian noise artificially added to the raw data.The result obtained with this data-set is shown in figure 13. Clearly, the vascular network was fully segmented. Even though no reference is available, a visual inspection indicates that the method was successful in this case.

3.1.2 Synthetic data 2

In these experiments a more complex synthetic data-set provided by MiguelGalarreta2012mastery; Miguel2013Vessels has been used. The data-set of a network of three-dimensional synthetic blood vessels generated by so called stochastic Lindenmayer systems (L-systems) relying on grammars that represent blood vessel architectures so as to produce nearly realistic vascular networks. Nine different sequences were used in our experiments, each of them with different morphological characteristics. Our method achieved maximum accuracies ranging from 98% to 100% according to the overlap Dice similarity index described in cls2009.

Figure 14 shows the results of experiments to assess the sensibility of our method to its parameters.Two of them were the sampling distance ($s$) and aperture angle ($\alpha$) of the conical sampling cloud described in section 2.2. Two other parameters related to the Linear Programming optimization described in section 2.4 have also been tested, corresponding to $D_{\textrm{radius}}$ from equation 16 and $D_{\textrm{max}}$ from equation 12. We used the overlap Dice similarity index to compare the volume found during the segmentation process and the one of the original image.

Figure 14 shows the average overlap measure computed over all nine sequences for varying values of $D_{\textrm{radius}}$. It is possible to notice that on average $D_{\textrm{radius}}$ must be kept small (from 3 to 5), otherwise bifurcations are not detected because toll costs are too high, causing the drop in accuracy. In figure 14 we make a similar study for the aperture angle $\alpha$ parameter. Interestingly, this parameters is related to anatomical characteristics of the vascular network (or part of it) to be segmented, which might spread more or less widely. This fine tuning is necessary because if the angle is too small bifurcations cannot be properly followed because they do not fall inside the cone; on the other hand, if the angle is too large, surrounding structures are included in the conical volume and spurious paths can be found.

Figure 14 presents the relationship between $D_{\textrm{radius}}$, the aperture angle $\alpha$ and the resulting average overlap measure. As the $\alpha$ value used gets higher, the optimum $D_{\textrm{radius}}$ also increases, so as to avoid that too many false bifurcations are found. Nonetheless, the results are very stable with a wide range of parameters configuration. Finally, the plot of figure 14 shows the relation between $D_{\textrm{max}}$ and the sampling distance $s$. Clearly, the optimum value of $D_{\textrm{max}}$ increases with the sampling distance $s$. This is not unexpected, since $C_{t}$ is related to the distance of the corresponding vessel points, which are clearly affected by the sampling distance. A sampling distance of 1-1.5 mm gave in our experiments good and stable results.

Figure 15 shows the outcome and provides a visual evaluation of its accuracy. Since the paths are composed by sampling points, the pink part was artificially filled, and this is the reason of the split observed in the bottom model.

3.1.3 Pulmonary Data

This dataset was provided by the Extraction of Airways from CT 2009 (EXACT09) challenge pulmonary and contains pulmonary vessels inside the lungs.

The goal of this challenge was to compare the results of various algorithms to extract the airway tree from CT scans using a common dataset and performance evaluation procedure. The challenge provides training and testing datasets, and there is also a reference data available for the pulmonary airways but not for the pulmonary vascular network, so, our evaluation in this case was qualitative.

Figure 16 shows the results produced by our method on this data set. It can be seen that most of the vascular network, including bifurcations, were found. Even though this kind of evaluation is not numerical, the results show the potential of the method, specially if one takes into consideration that a single seed point was used.

3.1.4 Coronary Data

In this experiment we tested our method for cardiac vessels segmentation. To evaluate our algorithm for such application, we used the coronary dataset provided for the MICCAI ”3D Segmentation in the Clinic: A Grand Challenge II” coronary. This database contains thirty-two cardiac CT datasets with reference data available for the four main coronary vessels. It is important to notice that the reference is composed by four single vessel segmentation, and therefore does not test the full potential of our method since they do not form a full vascular network. Nonetheless, it is a very interesting dataset and the available references can be used as guidance for visual assessment.

Figure 17 shows that the reference single vessel is among the branches segmented (in green) using our method, which segmented other two extra branches as well.

Considering the single vessel reference available for this dataset, an interesting effect is noticed. Since we use a single start point for segmenting the vascular network, a so-called ”blind effect” is observed. Even though the segmentation process segments the network and includes the reference single vessel, it also segments other branches, and therefore would be badly evaluated by the competition evaluation tool. If we restrict the algorithm to find a single path, then it is not possible to ensure a-priori which of the three detected branches would be chosen. Figure 18 depicts the effect.

3.1.5 Carotids Data

Carotids are the vessels that irrigate the brain, and therefore are very important for many medical conditions such as the cerebrovascular accidents, commonly known as strokes. We also tested our algorithm in segmenting the carotid vessels, using the dataset carotid created for the 3rd MICCAI Workshop in the series ”3D Segmentation in the Clinic: a Grand Challenge III”.

The reference data available for this challenge concerns the identification of the carotid bifurcation. Even though the carotid vascular network reference is not available, it is possible to assess visually the outcome, since these vessels are morphologically simple, usually having a single main bifurcation. Figure 19 shows that both vessels and the bifurcation are found.

3.1.6 Olfactory Projection Fibers (OPF)

The OPF dataset is actually not from a vascular system, but from a nervous fibers network of the olfactory system. This, of course, hinders the detection of the vascular network (the use of a more appropriate model for vesselness would be advisable), albeit good results were achieved. The dataset is available at the DIADEM (short for Digital Reconstruction of Axonal and Dendritic Morphology) challenge website diadem, which was a competition for evaluating algorithmic methods for automated neuronal tracing.

OPF are network-like structures, but the inner part, corresponding to the lumen in vessels, are not exactly homogeneous and therefore the Gaussian mixture model proposed for vesselness evaluation, delivers low likelihood values. Still, it is possible to take advantage of the structure and segment at least part of the network.

We used this non-vascular network dataset in our experiments mainly because of its reference data, which allowed us to evaluate our method quantitatively. For this, we used the evaluation tool proposed for the DIADEM challenge.

The obtained results are summarized in Table 1, which allows for a comparison of our method with other works. In the table we compare with turetkenneuro2011; turetkencvpr2012, but the interested reader can refer to the challenge website to check other methods results. The results obtained are comparable with the ones obtained by methods designed for segmenting nervous fibers networks. This is encouraging and we believe that a more suited metric for nervous fibers (instead of the proposed vesselness) would improve substantially the results, even though this is not in the scope of this paper. Figure 20 gives a visual feedback of the results obtained in one of the available datasets.

3.2 Topological Description

Our method delivers not only the vascular network segmentation but its topology as well. It is possible to detect vascular branches and identify connecting points, which are assumed to be bifurcations (even though the exact point of a bifurcation is not very accurately determined).

We have not evaluated nummerically the results for topology with the datasets used, but some visual feedback is given in figure 21. It is possible to visualize branches in different colours and therefore inspect visually the vascular network topology. The results are coherent with visual inspection.

4 Conclusions

This paper presents a method to segment full vascular networks through an iterative procedure that starts at a single vessel point provided by the user. A cloud of samples is defined within conical volume having at its apex the current vessel point. The decision whether or not a sample is a vessel point is based on a metric of how well the sample fits a vessel model. The vascular network is modeled by a directed graph. The final vascular network connecting the detected vessel points is obtained by solving a a flow problem using linear programming. This procedure finds branches of vessels iteratively until no new branch is found.

We tested our method using several datasets. An extensive performance comparison of our method with alternative approaches could not be done in the present work, due the lack of a public data set for full vascular network segmentation. All publicly available data sets just contain single vessel branches, or non vascular network structures, such as a nervous fiber network.

Synthetic data was correctly segmented even for datasets with many bifurcations and vessels with different diameters. The method also produced visually coherent vascular networks for different organs. A dataset of pulmonary CT images containing vascular networks with high capillarity, was reasonably segmented, specially if one considers that a single seed point was used. Datasets of carotid and coronary CT scans stressed the algorithm due to the very nature of their data with many different structures surrounding the vessels, some of them sharing vascular textures and densities. Nevertheless, the results obtained were visually coherent with an specialist expectancy. Finally, we tested our algorithm using an OPF dataset, which, despite not being a vascular dataset but a nervous fibers network, provided a reference with which we could evaluate numerically our method. Even though the vesselness measurement is not very suited for evaluating nervous fibers, good results were obtained and the olfactory fibers network was segmented properly. The results obtained show the potential of the proposed method to segment and extract the topology of different vascular networks.

Acknowledgements

The authors acknowledge Faperj, CNPq and CAPES for funding this research.