Local motion-compensated method for high-quality 3D coronary artery reconstruction

Bo Liu; Xiangzhi Bai; Fugen Zhou

doi:10.1364/BOE.7.005268

1. Introduction

Coronary artery disease (CAD) is a common disease and main cause of death worldwide [1]. C-arm based percutaneous coronary intervention (PCI) is a common procedure for the diagnosis and treatment of CAD, in which X-ray coronary angiography is one of the most commonly utilized method to assess CAD and is still considered as the gold standard in clinical procedures [1]. However, due to its projection characteristic, X-ray coronary angiography suffers from some clinical limitations, such as vessel overlap, loss of 3D spatial information. To better assist PCI procedures, many efforts have been devoted to perform 3D tomographic reconstruction of high-contrast coronary arteries using angiograms rotationally acquired on a C-arm. Compared with other 3D imaging modalities, such as cardiac CT [2], C-arm based 3D visualization has certain advantages. For example, it can be directly used for guiding PCI without the need of complicated registration process [3]. However, because of the vessel motion induced by respiratory and heart beating, accurate 3D reconstruction of coronary vessels still remains a challenge [1]. Poor reconstruction affects the accuracy of the quantitative coronary analysis (QCA) measurement, such as the length and diameter of the vessel branch or lesion. Inaccurate measurement may lead to incorrect selection of treatment apparatus and would reduce the effectiveness of the treatment consequently [1].

Cardiac gating and motion compensation are two primary ways to handle vessel motion for 3D reconstruction of coronary arteries. In cardiac-gated reconstruction for a selected cardiac phase p_r, a subset of angiograms is selected for reconstruction according to simultaneously recorded ECG signal or other surrogate signals. Nearest-neighbor (NN) and finite-window (FW) gating are two commonly used strategies for selecting angiograms [4]. For NN gating, only the angiogram closest to phase p_r is selected in each cardiac cycle. For FW gating, more angiograms around phase p_r are selected for reconstruction using some weighting functions, such as power of cosine function [5]. The gating window is a tradeoff between two different artifacts. As the commonly used reconstruction methods, such as the Feldkamp–Davis–Kress (FDK) method [6], are very sensitive to the angular coverage of the angiograms. Small gating window causes angular-undersampling problem and the reconstruction will be flooded with severe streak artifacts. Increasing gating window can ease the streak artifacts but will introduce more motion artifacts because of the morphological inconsistency of the coronary vessels in selected angiograms. In recent years, prior knowledge regularized iterative reconstruction has been studied to handle this dilemma [3, 7–12]. By imposing proper regularization during the iterative reconstruction process, the streak artifacts can be significantly suppressed even using less than 10 angiograms. Some of these methods have been investigated in clinical context and shown promising results [3].

Motion compensation is a commonly used technique in tomographic reconstruction of moving objects [1, 13]. For 3D coronary artery reconstruction, the motion of coronary artery is firstly estimated by some kinds of registration and then compensated in the reconstruction process using dynamic reconstruction methods [14, 15]. In this way, all angiograms could be used to perform the reconstruction at each phase without introducing severe motion artifacts. However, as the key to motion-compensated methods, accurate estimation of vessel motion is challenging and time-consuming [16]. Current methods either use simplified vessel motion model [17, 18] or pre-computed gated reconstruction (3D image or model) [19, 20] to ease the motion estimation problem. The accuracy of estimated motion is compromised or depends on the quality of gated reconstruction [20].

Besides the angular-undersampling problem, cardiac-gated reconstruction still suffers from the problem of residual motion. Because of the irregular heart motion or other reasons, there still exists residual vessel motion among selected angiograms which degrades the quality of cardiac-gated reconstruction. To handle this issue, strategies utilizing local motion compensation to improve the cardiac-gated reconstruction have been proposed. Typically, an initial cardiac-gated reconstruction is firstly performed and then the vessel motion between the reconstructed image and angiograms is estimated and compensated in the final reconstruction [7, 21–23]. Different from the above-mentioned motion-compensated reconstruction methods which try to estimate the vessel motion in the whole cardiac cycle, the local motion compensation methods only estimate the vessel motion between the reconstruction phase and the angiograms near that phase. As the magnitude of the local vessel motion is comparatively small, local motion estimation is easier.

Current local motion-compensated reconstruction methods typically modeled and estimated the vessel motion in the 2D projection domain by registering the angiograms to the reprojections of the initial reconstruction. The estimated transformation was applied to the angiograms and these transformed angiograms were used for final reconstruction. Hansis et al used 2D elastic transformation to model the vessel motion and proposed to estimate the parameters by matching the 2D vessel centerlines [7, 22]. This method requires accurate extraction of the centerline from the angiograms and the reprojections, which is a nontrivial task and commonly requires manual interaction [21]. Schwemmer et al modeled vessel motion using a multi-scale scheme of 2D affine and B-spline transformations, and utilized an intensity-based deformable registration method to estimate the parameters [21, 23]. As no complex pre-processing steps like centerline extraction was needed, this method could be fully automatic. However, it used FDK method for reconstruction and suffered the problem of angular-undersampling.

In this work, a new local motion-compensated cardiac-gated reconstruction method is proposed. We modeled the vessel motion in the 3D space and proposed a 3D/2D registration method to estimate the motion parameters. The initial reconstruction was performed using our previously developed regularized iterative reconstruction method with alternate reconstruction and segmentation strategy [12], which could handle the angular-undersampling problem well by imposing proper constraints. And the iterative reconstruction method was also extended with motion-compensation components and used to conduct the final reconstruction. The proposed method is discussed in detail in the next section. In Section 3, the performance of the proposed method is quantitatively evaluated. Discussion and conclusion are presented in Section 4.

2. Methods

As shown in Fig. 1, the proposed method consists of three main steps. Firstly, a subset of angiograms was selected using NN gating and used to perform the initial reconstruction. Secondly, the proposed 3D/2D registration method was utilized to estimate the vessel motion between the initial reconstruction and each angiogram. Finally, a regularized iterative reconstruction method with motion-compensation was proposed to reconstruct the final result.

Fig. 1 The flowchart of the proposed method.

Download Full Size | PDF

2.1 Initial reconstruction

NN gating results in angiograms which is more consistent with vessel morphology than FW gating. However, the number of chosen angiograms is so limited (typically 5-10) that it is not possible to reconstruct usable image using conventional methods, such as the FDK method. Previously, a regularized iterative reconstruction method was developed to perform NN cardiac-gated reconstruction [12]. Based on the observation of the sparsity of the image and the constraint of the non-negativity, this method try to solve the following problem using simultaneous threshold ART (START) update formula [8]:

\underset{X \in R^{N}}{\arg \min} {‖ X ‖}_{1} s . t . X \geq 0, A X = b

where

X \in R^{N}

is the image to be reconstructed;

b \in R^{M}

is the angiograms stored in a single vector;

A

is the system matrix which models the cone-beam imaging process of the angiograms;

{‖ X ‖}_{1}

is the

l_{1}

norm of the reconstructed image which is minimized to enforce the sparsity of the image [24]. A novel strategy of alternate reconstruction and segmentation was developed to suppress the non-vessel background of the angiograms which have great adverse effect on the reconstruction quality. During reconstruction, while the image intensity is iteratively updated, a contour is also evolved to segment the reconstructed vascular tree based on level set segmentation method. The segmented vascular tree is later used to reduce the angiogram background. As previously validated, this method could reconstruct adequate image using less than 10 angiograms and was used in this work to perform the initial reconstruction.

2.2 Motion estimation

Due to the residual vessel motion, the initial reconstruction exhibited motion artifacts which hindered its clinical usage. A 3D/2D registration method was proposed to estimate the vessel motion between the initial reconstruction and each angiogram. The estimated residual vessel motion will be compensated in the final reconstruction.

2.2.1 The registration method

The motion of coronary artery is complex and non-rigid. During cardiac contraction, the arteries move consistently toward the left, inferior and anterior. The right coronary tree has larger displacement and velocity than the left coronary tree, and the distal displacement is larger than the proximal displacement [25]. To accurately model the coronary motion, three-order B-spline representation was used in this work to model the motion. Given predefined 3D grid of control points of size $C_{x} \times C_{y} \times C_{z}$ , the transformed position of a point $X \in R^{3}$ is determined by the motion vector of its 4 × 4 × 4 neighboring control points [26]

T (X, D) = X + \sum_{l = 0}^{3} \sum_{m = 0}^{3} \sum_{n = 0}^{3} B_{l} (u) B_{m} (v) B_{n} (w) d_{i + l, j + m, k + n}

where

B_{0}

through

B_{3}

are the approximating third-order spline polynomials;

d_{i + l, j + m, k + n}

are the motion vectors of the neighboring control points;

D

is the parameter set to be determined and contains the motion vectors of all control points (

d_{i, j, k}

):

D = {d_{i, j, k} | i = 1 \dots C_{x}, j = 1 \dots C_{y}, k = 1 \dots C_{z}}

_.

Using the B-spline motion model, the aim of the 3D/2D registration is to determine the optimal motion parameters $D^{*}$ that make the reprojected vessel of the transformed 3D initial reconstruction $I (T (X, D))$ match each angiogram $A_{i}$ :

D^{*} = \arg \max_{D} M (A_{i}, P_{i} I (T (X, D)))

where

P_{i} I (T (X, D))

represents the reprojection of the transformed 3D initial reconstruction according to the imaging geometry of the i-th angiogram.

M

is the similarity metric that measures the consistence between the reprojection and the angiogram. As the gray value of angiogram and reprojection differs due to the error within the initial reconstruction and the reprojection process, we chose normalized correlation (NC) as the similarity metric which is more robust to intensity discrepancy than purely intensity-based metric like the sum of squared differences [21].

The deformable registration between the 3D image and 2D angiogram is an ill-posed problem. There are many transformations that could make the similarity increase. However, some of them may be wrong and implausible in reality. Therefore, according to the characteristics of the vessel motion, several constraints were imposed to obtain the accurate vessel motion. One of the constraints is the smoothness constraint $E_{s m o o t h}$ which ensures the vessel motion is smooth and continuous by restricting the second-order partial derivatives of the deformation [27]. The other constraint used is the volume preserving constraint $E_{v o l u m e}$ which ensures the volume of the vessel will not change much during the motion [27].

Combined the similarity measure and motion constraints, we arrived at the following registration optimization problem

\underset{D}{\arg \min} E (D) = - E_{s i m i l a r i t y} (D) + β_{1} E_{s m o o t h} (D) + β_{2} E_{v o l u m e} (D)

where

β_{1}

and

β_{2}

are the weighting factors for the smoothness and volume preserving constraints. As the gradient of the optimization function is available, we proposed to solve the optimization problem using Limited-memory LBFGS method [28] with Moré-Thuente line search method [29], which has been widely used for registration problems [30]. And a multi-resolution strategy is used to improve the efficiency and robustness of the optimization. The derivatives of the

E_{s m o o t h}

and

E_{v o l u m e}

with respect to the transformation parameters were calculated according to the reference [31]. Writing

E_{s i m i l a r i t y}

as

N C (A_{i}, A_{i}^{p}) = N C (A_{i}, P_{i} I (X^{'})) = N C (A_{i}, P_{i} I (T (X, D)))

its derivative with respect to a parameter d is derived as follows using the chain rule:

\frac{\partial N C}{\partial d} = \frac{\partial T}{\partial d} \frac{\partial I}{\partial X'} P_{i}^{T} \frac{\partial N C}{\partial A_{i}^{p}}

It was computed according to the following procedures in this work:

(1) transform the initial reconstruction $I$ using the current transformation and reproject the transformed image to obtain $A_{i}^{p}$ ;
(2) calculate the $\partial N C / \partial A_{i}^{p}$ which is also a 2D image as the angiogram $A_{i}$ ;
(3) backproject $\partial N C / \partial A_{i}^{p}$ to obtain the corresponding 3D image $G$ of the derivative value;
(4) calculate $\partial N C / \partial d$ as the sum of contribution of all sample points $\frac{\partial N C}{\partial d} = \sum_{X_{i} \in Ω} {(\frac{\partial N C}{\partial d})}_{X_{i}} = \sum_{X_{i} \in Ω} G_{X_{i}} (\frac{\partial T}{\partial d} �^{T} \frac{\partial I}{\partial X_{i}'}$ , where $G_{X_{i}}$ is the scalar value of the 3D gradient image $G$ at position $X_{i}$ ; $Ω$ is the set of all sample points; $\partial I / \partial X_{i}^{'}$ is the image gradient of $I$ at the transformed position $X_{i}^{'} = T (X_{i}, D)$ ; $\partial T / \partial d$ is the derivative of B-spline transform with respected to d at position $X_{i}$ .

2.2.2 Registration region of interest (ROI)

We desired to align the reconstructed 3D vasculature with the vessel in angiograms. As the vessel only occupy a small fraction of the reconstructed image, the registration could be sped up and stabilized by restricting the registration ROI to a small neighborhood around the reconstructed 3D vasculature. During registration, only voxels within the ROI were used for metric evaluation.

Using the alternate reconstruction and segmentation strategy, a binary segmentation of the reconstructed vasculature could be obtained [12]. Therefore, in this work, we defined the registration ROI to be the region

R O I = {x | d i s t (x, v e s s e l) \leq d}

where

d i s t (x, v e s s e l)

is the distance of a point to the surface of the segmented vessel and the distance of the points within the vessel is set to be zero. The threshold d determines the size of the ROI and should be chosen largely enough to accommodate the possible displacement of the vessel branch. It can be chosen by inspecting whether the reprojection of the ROI encloses the vessels in all angiograms and was fixed to be 6mm in the following experiments.

2.2 Motion-compensated reconstruction

In the final reconstruction, the regularized iterative reconstruction method used for initial reconstruction was extended to compensate the estimated vessel motion. In its initial form, the reconstructed image $X$ is iteratively updated according to each angiogram $b_{i}$ as follows:

χ^{k + 1} = χ^{k} + α V^{- 1} A^{T} W^{- 1} (b_{i} - A X^{k})

X^{k + 1} = H (χ^{k + 1})

where

χ

is an auxiliary image;

H

is the Heaviside function and

X

is obtained by setting the negative elements of

χ

to zero;

V

and

W

are diagonal matrix with the j-th diagonal elements as

\sum_{i = 1}^{m} A_{i j}

and

\max (C_{\min}, \sum_{i = 1}^{n} A_{j i} B i n (X_{i}^{k}))

respectively;

C_{\min}

is a constant and

B i n (X_{i}^{k})

is the binary segmentation of

X_{i}^{k}

using the threshold of zero;

A_{i j}

is the element of the system matrix

A

. As shown in Fig. 2, one iteration mainly contains four steps: reprojection (

A X^{k}

), calculation of data discrepancy (

b_{i} - A X^{k}

), backprojection (

V^{- 1} A^{T} W^{- 1} (b_{i} - A X^{k})

) and finally the update of the image.

Fig. 2 The main steps of the reconstruction iteration with motion compensation.

Download Full Size | PDF

When there are residual vessel motion between the reconstructed image $X$ and angiogram $b_{i}$ , the direct reprojection of $X$ mismatch the vessel in $b_{i}$ . $X$ needs to be transformed according to the estimated motion field $D_{i}$ before reprojection to calculate the accurate data discrepancy, which makes the formula for discrepancy calculation change to $b_{i} - A T (X^{k}, D_{i})$ . Meanwhile, the residual vessel motion will also cause the backprojection of the data discrepancy to mismatch the reconstructed image. The backprojection of the data discrepancy needs to be transformed according to the inverse of the estimated motion field $D_{i}$ to accurately update the reconstructed image. In this way, the formula for auxiliary image update (Eq. (8)) changes to

χ^{k + 1} = χ^{k} + α T (V^{- 1} A^{T} {\hat{W}}^{- 1} (b_{i} - A T (X^{k}, D_{i})), D_{_{i}}^{- 1})

where the j-th diagonal elements of

\hat{W}

changes to

\max (C_{\min}, \sum_{i = 1}^{n} A_{j i} B i n (T (X^{k}, D_{i})))

;

D_{_{i}}^{- 1}

is the inverse of the estimated motion field. Finally, the update formula of the iterative reconstruction with motion compensation is the combination of the Eq. (9) and Eq. (10). Figure 2 shows the main steps of the iteration with motion compensation. The inverse of the motion field is needed and calculated using a fixed-point approach [32].

As the segmentation of the reconstructed vessel could be obtained from initial reconstruction, the background of the angiograms could be suppressed before reconstruction and the segmentation part of the method was omitted.

3. Experiment and results

3.1 Experiment setup

The proposed method was evaluated on angiograms simulated using the 4D XCAT phantom [33]. This phantom contains anatomy information derived from real CT/MRI data sets, which allows the simulation of realistic data sets for cardiac C-arm CT. Projections were generated using cone-beam imaging geometry with source to detector distance (SDD) of 150 cm and source to axis distance (SAD) of 50 cm. The dimension of the simulated angiograms is 512 × 512 with an isotropic resolution of 0.5 mm. According to real C-arm protocols, the angiograms were simulated using the following setting by default. The C-arm rotated 220 degrees around the patient in 7 seconds. The heart-beating rate was set to be 80 beats/minute. The acquisition started from cardiac phase 20% at imaging frequency of 30 fps. In total, 210 angiograms were acquired which covered 9 cardiac cycles. Figure 3 (a) shows the 3D vessel morphology at different cardiac phases. Figure 3(b) shows the vessel centerlines used for quantitative evaluation. Figure 4 shows the cardiac phases of the simulated angiograms and the chosen angiograms using NN gating and FW gating (window width = 10%) for the cardiac phase 90%.

Fig. 3 (a) 3D vessel morphology at different cardiac phases: 0% (white), 40% (red), 70% (green), and 90% (blue); (b) the vessel centerlines used for quantitative evaluation.

Download Full Size | PDF

Fig. 4 The cardiac phase of the simulated angiograms and the chosen angiograms using NN gating (red circled) and FW gating (green circled, window width = 10%) for cardiac phase 90%.

Download Full Size | PDF

We tried to reconstruct the vessel at cardiac phase 90%. The residual motion was caused by the discrepancy between the phases of the selected angiograms and target phase, and its magnitude was doubled to simulate imperfect gating. The quality of the reconstruction was evaluated using the maximum mean overlap (MMO) and the relative radius error (RRE) metrics [12]. MMO evaluates the overall morphological accuracy of reconstructed vasculature and a larger MMO value means better reconstruction. RRE dedicates to the accuracy of vessel radius which is of great clinical values in PCI and a smaller RRE value means more accurate vessel size. We utilized the Vascular Modeling Toolkit (VMTK) [34] to extract the vessel centerline used for RRE calculation. The centerline used for evaluation is shown in Fig. 3(b).

The experiments were run on a personal laptop with an Intel i5-3210 CPU, 12 GB Memory, and a NVIDIA NVS 5400M Video Card. The reconstruction framework was implemented in MATLAB while some procedures, such as the deformable 3D/2D registration, were implemented using C/C + + and called from MATLAB. We also used graphic processing unit (GPU) to implement computation-extensive reprojection and backprojection procedures to improve reconstruction efficiency.

3.2 Parameter selection

The parameters $β_{1}$ and $β_{2}$ in Eq. (4) control the weights of the smoothness and volume preserving constraints. Through experiments, we found it worked well to set the two parameters to 1 and 0.1 respectively.

Three levels of registration resolution were adopted. The spacing of control points were set to 16 mm, 8 mm and 4 mm respectively. The registration iteration stopped at each resolution when the Moré-Thuente method [29] failed to find a suitable step or when the decrease of the objective function between two iterations was below 1e-2, 1e-3, 5e-4 for the first, second and third registration resolution or when the number of iterations reached 60.

The max number of iteration of the fixed-point approach for inverting deformation field was set to 20. The iteration also stopped when the mean and max error were below 1e-6 and 0.1.

For motion-compensated reconstruction, the max number of the iterations was set to 15.

3.3 Performance evaluation

3.3.1 Performance evaluation under different magnitude of residual motion

To evaluate the performance of proposed method under different magnitude of residual motion, three sequences of angiograms were simulated with increasing magnitude of residual vessel motion. The motion magnitude for the second and third sequences was two and three times the motion magnitude of the first sequence respectively. The motion magnitude will be denoted as 1 × , 2 × and 3 × in the following. The motion-compensated reconstruction was performed using angiograms selected using the NN gating (8 angiograms) and the FW gating (window width = 5%, 16 angiograms), and the results were compared. As more vessel motion exists in the FW gated angiograms, the quality of initial reconstruction using FW gating would be much worse than that using NN gating. Poor initial reconstruction would make the accurate registration very difficult. Because of this, for the case of motion-compensated reconstruction using FW gating (MC_fw), the initial reconstruction was performed using NN gating. The motion estimation and final reconstruction were performed for the angiograms chosen using FW gating. The same strategy was also used in other work [7].

Figure 5 shows the volume rendering of the initial reconstruction with different magnitude of residual motion. With larger motion magnitude, the motion artifacts were increased and the edges were getting more blurred, which made the branch thinner visually.

Fig. 5 The volume rendering of the initial and motion-compensated reconstructions of the proposed method (using NN gating) under different magnitudes of residual motion: (a) no residual motion; (b) 1 × ; (c) 2 × ; (d) 3 × . The top row of (b-d) show the initial reconstructions; the bottom row of (b-d) show the motion-compensated reconstructions.

Download Full Size | PDF

Figure 6 shows the registration iteration information for some angiograms. It can be observed that the registration can converge in finite number of iterations. For angiograms with large vessel motion, more iterations were needed to match the initial reconstruction to them. The NC between original angiograms and the reprojection of the initial reconstruction was found to increase from 0.87 ± 0.08, 0.74 ± 0.17, 0.66 ± 0.21 before registration to 0.96 ± 0.01, 0.95 ± 0.01, 0.94 ± 0.02 after registration for the three sequences. This demonstrated the efficacy of the proposed 3D/2D registration method. Figure 7 compares the angiograms and the reprojections of the initial reconstruction before and after the registration. As we can see, the proposed registration method could match the vessel morphology between initial reconstruction and the angiograms well.

Fig. 6 The registration iteration information for angiograms with different magnitude of residual motion. (a) 1 × ; (b) 3 × . The iterations in different registration resolution is represented using different markers.

Download Full Size | PDF

Fig. 7 Chessboard display of the angiogram and reprojection of the initial reconstruction for the data with 3 × motion before the registration (the first row) and after the registration (the second row).

Download Full Size | PDF

Figure 8 shows the quantitative metrics for different reconstructions. Overall, the reconstruction degraded with larger motion magnitude and the proposed method could improve the reconstruction quality effectively. Compared with the motion-compensated reconstruction using NN gating (MC_nn), the result of MC_fw was comparable for the data with 1 × motion but worse when motion magnitude is larger. One reason for this is that the regularized iterative reconstruction method could obtain high-quality reconstruction using NN gating if there is no residual motion. There is not much room to improve by using more angiograms. What is more important is that the proposed registration method could not perfectly align the initial reconstruction and angiograms. When motion magnitude is large, the registration error will be significant enough to degrade the reconstruction result. In the following, we only consider using NN gating for the proposed motion-compensated reconstruction method.

Fig. 8 Quantitative metrics for different reconstructions. Initial: the initial reconstruction; MC_nn and MC_fw: the motion-compensated reconstruction of the proposed method using NN gating and using a finite gating width of 10% respectively.

Download Full Size | PDF

Figure 5 compares the motion-compensated reconstruction with the initial reconstruction. As we can see, the motion artifacts were efficiently reduced and the edges became sharper than the initial reconstruction.

3.3.2 Analysis of the registration constraints

To analyze the influence of the smoothness and volume preserving constraints, we performed the motion-compensated reconstruction using different settings of constraint weights ( $β_{1}$ and $β_{2}$ ). The weight settings and the metrics for the reconstruction results are shown in Table 1. The NC between reprojections and origin angiograms were 0.94 ± 0.02 before registration. After registration, the NC was all 0.96 ± 0.01 for the three cases. When the motion magnitude is small (1 × ), the results of different reconstructions were comparable. However, with larger motion magnitude, the reconstruction without imposing constraints ( $β_{1}$ = 0 $β_{2}$ = 0) was poorer than the other two cases. The reason for this could be perceived by inspecting the registration result and the obtained motion field. Figure 9 shows these results for one angiogram. Because of the residual vessel motion, the vessel edge of the initial reconstruction was blurred which made the vessels in the reprojection thinner than the original angiograms. Without imposing the constraints, the vessels would be expanded to match the initial angiograms during registration and the obtained motion field was very irregular. As an example, please compare the vessel cross-sections indicated by the red arrow in Fig. 9(a)-9(d). Though the NC was increased by registration, the obtained vessel motion was inconsistent with the reality and the quality of the motion-compensated reconstruction was degraded as a consequence. In comparison, the obtained motion field using constraints was very smooth and regular. And the registration method was not sensitive to the constraint weights.

Table 1. Comparison of the metric values for reconstructions using different constraint weights. For comparison, the MMO and RRE for the initial reconstruction are also shown

View Table | View all tables in this article

Fig. 9 (a): one original angiograms; (b): the corresponding reprojection of the initial reconstruction; (c) and (d): the reprojection of the registration results with ( $β_{1}$ = 1 $β_{2}$ = 0.1) and without imposing constraints; (e) and (f): the obtained motion fields near the vessel surface with and without imposing constraints

Download Full Size | PDF

3.3.3 Performance evaluation for different angles

The quality of iterative reconstruction depends on the angles of the angiograms [35]. Besides, angiograms at different angles pose different situations to the registration method. Therefore, we evaluated the performance of the proposed method for different sets of angiograms at different angles.

According to the experiment setting, changing the starting cardiac phase will change the angles of the selected angiograms used for reconstruction. Three sequences of angiograms with 2 × motion were simulated by setting the starting cardiac phase to 0%, 20% and 40% and used for reconstruction. The metrics for the reconstruction results are shown in Table 2. As we can see, the metrics varied for different reconstruction but did not differ much. The proposed method could improve the reconstruction quality for all data and the motion-compensated reconstruction results had comparable MMO and RRE values.

Table 2. Comparison of the metric values for different reconstruction. Ideal: reconstruction using data without residual motion; Initial: initial reconstruction; MC: motion-compensated reconstruction

View Table | View all tables in this article

3.3.4 Comparison with a state-of-the-art method

For comparison, we also performed the reconstruction using the most recent state-of-the-art method proposed by Schwemmer [21] on three sequences of angiograms with 1 × , 2 × and 3 × motion using the default imaging setting. This method used 2D registration for motion estimation and FDK method for the reconstruction. Promising results have been obtained by this method [21, 23]. In the following, we will refer to this method as MC_2D&FDK. The algorithm parameters were tuned to get the best reconstruction result. Two iterations of motion-compensation were conducted and the gating width was selected as 40%.

Figure 10 compares the quantitative metrics for the reconstruction of MC_2D&FDK and the proposed method. Figure 11 shows the volume rendering the reconstruction results of the MC_2D&FDK. As we can see, the MC_2D&FDK effectively improved the reconstruction quality. A second iteration of motion compensation could further improve the quality especially when the residual motion magnitude was large.

Fig. 10 Quantitative metrics for different reconstructions. Initial_proposed and MC_proposed: the initial and the motion-compensated reconstruction of the proposed method respectively; Initial_2D&FDK: the initial reconstruction of the MC_2D&FDK method; MC_2D&FDK-FirstIter and MC_2D&FDK-SecondIter: the results of MC_2D&FDK after first and second iteration of motion compensation respectively

Download Full Size | PDF

Fig. 11 The volume rendering of the gated and motion-compensated reconstruction of the MC_2D&FDK method under different magnitude of residual motion: (a) no residual motion; (b) 1 × ; (c) 2 × ; (d) 3 × . The top row of (b-d) shows the initial reconstructions; the middle row of (b-d) shows the results after first iteration of motion-compensation; the bottom row of (b-d) shows the results after second iteration

Download Full Size | PDF

Though only one iteration of the proposed motion compensation method was performed, its reconstruction quality was better than MC_2D&FDK. As for the RRE, the results for two methods were comparable and the proposed method was slightly better. However, as for the MMO, the proposed method significantly outperformed MC_2D&FDK. This is because that MC_2D&FDK used FDK method for the reconstruction which is very sensitive to under-sampling problem. Though wider gating window (40%) and streak reduction strategy [5] were used, the reconstruction still suffered from severe streak artifacts, as shown in Fig. 11.

As stated in the next section, performing one iteration of motion compensation took almost the same time for both methods. A second iteration of motion compensation would approximately double the time. Therefore, the proposed method could obtain better reconstruction using much less time.

3.3.5 Reconstruction time

lists the reconstruction time of the proposed method and MC_2D&FDK. For the proposed method, registration of one angiogram took about 60 seconds and inverting the deformation field took about 30 seconds. It took about 12 minutes to process 8 angiograms using NN gating. The iterative reconstruction with motion-compensation took 6 minutes. In total, the whole reconstruction took about 18 minutes. In comparison, 2D registration for one angiogram took about 17 seconds. Using 40% gating width resulted in 73 angiograms and the total registration time was about 20 minutes. Plus the FDK reconstruction time, the total reconstruction time for the MC_2D&FDK method was about 21 minutes. The reconstruction time of the two methods was comparable for one iteration of motion-compensated reconstruction (Table 3).

Table 3. Reconstruction time (minutes)

View Table | View all tables in this article

4. Discussion and conclusion

The angular-undersampling and the residual motion are two main factors that degenerate the cardiac-gated reconstruction of coronary artery [7]. In this work, we propose a fully automatic iterative reconstruction method with local motion compensation to address these problems. This method employs a regularized iterative reconstruction method with alternate reconstruction and segmentation strategy to handle angular-undersampling and 3D/2D registration based motion estimation and compensation methods to handle residual motion. Experimental results show that the method could effectively estimate and compensate the residual vessel motion, and achieve high-quality reconstruction.

When the motion magnitude is large, there still is a quality gap between the motion-compensated reconstruction and the ideal reconstruction without residual motion. Therefore, it is commonly suggested to perform the cardiac-gated reconstruction at cardiac phase with slow vessel motion, such as the end-systolic or late-diastolic phases [7]. Performing multiple iterations of motion compensation can further improve the reconstruction quality [7, 21, 23]. Take the experiment in section 3.3.4 for example, after second iteration of the proposed method, the MMO increased from 0.88 to 0.90 and the RRE decreased from 7.0 to 5.6 for the data with 3 × motion. However, the reconstruction time would be too long for clinical application. Therefore, in this work, we did not evaluate multiple iterations of motion compensation in detail.

The choice of reconstruction method has a great effect on the final reconstruction quality. Firstly, the quality of initial reconstruction affects the accuracy of estimated vessel motion. Secondly, the quality of final reconstruction directly depended on the reconstruction method. Compared with the FDK method, the regularized iterative reconstruction method uses NN gated angiograms for reconstruction which have more consistent vessel morphology and thus introduce less motion artifacts to the result. And the streak artifacts induced by angular-undersampling could be successfully suppressed by imposing proper constraints. Therefore, using regularized iterative reconstruction method for motion-compensated reconstruction could possibly lead to better result.

Motion-compensated reconstruction consists of image registration and reconstruction steps, which is very time-consuming. Currently, the proposed method took about 18 minutes to finish the whole motion-compensated reconstruction procedure. Though the reconstruction time is comparable to other methods, it needs to be reduced to around or less than 3 minutes to be clinical usable [19]. In the next work, we try to reduce the reconstruction time in the following ways. Firstly, we will implement more procedures using the GPU. For example, the image transformation and the fixed-point approach for inverting deformation field are very suitable for parallel computing and can be implemented very efficiently using GPU. Secondly, current implementation utilized hardware disk for data sharing between different procedures. Reducing disk access could greatly improve reconstruction efficiency.

The 3D/2D registration method used NC as similarity metric to accommodate the gray value difference between the reprojection and the angiogram. Compared with NC, mutual information (MI) metric can handle nonlinear relation and is more robust to gray value difference, but is more computing extensive. Using MI as registration similarity metric may improve the registration accuracy and thus the reconstruction quality, which will be investigated in the future work.

Currently, the proposed method was evaluated on simulated data. The availability of the ground truth makes objective and quantitative evaluation possible. As the experiments were carefully designed to simulate the clinical conditions, the performance of the method can be verified. Further quantitative evaluation of the method using clinical data will be carried out in the next work.

Funding

Natural National Science Foundation of China (NSFC) (61601012, 61171005, 61271023), Program for New Century Excellent Talents in University (NCET-13-0020) and Fundamental Research Funds for the Central Universities (YWF-16- BJ-Y-28).

Acknowledgments

A short abstract (less than 300 words) of this work has been shown at the Annual Meeting of the American Association of Physicists in Medicine in 2015, Anaheim, CA, United States.

References and links

1. S. Çimen, A. Gooya, M. Grass, and A. F. Frangi, “Reconstruction of Coronary Arteries from X-ray Angiography: A Review,” Med. Image Anal. 32, 46–68 (2016). [CrossRef] [PubMed]

2. G. Cao, B. Liu, H. Gong, H. Yu, and G. Wang, “A Stationary-Sources and Rotating-Detectors Computed Tomography Architecture for Higher Temporal Resolution and Lower Radiation Dose,” IEEE Access 2, 1263–1271 (2014). [CrossRef]

3. A. M. Neubauer, J. A. Garcia, J. C. Messenger, E. Hansis, M. S. Kim, A. J. Klein, G. A. Schoonenberg, M. Grass, and J. D. Carroll, “Clinical feasibility of a fully automated 3D reconstruction of rotational coronary X-ray angiograms,” Circ. Cardiovasc. Interv. 3(1), 71–79 (2010). [CrossRef] [PubMed]

4. G. Schoonenberg, A. Neubauer, and M. Grass, “Three-Dimensional Coronary Visualization, Part 2: 3D Reconstruction,” Cardiol. Clin. 27(3), 453–465 (2009). [CrossRef] [PubMed]

5. C. Rohkohl, G. Lauritsch, A. Nottling, M. Prummer, and J. Hornegger, “C-arm CT: Reconstruction of dynamic high contrast objects applied to the coronary sinus,” in Nuclear Science Symposium Conference Record,2008. NSS '08. IEEE, 2008), 5113–5120. [CrossRef]

6. L. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical cone-beam algorithm,” J. Opt. Soc. Am. A 1(6), 612–619 (1984). [CrossRef]

7. E. Hansis, J. D. Carroll, D. Schäfer, O. Dössel, and M. Grass, “High-quality 3-D coronary artery imaging on an interventional C-arm x-ray system,” Med. Phys. 37(4), 1601–1609 (2010). [CrossRef] [PubMed]

8. E. Hansis, D. Schäfer, O. Dössel, and M. Grass, “Evaluation of iterative sparse object reconstruction from few projections for 3-D rotational coronary angiography,” IEEE Trans. Med. Imaging 27(11), 1548–1555 (2008). [CrossRef] [PubMed]

9. M. Li, H. Yang, and H. Kudo, “An accurate iterative reconstruction algorithm for sparse objects: application to 3D blood vessel reconstruction from a limited number of projections,” Phys. Med. Biol. 47(15), 2599–2609 (2002). [CrossRef] [PubMed]

10. H. Wu, C. Rohkohl, and J. Hornegger, “Total Variation Regularization Method for 3D Rotational Coronary Angiography,” in Bildverarbeitung für die Medizin 2011, H. Handels, J. Ehrhardt, T. M. Deserno, H.-P. Meinzer, and T. Tolxdorff, eds. (Springer Berlin Heidelberg, 2011), pp. 434–438.

11. H. Yining, X. Lizhe, J. C. Nunes, J. J. Bellanger, M. Bedossa, and C. Toumoulin, “ECG gated tomographic reconstruction for 3-D rotational coronary angiography,” in Engineering in Medicine and Biology Society (EMBC),2010Annual International Conference of the IEEE, 2010), 3614–3617. [CrossRef]

12. B. Liu, F. Zhou, and X. Bai, “Improved C-arm cardiac cone beam CT based on alternate reconstruction and segmentation,” Biomed. Signal Process. Control 13, 113–122 (2014). [CrossRef]

13. M. F. Kraus, B. Potsaid, M. A. Mayer, R. Bock, B. Baumann, J. J. Liu, J. Hornegger, and J. G. Fujimoto, “Motion correction in optical coherence tomography volumes on a per A-scan basis using orthogonal scan patterns,” Biomed. Opt. Express 3(6), 1182–1199 (2012). [CrossRef] [PubMed]

14. D. Schäfer, J. Borgert, V. Rasche, and M. Grass, “Motion-compensated and gated cone beam filtered back-projection for 3-D rotational X-ray angiography,” IEEE Trans. Med. Imaging 25(7), 898–906 (2006). [CrossRef] [PubMed]

15. T. Pengpan, W. Qiu, N. D. Smith, and M. Soleimani, “Cone Beam CT using motion-compensated algebraic reconstruction methods with limited data,” Comput. Methods Programs Biomed. 105(3), 246–256 (2012). [CrossRef] [PubMed]

16. E. Hansis, H. Schomberg, K. Erhard, O. Dössel, and M. Grass, “Four-dimensional cardiac reconstruction from rotational x-ray sequences: first results for 4D coronary angiography,” Proc. SPIE 7258, 72580B (2009).

17. A. Keil, J. Vogel, G. Lauritsch, and N. Navab, “Dynamic cone beam reconstruction using a new level set formulation,” Med Image Comput Comput Assist Interv 12(2), 389–397 (2009). [PubMed]

18. C. Rohkohl, G. Lauritsch, L. Biller, and J. Hornegger, “ECG-Gated Interventional Cardiac Reconstruction for Non-periodic Motion,” in Medical Image Computing and Computer-Assisted Intervention – MICCAI 2010 (Springer Berlin / Heidelberg, 2010), pp. 151–158.

19. C. Rohkohl, G. Lauritsch, L. Biller, M. Prümmer, J. Boese, and J. Hornegger, “Interventional 4D motion estimation and reconstruction of cardiac vasculature without motion periodicity assumption,” Med. Image Anal. 14(5), 687–694 (2010). [CrossRef] [PubMed]

20. K. Müller, A. K. Maier, C. Schwemmer, G. Lauritsch, S. De Buck, J. Y. Wielandts, J. Hornegger, and R. Fahrig, “Image artefact propagation in motion estimation and reconstruction in interventional cardiac C-arm CT,” Phys. Med. Biol. 59(12), 3121–3138 (2014). [CrossRef] [PubMed]

21. C. Schwemmer, C. Rohkohl, G. Lauritsch, K. Müller, and J. Hornegger, “Residual motion compensation in ECG-gated interventional cardiac vasculature reconstruction,” Phys. Med. Biol. 58(11), 3717–3737 (2013). [CrossRef] [PubMed]

22. E. Hansis, D. Schäfer, O. Dössel, and M. Grass, “Projection-based motion compensation for gated coronary artery reconstruction from rotational x-ray angiograms,” Phys. Med. Biol. 53(14), 3807–3820 (2008). [CrossRef] [PubMed]

23. C. Schwemmer, G. Lauritsch, A. Kleinfeld, C. Rohkohl, K. Müller, and J. Hornegger, “Clinical Data Evaluation of C-arm-based Motion Compensated Coronary Artery Reconstruction,” in The third international conference on image formation in X-ray computed tomography(2014), pp. 60–63.

24. H. Guo, J. Yu, X. He, Y. Hou, F. Dong, and S. Zhang, “Improved sparse reconstruction for fluorescence molecular tomography with L1/2 regularization,” Biomed. Opt. Express 6(5), 1648–1664 (2015). [CrossRef] [PubMed]

25. G. Shechter, J. R. Resar, and E. R. McVeigh, “Displacement and velocity of the coronary arteries: cardiac and respiratory motion,” IEEE Trans. Med. Imaging 25(3), 369–375 (2006). [CrossRef] [PubMed]

26. D. Rueckert, L. I. Sonoda, C. Hayes, D. L. Hill, M. O. Leach, and D. J. Hawkes, “Nonrigid registration using free-form deformations: application to breast MR images,” IEEE Trans. Med. Imaging 18(8), 712–721 (1999). [CrossRef] [PubMed]

27. T. Rohlfing, C. R. Maurer Jr, D. A. Bluemke, and M. A. Jacobs, “Volume-preserving nonrigid registration of MR breast images using free-form deformation with an incompressibility constraint,” IEEE Trans. Med. Imaging 22(6), 730–741 (2003). [CrossRef] [PubMed]

28. D. C. Liu and J. Nocedal, “On the limited memory BFGS method for large scale optimization,” Math. Program. 45(1-3), 503–528 (1989). [CrossRef]

29. J. J. Moré and D. J. Thuente, “Line search algorithms with guaranteed sufficient decrease,” ACM Trans. Math. Softw. 20(3), 286–307 (1994). [CrossRef]

30. F. Yang, M. Ding, X. Zhang, W. Hou, and C. Zhong, “Non-rigid multi-modal medical image registration by combining L-BFGS-B with cat swarm optimization,” Inf. Sci. 316, 440–456 (2015). [CrossRef]

31. M. Staring and S. Klein, “Itk:Transforms supporting spatial derivatives,” Insight J. (2010).

32. M. Chen, W. Lu, Q. Chen, K. J. Ruchala, and G. H. Olivera, “A simple fixed-point approach to invert a deformation field,” Med. Phys. 35(1), 81–88 (2008). [CrossRef] [PubMed]

33. W. P. Segars, G. Sturgeon, S. Mendonca, J. Grimes, and B. M. Tsui, “4D XCAT phantom for multimodality imaging research,” Med. Phys. 37(9), 4902–4915 (2010). [CrossRef] [PubMed]

34. L. Antiga, M. Piccinelli, L. Botti, B. Ene-Iordache, A. Remuzzi, and D. A. Steinman, “An image-based modeling framework for patient-specific computational hemodynamics,” Med. Biol. Eng. Comput. 46(11), 1097–1112 (2008). [CrossRef] [PubMed]

35. M. Li, H. Kudo, J. Hu, and R. H. Johnson, “Improved iterative algorithm for sparse object reconstruction and its performance evaluation with micro-CT data,” IEEE Trans. Nucl. Sci. 51(3), 659–666 (2004). [CrossRef]

Local motion-compensated method for high-quality 3D coronary artery reconstruction

Abstract

1. Introduction

2. Methods

2.1 Initial reconstruction

2.2 Motion estimation

2.2.1 The registration method

2.2.2 Registration region of interest (ROI)

2.2 Motion-compensated reconstruction

3. Experiment and results

3.1 Experiment setup

3.2 Parameter selection

3.3 Performance evaluation

3.3.1 Performance evaluation under different magnitude of residual motion

3.3.2 Analysis of the registration constraints

3.3.3 Performance evaluation for different angles

3.3.4 Comparison with a state-of-the-art method

3.3.5 Reconstruction time

4. Discussion and conclusion

Funding

Acknowledgments

References and links

Cited By

Figures (11)

Tables (3)

Equations (10)

Biomedical Optics Express

Motion Magnitude	Metric	Initial reconstruction	$β_{1}$ = 0 $β_{2}$ = 0	$β_{1}$ = 1 $β_{2}$ = 0.1	$β_{1}$ = 2 $β_{2}$ = 0.2
1 ×	MMO	0.95	0.96	0.96	0.96
1 ×	RRE	6.8	3.7	3.8	3.9
2 ×	MMO	0.87	0.90	0.93	0.93
2 ×	RRE	15.9	6.6	4.4	4.9
3 ×	MMO	0.79	0.85	0.89	0.89
3 ×	RRE	23.4	8.2	7.0	7.0

Metric	Reconstruction	0%	20%	40%
MMO	Ideal	0.98	0.98	0.98
	Initial	0.88	0.87	0.88
	MC	0.94	0.93	0.94
RRE	Ideal	3.2	3.0	3.4
	Initial	15.6	15.9	15.7
	MC	4.1	4.4	4.2