## 1.

## Introduction

Vessel centerline has been used in computing edge gradients and searching for border positions,^{1} to derive video-densitometric profiles,^{2} for measurement of vessel diameter,^{3} for calculation of lesion symmetry,^{4} and as the basis for three-dimensional (3-D) reconstruction of vessel segments or of the entire arterial tree.^{3}^{,}^{5}

Many techniques in the literature propose to detect vessel centerline by thinning algorithm.^{6}^{–}^{8} These algorithms are notoriously sensitive to noise and can result in ambiguous results at bifurcations. The results will lead to inaccuracies in width, vessel direction, and touristy measurements. Several articles have introduced tracking-based techniques for extraction of the vessel centerlines.^{9}^{–}^{12} In these approaches, vessel centerlines are extracted segment-by-segment while the vessel boundary associated with each segment is delineated. These algorithms require user intervention to determine the seed point. The axial tangent of the seed is then generated. This tangent helps to define a cross-sectional plane (normal plane) by shifting forward the current one along the tangential direction, on which the vessel is segmented and the next axial point is calculated as the centroid of the planar segmentation. This process keeps iterating until a termination criterion is met or the user preempts the tracing.

In most images, however, noise and low-contrast pose significant challenges to the extraction of vessel centerline. For example, the central intensity of some vessels differs from the background by as little as 3 gray levels, yet the background noise standard deviation is 2.3 gray levels. Narrow vessels (in diameter $<3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{voxels}$) often have the lowest contrast. In these cases, a disconnected vessel centerline may be obtained and, more importantly, the aforementioned methods^{9}^{–}^{12} may only produce an incomplete extraction of the vessel centerline.

Some authors attempt to overcome these problems by improving vessel tracking.^{13}^{–}^{17} For example, Manniesing et al.^{15} presented an approach for vessel centerline tracking based on surface evolution. In this work, the main idea is to guide the evolution of the surface by analyzing its skeleton topology during evolution. Friman et al.^{16} proposed a multiple hypothesis template tracking approach for the segmentation of small 3-D vessel structures. This work also contributes a novel mathematical vessel template model, with which an accurate vessel centerline extraction is obtained. Wong and Chung^{17} proposed a probabilistic framework for vessel centerline tracing. In order to obtain a satisfactory vessel centerline extraction with the above methods, user guidance on centerline tracing is necessary. Otherwise, these methods will unavoidably lead to the extraction of disconnected centerline or a failure to extract complete centerline.

The main purpose of the present paper is to design an algorithm, which is robust and accurate and requires no human interaction for extraction of vessel centerline. The vessel centerline extracted by the proposed method is continuous and complete. Our methods consisted of the main three steps: (1) enhancing small vessel by processing with a set of line detection filters, corresponding to the 13 orientations; (2) extracting vessel centerline segment candidates by a thinning algorithm; and (3) grouping and selecting vessel centerline segmentations by a global optimization algorithm. We validate the method quantitatively on a number of synthetic data sets and real computed tomography (CT) angiography data sets.

## 2.

## Methods

Because of noise and low-contrast, disconnections occur in vessel centerline detection. In order to obtain complete vessel centerline, we will focus on these centerline disconnections. The tracing methods may unavoidably lead to the extraction of disconnected centerlines if no user guidance is provided, we will locate disconnections firstly, and then show how the global optimization algorithm may be used in grouping them to select optimal connections.

## 2.1.

### Extracting Vessel Centerline Segment Candidates

In this section, we focus on extracting vessel centerline segment candidates for subsequent processes. Starting by enhancing small vessel, we then employ histogram-based threshold method to extract vessel region. Finally, we introduce a thinning algorithm to obtain vessel centerline segment candidates.

## 2.1.1.

#### Preprocessing for small vessel enhancement

In the step of vessel enhancement, the vessels are enhanced and the noises are suppressed. Multiple scale filtering based on the analysis of the Hessian matrix and its eigenvalues have received a large amount of attention in the vessel enhancement.^{18} In addition, Gabor filters^{19} are applied to images at different scales in order to account for vessels of different widths. However, some small and weak vessels still cannot be detected because they are not successfully enhanced at any of the multiple scales. Further, how to extract accurate vessel width has not been sufficiently discussed in these multiscale schemes.

To improve the discrimination between the small and weak vessels and the background noise, the image is processed, corresponding to the 13 orientations. Figure 1(a) shows six kinds of major directions in 3-D images on a cubic grid; the 13 orientations are denoted by SD, ED, ND, WD, ES, SW, S- (South direction), D- (Down direction), E- (East direction), USW, USE, UNE, and UNW [see Fig. 1(b)]. The set of convolution kernels used in this operation is shown in Fig. 2. For each voxel, the highest filter response is kept and added to the image. This resulting image is the base of all subsequent operations used for the detection of vessel centerlines.

## 2.1.2.

#### Extraction of vessel region

In order to obtain the vessel region, we employ the Parzen window estimate^{20} for automatic threshold selection. The Parzen window integrates image histogram information with the explicit spatial information about voxels of different gray levels to estimate the threshold value.

Assume $f(x,y,z)$ be the gray value of the voxel located at the point ($x,y,z$). In an image $\{f(\mathrm{x},\mathrm{y},\mathrm{z})\}$ of size $N=m\times n\times k$ with $L$ gray levels, let ${w}_{i}=\{(x,y,z)|f(x,y,z)=i\}$, $i=1,2,\dots ,L$, and ${C}_{i}$ denote the number of points in ${w}_{i}$. According to the Parzen window technique, for the point space ${w}_{i}$, its probability density function $p(x,y,z,{w}_{i})$ can be approximated as the following,

## Eq. (1)

$$p(x,y,z,{w}_{i})=p({w}_{i})p(x,y,z|{w}_{i})\phantom{\rule{0ex}{0ex}}=\frac{{C}_{i}}{N}\xb7\frac{1}{{C}_{i}}\sum _{j=1}^{{C}_{i}}\frac{1}{V{c}_{i}}\varphi (x,y,z;{x}_{j},{y}_{j},{z}_{j};{h}_{{C}_{i}}^{3}),$$## Eq. (2)

$$\varphi (x,y,z;{x}_{j},{y}_{j},{z}_{j};{h}_{{C}_{i}}^{3})\phantom{\rule{0ex}{0ex}}=\frac{1}{3\pi}\mathrm{exp}[-\frac{{(x-{x}_{j})}^{2}+{(y-{y}_{j})}^{2}+{(z-{z}_{j})}^{2}}{3{h}_{{C}_{i}}^{3}}].$$Let $t$ be a threshold which separates all gray levels into two classes, the background class and the vessel class. These two classes constitute the corresponding probability spaces and their probability density function $p(x,y,z)$ and $r(x,y,z)$ can be formulated as follows using the Parzen window estimate,

## Eq. (3)

$$p(x,y,z)=\frac{1}{N}\sum _{i=1}^{t}\sum _{j=1}^{{C}_{i}}\frac{1}{V{c}_{i}}\varphi (x,y,z;{x}_{j},{y}_{j},{z}_{j};{h}_{{C}_{i}}^{3})\phantom{\rule{0ex}{0ex}}r(x,y,z)=\frac{1}{N}\sum _{i=t+1}^{L}\sum _{j=1}^{{C}_{i}}\frac{1}{V{c}_{i}}\varphi (x,y,z;{x}_{j},{y}_{j},{z}_{j};{h}_{{C}_{i}}^{3})\mathrm{.}$$In general, a thresholding method determines the optimal gray value ${t}_{\text{opt}}$ of $t$ based on a certain criterion function. In order to achieve this goal, we should choose ${t}_{\text{opt}}$ such that the class $p$ and the class $r$ can be well separated. Here, we adopt the following criterion function,

## Eq. (4)

$${t}_{\text{opt}}=\underset{t}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{max}}\iiint {(p(x,y,z)-r(x,y,z))}^{2}\mathrm{d}x\mathrm{d}y\mathrm{d}z\mathrm{.}$$Using this threshold ${t}_{\text{opt}}$, each voxel of the image is divided into vessel region and background region: the voxels having an intensity that exceeds the threshold are referred as vessel region (see Fig. 3), and the remainder as background region.

## 2.1.3.

#### Extracting the segment candidates of the vessel centerline

In this section, we extract the segment candidates of the vessel centerlines by thinning 3-D binary images obtained in Sec. 2.1.2. We give a brief description of the thinning algorithm, and a detailed description of this algorithm has been published previously.^{6}

In the thinning algorithm, we use directional strategy in which each iteration step contains 12 successive parallel reduction operations according to the 12 directions. Thinning algorithm deletes border points of binary image in 12 iteration steps, and the process is repeated until no border points can be deleted. The ordered list of the deletion directions is proposed: (US, NE, WD, ES, UW, ND, SW, UN, ED, NW, UE, SD), as shown in Fig. 4.

Deletable points in a subiteration are given by a set of $3\times 3\times 3$ matching templates. A black point is deletable if at least one template in the set of templates matches it. Templates are usually described by three kinds of elements, “filled circle” (black), “open circle” (white), and “.” (“don’t care”), where “don’t care” matches either a black or white point in a given picture. In order to reduce the number of masks, we use additional notations. The first subiteration assigned to the deletion direction US can delete certain U- or S-border points; the second subiteration associated with the deletion direction NE attempt to delete N- or E-border points, and so on. The set of templates ${T}_{\mathrm{US}}$ (containing 16 templates) is presented in Fig. 5. The deletable points of the other subiterations can be obtained by proper rotations and/or reflections of the templates in Fig. 5.

A black point $x$ is deletable if it can be deleted by at least one template in ${T}_{\mathrm{US}}$; $x$ is nondeletable otherwise. The thinning algorithm goes around the object to be thinned according to the 12 deletion directions. As shown in Fig. 6, the thinning algorithm directly extracts “skeleton (vessel centerline)”.

## 2.2.

### Grouping and Selecting Candidate Segments

As shown in Fig. 6, vessel centerlines are composed of a good number of candidate segments. Let us suppose that there are $n$ candidate segments, that means the number of possible connections is $N={C}_{n}^{2}$. Thus, the complexity of our algorithm with ${C}_{n}^{2}$ possible connections is $T(n)=\phantom{\rule{0ex}{0ex}}O({n}^{3})$. Since the use of only simple image processing techniques will unavoidably lead to the extraction of false segments or a failure to extract complete segments, we present the following scheme for grouping and selecting the candidate segments.

In Fig. 7(a), there are four segments, two of which cross each other. At an intersection point trace segments tend to split, as shown in Fig. 7(b). In this case, there are four possibilities for connecting two segments. If we follow only one segment, we select one possibility, which is the most likely connection [Fig. 7(b), left]. However, when we select two connections simultaneously, other possibilities may be better [Fig. 7(b), right]. In general, such situations may occur in several places. Therefore, in order to select globally optimal connections, it is necessary to consider all the possibilities for connections simultaneously.

We formulate the problem as a combinatorial optimization problem which can be solved by a Hopfield-style network. For possible connections $N$ in an image, let $p=({p}_{1},{p}_{2},\dots ,{p}_{N})$ ($0\le {p}_{i}\le 1$) represent the probabilities that possible connections are true. Let $R=\{{r}_{ij}\}$ be the symmetric matrix that expresses the compatibility relations between the possible connections. Let ${c}_{i}$ represent the amount of certainty of each possible connection based only on local measurements. We find $p$ minimizing

Dynamic programming and gradient-descent methods are two widely used techniques in the optimization problem. Dynamic programming prescribes a search which tracks backward from the final step, retaining in memory all suboptimal paths from any given point to the finish, until the starting point is reached. The result of this is that the procedure is too computationally expensive for the large-scale problem. As mentioned above, the number of possible connections is extremely huge. Hereto, we use a gradient descent optimizer. The main advantage of this optimizer compared with the dynamic programming optimizer is that it causes a significant reduction in computation time. Let $q=({q}_{1},{q}_{2},\dots ,{q}_{n})$ be the gradient vector, where ${q}_{i}=\phantom{\rule{0ex}{0ex}}-\frac{\partial E(p)}{\partial {p}_{i}}$. By the restriction of $0\le {p}_{i}\le 1$, the update of $p$ can be expressed as follows,

where $m$ is the number of iterations, $\alpha $ is the step size, $|\xb7|$ represents the ${L}_{2}$ norm, and $g(x)$ is defined as follows:## Eq. (7)

$$g(x)=\{\begin{array}{lll}0,& \text{if}& x\le 0\\ x,& \text{if}& 0<x<1\\ 1,& \text{if}& x\ge 1\end{array}.$$The update Eq. (6) is iterated until $|{p}^{(m+1)}-{P}^{(m)}|\le \delta $ is satisfied.

We then select connections satisfying ${p}_{i}=1$ as the best connections. In the following, we describe the method for the detection of possible connections and determination of the compatibility relations ${r}_{ij}$ and certainty ${c}_{i}$.

We detect possible connections based on the smoothness and relative distance measures between two endpoints. The smoothness measure, $g$, is defined as

## Eq. (8)

$$g=({g}_{1}+{g}_{2})/2\phantom{\rule{0ex}{0ex}}{g}_{1}=d\xb7\mathrm{tan}\text{\hspace{0.17em}}{\theta}_{a}\phantom{\rule{0ex}{0ex}}{g}_{2}=d\xb7\mathrm{tan}\text{\hspace{0.17em}}{\theta}_{b},$$## Eq. (10)

$$g<{T}_{g}\phantom{\rule[-0.0ex]{1em}{0.0ex}}\text{and}\phantom{\rule[-0.0ex]{1em}{0.0ex}}d/({l}_{a}+{l}_{b})<{T}_{d},$$We embed into Eq. (5) the constraints which are used to select the optimum connection pattern among the detected possible connections. First, we assume that segments do not branch off. Therefore, if there are more than two possible connections at the same end point, they are incompatible. This relation can be represented as

## Eq. (11)

$${r}_{ij}=\{\begin{array}{cc}-\epsilon ,& \mathrm{if}\text{\hspace{0.17em}}\text{\hspace{0.17em}}i\text{'}\mathrm{th}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{segments are incompatible}\\ 0,& \text{else}\end{array},$$## Eq. (12)

$${c}_{i}={g}_{i}+\lambda {e}^{-{p}_{i}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{ln}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{p}_{i}},$$The iteration is terminated when any new connection is not selected after the grouping process. Furthermore, the gaps of the elected connections are filled using the B-spline.

Finally, the selected and grouped curves are identified as vessel centerlines (Fig. 9), and the remaining components (points and segments) are discarded as noise.

## 3.

## Results

We validate the proposed method quantitatively on a number of synthetic data sets, liver artery from abdominal CT, and pulmonary vessel from chest CT. In addition, a comparison is made with two state-of-the-art centerline detection methods using the liver artery from abdominal CT.

## 3.1.

### Validation of the Algorithm on Synthetic Data

An important consideration in the evaluation of a vessel centerline extraction algorithm is its validation. For methods which have been designed to work on segmented tubular structures, such as ours, quantitative validation on medical images is a challenge because a ground truth centerline is typically not defined. Thus, we have chosen to evaluate our approach quantitatively by creating synthetic tubular structures with known centerlines.

We created several ground truth centerlines by either using a fixed shape or by drawing them by hand. For each centerline, we then centered a sphere at one location and then translated the sphere along it. The volume swept by the moving sphere thus defined a synthetic tubular structure with a known ground truth centerline. The specific examples we created were:

• Helix: A helical centerline with a sphere of constant radius designed to test the method when both curvature and torsion are present in the anatomical structure.

• Spiral: A 3-D spiral with small gaps simulating low-contrast areas is generated to test the robustness of the method.

• Tree: A hand drawn tree-like centerline with a sphere of varying radius designed to test the ability of the method to handle complex branching structures.

Figures 10 and 11 show the results of applying our algorithm to the original and noisy versions of the objects with known centerlines, again with the same color coding. The figures indicate that the overlap between the ground truth and computed centerlines is very strong for every synthetic object considered. A quantitative analysis was carried out by computing the amount of overlap, the average distance and the maximum distance between points on the ground truth and computed centerlines. The results are shown in Table 1 for the original objects and their noisy versions.

## Table 1

Validation of the centerline extraction algorithm for the synthetic objects.

Name | Overlap (%) | Mean distance (mm) | Maximum distance(mm) |
---|---|---|---|

Helix structure | 100 | 0.03 | 0.76 |

Noisy helix structure | 100 | 0.03 | 0.78 |

Spiral structure | 100 | 0.02 | 1.02 |

Noisy spiral structure | 100 | 0.02 | 0.97 |

Tree structure | 98.2 | 0.09 | 1.32 |

Noisy tree structure | 97.3 | 0.10 | 1.57 |

The amount of overlap is consistently above 98% except for the noisy tree which has an overlap of 97.3%. The maximum distance is within 1.57 mm. However, the average distance is never more than 0.1 mm and in most cases is far less. These results provide some insight into how the method would behave for real medical data. First, the helix tube example shows that the method can handle structures with curvature. Second, the spiral phantom displays the ability of the approach to handle small vessel with small gaps. Third, the method can handle complex tree structures with small branches. Overall, our results provide quantitative evidence that the proposed algorithms perform well.

## 3.2.

### Extracting the Centerlines of Liver Artery from Abdominal CT

In liver surgery planning, the liver arteries play an important role in determining the optimal resections and in predicting blood supply to remaining liver tissue. The main problem in segmenting the liver arteries is that they are small and very poorly contrasted in the CT data, see Fig. 12. Another complicating factor is that the arteries run next to the larger portal veins, which increase the risk of segmentation leakage.

Eight abdominal CT angiography data sets were acquired using 64-row and dual-source CT scanners. The acquired image volumes are typically of the size $512\times 512\times \phantom{\rule{0ex}{0ex}}182\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{voxels}$ with a pixel resolution of about $0.68\times 0.68\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{mm}$ and a slice thickness of 1.5 mm. The liver regions in the images were roughly hand-segmented by a radiology specialist and used as a mask. In our study, the manually marked center points of the vessels by two radiologists were used as the gold standard for quantitative evaluation of the extraction accuracy of the vessel centerline. The radiologists provided gold standards by manually tracking the vessel tree and marking the center of the vessels using a graphical user interface (GUI) developed in our laboratory. On the GUI, the sagittal view, axial view, and coronal view of the CT images corresponding to the region where a vessel is being tracked are displayed on a monitor. The GUI has functions allowing the user to cine-page through the CT slices, scroll in and out of individual vessels, adjust window setting, and zoom to improve visualization. The user can manually track the vessel trees by marking the vessel center points in any one of the three views at each vessel branch and the center point location will automatically propagate to all three views. Using the GUI, the radiologists marked the centerline for each branch, and labeled the anatomical level of the arteries. After marking center points of the arteries, the reconstruction of the discrete center points was performed using B-spline interpolation.

In this section, a comparison is made with two state-of-the-art centerline detection methods. We briefly review the concept of both algorithms.

• Wörz method: Wörz and Rohr

^{21}utilized an incremental process that starts from a given point of the vessel and proceeds along the tangential direction. Along the tangential direction a prediction is placed, then corrects the prediction by searching for the point of maximum vesselness in an orthogonal plane. Because of the variations of vessel they introduced an adaptive model at the new position to improve the prediction. The tracking process starts from a seed point. If the tracking process terminates prematurely, further seed points will be added.• Agam method: Agam et al.

^{22}used a set of filters to extract centerlines, which is based on the eigenvalues of a correlation matrix of regularized gradient vectors. In this method, vessels are modeled by tubular segments where each tubular segment is approximately a cylinder with a Gaussian intensity profile at the section plane. Based on this model, the centerline is in a direction which is orthogonal to all the gradient vectors inside the vessel.

The vessel centerlines of the liver images included in the database were again segmented using the Wörz method, Agam method, and our method. Segmentation results of hepatic artery for three CT liver volumes are shown in Fig. 13. As shown in Fig. 13, the Wörz method and Agam method failed to extract small vessels in low contrast areas (rectangular regions). For the Wörz method, the tracking process terminated at the low contrast area, more seed points should be given for further process. Even when four seed points were given, the Wörz method still failed to extract the completed centerline in case 3. For the Agam method, the correlation-based filters can detect both small and large vessels due to the use of multiple sets of eigenvalues of the correlation matrix. However, the Agam method also failed for some small vessels. Averaged across all CT data sets, the overlap with the human expert extractions was used for accuracy evaluation. Accuracy results for the comparison are presented in Table 2. The results show that our method achieved the highest segmentation accuracy in comparison with these methods.

## Table 2

The overlap (%) between manual extraction and the computed extraction by three different methods on the eight liver artery data sets.

Methods | Data set numbers | |||||||
---|---|---|---|---|---|---|---|---|

1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | |

Wörz method | 75.1 | 87.1 | 81.2 | 81.1 | 90.1 | 84.3 | 73.2 | 90.7 |

Agam method | 74.0 | 83.5 | 88.5 | 83.7 | 89.9 | 83.7 | 79.3 | 87.0 |

Our method | 96.3 | 97.2 | 92.8 | 98.5 | 95.7 | 98.6 | 92.7 | 96.8 |

## 3.3.

### Extracting the Centerlines of Pulmonary Vessel from Chest CT

Many of the published vessel segmentation and tracking methods provided accurate results in 2-D or 3-D images for vascular structures in the retina, brain, heart, etc. However, few studies have been conducted for segmentation, tracking, and reconstruction of the pulmonary vessel tree on CT images because the pulmonary vessels are more complicated compared to the vessels in other parts of the body in several aspects: widely distributed CT values, large variations of vessel sizes, and the complicated branching structures.

We used three chest CT angiography data sets taken by 64-row and dual-source CT scanners. The volume size is $512\times 512\times 265$ with an in-slice pixel resolution of $0.68\times 0.68\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{mm}}^{2}$ and slice thickness was 1.5 mm. In our vessel segmentation scheme, the lung region is first extracted to reduce computation time and avoid the effects caused by other nonvascular structures such as ribs, motion artifacts of the heart, and the edges of the heart and chest walls. An automatic segmentation method^{23} was applied to the volume of the scan within the patient body to extract lung regions.

For the three clinical cases, there is no “ground truth” to prove the presence or absence of the vessels and their sizes. In our study, the manually marked center points of the vessels by three radiologists were used as the gold standard for quantitative evaluation of the segmentation accuracy of the vessels. A total of 2655, 2782, and 2853 vessel center points were marked by radiologists, respectively, for cases 1, 2, and 3. It is a demanding task for radiologists to mark the center points manually along small peripheral vessels on images displayed on a computer monitor. The radiologists decided to stop marking the center points on the peripheral vessels at certain levels, leaving some vessels not marked. Therefore, a visual inspection was performed to evaluate those vessels that were not marked by radiologists. Figure 14(b) shows an example of the manually marked vessel center points (red points) superimposed on the computer segmented vessel centerline. Figure 14(c) shows another view from the automatically segmented vessel centerlines shown in Fig. 14(b). It can be seen that many of the small vessels not marked by radiologists were extracted by our segmentation method. These small vessels were estimated to be smaller than 2 mm in diameter by manual measurement using an electronic ruler on the computer GUI. The accuracy of the vessel centerline segmentation was evaluated using three data sets. Averaged across three data sets, average distance to the ground-truth centerlines, was 0.32 mm.

## 4.

## Discussions

A fully automatic framework for vessel centerline extraction based on a multistep scheme has been presented. Our approach includes three major steps: (1) enhancing thin vessel, (2) extracting vessel centerline segment candidates, and (3) grouping and selecting vessel centerline segment. Not all of these steps are always required. For the thick vessels (in diameter $>2\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{voxels}$), vessel enhancement step may be omitted.

This approach has five main advantages. First, compared to existing vessel centerline tracking approaches, it does not need user intervention; and an automatic vessel centerline extracting is performed by using a global optimization algorithm. Second, with a filter of 13 orientations, the approach can detect low-contrast and narrow vessel. Third, some methods often fail to connect segments belonging to one vessel in low contrast areas, our approach can fill gaps in blood vessels. Fourth, our framework avoids incomplete and disconnected extraction of the vessel centerline. And fifth, the approach is robust to noise and can produce vessel centerline extractions that have level of variability similar to those segmentations obtained from human raters.

Quantitative evaluation of vessel centerline accuracy is a challenging problem because there is no ground truth for the vessel centerline extraction for clinical cases. Although a tubular phantom can be constructed with known vessels, it is difficult to construct a realistic phantom with a complete arterial tree mimicking complicated vascular structures. In a recent study by Shikata et al.^{24} about 2000 points were manually placed along the pulmonary vessels in each case as true vessel points for quantitative evaluation of a pulmonary vascular tree segmentation method. This method will provide an estimate of the fraction of the vessel tree segmented by the computer relative to that marked manually if the representative points are placed with reasonably balanced distribution in the entire pulmonary vascular tree.

It is a demanding task for radiologists to mark the center points manually along small peripheral vessels on images displayed on a computer monitor. The radiologists decided to stop marking the center points on the peripheral vessels at certain levels, leaving some vessels not marked. Therefore, a visual inspection was performed to evaluate those vessels that were not marked by radiologists (Fig. 14).

In our experiment using the chest CT data sets, three patient cases were used. Although the number of cases is small, the total number of center points that mark the entire pulmonary vessel tree down to very small vessels, for example, the segmental and sub segmental levels of arteries is quite large (more than 8290 for three cases). The evaluation including quantitative comparison and visual inspection on these three cases allows us to obtain a reasonable assessment of the vessel centerline extraction performance in this preliminary study.

The computation time has also been considered, which might be critical in some applications. In the present study, the methods have been programmed using MATLAB 6.5 in an Intel(R) Core(TM) i5-4440 computer, 3.2 GHz, 16 G RAM. The processing time of the proposed method was about 8 min for one liver artery system and about 72 min for one chest CT angiography data.

## 5.

## Conclusion

We have successfully developed an automatic method for the 3-D vessel centerline extraction in CT images. Experimental results on synthetic and clinical data sets have shown that our method can extract complete and continuous centerline in the problematic regions that contain small and low-contrast vessels, and gaps between vessels. The extraction algorithm has high robustness toward noise and can produce vessel centerline extractions that have level of variability similar to those obtained manually. This study suggests that the proposed framework can be good enough to replace the time-consuming and tedious slice-by-slice manual extraction approach.

In the future work, we plan to get a complete segmentation of the vessels. These structures need to be filled starting from the detected centerlines.

## Acknowledgments

The authors would like to thank the anonymous reviewers for constructive suggestions. This work is supported by the National Natural Science Foundation of China under Grant No. 61073124; the National Basic Research Program of China (973) under Grant 2011CB707701.

## References

## Biography

**Yuanzhi Cheng** received his PhD degree from Osaka University, Osaka, Japan, in 2007. He is currently an associate professor in the School of Computer Science and Technology, Harbin Institute of Technology. He has published more than 30 papers in scientific journals. His research interests are concentrated on pattern recognition, medical image processing, and computer-assisted surgical system.

**Yadong Wang** is currently a professor in the School of Computer Science and Technology, Harbin Institute of Technology. His research interests include pattern recognition, machine learning and bioinformatics. He has published more than 60 papers in scientific journals.

**Shinichi Tamura** received his PhD degrees in electrical engineering from Osaka University, Osaka, Japan, in 1971. He is currently a professor in the Center for Advanced Medical Engineering and Informatics, Osaka University. He has published more than 260 papers in scientific journals and received several awards from journals, including *Pattern Recognition* and *Investigative Radiology*. His current research activities include works in the field of medical image analysis and its applications. He is a fellow of IEEE.