Document Detail

Automatic centerline extraction of coronary arteries in coronary computed tomographic angiography.
Jump to Full Text
MedLine Citation:
PMID:  21637981     Owner:  NLM     Status:  MEDLINE    
Abstract/OtherAbstract:
Coronary computed tomographic angiography (CCTA) is a non-invasive imaging modality for the visualization of the heart and coronary arteries. To fully exploit the potential of the CCTA datasets and apply it in clinical practice, an automated coronary artery extraction approach is needed. The purpose of this paper is to present and validate a fully automatic centerline extraction algorithm for coronary arteries in CCTA images. The algorithm is based on an improved version of Frangi's vesselness filter which removes unwanted step-edge responses at the boundaries of the cardiac chambers. Building upon this new vesselness filter, the coronary artery extraction pipeline extracts the centerlines of main branches as well as side-branches automatically. This algorithm was first evaluated with a standardized evaluation framework named Rotterdam Coronary Artery Algorithm Evaluation Framework used in the MICCAI Coronary Artery Tracking challenge 2008 (CAT08). It includes 128 reference centerlines which were manually delineated. The average overlap and accuracy measures of our method were 93.7% and 0.30 mm, respectively, which ranked at the 1st and 3rd place compared to five other automatic methods presented in the CAT08. Secondly, in 50 clinical datasets, a total of 100 reference centerlines were generated from lumen contours in the transversal planes which were manually corrected by an expert from the cardiology department. In this evaluation, the average overlap and accuracy were 96.1% and 0.33 mm, respectively. The entire processing time for one dataset is less than 2 min on a standard desktop computer. In conclusion, our newly developed automatic approach can extract coronary arteries in CCTA images with excellent performances in extraction ability and accuracy.
Authors:
Guanyu Yang; Pieter Kitslaar; Michel Frenay; Alexander Broersen; Mark J Boogers; Jeroen J Bax; Johan H C Reiber; Jouke Dijkstra
Publication Detail:
Type:  Journal Article; Research Support, Non-U.S. Gov't; Validation Studies     Date:  2011-06-03
Journal Detail:
Title:  The international journal of cardiovascular imaging     Volume:  28     ISSN:  1875-8312     ISO Abbreviation:  Int J Cardiovasc Imaging     Publication Date:  2012 Apr 
Date Detail:
Created Date:  2012-05-28     Completed Date:  2012-09-19     Revised Date:  2013-06-28    
Medline Journal Info:
Nlm Unique ID:  100969716     Medline TA:  Int J Cardiovasc Imaging     Country:  United States    
Other Details:
Languages:  eng     Pagination:  921-33     Citation Subset:  IM    
Affiliation:
Division of Image Processing, Department of Radiology, Leiden University Medical Center, Leiden, The Netherlands. gyyang1980@gmail.com
Export Citation:
APA/MLA Format     Download EndNote     Download BibTex
MeSH Terms
Descriptor/Qualifier:
Algorithms*
Automation
Coronary Angiography / methods*
Coronary Artery Disease / pathology,  radiography*
Coronary Vessels / pathology*
Discriminant Analysis
Humans
Netherlands
Predictive Value of Tests
Radiographic Image Interpretation, Computer-Assisted / methods*
Registries
Reproducibility of Results
Tomography, X-Ray Computed / methods*
Comments/Corrections

From MEDLINE®/PubMed®, a database of the U.S. National Library of Medicine

Full Text
Journal Information
Journal ID (nlm-ta): Int J Cardiovasc Imaging
Journal ID (iso-abbrev): Int J Cardiovasc Imaging
ISSN: 1569-5794
ISSN: 1875-8312
Publisher: Springer Netherlands, Dordrecht
Article Information
Download PDF
© The Author(s) 2011
Received Day: 3 Month: 2 Year: 2011
Accepted Day: 19 Month: 5 Year: 2011
Electronic publication date: Day: 3 Month: 6 Year: 2011
pmc-release publication date: Day: 3 Month: 6 Year: 2011
Print publication date: Month: 4 Year: 2012
Volume: 28 Issue: 4
First Page: 921 Last Page: 933
ID: 3360862
PubMed Id: 21637981
Publisher Id: 9894
DOI: 10.1007/s10554-011-9894-2

Automatic centerline extraction of coronary arteries in coronary computed tomographic angiography
Guanyu Yang1 Address: gyyang1980@gmail.com
Pieter Kitslaar12
Michel Frenay1
Alexander Broersen1
Mark J. Boogers34
Jeroen J. Bax3
Johan H. C. Reiber1
Jouke Dijkstra1
1Division of Image Processing, Department of Radiology, Leiden University Medical Center, Leiden, The Netherlands
2Medis Medical Imaging Systems B.V., Leiden, The Netherlands
3Department of Cardiology, Leiden University Medical Center, Leiden, The Netherlands
4Interuniversity Cardiology Institute of the Netherlands, Utrecht, The Netherlands

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 [812]. 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.


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.

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 NS 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 [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ |{I(\user2{p})} - I(\user2{x})\left| {/I(\user2{x})} \right.> \varepsilon $$\end{document}] , 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.

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

[Formula ID: Equ1]
[\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ G(\user2{x}) = \prod\limits_{j = 1}^{M} {F(r_{j} )||\user2{v}_{j} ||} $$\end{document}]
in which, j is the group index. rj equals to Nj/NS (Nj is the number of rays in the j-th group). vj is the average vector of the normalized direction vectors of the rays in the j-th group and ||·|| is the Euclidean norm. The term ||vj|| 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,
[Formula ID: Equ2]
[\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ F(r) = 1 - \frac{1}{{1 + e^{ - \rho (r - \gamma )} }}, \, r \in (0,1] $$\end{document}]
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 γ·NS. 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 VF(x) and the optimal scale σ(x) are above thresholds TF and σF respectively. Thus, the improved vesselness response VI (x) can be written as

[Formula ID: Equ3]
[\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ 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. $$\end{document}]
TF 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 TC or dC respectively. In this paper, TC and dC 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].


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 (VI (x) < TI) 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 c2 and c5 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,

[Formula ID: Equ4]
[\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ 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 $$\end{document}]
in which tSE is the arrival time of the wave propagation from the start to the end points. dS and dE are the normalized tangential directions at the start and end points, dSE is the normalized direction of the found connection. Finally, [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ t_{SE}^{\max } $$\end{document}] is the maximum tSE 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.


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


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 NS = 429, ρ = 40, γ = 0.1, TF = 50, σF = 2, TI = 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.


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.


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. 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. Based on this path line, lumen contours in transversal planes along the vessel were initialized by a dedicated lumen contour detector [12].
  3. An expert from the cardiology department manually corrected the lumen contours if required.
  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.

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.

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.


This work is supported by the Dutch Technology Foundation STW (Utrecht, the Netherlands), applied science division of NWO, and the Technology Program of the Ministry of Economic Affairs, grant No. 10084. Mark J. Boogers is supported by the Dutch Heart Foundation, grant No. 2006T102.

Conflict of interest

None.

Open Access

This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.


References
1.. Mark DB,Berman DS,Budoff MJ,Carr JJ,Gerber TC,Hecht HS,Hlatky MA,Hodgson JM,Lauer MS,Miller JM,Morin RL,Mukherjee D,Poon M,Rubin GD,Schwartz RS. ACCF/ACR/AHA/NASCI/SAIP/SCAI/SCCT 2010 expert consensus document on coronary computed tomographic angiography: a report of the American College of Cardiology Foundation Task Force on Expert Consensus DocumentsJ Am Coll CardiolYear: 201055232663269910.1016/j.jacc.2009.11.01320513611
2.. Mowatt G,Cook J,Hillis G,Walker S,Fraser C,Jia X,Waugh N. 64-Slice computed tomography angiography in the diagnosis and assessment of coronary artery disease: systematic review and meta-analysisHeartYear: 2008941386139310.1136/hrt.2008.14529218669550
3.. Raff GL,Abidov A,Achenbach S,Berman DS,Boxt LM,Budoff MJ,Cheng V,DeFrance T,Hellinger JC,Karlsberg RP,Society of Cardiovascular Computed TomographySCCT guidelines for the interpretation and reporting of coronary computed tomographic angiographyJ Cardiovasc Comput TomogrYear: 20093212213610.1016/j.jcct.2009.01.00119272853
4.. Wesarg S,Khan MF,Firle EA. Localizing calcifications in cardiac CT data sets using a new vessel segmentation approachJ Digit ImagingYear: 200619324925710.1007/s10278-006-9947-616741661
5.. Giessen AG,Schaap M,Gijsen FJH,Groen HC,Walsum T,Mollet NR,Dijkstra J,Vosse FN,Niessen WJ,Feyter PJ,Steen AFW,Wentzel JJ. 3D fusion of intravascular ultrasound and coronary computed tomography for in vivo wall shear stress analysis: a feasibility studyInt J Cardiovasc ImagingYear: 201026778179610.1007/s10554-009-9546-y19946749
6.. Marquering H,Dijkstra J,Besnehard Q,Duthé J,Schuijf J,Bax J,Reiber JHC. Coronary CT angiography: IVUS image fusion for quantitative plaque and stenosis analysesProc SPIEYear: 2008691869181G10.1117/12.772582
7.. Lesage D,Angelini ED,Bloch I,Funka-Lea G. A review of 3D vessel lumen segmentation techniques: models, features and extraction schemesMed Image AnalYear: 200913681984510.1016/j.media.2009.07.01119818675
8.. Friman O,Hindennach M,Kühnel C,Peitgen HO. Multiple hypothesis template tracking of small 3D vessel structuresMed Image AnalYear: 200914216017110.1016/j.media.2009.12.00320060770
9.. Szymczak A,Tannenbaum A,Mischaikow K. Coronary vessel cores from 3D imagery: a topological approachMed Image AnalYear: 200610450551310.1016/j.media.2006.05.002
10.. Metz CT,Schaap M,Weustink AC,Mollet NR,Walsum T,Niessen NJ. Coronary centerline extraction from CT coronary angiography images using a minimum cost path approachMed PhysYear: 200936125568557910.1118/1.325407720095269
11.. Lesage D,Angelini E,Bloch I,Funka-Lea G. Bayesian maximal paths for coronary artery segmentation from 3D CT angiogramsProc Med Image Comput Assist Interv (MICCAI)Year: 20095761222229
12.. Marquering HA,Dijkstra J,Koning PJ,Stoel BC,Reiber JHC. Towards quantitative analysis of coronary CTAInt J Cardiovasc ImagingYear: 2005211738410.1007/s10554-004-5341-y15915942
13.. Bouraoui B, Ronse C, Baruthio J, Passat N, Germain P (2008) Fully automatic 3D segmentation of coronary arteries based on mathematical morphology. In: IEEE international symposium biomedical imaging: from nano to macro (ISBI) 1059–1062
14.. Zambal S, Hladuvka J, Kanitsar A, Buhler K (2008) Shape and appearance models for automatic coronary artery tracking. The Midas Journal—2008 MICCAI Workshop Grand Challenge Coronary Artery Tracking. http://hdl.handle.net/10380/1420. Accessed 9 Dec 2010
15.. Tek H, Gulsun M, Laguitton S, Grady L, Lesage D, Funka-Lea G (2008) Automatic coronary tree modeling. The Midas Journal—2008 MICCAI Workshop Grand Challenge Coronary Artery Tracking. http://hdl.handle.net/10380/1426. Accessed 9 Dec 2010
16.. Wang C, Smedby Ö (2008) An automatic seeding method for coronary artery segmentation and skeletonization in CTA. The Midas Journal—2008 MICCAI Workshop Grand Challenge Coronary Artery Tracking. http://hdl.handle.net/10380/1434. Accessed 9 Dec 2010
17.. Bauer C, Bischof H (2008) Edge based tube detection for coronary artery centerline extraction. The Midas Journal—2008 MICCAI Workshop Grand Challenge Coronary Artery Tracking. http://hdl.handle.net/10380/1403. Accessed 9 Dec 2010
18.. Xu C,Prince JL. Snakes, shapes, and gradient vector flowIEEE Trans Image ProcYear: 19987335936910.1109/83.661186
19.. Frangi AF, Niessen WJ, Vincken KL, Viergever MA (1998) Multiscale vessel enhancement filtering. In: Proceedings of medical image comput assisted intervention (MICCAI). Lecture notes in computer science 1496: 130–137
20.. Kitslaar PH, Frenay M, Oost E, Dijkstra J, Stoel B, Reiber JHC (2008) Connected component and morphology based extraction of arterial centerlines of the heart (CocomoBeach). The Midas Journal—2008 MICCAI Workshop Grand Challenge Coronary Artery Tracking. http://hdl.handle.net/10380/1460. Accessed 9 Dec 2010
21.. Schaap M,Metz C,et al. Standardized evaluation methodology and reference database for evaluating coronary artery centerline extraction algorithmsMed Image AnalYear: 200913570171410.1016/j.media.2009.06.00319632885
22.. Sofka M,Stewart CV. Retinal vessel centerline extraction using multiscale matched filters, confidence and edge measureIEEE Trans Med ImagingYear: 200625121531154610.1109/TMI.2006.88419017167990
23.. Boogers MJ,Schuijf JD,et al. Automated quantification of stenosis severity on 64-slice CT: a comparison with quantitative coronary angiographyJ Am Coll Cardiol ImgYear: 20103699709
24.. Kitslaar PH,Marquering H,Jukema W,Koning G,Nieber M,Vossepoel A,Bax J,Reiber JHC. Automated determination of optimal angiographic viewing angles for coronary artery bifurcations from CTA dataProc SPIEYear: 2008691869181J10.1117/12.770255

Figures

[Figure ID: Fig1]
Fig. 1 

Processing pipeline of the coronary tree extraction algorithm



[Figure ID: Fig2]
Fig. 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]



[Figure ID: Fig3]
Fig. 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



[Figure ID: Fig4]
Fig. 4 

Illustration of ostia detection step. Two voxels with the maxima vesselness responses, OL and OR, 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



[Figure ID: Fig5]
Fig. 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



[Figure ID: Fig6]
Fig. 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)



[Figure ID: Fig7]
Fig. 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



[Figure ID: Fig8]
Fig. 8 

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



[Figure ID: Fig9]
Fig. 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



[Figure ID: Fig10]
Fig. 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)



Tables
[TableWrap ID: Tab1] Table 1 

Comparison to five automatic extraction algorithms evaluated in the CAT08 [21]


OV (%) OF (%) OT (%) Acc. (mm) Time
DepthFirstModelFit 84.7 65.3 87.0 0.28 4–8 min
AutoCoronaryTree 84.7 59.5 86.2 0.34 <30 s
GVFTube’n’linkage 92.7 71.9 95.3 0.37 10 min
CocomoBeach 78.8 64.4 81.2 0.29 70 s
VirtualContrast 75.6 56.1 78.7 0.39 5 min
Our method 93.7 74.2 95.9 0.30 2 min


Article Categories:
  • Original Paper

Keywords: Keywords Coronary computed tomographic angiography (CCTA), Coronary artery, Centerline, Automatic extraction.

Previous Document:  Optimization of energy level for coronary angiography with dual-energy and dual-source computed tomo...
Next Document:  Right ventricular function declines after cardiac surgery in adult patients with congenital heart di...