Introduction

Coronary computed tomographic angiography (CCTA) is a non-invasive imaging modality which can visualize the heart including the coronary arteries. It can provide not only anatomical information on the coronary arteries, but also pathological information such as presence and extent of calcifications, plaque burden, stenoses and occlusions, which are useful in the diagnosis and treatment of coronary artery disease (CAD) [1, 2]. As each CCTA dataset is built up by a stack of 2D images, typically more than 200 slices, post-processing tools are required to facilitate the image interpretation by cardiologists and radiologists. For example, the multi-planar reformatted (MPR) image in transversal direction along the vessel allows observing the plaque morphology and its effect on the vessel lumen diameter. A curved multi-planar reformatted (curved MPR) image can display a tortuous branch in a single 2D image which enables the examination of a lesion and calcified spots along the course of the branch [3]. Since those post-processing tools rely on the availability of accurate centerlines, a robust coronary artery extraction algorithm is an essential part of these processing tools. Furthermore, centerline extraction of the coronary arteries is also one of the prerequisite steps in subsequent analysis steps, such as detection of lesions [4], vascular shear stress estimation [5] and image fusion [6].

A significant amount of research has been done on vessel segmentation in medical images, including coronary artery extraction [7]. However, in most solutions of coronary artery extraction in CCTA images, user interactions were required. For instance, initializing the vessel of interest by means of the manual definition of the start and end points, or adding intermediate points to bridge gaps [8–12]. These manipulations require anatomical and pathological knowledge and increase the processing time. If these tools are to be used in clinical practice, further automation towards a robust solution is required.

Several approaches were dedicated to the automatic centerline extraction of the coronary tree in CCTA images, which were mainly based on morphological operators [13], model-fitting [14], medialness filter [15] and fuzzy connectedness [16]. These methods all had some difficulties to extract the distal parts of coronary arteries or the segments with narrowings, calcifications and motion artifacts. An inaccurate extraction or a premature termination of the algorithm may occur if, for example, the intensity distribution or morphology differs from the normal vessel. Bauer et al. [17] presented another method which combined the gradient vector flow (GVF) [18] with Frangi’s vesselness measure [19] to avoid isotropic diffusion in Gaussian scale space. It obtained a better vessel enhancement especially when the structures near a coronary artery, such as cardiac chambers and calcifications, had the same or higher gray level. However, the extracted centerlines only had voxel-accuracy at most and the accuracy of centerlines was decreased at vessel junctions because of the small gaps generated by vesselness filter at vessel junctions.

In our previous work [20], an automatic method was presented based on connected component analysis and wave propagation. The evaluation results showed that the intensity-based thresholding used in this method should be improved for better performance of extraction ability.

Therefore, in this paper, we present a novel fully automatic coronary artery extraction method for CCTA images mainly relying on an improved Frangi’s vesselness filter. The new vesselness filter adopts a discriminator based on local geometric features to decrease the false positive responses acquired in Frangi’s vesselness filter. Automatic ostium detection, branch searching and centerline refinement steps followed by the extraction of the centerlines of the whole coronary tree. This method was validated in two ways: first using the Rotterdam Coronary Artery Algorithm (RCAA) Evaluation Framework [21], and secondly by comparing our results with reference data obtained from 50 clinical CCTA datasets with various vessel pathologies and image qualities.

Methods

According to the processing pipeline displayed in Fig. 1, our automatic coronary tree extraction algorithm can be divided into several successive steps, i.e. pre-processing, improved Frangi’s vesselness filtering, centerline extraction, branch searching and centerline refinement. These steps will be described in the next paragraphs.

Fig. 1
figure 1

Processing pipeline of the coronary tree extraction algorithm

Full size image

Pre-processing

In order to achieve better computational efficiency, CCTA images are interpolated linearly and down-sampled to acquire isotropic volumes with the axial image size equal to 256 × 256 pixels. Pulmonary vessels are removed by a morphological closing operator. A spherical kernel with a radius of 8 voxels is applied to the pulmonary region with image intensities lower than −224HU. A threshold of 676HU as adopted in [8] is used to replace the extreme high image intensity caused by the presence of calcifications, stent or pacemaker, etc.

Improved Frangi’s vesselness filter

Frangi’s vesselness filter [19] measures the similarity to a tubular structure with a Hessian matrix, which has been widely used for the vessel enhancement in medical images. However, as reported in [22], a major drawback of this filter is that a step-edge response is obtained at the boundary of a non-vessel structure. This step-edge response becomes more pronounced in CCTA images, since cardiac chambers have the same intensity as coronary arteries, which can easily lead to a false extraction result. For example, in Fig. 2a, the boundary of the right atrium was extracted erroneously with Frangi’s vesselness filter in which a threshold of 80 was used to remove low filter responses. Increasing this threshold could alleviate this false extraction. Unfortunately, this came at the expense of the extraction of distal parts and side branches as shown in Fig. 2b. In addition, this threshold could not be fixed for different datasets as the value of vesselness response depends on the contrast of coronary arteries. In Fig. 2c, when our improved Frangi’s vesselness filter was used, the erroneous extraction was removed and meanwhile the coronary tree could be extracted completely.

Fig. 2
figure 2

Coronary tree extraction results by Frangi’s vesselness filter with the thresholds of 80 (a) and 140 (b) respectively, and with our improved Frangi’s vesselness filter (c). The aorta was extracted by the aorta detection algorithm. The red circle marks the false extraction of the right atrium. Increasing the threshold could alleviate the false extraction. But the distal parts and side branches were not extracted at the threshold of 140 as shown in (b). The result of Frangi’s vesselness filter was rescaled to [0, 1024]

Full size image

We present a novel method to discriminate between the false step-edge responses and the positive vessel responses in Frangi’s vesselness measurement by adding local geometric features. These features are obtained by performing ray-casting in a local sphere S at a voxel x with radius R. A total number of N S rays are uniformly distributed cast from the center x to the sphere’s surface. For each ray, the casting stops when it reaches the boundary of the sphere or when it encounters the boundary of a local structure. This boundary is measured by comparing the image intensity at the current position p on the ray and at the spherical center x. When \( |{I(\user2{p})} - I(\user2{x})\left| {/I(\user2{x})} \right.> \varepsilon \), this indicates a local boundary (ε is the threshold that controls the sensitivity of local boundary detection). After casting in all directions, only the rays whose lengths are longer than 0.8 R are preserved. Figure 3a shows ray-casting results in 4 different synthetic shapes; a tube, a tube junction, a plane and a sphere. In the tube and the tube junction, the rays are oriented towards the tube directions. While, in the plane and the sphere, the rays are scattered without giving any explicit direction. Figure 3b illustrates the ray-casting results at four points in a CCTA image. Three points are located in the coronary arteries and one at the boundary of cardiac cavity. The rays displayed in the image show that points A and B are a trifurcation and a bifurcation, respectively. Point C is located in a vessel where the rays give two opposite directions. Point D is at the boundary of a large structure and the rays do not have any clear direction. The distribution of these rays describes the geometric features of the local shape, which is then measured to discriminate between the false step-edge responses of cardiac cavities and the desired responses of coronary arteries.

Fig. 3
figure 3

a Ray-casting results in 4 synthetic datasets. b Ray-casting results in a CCTA image. Points A, B, C are located in the coronary arteries and point D is located at the boundary of cardiac cavity. c shows connected component labeling results at points A and D with different colors. Each preserved ray was mapped to a point on the spherical surface. Gray lines represent the neighborhood relationship of the points on the spherical surface. Lines from sphere center are the average directions of each group

Full size image

To quantify the distribution of the rays, connected component labeling on the sphere is performed to cluster the rays oriented in the same local direction into M groups. Figure 3c shows the results of connected component labeling at points A and D from Fig. 3b. The distribution of the rays at voxel x is then measured as follows

$$ G(\user2{x}) = \prod\limits_{j = 1}^{M} {F(r_{j} )||\user2{v}_{j} ||} $$
(1)

in which, j is the group index. r j equals to N j /N S (N j is the number of rays in the j-th group). v j is the average vector of the normalized direction vectors of the rays in the j-th group and ||·|| is the Euclidean norm. The term ||v j || is introduced into G(x) to decrease the response when the rays in one group located in a plane-like structure (bottom-left figure in Fig. 3a). The function F(r) is designed to penalize groups with a large number of rays and is based on a sigmoid function. This function is defined as,

$$ F(r) = 1 - \frac{1}{{1 + e^{ - \rho (r - \gamma )} }}, \, r \in (0,1] $$
(2)

in which ρ controls the sharpness of the curve at the threshold of γ (γ ∈ (0, 1]). Since the voxels located in the cardiac cavities have relatively large values of r, γ is set to be a relatively small value to let F(r) approach 0.0 rapidly when the number of rays in one group is larger than γ·N S . Thus, G(x) can discriminate the voxels located in the coronary arteries (e.g. points A, B and C in Fig. 3b) and the large structure such as cardiac chamber (e.g. point D in Fig. 3b). For example, if we set ρ = 40 and γ = 0.1, the value of G(x) at points A, B, C and D in Fig. 3b were 0.71, 0.85, 0.94 and 2.2 × 10−7 respectively. This geometric measure serves as a weighting factor to derive our novel vesselness measurement by multiplying it with Frangi’s vesselness result.

Frangi’s vesselness filter is calculated at several scales to adapt for different vessel sizes, so the maximum response at the optimal scale is kept. As the step-edge responses of heart chambers are obtained at large scales in the multi-scale strategy, the calculation of geometric feature measurement G(x) only needs to be performed when the Frangi’s vesselness response V F (x) and the optimal scale σ(x) are above thresholds T F and σ F respectively. Thus, the improved vesselness response V I (x) can be written as

$$ V_{I} (\user2{x}) = \left\{ \begin{array}{l} G(\user2{x})V_{F} (\user2{x}), \, V_{F} (\user2{x}) \ge T_{F} , \, \sigma (\user2{x}) \ge \sigma_{F} \, \hfill \\ V_{F} (\user2{x}), \, V_{F} (\user2{x}) \ge T_{F} , \, \sigma (\user2{x}) < \sigma_{F} \hfill \\ 0, \, V_{F} (\user2{x}) < T_{F} \hfill \\ \end{array} \right. $$
(3)

T F is used to ignore very small responses to improve computational efficiency. σ F only includes the large scales used in the Frangi’s vesselness filter.

Aorta and ostium detection

In order to find the coronary tree correctly, its origins from the ascending aorta should be detected. Therefore, the ascending aorta is segmented first. It is identified based on the Hough circle transform in the first several axial slices. Then the estimated center and radius are used for the initialization of the aorta detection in the next axial section. Aorta detection is performed slice by slice until the relative difference of the average intensity inside the circle or the distance of the circle centers between two successive slices is larger than the threshold T C or d C respectively. In this paper, T C and d C were fixed at 5% and 7 mm respectively. These detected circles define a region of interest in which the aorta is segmented by intensity-based region growing in each 2D slice. Next, a 3D morphological opening operator using a spherical kernel with the radius of 6 voxels is applied to remove possible leaks in 2D slices and to generate the 3D segmentation of the aorta. Once the aorta region is obtained, two voxels with the maximum vesselness responses are detected around the aorta region to initialize the ostia of the left and right coronary trees. The angle θ defined in the polar coordinate system shown in Fig. 4 is used to distinguish between the ostia of the left and right coronary trees. Based on the normal anatomy of the coronary arteries, the range of θ for the left ostium is θ L  ∈ [−90∘, 45∘] and for the right ostium is θ R  ∈ (45∘, 135∘].

Fig. 4
figure 4

Illustration of ostia detection step. Two voxels with the maxima vesselness responses, O L and O R , are detected to be the ostia of the left and right coronary trees respectively. The angular search regions for the left and right ostium are θ L  ∈ [−90∘, 45∘] and θ R  ∈ (45∘, 135∘] respectively. θ = 90° indicates the anterior side of the patient

Full size image

Centerline extraction

As a first step in the centerline extraction, a binary vessel enhanced image is created by removing those voxels with very low responses (V I (x) < T I ) from the vesselness result. Connected component labeling follows to remove small components which only consist of a few voxels. Then, skeletonization by successively eroding border voxels is done to generate the central axes of the structures in the binary image. The tree-like skeletons originating from the positions closest to the detected ostium points are extracted as the ‘initial tree’. The ostium points are then updated to the beginning points of these extracted tree-like skeletons.

Branch searching

In CCTA images, the intensity along a coronary artery may vary because of artifacts and the presence of stenoses and calcifications, which can result in gaps in the binary vessel enhanced image. In order to bridge these gaps, branch searching is applied to connect the separated distal parts to the ‘initial tree’.

Figure 5 illustrates the strategy of branch searching. The search starts from the points with large curvature changes (e.g. point B) and the end points (e.g. points A and C) on the ‘initial tree’. If an unconnected branch (Branch I to IV) is found within a rectangular searching region, a connection (dot lines) with this branch is found by a wave propagation algorithm [12]. It could be that connections to the same unconnected branch are found from different starting points (for example connections c 2 and c 5 with Branch IV in Fig. 5). In this case, only the connection with the minimum cost is selected as the optimal connection. The connection cost is calculated by combining the arrival time of the wave propagation and the angle differences between the parent and child branches. It is defined as,

$$ C(t_{SE} ,\user2{d}_{S} ,\user2{d}_{E} ,\user2{d}_{SE} ) = t_{SE} /t_{SE}^{\max } + (\left\| {\user2{d}_{S} - \user2{d}_{E} } \right\| + \left\| {\user2{d}_{S} - \user2{d}_{SE} } \right\| + \left\| {\user2{d}_{E} - \user2{d}_{SE} } \right\|)/6 $$
(4)

in which t SE is the arrival time of the wave propagation from the start to the end points. d S and d E are the normalized tangential directions at the start and end points, d SE is the normalized direction of the found connection. Finally, \( t_{SE}^{\max } \) is the maximum t SE of all the found connections to the same unconnected branch. ||·|| is the Euclidean norm. The branch extension continues until no more parts are found or when the number of found parts reaches a user-defined maximum.

Fig. 5
figure 5

Illustration of the branch searching. The bold lines represent the extracted ‘initial tree’. The thin lines are unconnected branches, and the dotted lines are connections between the ‘initial tree’ and the unconnected branches found within the rectangular searching regions

Full size image

Centerline refinement

Although the down-sampled images used in the previous steps allow our method to be computational efficiency, they decrease the centerline accuracy. In addition, the intensities of calcifications are changed to avoid extreme high image intensities in the pre-processing step, and in consequence, centerlines obtained by skeletonization within these calcified regions may not be the real lumen centerlines. For these reasons, a centerline refinement step using the original images is included to improve the accuracy of centerlines.

Straightened multi-planar reformatted (MPR) images are constructed from the initial centerlines for the refinement. Lumen contours are automatically detected in 4 longitudinal cuts (Fig. 6a). The lumen contour detector is based on a minimum cost path method [12] and excludes calcified regions to provide a better delineation of vessel lumen. The calcified regions are avoided by detecting large deviations from the expected lumen intensity based on the average intensity in the center of the stretched MPR images. Finally, the refined centerline is generated as the average of these detected lumen contours (Fig. 6b).

Fig. 6
figure 6

a Lumen contour detection in 4 longitudinal cuts through a straighten MPR image. b Refined centerline (red) is the average of the detected lumen contours transformed back to the original 3D space (white)

Full size image

Results

The whole coronary tree extraction algorithm was implemented in C++ within the MeVisLab platform (http://www.mevislab.de). The method was first evaluated with the RCAA evaluation framework used in the MICCAI Coronary Artery Tracking challenge 2008 (CAT08) [21], which also allows us to compare our algorithm with other published algorithms. In this framework, 32 CCTA datasets were acquired by two CT scanners (Siemens Sensation 64 and Siemens Somatom Definition) and graded according to their image qualities and calcium scores. In each dataset, four coronary branches (LAD, LCx, RCA and one large side-branch) were manually annotated by three experts. The reference centerlines of 8 datasets (No. 0–7) are publicly available as training datasets; leaving the other 24 datasets (No.8–31) for testing. The evaluation results in this paper were obtained by uploading our extraction results through the website of the RCAA evaluation framework (http://coronary.bigr.nl).

The parameters α, β and c in Frangi’s vesselness filter [19] were set to be 0.5, 0.5 and 300 respectively. The scale range for multi-scale strategy in Frangi’s vesselness filter [19] was from 1 to 3 voxels with 1 voxel step. The results of Frangi’s vesselness filter were normalized to [0, 1024]. The raycasting parameters R and ε varied with the optimal scale σ F (x): R(x) = 3(σ F (x) + 1), ε = 0.10 if σ F (x) = 2 and ε = 0.15 if σ F (x) = 3. The other parameters used were N S  = 429, ρ = 40, γ = 0.1, T F  = 50, σ F  = 2, T I  = 20. These parameters were determined by comparing the extracted centerlines with the reference standard of the training datasets to obtain the highest scores in the RCAA evaluation framework. Later, these parameters were fixed for all the other datasets. In all of the datasets, the aortas and ostium points were detected correctly. Figure 7 displays the extraction results in 2 datasets (No. 11 and 21). Aortas were segmented first and cropped automatically based on cutting planes defined by two detected ostium points and aorta center. As the results show, most coronary arteries were extracted as the ‘initial tree’ (red solid lines) and then the distal parts where added after branch searching (yellow dashed lines). In dataset 11 (left figure in Fig. 7) the separated LAD artery and its side-branches were successfully reconnected to the ‘initial tree’ by branch searching.

Fig. 7
figure 7

Coronary tree extraction results in 2 datasets from the RCAA evaluation framework (No. 11 and 21). Red solid lines represent the ‘initial tree’ and yellow dashed lines are branches which have been found in the branch searching step. The aortas were segmented automatically in the aorta detection step. The blue points indicate the ostia of the left and right coronary trees

Full size image

RCAA evaluation

Four centerlines for evaluation were selected automatically by reference points provided by the RCAA evaluation framework. The algorithm was evaluated by comparing the extracted centerlines with the reference standard in two aspects: extraction ability and accuracy. The extraction ability is measured by three overlap values, which are overall overlap (OV), overlap until first error (OF) and overlap with the clinically relevant part of the vessel (OT). The accuracy is measured by the average distance to the reference centerlines. More details about these measurements can be found in [21].

The evaluation results of the testing datasets (No.8–31) are summarized in Fig. 8. The average overlap measurements are OV = 93.7%, OF = 74.2% and OT = 95.9%. The high OV and OT values demonstrate good extraction capability of our method. The lower values of OF obtained in some datasets are caused by severe calcifications or artifacts in the proximal parts, where vessel lumen sometimes is difficult to delineate in the refinement step. If we only consider the datasets with good image quality and low calcium score (No. 9, 16, 21, 22, 28 and 30), the average OF is 92.4% which approximates the average OV of 95.3%. The overall average centerline accuracy is 0.30 mm which is smaller than the mean voxel size of the datasets, 0.32 × 0.32 × 0.4 mm3. Table 1 shows the comparison with the 5 automatic methods presented at the CAT08 workshop. The results show that overlap measurements of our method are the highest. Our accuracy is ranked at the 3rd position, and only has a slightly difference with the best one. Compared with our previous method ‘CocomoBeach’ (in Table 1), the extraction ability has improved significantly.

Fig. 8
figure 8

Results of the average overlap (top) and accuracy (bottom) from the RCAA evaluation framework per dataset

Full size image
Table 1 Comparison to five automatic extraction algorithms evaluated in the CAT08 [21]
Full size table

Clinical evaluation

In the secondary evaluation phase, additional analyses were performed with another 50 clinical datasets to test our method with various vessel pathologies and image qualities. These datasets originate from patients who had sequentially undergone CCTA imaging and conventional invasive coronary. Patients were derived from ongoing clinical registry if they met the following selection criteria: (1) presence of both CCTA and invasive coronary angiography, (2) diagnostic image quality of both CCTA and invasive coronary angiography, (3) presence of sinus rhythm and (4) absence of atrial fibrillation, renal dysfunction (glomerular filtration rate <30 ml/min), documented iodine-containing contrast allergy and pregnancy. CCTA was performed clinically for non-invasive assessment of known or suspected coronary atherosclerosis. Subsequently, patients were referred for invasive coronary angiography because of imaging results or clinical presentation to provide further information on the extent and severity of CAD. Ten images were acquired by a 320-slice CT scanner (Toshiba Aquilion One) and the other 40 were acquired by a 64-slice CT scanner (Toshiba Aquilion 64). The average voxel size is 0.37 × 0.37 × 0.34 mm3. In each image, the reference centerlines were obtained with the following steps:

  1. (1)

    An initial path line of a vessel with clinical interest was extracted between user-defined start and end points using the wave propagation method [12].

  2. (2)

    Based on this path line, lumen contours in transversal planes along the vessel were initialized by a dedicated lumen contour detector [12].

  3. (3)

    An expert from the cardiology department manually corrected the lumen contours if required.

  4. (4)

    The reference centerline was defined as the linear interpolation of the centers of gravity of lumen contours in successive transversal planes along the vessel.

In total 100 reference centerlines were generated with an average length of 113 mm. Quantitative analysis on the generated lumen contours [23] showed that the average area stenosis in these vessels is 44% (±20%) and the average diameter stenosis is 27% (±14%). The maximum lumen area stenosis is 92% and the maximum diameter stenosis is 73%. After applying our new coronary artery extraction algorithm, the overlap and accuracy were measured in the same way as in the framework of the RCAA.

Since the reference centerlines were limited to the proximal parts with clinical interest, any additional distal parts extracted by our new method were discarded in the evaluation. For example, Fig. 9 illustrates one coronary artery extraction result in this evaluation phase. The centerline in green was automatically selected from the extracted coronary tree and compared with the reference which was defined by the lumen contours in red. The distal part of the green centerline which was longer than the reference centerline was discarded in the evaluation. The centerline is correctly extracted when the green line travels within the region defined by the red lumen contours.

Fig. 9
figure 9

A centerline extraction result of one dataset used in the secondary evaluation phase. The blue lines are extracted centerlines of the coronary arteries. The LCx in green is selected to compare with the reference centerline which is defined by the centers of gravity of the lumen contours in red

Full size image

Our coronary extraction algorithm was applied to these 50 datasets with the same parameters as used for the datasets of the RCAA evaluation framework. The aorta in each dataset was detected successfully; while ostium detection failed in four datasets because of the anomalous origins of coronary trees. Two user-defined points to locate the ostia were used to extract the coronary tree in these four datasets. The average centerline overlap and accuracy of these 100 vessels is 96.1% and 0.33 mm, respectively. The high overlap value demonstrates the centerlines extracted by our automatic method are in high accordance with the reference centerlines. The accuracy is still smaller than the average voxel size as shown in the RCAA evaluation framework.

With the datasets in these two evaluation phases, different pathologies and image qualities were evaluated intensively. In Fig. 10, curved MPR images display the centerlines extracted successfully in the coronary arteries with lesions, severe calcifications, stents and datasets with significant image noise.

Fig. 10
figure 10

Examples of four curved MPR images generated from centerlines (in red) extracted successfully in the coronary arteries with a lesion (a), severe calcifications (b), a stent (c) and a dataset with significant image noise (d)

Full size image

It is worth mentioning that our method has a competitive computation time which is less than 2 min for one dataset on a computer with CPU Core2 Q9550 and 4G RAM using four CPU threads in the lung vessel removal and vesselness filter computation.

Conclusions and discussions

In this paper, we have presented a new method for fully automatic coronary tree extraction which is needed for clinical practice and related research of CCTA images. We first presented an improvement to Frangi’s vesselness filter by a response discriminator which can decrease the step-edge responses at cardiac cavities based on local geometric features. The presented extraction scheme based on this improved vesselness filter allows coronary arteries to be extracted with a fixed set of parameters in all of the 82 evaluated CCTA images.

Quantitative evaluations were performed in two phases. The first phase with the RCAA evaluation framework demonstrated good extraction ability and accuracy of our method in comparison with the other automatic methods. Worth mentioning is that the 95.9% overlap in clinical relevant parts approximates the highest measurement, 98.7%, obtained by an interactive method in the CAT08. The secondary phase with additional clinical datasets gave similar results, and at the same time showed that our method could be a reliable alternative of the user-interactive method in clinical practice. In addition, the robustness of our method was evaluated as well, since the images acquired by different CT scanners (Siemens and Toshiba scanners) and different image acquisition settings, such as contrast injection protocols, scanning parameters and reconstruction kernels, etc.

Benefiting from the improvements of Frangi’s vesselness filter and newly developed branch searching step, this new centerline extraction method has a higher extraction ability than our previous method [20], e.g. 93.7% vs. 78.8% for the overall overlap. The other steps except centerline refinement in this new method were also re-implemented to obtain better performance. For example, in the RCAA evaluation framework, aorta detection was successful in all the datasets with our new method while it failed for dataset No. 15 and 16 with the previous method.

The response discriminator presented is based on ray-casting which in general is considered to be a time-consuming operation. However, it is only performed for those voxels where Frangi’s vesselness response and the optimal scale are both larger than their respective thresholds. According to our statistics, only 5% of the voxels in a CCTA image are candidates for this discriminator. Furthermore the calculation of vesselness filter was performed in parallel by a multi-threading technique. As a consequence, the whole processing pipeline takes 2 min for one dataset on a standard desktop computer. The computational efficiency can be improved even further by a GPU-accelerated strategy.

As emphasized, a major advantage of our approach is that user interaction is unnecessary in most of the datasets. For example, ostium points were defined manually in only 4 of the 82 datasets. Benefiting from this automatic pipeline, the extraction can be performed as an offline pre-processing program on a server before cardiologists or radiologists begin their diagnosis. Even in those cases where the ostium detection failed, this pre-processing strategy can greatly reduce the time for coronary tree extraction because the most time-consuming part, i.e. vesselness filter, is already calculated beforehand. In these cases the coronary tree can be extracted almost in real-time after some simple mouse clicks through a user-friendly interface on the client computer.

Another advantage of our method is that the whole coronary tree can be retrieved simultaneously including the main branches and side branches. Using a user-interactive method to extract such a complete coronary tree would require a lot of user interactions. The anatomical and pathological information provided by this coronary tree is helpful in the diagnosis of CAD or the planning of cardiac interventions. Furthermore, it can reduce the processing time or improve the results of other related research work, such as vascular shear stress estimation [5], image fusion between CCTA and IVUS datasets [6] and optimal angle estimation for X-ray angiography [24].

A remark can be made that the semi-automatic method used for generating reference centerlines in the second evaluation phase resembled the centerline refinement in our centerline extraction algorithm. This might introduce a possible bias into the reference centerlines and the results. In addition, the evaluation in the second phase was limited to the parts of coronary arteries with clinical interest. Thus, the distal part extracted by our automatic method was not evaluated. However, it does show that our new fully automated method performs at least as good as our old semi-automated method.

With respect to the centerline refinement step it is seen that the accuracy of the centerlines is in general improved by using the detected lumen contours. However, sometimes the boundary of vessel lumen cannot be detected precisely especially in the region with lesions or severe calcifications. In these situations, the measurements of centerline accuracy, and sometimes even the overlap measurements, decreases. Future improvements to the contour detection could alleviate this problem.

In conclusion, we present a fully automatic centerline extraction algorithm for coronary arteries in CCTA image which is mainly based on an improved Frangi’s vesselness filter. Quantitative evaluations show that our method is able to extract the coronary arteries with high overlap and accuracy measurements. This automatic extraction algorithm has promising potential for both the clinical practice and the related research work.