4.1. Cardiac CTA data

The CTA data was acquired in the Erasmus MC, University Medical Center Rotterdam, The Netherlands. Thirty-two datasets were randomly selected from a series of patients who underwent a cardiac CTA examination between June 2005 and June 2006. Twenty datasets were acquired with a 64-slice CT scanner and twelve datasets with a dual-source CT scanner (Sensation 64 and Somatom Definition, Siemens Medical Solutions, Forchheim, Germany).

A tube voltage of 120 kV was used for both scanners. All datasets were acquired with ECG-pulsing [49]. The maximum current (625 mA for the dual-source scanner and 900 mA for the 64-slice scanner) was used in the window from 25% to 70% of the RR-interval and outside this window the tube current was reduced to 20% of the maximum current.

Both scanners operated with a detector width of 0.6 mm. The image data was acquired with a table feed of 3.8 mm per rotation (64-slice datasets) or 3.8 mm to 10 mm, individually adapted to the patient’s heart rate (dual-source datasets).

Diastolic reconstructions were used, with reconstruction intervals varying from 250 ms to 400 ms before the R-peak. Three datasets were reconstructed using a sharp (B46f) kernel, all others were reconstructed using a medium-to-smooth (B30f) kernel. The mean voxel size of the datasets is 0.32 × 0.32 × 0.4mm³.

4.2. Reference standard

In this work we define the centerline of a coronary artery in a CTA scan as the curve that passes through the center of gravity of the lumen in each cross-section. We define the start point of a centerline as the center of the coronary ostium (i.e. the point where the coronary artery originates from the aorta), and the end point as the most distal point where the artery is still distinguishable from the background. The centerline is smoothly interpolated if the artery is partly indistinguishable from the background, e.g. in case of a total occlusion or imaging artifacts.

This definition was used by three trained observers to annotate centerlines in the selected cardiac CTA datasets. Four vessels were selected for annotation by one of the observers in all 32 datasets, yielding 32 × 4 = 128 selected vessels. The first three vessels were always the right coronary artery (RCA), left anterior descending artery (LAD), and left circumflex artery (LCX). The fourth vessel was selected from the large side branches of these main coronary arteries and the selection was as follows: first diagonal branch (14x), second diagonal branch (6x), optional diagonal coronary artery (6x), first obtuse marginal branch (2x), posterior descending artery (2x), and acute marginal artery (2x). This observer annotated for all the four selected vessels points close to the selected vessels. These points (denoted with ‘point A’) unambiguously define the vessels, i.e. the vessel of interest is the vessel closest to the point and no side-branches can be observed after this point.

After the annotation of these 128 points, the three observers used these points to independently annotate the centerlines of the same four vessels in the 32 datasets. The observers also specified the radius of the lumen at least every 5 mm, where the radius was chosen such that the enclosed area of the annotated circle matched the area of the lumen. The radius was specified after the complete central lumen line was annotated (see Figure 4).

Open in a separate window

Figure 4

An example of the annotations of the three observers in black and the resulting reference standard in white. The crosses indicate the centers and the circles indicate the radii.

The paths of the three observers were combined to one centerline per vessel using a Mean Shift algorithm for open curves: The centerlines are averaged while taking into account the possibly spatially varying accuracy of the observers by iteratively estimating the reference standard and the accuracy of the observers. Each point of the resulting reference standard is a weighted average of the neighboring observer centerline points, with weights corresponding to the locally estimated accuracy of the observers [50].

After creating this first weighted average, a consensus centerline was created with the following procedure: The observers compared their centerlines with the average centerline to detect and subsequently correct any possible annotation errors. This comparison was performed utilizing curved planar reformatted images displaying the annotated centerline color-coded with the distance to the reference standard and vice-versa (see Figure 2). The three observers needed in total approximately 300 hours for the complete annotation and correction process.

Open in a separate window

Figure 2

An example of one of the color-coded curved planar reformatted images used to detect possible annotation errors.

After the correction step the centerlines were used to create the reference standard, using the same Mean Shift algorithm. Note that the uncorrected centerlines were used to calculate the inter-observer variability and agreement measures (see section 4.5).

The points where for the first time the centerlines of two observers lie within the radius of the reference standard when traversing over this centerline from respectively the start to the end or vice versa were selected as the start- and end point of the reference standard. Because the observers used the abovementioned centerline definition it is assumed that the resulting start points of the reference standard centerlines lie within the coronary ostium.

The corrected centerlines contained on average 44 points and the average distance between two successive annotated points was 3.1 mm. The 128 resulting reference standard centerlines were on average 138 mm (std. dev. 41 mm, min. 34 mm, max. 249 mm) long.

The radius of the reference standard was based on the radii annotated by the observers and a point-to-point correspondence between the reference standard and the three annotated centerlines. The reference standard centerline and the corrected observer centerlines were first resampled equidistantly using a sampling distance of 0.03 mm. Dijkstra’s graph searching algorithm was then used to associate each point on the reference standard with one or more points on each annotated centerline and vice versa. Using this correspondence, the radius at each point of the reference standard was determined by averaging the radius of all the connected points on the three annotated centerlines (see also Figure 3 and Figure 4). An example of annotated data with corresponding reference standard is shown in Figure 1. Details about the connectivity algorithm are given in section 4.3.

Open in a separate window

Figure 1

An example of the data with corresponding reference standard. Top-left: axial view of data. Top-right: coronal view. Bottom-left: sagittal view. Bottom-right: a 3D rendering of the reference standard.

Open in a separate window

Figure 3

An illustrative example of the Mean Shift algorithm showing the annotations of the three observers as a thin black line, the resulting average as a thick black line, and the correspondence that are used during the last Mean Shift iteration in light-gray.

4.3. Correspondence between centerlines

All the evaluation measures are based on a point-to-point correspondence between the reference standard and the evaluated centerline. This section explains the mechanism for determining this correspondence.

Before the correspondence is determined the centerlines are first sampled equidistantly using a sampling distance of 0.03 mm, enabling an accurate comparison. The evaluated centerline is then clipped with a disc that is positioned at the start of the reference standard centerline (i.e. in or very close to the coronary ostium). The centerlines are clipped because we define the start point of a coronary centerline at the coronary ostium and because for a variety of applications the centerline can start somewhere in the aorta. The radius of the disc is twice the annotated vessel radius and the disc normal is the tangential direction at the beginning of the reference standard centerline. Every point before the first intersection of a centerline and this disc is not taken into account during evaluation.

The correspondence is then determined by finding the minimum of the sum of the Euclidean lengths of all point-point connections that are connecting the two centerlines over all valid correspondences. A valid correspondence for centerline I, consisting of an ordered set of points p_i (0 ≤ i < n, p₀ is the most proximal point of the centerline), and centerline II, consisting of an ordered set of points q_j (0 ≤ j < m, q₀ is the most proximal point of the centerline), is defined as the ordered set of connections C = {c₀, …, c_n+m−1}, where c_k is a tuple [p_a, q_b] that represents a connection from p_a to q_b, which satisfies the following conditions:

• The first connection c₀ connects the start points: c₀ = [p₀, q₀].

• The last connection c_n+m−1 connects the end points: c_n+m−1 = [p_n−1, q_m−1].

• If connection c_k = [p_a, q_b] then connection c_k+1 equals either [p_a+1, q_b] or [p_a, q_b+1].

These conditions guarantee that each point of centerline I is connected to at least one point of centerline II and vice versa.

Dijkstra’s graph search algorithm is used on a matrix with connection lengths to determine the minimal Euclidean length correspondence. See Figure 3 for an example of a resulting correspondence.

4.4. Evaluation measures

Coronary artery centerline extraction may be used for different applications, and thus different evaluation measures may apply. We account for this by employing a number of evaluation measures. With these measures we discern between extraction capability and extraction accuracy. Accuracy can only be evaluated when extraction succeeded; in case of a tracking failure the magnitude of the distance to the reference centerline is no longer relevant and should not be included in the accuracy measure.

4.4.1. Definition of true positive, false positive and false negative points

All the evaluation measures are based on a labeling of points on the centerlines as true positive, false negative or false positive. This labeling, in its turn, is based on a correspondence between the points of the reference standard centerline and the points of the centerline to be evaluated. The correspondence is determined with the algorithm explained in section 4.3.

A point of the reference standard is marked as true positive TPR_ov if the distance to at least one of the connected points on the evaluated centerline is less than the annotated radius and false negative FN_ov otherwise.

A point on the centerline to be evaluated is marked as true positive TPM_ov if there is at least one connected point on the reference standard at a distance less than the radius defined at that reference point, and it is marked as false positive FP_ov otherwise. With · we denote the cardinality of a set of points, e.g. TPR_ov denotes the number of reference points marked true positive. See also Figure 5 for a schematic explanation of these terms and the terms mentioned in the next section.

Open in a separate window

Figure 5

An illustration of the terms used in the evaluation measures (see section 4.4). The reference standard with annotated radius is depicted in gray. The terms on top of the figure are assigned to points on the centerline found by the evaluated method. The terms below the reference standard line are assigned to points on the reference standard.

4.5. Observer performance and scores

Each of the evaluation measures is related to the performance of the observers by a relative score. A score of 100 points implies that the result of the method is perfect, 50 points implies that the performance of the method is similar to the performance of the observers, and 0 points implies a complete failure. This section explains how the observer performance is quantified for each of the four evaluation measures and how scores are created from the evaluation measures by relating the measures to the observer performance.

4.5.1. Overlap measures

The inter-observer agreement for the overlap measures is calculated by comparing the uncorrected paths with the reference standard. The three overlap measures (OV, OF, OT) were calculated for each uncorrected path and the true positives, false positives and false negatives for each observer were combined into inter-observer agreement measures per centerline as follows:

$$\begin{array}{cc}\hfill {\mathrm{OV}}_{\mathrm{ag}}& =\frac{\sum (\Vert {\mathrm{TPR}}_{\mathrm{ov}}^i\Vert +\Vert {\mathrm{TPM}}_{\mathrm{ov}}^i\Vert )}{\sum (\Vert {\mathrm{TPR}}_{\mathrm{ov}}^i\Vert +\Vert {\mathrm{TPM}}_{\mathrm{a}}^i\Vert +\Vert {\mathrm{FP}}_{\mathrm{ov}}^i\Vert +\Vert {\mathrm{FN}}_{\mathrm{ov}}^i\Vert )}\hfill \\ \hfill {\mathrm{OF}}_{\mathrm{ag}}& =\frac{\sum \Vert {\mathrm{TPR}}_{\mathrm{of}}^i\Vert }{\sum (\Vert {\mathrm{TPR}}_{\mathrm{of}}^i\Vert +\Vert {\mathrm{FN}}_{\mathrm{of}}^i\Vert )}\hfill \\ \hfill {\mathrm{OT}}_{\mathrm{ag}}& =\frac{\sum (\Vert {\mathrm{TPR}}_{\mathrm{ot}}^i\Vert +\Vert {\mathrm{TPM}}_{\mathrm{ot}}^i\Vert )}{\sum (\Vert {\mathrm{TPR}}_{\mathrm{ot}}^i\Vert +\Vert {\mathrm{TPM}}_{\mathrm{ot}}^i\Vert +\Vert {\mathrm{FP}}_{\mathrm{ot}}^i\Vert +\Vert {\mathrm{FN}}_{\mathrm{ot}}^i\Vert )},\hfill \end{array}$$

where i = {0, 1, 2} indicates the observer.

After calculating the inter-observer agreement measures, the performance of the method is scored. For methods that perform better than the observers the OV, OF, and OT measures are converted to scores by linearly interpolating between 100 and 50 points, respectively corresponding to an overlap of 1.0 and an overlap similar to the inter-observer agreement value. If the method performs worse than the inter-observer agreement the score is obtained by linearly interpolating between 50 and 0 points, with 0 points corresponding to an overlap of 0.0:

$${\text{Score}}_{\mathrm{O}}=\{\begin{array}{cc}({\mathrm{O}}_{\mathrm{m}}∕{\mathrm{O}}_{\mathrm{ag}})\ast 50\hfill & {\mathrm{O}}_{\mathrm{m}}\le {\mathrm{O}}_{\mathrm{ag}}\hfill \\ 50+50\ast \frac{{\mathrm{O}}_{\mathrm{m}}-{\mathrm{O}}_{\mathrm{ag}}}{1-{\mathrm{O}}_{\mathrm{ag}}}\hfill & {\mathrm{O}}_{\mathrm{m}}>{\mathrm{O}}_{\mathrm{ag}},\hfill \end{array}\phantom{\}}$$

where O_m and O_ag define the OV, OF, or OT performance of respectively the method and the observer. An example of this conversion is shown in Figure 6(a).

Open in a separate window

Figure 6

Figure (a) shows an example of how overlap measures are transformed into scores. Figure (b) shows this transformation for the accuracy measures.

5.1. Category 1: automatic extraction

Automatic extraction methods find the centerlines of coronary arteries without user interaction. In order to evaluate the performance of automatic coronary artery centerline extraction, two points per vessel are provided to extract the coronary artery of interest:

• Point A: a point inside the distal part of the vessel; this point unambiguously defines the vessel to be tracked;

• Point B: a point approximately 3 cm (measured along the centerline) distal of the start point of the centerline.

Point A should be used for selecting the appropriate centerline. If the automatic extraction result does not contain centerlines near point A, point B can be used. Point A and B are only meant for selecting the right centerline and it is not allowed to use them as input for the extraction algorithm.

5.2. Category 2: extraction with minimal user interaction

Extraction methods with minimal user interaction are allowed to use one point per vessel as input for the algorithm. This can be either one of the following points:

• Point A or B, as defined above;

• Point S: the start point of the centerline;

• Point E: the end point of the centerline;

• Point U: any manually defined point.

Points A, B, S and E are provided with the data. Furthermore, in case the method obtains a vessel tree from the initial point, point A or B may be used after the centerline determination to select the appropriate centerline.

7.1. Fully-automatic methods

• AutoCoronaryTree [54, 55]: The full centerline tree of the coronary arteries is extracted via a multi-scale medialness-based vessel tree extraction algorithm which starts a tracking process from the ostia locations until all coronary branches are reached.

• CocomoBeach [56]: This method starts by segmenting the ascending aorta and the heart. Candidate coronary regions are obtained using connected component analysis and the masking of large structures. Using these components a region growing scheme, starting in the aorta, segments the complete tree. Finally, centerlines within the pre-segmented tree are obtained using the WaveProp [10] method.

• DepthFirstModelFit [57]: Coronary artery centerline extraction is accomplished by fitting models of shape and appearance. A large-scale model of the complete heart in combination with symmetry features is used for detecting coronary artery seeds. To fully extract the coronary artery tree, two small-scale cylinder-like models are matched via depth-first search.

• GVFTube’n’Linkage [58]: This method uses a Gradient Vector Flow [59] based tube detection procedure for identification of vessels surrounded by arbitrary tissues [60, 61]. Vessel centerlines are extracted using ridge-traversal and linked to form complete tree structures. For selection of coronary arteries gray value information and centerline length are used.

• VirtualContrast [62]: This method segments the coronary arteries based on the connectivity of the contrast agent in the vessel lumen, using a competing fuzzy connectedness tree algorithm [34]. Automatic rib cage removal and ascending aorta tracing are included to initialize the segmentation. Centerline extraction is based on the skeletonization of the tree structure.

7.2. Semi-automatic methods

• AxialSymmetry [63]: This method finds a minimum cost path connecting the aorta to a user supplied distal endpoint. Firstly, the aorta surface is extracted. Then, a two-stage Hough-like election scheme detects the high axial symmetry points in the image. Via these, a sparse graph is constructed. This graph is used to determine the optimal path connecting the user supplied seed point and the aorta.

• CoronaryTreeMorphoRec [64]: This method generates the coronary tree iteratively from point S. Pre-processing steps are performed in order to segment the aorta, remove unwanted structures in the background and detect calcium. Centerline points are chosen in each iteration depending on the previous vessel direction and a local gray scale morphological 3D reconstruction.

• KnowledgeBasedMinPath [65]: For each voxel, the probability of belonging to a coronary vessel is estimated from a feature space and a vesselness measure is used to obtain a cost function. The vessel starting point is obtained automatically, while the end point is provided by the user. Finally, the centerline is obtained as the minimal cost path between both points.

8.1. Results categorized on image quality, calcium score and vessel type

Separate rankings are made for each group of datasets with corresponding image quality and calcium rating to determine if the image quality or the amount of calcium has influence on the rankings. Separate rankings are also made for each of the four vessel types. These rankings are presented in Table 8. It can be seen that some of the methods perform relatively worse when the image quality is poor or an extensive amount of calcium is present (e.g. CocomoBeach [56] and DepthFirstModelFit [57]) and vice versa (e.g. KnowledgeBasedMinPath [65] and VirtualContrast [62]).

Table 8 also shows that on average the automatic methods perform relatively worse for datasets with poor image quality (i.e. the ranks of the automatic methods in the P-column are on average higher compared to the ranks in the M- and G-column). This is also true for the extraction of the LCX centerlines. Both effects can possibly be explained by the fact that centerline extraction from poor image quality datasets and centerline extraction of the (on average relatively thinner) LCX is more difficult to automate.

8.2. Algorithm performance with respect to distance from the ostium

For a number of coronary artery centerline extraction applications it is not important to extract the whole coronary artery; only extraction up to a certain distance from the coronary ostium is required (see e.g. [73, 74]).

In order to evaluate the performance of the methods with respect to the distance from the ostium, charts are generated that demonstrate the average performance over all 96 evaluated centerlines for each of the methods at a specific distance from the ostium (measured along the reference standard). Figure 7(a) shows these results for the automatic methods, Figure 7(b) shows the results for the methods with minimal user-interaction, and Figure 7(c) shows the results for the semi-automatic methods.

Open in a separate window

Figure 7

The algorithm performance of each method with respect to the distance from the ostium averaged over all 96 evaluated vessels over the first 175mm (only 10% of the vessels were longer than 175mm). Overlap: the fraction of points on the reference standard marked as true positive. Accuracy: the average distance to the centerline if the point is marked true positive. Each of the three graphs shows in light-gray the results of all the thirteen evaluated methods and in color the results of the respective algorithm category. The graphs also show in black the average accuracy and overlap for all thirteen evaluated methods.

The graphs show that all the evaluated methods are better able to extract the proximal part of the coronaries than the more distal part of the vessels.

They also show that after approximately 5 cm the accuracy of almost all the methods is relatively constant. Furthermore, the graphs again demonstrate the fact that the automatic methods are on average more accurate than the semi-automatic or interactive methods.

Michiel Schaap, Coert Metz, Theo van Walsum and Wiro Niessen are supported by the Stichting voor de Technische Wetenschappen (STW) of The Netherlands Organization for Scientific Research (NWO). Christian Bauer was supported by the Austrian Science Fund (FWF) under the doctoral program Confluence of Vision and Graphics W1209. The work of Hrvoje Bogunović has been partially funded by the Industrial and Technological Development Centre (CDTI) under the CENIT Programme (CDTEAM Project), the EC @neurIST (IST-FP6-2004-027703) project and the Networking Research Center on Bioengineering, Biomaterials and Nanomedicine (CIBER-BBN). Michel Frenay was supported by a research grant from Bio-Imaging Technologies. Marcela Hernández Hoyos and Maciej Orkisz were supported by ECOS-Nord #C07M04 and Region Rhone-Alpes PP3/I3M of cluster ISLE. Pieter Kitslaar was supported by innovation grant IS044070(ADVANCE) from the Dutch Ministry of Economic Affairs. Karl Krissian is funded by the Spanish government and the University of Las Palmas of Gran Canaria under the Ramon y Cajal program and he would like to acknowledge J. Pozo, M. Villa-Uriol, and A. Frangi for their contribution to this work. Luengo-Oroz and Carlos Castro would like to acknowledge Maria Jesus Ledesma-Carbayo for her supervision and the work of Luengo-Oroz and Carlos Castro was supported by the research projects TIN2007-68048-C02-01 from the Spanish Ministry of Education and Science and the CDTEAMproject from the Spanish Ministry of Industry (CDTI). Simon K. Warfield was supported in part by NIH grants R01 RR021885, R01 GM074068 and R01 EB008015. The work of Zambal et al. was in part funded by Agfa HealthCare and in part by the Austrian Kplus funding program.