Towards Accurate and Complete Registration of Coronary Arteries in CTA Images

Abstract

Coronary computed tomography angiography (CCTA) has been widely used nowadays. By combining multiple intra-subject CCTA images from different dates or different phases, cardiologists can monitor the disease progress, and researchers can explore the rules of coronary artery motion and changes within a cardiac cycle. For direct comparison and high efficiency, alignment of arteries is necessary. In this paper, we propose an automated method for accurate and complete registration of coronary arteries. Our method includes bifurcation matching, segment registration, and a novel approach to further improve the completeness of registration by combining the previous results and a level set algorithm. Our method is evaluated using 36 CCTA image pairs captured at different dates or different phases. The average distance error is $0.044 \pm 0.008$ mm and the average correct rate of registration is 90.7%.

1 Introduction

Coronary artery disease (CAD) is one of the leading causes for death in the whole world. Coronary computed tomography angiography (CCTA) examinations are growing in use as a non-invasive method for detecting narrowing, calcifications and fatty deposits in coronary arteries, with the ability to display anatomical details. By combining multiple CCTA images of the same patient at different dates (e.g. initial and follow-up visits), cardiologists can monitor CAD progress. Furthermore, based on the combination of CCTA images at different phases within a cardiac cycle (e.g. end-systole and end-diastole), many valuable medical researches can be done, such as exploring the rules of coronary artery motion and lumen changes. Coronary arteries in intra-subject CCTA images, even within the same cardiac cycle, vary in shape and posture due to factors such different relative positions, cardiac motions and respiratory motions. Therefore, comparative observation and analysis for multiple CCTA images is time-consuming, and an automated registration method of coronary arteries is beneficial.

Medical image registration has been studied a lot in recent years, but most of them are on large organs such as brains and lungs [7]. Existing methods of coronary artery registration take segmented coronary arteries or extracted artery centerlines as objects, rather than whole CCTA images, because the tree-like arteries with small lumen and wide spatial distribution will be severely disturbed by surrounding structures. In [3], point sets sampled on artery centerlines are registered based on the coherent point drift (CPD) algorithm and depending on the temporal consistency of multiple phases covering the whole cardiac cycle. The approach in [6] is limited to register straightened centerlines of left anterior descending artery (LAD), left circumflex artery (LCX) and right coronary artery (RCA) separately. The methods in [3, 6] are validated on only 3 and 5 real cases. [4] uses a bigger dataset of 26 cases, whose approach consists of an affine transform and a cubic B-spline transform. In [1], the coronary artery registration is converted to aligning nearby cardiac surface landmarks, so the accuracy highly relies on the performance of heart segmentation and landmark detection. Most studies [1, 3, 4] focus on cases of multiple phases, only [6] considers images from two dates. The registration accuracy reported in these work is approximately 1–2 mm, measured by average distance of centerlines or manual landmarks. A main challenge for registration is topological difference such as missing branches in segmentation results, which may undermine some methods. Another limitation of these work is that only the common parts of arteries are considered. This is undesirable because missed arteries are usually pathological and important.

In this paper, we propose an automated method to register coronary arterial trees extracted from CCTA images, which can deal with two images from different dates or phases, and can also be extended to register multiple images. The method focuses on both accuracy and completeness of registration. The arterial trees used as input can be obtained from any centerline extraction method. Bifurcation matching is the first and fundamental step, which determines the correspondence of centerline segments. And then correctly matching the sampling points on each segment pair contributes largely to accuracy. When input arteries have different topologies, level set segmentation is used to find missing arteries in one image with the guidance of unmatched arteries in the other image. The registration process is then repeated using the new segmented arteries. The performance is assessed in terms of accuracy and completeness, on a dataset of 36 pairs of CCTA images from 36 patients, among which 9 pairs are acquired with an interval of 3–7 years and 27 pairs are from end-systole and end-diastole.

2 Methods

In this work, the coronary arteries are segmented using a method combining thresholding and region growing. Then the centerlines used as the input of registration are extracted using minimal cost path algorithm. As seen in Fig. 1, the registration is composed of three main steps, including matching coronary bifurcations, matching sampling points on the centerlines, whose results are utilized to compute the parameters of a Thin Plate Spline (TPS) model [2], and further improving the completeness of segmentation and registration.

Fig. 1.

The flow chart of registration. Step1 is bifurcation matching. Step2 is segment registration. Step3 is further registration. The source and aligned centerlines are shown in blue, the target centerline in green, the matched parts after step2/step3 in red/cyan. The matched/unmatched bifurcations are displayed by black lines/circles.

2.1 Bifurcation Matching

In the first step, the bifurcations on target and source centerlines undergo a point matching process, with properties of location and direction. The direction is defined as the tangent direction of branch before the bifurcation and pointing from the proximal to distal. The matching method takes both similarity and compatibility of pairs into account.

Consider two bifurcation sets: set $B_\mathrm {s}$ with $n_\mathrm {s}$ bifurcations on source centerline and set $B_\mathrm {t}$ with $n_\mathrm {t}$ bifurcations on target centerline. Based on the 3D scale-invariant feature transform (SIFT) [8] descriptors, for each point $b_\mathrm {t}\in B_\mathrm {t}$, we find a point $b_\mathrm {s}\in B_\mathrm {s}$ with the shortest Euclidean distance from $b_\mathrm {t}$, thus we have $n_\mathrm {t}$ point pairs. Then a correspondence graph G(V, E, M) is constructed, in which the node set V represents bifurcation pairs, the edge set E consists of relations between pairs, and M is an attribute matrix. The attribute value M(i, j) assigned to edge $e_{ij}\in E (i\ne j)$ reflects the compatibility of ith and jth pairs, and M(i, i) indicates the similarity of ith pair. Since correct pairs tend to establish links among each other, the matching problem can be formulated as a problem of detecting the main strongly connected cluster in graph G, which can be solved by eigenvector technique [5].

M is related to three geometric parameters. As illustrated in Fig. 2, d is the Euclidean distance between bifurcations, $\theta$ is the relative angle of bifurcation directions, and $\varphi$ is the angle of one bifurcation measured when taking the other as the polar axis. M is defined as:

$$\begin{aligned}&M(i,j)={\left\{ \begin{array}{ll} Z(i,j),&{}i\ne j\\ S(i,j),&{}i=j \end{array}\right. } ,\end{aligned}$$

(1)

$$\begin{aligned}&Z(i,j)=\frac{ratio(i,j)-TH_{\mathrm {dist}}}{1-TH_{\mathrm {dist}}}\times \left( 1-\frac{|\theta _i-\theta _j|}{TH_{\theta }} \right) \times \left( 1-\frac{|\varphi _i-\varphi _j|}{TH_{\varphi }} \right) ,\end{aligned}$$

(2)

$$\begin{aligned}&ratio(i,j)=\frac{\min \{d_i,d_j\}}{\max \{d_i,d_j\}} ,\end{aligned}$$

(3)

$$\begin{aligned}&S(i,j)={\left\{ \begin{array}{ll} 1-\frac{\Vert Des_i^\mathrm {s}-Des_i^\mathrm {t}\Vert }{\Vert Des_i^\mathrm {s}\Vert +\Vert Des_i^\mathrm {t}\Vert },&{}i=j\\ 0,&{}i\ne j \end{array}\right. } , \end{aligned}$$

(4)

where $TH_{\mathrm {dist}}$, $TH_{\theta }$ and $TH_{\varphi }$ are thresholds of the minimum ratio(i, j), the maximum $|\theta _i-\theta _j|$ and the maximum $|\varphi _i-\varphi _j|$, respectively. And Z(i, j) is small when any threshold is approached, representing bad compatibility. $Des_i^\mathrm {s}$ and $Des_i^\mathrm {t}$ are SIFT descriptors of ith pair on the source and target centerlines. Here, the SIFT descriptor is composed of $4\times 4\times 4$ orientation histograms each with 12 bins, thus having 768 elements.

Fig. 2.

Illustration of three parameters involved in the attribute matrix M. $m_{i1}$ and $m_{i2}$ are two bifurcations with direction.

Here we define a segment as a tree section between two bifurcations or between a bifurcation and an end point (or root point). The input trees may be different in topology due to incomplete segmentation. As illustrated in Fig. 3(a), centerline $C_\mathrm {s}$ has bifurcation $b_2$ while $C_\mathrm {t}$ has not. Our algorithm removes $b_2$ and related segment $c_5$ as well as merging $c_3$ and $c_4$, and then $C_\mathrm {s}$ is updated with new topology.

Note that which segment to be removed is decided by the continuity and smoothness between $c_5$, $c_4$ and $c_3$. Figure 3(b) shows a real case in which three segments will be removed one by one from the source centerline. All the unmatched bifurcations are processed like this before segment matching, in which segments of two centerlines were considered matched if they share two bifurcation pairs or have one bifurcation pair and similar directions. The unmatched segments will skip the step of segment registration, and undergo the TPS transform directly.

Fig. 3.

Centerlines with different topologies. In (b) the segments to be removed are shown in red.

2.2 Segment Registration

For each segment pair determined as above, refined registration is applied. Segment $Seg_\mathrm {s}$ on source centerline and segment $Seg_\mathrm {t}$ on target centerline are discretized by evenly sampling points, with intervals of $\triangle _\mathrm {s}$ and $\triangle _\mathrm {t}$, satisfying $\triangle _\mathrm {t}=K\triangle _\mathrm {s}$ ($K=10$ in our work). Suppose that $Seg_\mathrm {s}$ and $Seg_\mathrm {t}$ have $N_\mathrm {s}$ and $N_\mathrm {t}$ points respectively, we define a deformation model by a series of status:

$$\begin{aligned} \mathcal {S} := \left\{ S(i)= j | i\in \{0,1,\dots ,N_\mathrm {t}\}, j\in \{0,1,\dots ,N_\mathrm {s}\}\right\} , \end{aligned}$$

(5)

where S(i) is the status of ith point on $Seg_\mathrm {t}$. $S(i)=j$ means that ith point on $Seg_\mathrm {t}$ corresponds to jth point on $Seg_\mathrm {s}$. Here $i\in \{0,1,\dots ,N_\mathrm {t}\}$ can be seen as a label assigned to points on $Seg_\mathrm {s}$. Note that some labels may be omitted when $Seg_\mathrm {s}$ is shorter than $Seg_\mathrm {t}$.

To determine the parameters (a series of index numbers on $Seg_\mathrm {s}$) of model $\mathcal {S}$, an objective function is defined and maximized. The function consists of image similarity and geometric similarity. We define the image similarity as:

$$\begin{aligned} Sim1(i) = 1-\frac{\left\| Des_i^\mathrm {t}-Des_{S(i)}^\mathrm {s}\right\| }{\left\| Des_i^\mathrm {t} \right\| + \left\| Des_{S(i)}^\mathrm {s}\right\| } , \end{aligned}$$

(6)

where $Des_i^\mathrm {t}$ and $Des_{S(i)}^\mathrm {s}$ are SIFT descriptors (computed as the same way of Sect. 2.1) of ith point on $Seg_\mathrm {t}$ and S(i)th point on $Seg_\mathrm {s}$ respectively. Geometric constraints are imposed to avoid excessive stretching or centerlines folding, and thus geometric similarity for ith point on $Seg_\mathrm {t}$ is defined as:

$$\begin{aligned} {}{} Sim2(i)={\left\{ \begin{array}{ll} 1-\left| \log \left( \frac{S(i)-S(i-1)}{K}\right) \right| ,&{}S(i)\ge S(i-1)\\ -\infty ,&{}S(i)< S(i-1) \end{array}\right. } , \end{aligned}$$

(7)

where $K=\triangle _\mathrm {t}/\triangle _\mathrm {s}$ is 10 as mentioned above. We use logarithm to limit the difference in spacing between corresponding points on two segments, due to its better performance than piecewise linear functions in experiments.

The global objective function is as follows:

$$\begin{aligned} J\left( \mathcal {S}\right) =\sum _{i=1}^{N_\mathrm {t}} \omega _1 Sim1(i)+\omega _2 Sim2(i) . \end{aligned}$$

(8)

Since each term of $J(\mathcal {S})$ related to i depends only on status S(i) and status $S(i-1)$, which satisfies Markov property, J can be viewed as a Hidden Markov Model and optimized using Viterbi algorithm [6]. Matching points from all the segment pairs together form the control points of TPS transform, which is used to align the coronary arterial trees.

2.3 Further Segmentation and Registration

After registering common parts of input arterial trees, the third step is to further improve the completeness of segmentation and registration. Assuming that a branch mask $BR_\mathrm {s}$ is part of the input segmented arterial tree $T_\mathrm {s}$ but has no correspondence on tree $T_\mathrm {t}$, we attempt to segment the missing (potential) branch on tree $T_\mathrm {t}$. Firstly, $BR_\mathrm {s}$ is transformed using the TPS model from the previous registration, resulting in a mask $BR_\mathrm {a}$ in the other image. We define two volumes of interest ($V\mathrm {s}$ and $V_\mathrm {t}$) surrounding $BR_\mathrm {s}$ and $BR_\mathrm {a}$ with a 10-pixel margin in each direction parallel to coordinate axes. A B-spline transform model is computed to register $V\mathrm {s}$ and $V_\mathrm {t}$, which is then applied to $BR_\mathrm {s}$, resulting in $BR_\mathrm {t}^0$. Next, a level set algorithm is implemented initialized by $BR_\mathrm {t}^0$, giving the segmentation result $BR_\mathrm {t}$. The level set algorithm implemented is the very efficient sparse field method proposed by Whitaker [9] and the well-known Chan-Vese energy is minimized. Some branches absent in original segmentation can be obtained in this way, because the prior information derived from the other image contributes to overcoming the difficulties such as fuzzy edges and false connection to surrounding structures. Some examples are shown in Fig. 4. However, the method may fail in regions where the image quality is too bad and even experts cannot distinguish the vessels.

Fig. 4.

The original segmentation (blue) and branches segmented by level set (red).

After attempts of further segmentation, new centerlines are extracted from arteries, and the registration process is performed again. The completeness of registration is improved because larger parts of trees can be registered.

3 Experiments and Results

We validated our method using 36 CCTA image pairs from 36 patients, among which 9 pairs (Group 1) are acquired with an interval of 3–7 years and 27 pairs (Group 2) are from end-systole and end-diastole. The method is evaluated from aspects of accuracy and completeness. The accuracy is assessed by the average Euclidian distance between the aligned centerline $C_\mathrm {a}$ and target centerline $C_\mathrm {t}$:

$$\begin{aligned} D\left( C_\mathrm {a},C_\mathrm {t}\right) = \frac{1}{2N_\mathrm {a}}&\sum _{i=1}^{N_\mathrm {a}}D\left( c_\mathrm {a}^i,C_\mathrm {t}\right) +\frac{1}{2N_\mathrm {t}}\sum _{i=1}^{N_\mathrm {t}}D\left( c_\mathrm {t}^i,C_\mathrm {a}\right) \!\!,\end{aligned}$$

(9)

$$\begin{aligned} D\left( c_\mathrm {a}^i,C_\mathrm {t}\right)&=\min \left\{ d\left( c_\mathrm {a}^i,t\right) :t\in C_\mathrm {t}\right\} \!\!,\end{aligned}$$

(10)

$$\begin{aligned} D\left( c_\mathrm {t}^i,C_\mathrm {a}\right)&=\min \left\{ d\left( c_\mathrm {t}^i,a\right) :a\in C_\mathrm {a}\right\} \!\!, \end{aligned}$$

(11)

where $c_\mathrm {a}^i\in C_\mathrm {a}$, $c_\mathrm {t}^i\in C_\mathrm {t}$, $N_\mathrm {a}$ and $N_\mathrm {t}$ are the numbers of dense sampling points on $C_\mathrm {a}$ and $C_\mathrm {t}$, and d means Euclidean distance. The completeness is measured by the correct rate of registration, defined as $R^\mathrm {reg} = \left( N^\mathrm {reg}_\mathrm {a}/N_\mathrm {a} + N^\mathrm {reg}_\mathrm {t}/N_\mathrm {t} \right) /2$, where $N^\mathrm {reg}_\mathrm {a}$ and $N^\mathrm {reg}_\mathrm {t}$ are the numbers of correctly matched points on $C_\mathrm {a}$ and $C_\mathrm {t}$, respectively. Here a point on one centerline is considered correctly matched when the closest point on the other centerline is within 0.5 mm. The preliminary and final results (before and after further registration) of index $D(C_\mathrm {a},C_\mathrm {t})$ and $R^\mathrm {reg}$ on the dataset are displayed in Table 1. The further registration improves the correct rate of registration. The slight improvement in distance is not surprising as additional registered arteries are usually more noisy.

Table 1. $E_\mathrm {dist}$ represents the average and standard deviation of $D(C_\mathrm {a},C_\mathrm {t})$. $E_\mathrm {reg}$ indicates the average of $R^\mathrm {reg}$. Superscript 1 and 2 indicate preliminary and final results.

Figure 5(a) and (b) demonstrate the centerlines after preliminary and further registration respectively in two cases. A pair of further registered arteries in case 1 is indicated by arrows in Fig. 5(c), using the presentation method of straighten curved planar reformation (SCPR). The artery segment in the right image acquired in 2012 is missed in initial segmentation due to the serious calcification at the bifurcation, but it is important for reflecting the CAD progress.

Fig. 5.

(a)(b) The preliminary (left) and further (right) registration results. The preliminary/further matched parts are shown in red/cyan, unmatched parts in blue/green for two trees. (c) A pair of further registered arteries in case 1 is indicated by arrows.

4 Discussion and Conclusion

Our work has two main contributions in registration of coronary arteries. One is the hierarchical strategy for registration, consisting of bifurcation matching and segment registration. Bifurcation matching as the fundamental step properly handles the problem of different topologies, and segment registration contributes to the final accuracy. The other contribution is the improvement in completeness of registration, which is ignored by other related studies. In view of the fact that registration and segmentation are closely coupled problems, a level set method is used to further segment the difficult arteries missed by the initial segmentation algorithm, based on the initial registration results and prior information from the other image where the corresponding arteries are present. Moreover, more complete segmentation results in turn improve the completeness of registration.

Although the control points for TPS are totally on centerlines, the transform can be applied to whole arteries. The registration results can not only assist cardiologists in monitoring CAD progress, but also contribute to exploring the rules of coronary artery motion and changes within a cardiac cycle.

In conclusion, the method proposed for registration of coronary arteries in CCTA images performs well in aspects of accuracy and completeness, which is validated using 36 pairs of intra-subject CCTAs images captured at different dates or different phases. In future work, we will assess the method on larger datasets, and improve the algorithm for registering noisy arteries.