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.¹⁸ In addition, Gabor filters¹⁹ 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.

Fig. 1

(a) The six major directions. (b) The proposed 13 orientations for thin vessels enhancement. The 13 directions correspond to the 13 neighbors of a point on a face-centered cubic grid.

Fig. 2

Set of one-voxel width line detector filters used for thin vessel enhancement.

2.1.2.

Extraction of vessel region

In order to obtain the vessel region, we employ the Parzen window estimate²⁰ 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)

$$\begin{gathered} p(x,y,z,w_i)=p(w_i)p(x,y,z|w_i) \\ =\frac{C_i}{N}·\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), \end{gathered}$$

where $h$ $C_i$ denotes the window width of the $i$’th gray-level set, ($x_j,y_j,z_j$) denotes the $j$’th voxel in $w_i$, ${\mathrm{VC}}_i$ denotes the volume of the corresponding cube with $h$ $C_i$ as its length of each edge, ${\mathrm{VC}}_i=h_{Ci}^3;p(w_i)$ can be approximated by $C_i/N$, $1\le i\le L$, and $\varphi (•)$ is a kernel function. In this article, the following Gaussian function for $\varphi (•)$ is taken,

Eq. (2)

$$\varphi (x,y,z;x_j,y_j,z_j;h_{C_i}^3)=\frac{1}{3\pi }\exp [-\frac{{(x-x_j)}^2+{(y-y_j)}^2+{(z-z_j)}^2}{3h_{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)

$$\begin{gathered} 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) \\ 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{.} \end{gathered}$$

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}{\arg \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.

Fig. 3

The result of the liver artery vessel region extracted by a Parzen window method.

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

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.

Fig. 4

The 12 deletion directions of 12-subiteration thinning algorithm.

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.

Fig. 5

The set of templates $T_{\mathrm{US}}$ assigned to the first subiteration of the thinning algorithm. Notations: every position marked “filled circle” matches a black point; every position marked “open circle” matches a white point; every “.” (“don’t care”) matches either a black or a white point; at least one position marked “x” matches a black point; at least one position marked “Y” matches a black point; at least one position marked “v” matches a white point; at least one position marked “w” matches a white point; two positions marked “z” match different points (one of them matches a black point and the other one matches a white one).

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)”.

Fig. 6

The liver artery centerlines obtained by thinning three-dimensional binary images. Vessel centerlines are composed of a good number of curve segments.

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)=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.

Fig. 7

(a) Magnification of the rectangular region at the lower right corner of Fig. 6. (b) Graph showing four endpoints (A, B, C, and D) for possible connections. Left: One possibility, which is the most likely connection. Right: Two possibilities, which is globally better.

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=-\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,

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

Fig. 8

Smoothness and relative distance measures for detection of possible connections and determination of certainty $c_i$.

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

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.

Fig. 9

Final result of the extracted liver artery centerlines. Parameter used: $T_g=15$, $T_d=1.0$, $\lambda =0.001$, $\epsilon =–1$, and $\delta =0.0005$.

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:

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.

Fig. 10

Left column: Synthetic objects with centerlines of fixed shape. Right column: The corresponding ground truth and computed centerlines. For each object, the overlapping points between the ground truth centerline and the computed one are shown in black, the points belonging to the ground truth centerline but not the computed one are shown in green, and the points belonging to the computed centerline but not the ground truth centerline are shown in red. (a) and (b) A helix structure. (c) and (d) A spiral structure. (e) and (f) A tree structure. A vessel tree shows the ability to detect bifurcations. The vessel branches are 1 mm in diameter.

Fig. 11

Left column: Synthetic objects with centerlines of fixed shape. Boundary noise was added randomly. Right column: The corresponding ground truth and computed centerlines. For each object, the overlapping points between the ground truth centerline and the computed one are shown in black, the points belonging to the ground truth centerline but not the computed one are shown in green, and the points belonging to the computed centerline but not the ground truth centerline are shown in red. (a) and (b) A helix structure. (c) and (d) A spiral structure. (e) and (f) A tree structure. A vessel tree shows the ability to detect bifurcations. The vessel branches are 1 mm in diameter.

Table 1

Validation of the centerline extraction algorithm for the synthetic objects.

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.

Fig. 12

Slices from the arterial phase of a CT angiography data set. The main problem in segmenting the liver arteries (arrows) is that they are small or very poorly contrasted in the CT data.

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 182\mathrm{voxels}$ with a pixel resolution of about $0.68\times 0.68\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.

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.

Fig. 13

Liver artery centerline segmented with three different methods. The first row, the second row and the third row illustrate the three different liver artery systems, respectively. First column shows the Wörz segmentation, the second column shows the Agam segmentation, and the third column shows our segmentation. In figure, the rectangular regions indicate several obvious failures; the black dots give the initial seed points for Wörz method; the blue lines show the centerlines of the manual segmentation; the red lines are the computed centerlines.

Table 2

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

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{\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²³ 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.

Fig. 14

Comparison of radiologist’s marked vessel center points (gold standard vessel center points) with the automatically segmented vessel centerlines. To facilitate visualization, we give 700 marked vessel center points. (a) Slices of a chest CT angiography data set. (b) One view from the automatically segmented vessel centerlines within the left and right lungs superimposed with the gold standard vessel center points (red points). (c) Another view from the automatically segmented vessel centerlines within the left and right lungs superimposed with the gold standard vessel center points (red points).