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

06/08/2017 ∙ by Dario Augusto Borges Oliveira, et al. ∙ 1

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.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 17

page 21

page 24

page 25

page 26

page 27

page 28

page 30

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 -minimum spanning tree (-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.

Figure 1: Modeling a vascular network in a graph using a single direction. On the left side there is a model representing a vascular network with an identified green branch; on the right side its respective graph with green edges indicating the green branch path, and dashed edges representing other path possibilites.

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.

Figure 2: Issues of modeling a vascular network as a graph using a single direction (some edges were hidden to improve understanding). The model on the left side has its green and blue branches correctly modelled in the graph represented on the right side. The yellow path, though, is incorrectly represented - as the red line points out - due to directed graph representation issues.

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.

Figure 3: Definition of graphs using the local vessel direction of the vascular network. The use of vessel-driven sampling models allows the corrected use of directed graph models, overcoming the representation issues showed in Figure 2.

Our contribution is three-fold:

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

  • 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.

  • 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.

Figure 4: Flowchart of the proposed method.

The full procedure involves the following steps:

  1. 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. Sampling model creation. Here a cloud of sample points is created within a conical neighborhood with apex at a given seed point.

  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. 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. New seeds definition. In this step new seed points are created from the vascular network branches found in the previous step.

  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

as its orthonormal basis, whereby the unitary vector

is aligned with the vessel central axis.

Figure 5: Conical sampling model structured by means of a set of spherical calottes.

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 , and

  • the height .

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

Figure 6: The values of and the sampling distance define the sampling circles at a given layer. The radius and the sampling distance define the number of sample points in each circumference .
(1)

Over each calotte , 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 (actually a value close to it).

(2)

Notice that the number of circumferences determined over a calotte increases with the calotte radius. The radius of the circumference over the calotte (see Figure 6) is given by equation 3:

(3)

Notice that for , 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 . Thus, the number of sample points uniformly distributed along the circumference of the calotte is given by equation 4:

(4)

Hence, the actual arc length determined by adjacent sample points along the circumference of the calotte will be given by equation 5:

(5)

Recall that the unitary vector is aligned with the conic central axis, and form an orthonormal basis. It can be easily demonstrated that the coordinate vector of the sample point over the circumference of the calotte can be given by equation 6:

(6)

where and 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 having as its third row the transpose of the unitary vector 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 are arbitrary, provided that is unitary. Hence, we define as:

(7)

where is an arbitrary unitary vector orthogonal to D and .

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 of sample points associated to a seed located at , whose vessel section is oriented in the direction defined by the unitary vector D is given by:

(8)

where is given by equation 7 and 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.

  • 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.

  • 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.

Figure 7: Gaussian mixture for cylinder fitting. The red cylinder fits the vessel at a given point, and the blue volume models the neighbouring area.

In this work, we propose a vesselness measurement consisting of two concentric cylinders and , comprising respectively a vessel section and its outer neighborhood, as shown in figure 7. Let and , be the radii of and , respectively, the common central point and the unitary vector representing the orientation of their common central axis. Let also be such that the volume of is equal to the volume of - . 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

and , where stands for the intensity of the voxel at coordinate

. 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

and - .

The vesselness is then given by:

(9)

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 , the radius and the direction 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:

  1. 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. 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. 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. 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. 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.

Figure 8: Vessel point detections are structured as a directed graph that allows finding vascular network branches, as shown in different colours.

Let be a set of vessel detections comprehended by the conical sampling volume associated to the current seed (recall Section 2.3) whereby represented by its 3D coordinates. A vessel branch that emanates from the current seed brings about a set of vessel detections , laying on a sequence of layers , where is a natural number and , as seen in Figure 5. Let be an arbitrary set of vessel branches. Then, the set of vessel branches composing the vascular network can be defined as the set of vessel branches that best explains the detections in . Thus,

can be found by maximizing the posterior probability of

given the set of vessel detections .

Assuming that the detections are , the objective function is expressed as:

(10)

where is the likelihood of .

In order to reduce the space of , 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 as well, leads to the decomposition:

(11)

for each vessel , represented by an ordered chain. and is the probability that a vessel branch starts and ends at detection . is the probability that is followed by 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 be a directed network with costs associated to every edge . An example of such a network is shown in Fig. 9. It contains two special nodes, the source and the sink . In our proposition, all flow that goes through the graph starts at the node and ends at the node, and each flow represents a vessel candidate . Each vessel detection is represented by two nodes, the beginning node and the end node (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 and connected by an edge with cost associated to the vesselness measurement for the sampled point. Also, a pair of nodes (, ) representing a vessel point is connected to the vessel points laying in the next layer, say (, ), by so called transition edges with cost related to the Euclidean distance between points and . In this model the optimal path corresponds to the minimum cost, which combines vesselness and the distance between nodes along the path.

Figure 9: Proposed graph structure. First layer of the cone consists of node 1, second layer by nodes 2,3,4 and third layer by nodes 5,6. Each pair of begin/end represents a node referenced by the index associated.

Below we detail four of the five different types of edges present in the graphical model. The edges associated to are going to be explained later.

Transition edges. Transition edges connect end nodes over a given layer to beginning nodes of the next layers (orange edges in Fig. 9), with cost and flow , if and belong to , and 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 :

(12)

Detection edges. These edges (plotted in green in Fig. 9) connect the beginning node and end node , with flow if belongs to , and 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 with two nodes and a detection edge in between with negative cost . The higher the likelihood of a detection 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 to all end nodes , with cost and flow . Similarly, exit edges connect the start node with sink , with cost , to ensure that a trajectory ends at a detection with high probability, and flow . The flows are 1 if the trajectory starts/ends at . 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:

(13)

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 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 indirectly: the set of vessels of the optimal vascular network corresponds to the set of paths of minimum flow in . All the related symbols are intrinsically related to the graph shown in Figure 9.

(14)

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 in its linearly relaxed form.

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

    (15)

By now we have fully defined a linear program for finding and consequently . 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 , which finds the global solution for each vessel sequentially and also identifies the network topology.

The first vessel 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 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 , in order to be included in . 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:

(16)

where stands for a weight, and stands for a maximum penalty distance in millimeters.

For all vessel detections which are at a distance or less to any vessel detection found thus far in vessel branch , we compute a corresponding , which will be included to Eq. 14 as a positive penalty. We privilege this way the detection of branches other than , since too similar paths will be penalized through the summation of toll costs. We found out empirically that W 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 to be the starting point of a new vessel branch, as shown in Fig. 10. Formally this is done by allowing 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. 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.

Figure 10: Proposed optimization method example of synthetic images. (a) Initial point marked in orange. First path found is the one with minimum cost (red). In the next step, sources will be added along the path (black arrows). (b,c,d) Paths found iteratively, sources of each path marked by colored arrows. (e) 3D view of the vascular network.

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:

  1. 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. 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.

Figure 11: Post processing impose some anatomical constraints for a vascular network. The red ’X’ shows the branches which would be eliminated following each rule (a) and (b) described in section 2.4.2.

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.

Figure 12: Definition of the seeds for the next iteration.

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.

  while  do
     1. Get a seed from and find the vessel direction
     2. Compute a sampling cloud of points as shown in section 2.2
     3. Build the graph from the sample points
     4. Compute vesselness measurements at sample points given Eq. 9
     while  do
         5. Find vessel with minimum cost
         6. Compute the toll charges (Eq. 16) and new flow conditions.
     end while
     6. Define new seeds
  end while
Algorithm 1 Iterative vessel network tracking

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.

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

  • 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.

Figure 13: Synthetic data segmentation using sinusoidal shaped vascular-like structures. On the left side the input data; on the right side the segmented vascular network in green.

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 () and aperture angle () 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 from equation 16 and 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: Parameter study of in mm, in mm, aperture angle in radians and sampling distance in mm (see sections 2.2 and 2.4). In each graph we see the performance achieved using different pairs of parameters values.

Figure 14 shows the average overlap measure computed over all nine sequences for varying values of . It is possible to notice that on average 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 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 , the aperture angle and the resulting average overlap measure. As the value used gets higher, the optimum 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 and the sampling distance . Clearly, the optimum value of increases with the sampling distance . This is not unexpected, since 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.

Figure 15: Synthetic data segmentation using L-Systems. It is possible to see the reference in yellow and the results achieved in pink.

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.

Figure 16: Vascular segmentation for pulmonary real dataset. We see the results in the axial, coronal and sagittal views and the 3D model generated.

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.

Figure 17: Vascular segmentation for coronary real dataset. We see the results in the axial, coronal and sagittal views and the 3D model generated.

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.

Figure 18: This figure shows what we called the blind effect using an axial view and the 3D model generated. It is possible to visually understand that since our algorithm does not take into account a specific vessel end point to follow a vessel path, it will not necessarily find a desired vessel, but segment all the vessels connected to the given start point.

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.

Figure 19: Vascular segmentation for carothids real dataset. We see the results in the axial, coronal and sagittal views and the 3D model generated.

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.

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
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.
Figure 20: Segmentation results for OPF dataset. We see the results in the axial, coronal and sagittal views and the 3D model generated.

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.

Figure 21: Extracted topology for synthetic data, where each color represents a different branch. In each row a different synthetic model is showed, the first a 2D synthetic model, and the second a 3D model, as referred previously.

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.

References