Point-Cloud Method for Automated 3D Coronary Tree Reconstruction From Multiple Non-Simultaneous Angiographic Projections

X-ray angiography is the most commonly used imaging modality for the detection of coronary stenoses due to its high spatial and temporal resolution of lumen contour and its utility to guide coronary interventions in real time. However, the high inter- and intra-observer variability in interpreting the geometry of 3D vascular structure based on multiple 2D image projections is a limitation in the accurate determination of lesion severity. This could be addressed by the 3D reconstruction of the coronary arterial (CA) tree. The automated reconstruction of 3D CA tree from 2D projections is challenging due to the existence of several imaging artifacts, such as vessel overlap, foreshortening, and most importantly respiratory and cardiac motion. Along with these artifacts, the acquisition geometry introduces the possibility of generating false vessel segments in the reconstruction. Our approach aims to reduce the motion artifacts in angiographic projections by developing a new method for rigid and non-rigid motion correction. A novel point-cloud based approach is subsequently introduced for reconstruction of 3D vessel centerlines by iteratively minimizing the reconstruction error. The performance of the proposed 3D reconstruction is evaluated using angiographic projections from 45 patients, producing average reprojection errors of ${0.092} \pm {0.055}$ mm and ${0.910} \pm {0.352}$ mm for 3D centerlines reconstruction, when co-registered with the parent vessels on projection planes that were/were not used to derive the 3D reconstruction, respectively. A comparison of the reconstructed 3D lumen surface with optical coherence tomography (OCT) measurements has been performed, showing no statistically significant difference in the luminal cross-sections reconstructed with our method, compared to OCT.


I. INTRODUCTION
I NVASIVE Coronary Angiography (ICA) is the most commonly used imaging modality for the detection of coronary stenoses. Its advantages include simplicity, high spatial and temporal resolution of lumen structure, and most importantly, its utility to guide coronary interventions in real time [1]. However, despite these clinical advantages, X-ray angiograms pose several challenges, especially in relation to visualizing lesion adequately and judging lesion severity. The 2D projections of 3D vascular structure in different image planes produce vessel overlap and foreshortening and hence, make it difficult for the cardiologists to interpret the geometry. This leads to high inter-and intra-observer variability in understanding the global anatomical structure of the 3D vascular tree and, in turn, affects the accuracy of the estimation of lesion severity and stent size [2]. The interpretation gets further complicated due to the existence of several motion artifacts, including cardiac, respiratory, and patient or device movement. In addition, potential adverse effects of higher amount of radiographic contrast agent and exposure to X-rays limit the number of image acquisitions. To overcome these inherent limitations of ICA, 3D reconstruction of coronary arterial (CA) trees from a limited number of 2D X-ray projections has been attempted by a number of research groups [3]- [6].
Several 3D CA tree reconstruction methods have been developed throughout the literature that deal with problems such as vessel overlap, foreshortening, suboptimal projection angles, and tortuosity and eccentricity [3], [7], [8]. These methods either try to generate a model of the 3D structure This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ of the coronary arteries [7], or reconstruct the 3D volume of X-ray attenuation coefficients using tomography [8]. In order to obtain the acquisition geometry, some of the reconstruction methods rely on a prior calibration step [9]. Alternatively, some reconstruction methods prefer non-calibrated data in order to reduce possible table movements during image acquisition and noise in calibrated parameters [6], [10]. These methods usually estimate geometry parameters before computing the reconstruction or simultaneously estimate the geometry parameters and the optimal reconstruction. In cardiac interventions, single-plane systems are usually preferred over biplane ones due to their lower cost, ease of use, and the possibility of acquiring multiple and individually tailored image projections. However, single-plane systems generate non-simultaneous angiographic projections, which are prone to several motion artifacts, including mainly patient or imaging device related movements, respiratory motion, and cardiac motion. Although the patient or imaging device related movements can be minimized by following a careful protocol during image acquisition [6], [11], [12], this is hard to achieve for respiratory motion and not possible for heart beating motion. Even if short acquisitions can be completed during a breath-hold, potential misregistration artifacts between breath-holds remain, affecting geometry conditions [13]. Additionally, the patients are often not able to follow the breathhold protocol.
To overcome the motion artifacts due to cardiac movements, retrospective gating is commonly utilized, where the image frames are selected from the same cardiac phase. This is usually achieved using the ECG signal acquired simultaneously with the angiograms. Schwemmer et al. [14] developed an algorithm for residual motion correction, based on deformable 2D-2D registration between forward projection of initial, ECG-gated 3D reconstruction and the original 2D projection. Rohkohl et al. [15] used a time-continuous B-spline motion field model, parameterized by the acquisition time, for motion estimation and 4D vascular reconstruction and validated its performance on a left coronary artery (LCA), a right coronary artery (RCA), and a coronary sinus. Unberath et al. [16] estimated the respiratory motion by optimizing epipolar consistency among all images in projection domain. Their evaluation over two numerical phantoms and three clinical cases produced average reprojection errors of 0.78 ± 0.31 mm and 1.60 ± 0.233 mm, respectively, among key points. In another work, Unberath et al. [17] developed an approach for prior-free respiratory motion estimation by optimizing epipolar consistency conditions and a task-based auto-focus measure. Their experimental analysis over a phantom and three clinical images achieved average reprojection errors among key points 1.17 ± 0.05 mm and 1.91 ± 0.393 mm, respectively.
In most of the existing reconstruction methods, the frames are selected from the end of the diastolic phase when the heart motion is minimal, such that no cardiac motion can be assumed between image projections [18], [19]. However, in a single-plane rotational X-ray system, this assumption does not hold due to non-simultaneous acquisition. In [20], 3D reconstruction was attempted as an energy minimization problem, solved using α-expansion moves of graph cuts; while a snake-based semi-automated reconstruction method was proposed in [21]. Point correspondences between arteries in monoplane image pairs were refined in [22] by discarding outliers using RANSAC and filtering the remaining inliers with a geometric curvature constraint. A perspective projection model-based 3D CA reconstruction from two uncalibrated monoplane angiograms was developed in [23], introducing a new model for calculation of contour points in the vascular surface. An external force back-projective model was integrated into the deformable model framework in [10].
Recently, a nonuniform rational basis splines (NURBS) based CA reconstruction method was proposed where the 3D centerline is reconstructed as the intersection of surfaces from corresponding branches [4]. The point correspondences are identified at uniform distances on the projected planes and the 3D luminal contours are generated using NURBS contours over the computed 3D boundary. The final 3D surface is generated as an interpolation across control points on the 3D contours. Another NURBS-based 3D CA tree reconstruction from two X-ray projections was recently introduced in [12] with breath-hold image acquisition.
Blondel et al. [24] used 3 − 7 angiographic projections for identifying correspondences between 2D centerlines. Based on N projections, the approach generated N(N − 1) reconstructions and fused them in a single set. Although the approach considered the effects of respiratory motions between acquisitions, it did not consider non-rigid deformations and identified point correspondences based on epipolar constraints. This could result in falsely reconstructed vessel segments and discontinuities along centerlines. More recently, Unberath et al. [25] applied the same approach of [24], combined with graph cuts. After merging N(N −1) reconstructions from N projections, the outliers were removed based on reprojection errors on additional reference projections. From the final set of 3D points, the vascular tree was reconstructed based on minimum spanning tree. Jandt et al. [26] used multiplicative combination of back-projected 2D vesselness responses to compute 3D vesselness response for every voxel. Fast marching was applied for vessel segmentation and backtracing for centerlines extraction. Since the approach did not distinguish between point or line correspondences initially, it extracted a large number of vessel segments, with only a few of them representing the existing vessels. The reconstructed segments were assigned a significance score, based on the length and 3D location, for final 3D centerlines reconstruction.
Although the problem of reconstructing 3D coronary vessel trees has been investigated through the last two decades, several key issues remain with the existing methods. Most of the algorithms require more than two projections. Some of the algorithms assume non-existence of motion artifacts and require breath-hold acquisition and no patient movement in single-plane X-ray systems [6], which makes them often unsuitable during cardiac interventions. Some algorithms are only applicable on biplane angiograms and do not involve any motion correction or geometry calibration steps [27]- [30]. Many existing 3D reconstruction methods provide only qualitative evaluation of their performance [11], [28] or evaluate with respect to physiological information [4]. Although some methods have quantitatively evaluated their reconstruction performance over synthetic or clinical images, their average reprojection error is often high enough to introduce significant distortions or discontinuities along vessel geometry [10], [22].
Among all the artifacts that degrade the reconstruction of 3D CA trees, the possibility of generating falsely reconstructed vessels is the most challenging. The geometry of the 3D CA tree depends on its skeleton, reconstructed in 3D space based on the intersections between projection lines joining the X-ray sources and the corresponding points on 2D centerlines in the projected image planes [4], [10], [12]. These corresponding points can be identified by finding the nearest projection line that joins the X-ray source to a point on the 2D centerline in another projection plane [12], [31]. In this way, the 3D reconstruction can be seen as one of calculating the intersection curves between 3D projection surfaces corresponding to centerlines in projection images. If the centerlines are straight, this is straightforward. However, coronary vessels often create C-shapes along their length, generating multiple intersecting points (or curves) between projection surfaces. Only one of the mathematical solutions of the rays-to-rays (or sheet-to-sheet) problem represents the true vessel. In Fig. 1, an example of the existence of falsely reconstructed vessels is depicted for 3D reconstruction of the RCA. From the figure, it can be observed that the projection surfaces intersect at more than one 3D location, only one of which is the true one. Since the RCA produces a C-shaped structure in all of its angiographic projections, most pronounced in the left anterior oblique (LAO)-straight and anterior-posterior (AP) cranial planes, this problem is common.
The purpose of our work is to generate a 3D coronary arterial tree from multiple angiographic projections without the need for any specific image acquisition protocol. The proposed method first addresses the issue of motion artifacts in angiographic projections. It estimates an optimal rigid transformation from each of the angiographic acquisitions to adjust the relative rigid motion, mainly due to respiration and patient or device movements. The remaining non-rigid distortion at end-diastole, mostly due to cardiac motion, is modelled by a radial basis function based warping of the 2D vessel centerlines. Next, the method addresses the problem of falsely reconstructed vessels in 3D CA reconstruction. A novel point-cloud based approach is introduced for the reconstruction of 3D centerlines that minimizes the reconstruction error (the sum of 2D reprojection errors from all projection planes) and produces the optimally reconstructed 3D vessel skeleton. The performance of the proposed 3D reconstruction technique has been qualitatively and quantitatively evaluated for the reconstruction of 109 vascular trees, including left anterior descending (LAD), left circumflex (LCx), and RCA, from 45 patients admitted to hospital for suspected coronary stenosis. The proposed rigid motion correction algorithm results in average reprojection error of 0.448 ± 0.396 mm, while the proposed 3D centerlines reconstruction algorithm generates an average reprojection error of 0.092 ± 0.055 mm. The performance has been further evaluated by reprojecting the generated 3D vasculatures on additional projection planes, not used for reconstruction, and it produced average reprojection error of 0.910 ± 0.352 mm in 16 cases, available from our set of 45 patients. We have also validated the reconstructed 3D lumen surfaces using OCT imaging, available in 5 patients from our set of 45. The results show that there is no statistically significant difference in the reconstructed 3D luminal crosssections when compared to OCT.
The structure of the rest of this paper is as follows: The proposed rigid and non-rigid motion correction techniques are discussed in Section II. Section III mathematically proves the existence of falsely reconstructed vessel segments and introduces a novel point-cloud based approach for optimal reconstruction. Section IV validates the proposed 3D CA tree reconstruction algorithm, along with the proposed motion correction, over 45 patients with suspected coronary diseases. Discussion and concluding remarks are provided in Section V.

II. PROPOSED MOTION CORRECTION
The first step of our pipeline aims at correcting motion between angiographic projections. Using the simultaneously acquired ECG signals, we first select the end-diastolic frames as corresponding frames from angiographic sequences. An end-diastolic frame is identified by its amplitude lying within Q and R waves and closest to the R wave. Before applying the proposed motion correction and 3D centerline reconstruction algorithms, the 2D images are preprocessed for the extraction of relevant image features. The input images are first processed using a classic Hessian-based multiscale filter, which identifies tubular structures and removes background anatomical structures [32]. The resulting vesselness image is then applied for automatically segmenting vessels from the background using the active contour region-growing algorithm by Chan and Vese [33] (maximum number of iterations 200 and weight of smoothing parameter 0.2, where these parameters were tuned by exploring typical acquired clinical images). Finally, the 2D vessel centerlines are extracted from segmented vessels, along with the corresponding boundary points, using a fast-marching level set method [34].

A. Rigid Motion Correction in Object Domain
One of the major issues in 3D reconstruction of CA trees is that the acquired angiogram images from different views are misaligned. Even after selecting the end-diastole frames for non-simultaneous image acquisition, the selection of image frame from a discrete set retains certain amount of temporal misalignment, leading to deformations in the coronary vessels. Rigid movement of the heart, and hence, the vessels, can mostly be attributed to the respiratory movement of the lung. Patient or device-related movements during or between image acquisitions also cause similar rigid movement.
In order to remove the effects of the rigid motion, the proposed approach aims at estimating the optimal rigid transformations for all angiographic image acquisitions, which minimize the orthogonal distance between corresponding landmarks in object domain. The proposed rigid motion correction procedure requires the identification of a few (preferably 4−6) corresponding landmark points (in our current implementation, vessel bifurcation points) from all projection planes.
Let us assume the rigid-motion corrected 3D locations of r landmark points are given by B i for i = 1, · · · , r and the original location of landmark B i on the j th projection plane is b i j for j = 1, · · · , s. Our aim is to estimate B i , given known b i j , j = 1, · · · , s, ∀i = 1, · · · , r . (Note that, throughout the paper, 3D points on projected planes are represented by lower-case letters, while other 3D points are denoted by an upper-case letter.) Let us also assume the locations of the X-ray source, a representative point on the projection plane j and the normal to the projection plane are defined as F j , M j , and N j , j = 1, · · · , s, respectively. Since the motion corrected 3D locations B i are not known, these 3D landmark locations are initially estimated as the nearest point to all 3D projection lines where Throughout the rest of the paper, we refer to any 3D point minimizing the sum of orthogonal distances between lines or planes as their nearest orthogonal point. Let us also define B i j as the nearest point on − −− → F j b i j from B i . B i j is the location of the rigid-motion affected 3D landmark B i during j th image acquisition. Our objective is to estimate the optimal translation t j and rotation R j for each acquisition j = 1, · · · , s, so that the 3D landmarks B i j , i = 1, · · · , r match with the corresponding 3D landmarks B i , i.e., the effects of rigid motion during acquisitions is minimized: The optimization in Eq. (2) is solved using Horn's quaternion-based method [35]. Since the original 3D locations of the rigid motion corrected landmark points are not known, an iterative algorithm is developed, where, in each iteration, the 3D landmarks are generated using Eq. (1) and the optimal rigid transformations are estimated using Eq. (2). In each iteration, the sources and landmarks on projection planes are updated based on the rigid transformation (3)  The procedure is iterated until the nearest 3D landmarks for each image acquisition coincide or the difference between updated 3D landmarks in two successive iterations is smaller than a pre-defined threshold (we used 0.0001 mm). An example of the result of this algorithm for two projections (AP-cranial and RAO-cranial) of the LCA of Patient 2 is depicted in Fig. 2. The pseudo-code of the proposed algorithm is provided in the Appendix file [36].

B. Correction of Non-Rigid Motion Artifacts
Although the proposed rigid motion correction can remove most of the motion artifacts from angiogram images, non-rigid motion may still remain due mostly to cardiac contraction and relaxation. In order to model this non-rigid motion, our method employs the radial basis function based image warping technique, proposed by Bookstein [37].
The objective is to identify the deformation function f that transforms the manually selected landmark locations b i j in a 2D image plane to the back-projections b R i j of the rigid motion corrected 3D landmarks B i , i = 1, · · · , r on the same plane. Since B i s are the nearest 3D orthogonal points to the manually selected landmark locations, their reprojections on image planes represent the corresponding landmarks without any motion artifact (including both rigid and nonrigid movements). The assumption is that the deviation from the reprojections of rigid motion corrected 3D landmarks on projection planes to the corresponding manually selected landmarks is due to the underlying non-rigid deformation of vessels, and can thus be applied to correct the vessel centerlines.
For notational simplicity, let So, the warped 2D centerline, removing non-rigid motions at the j th projection plane, is given by where c j is the 2D centerline on the rigid-motion corrected projection plane j and Applying the warping function, defined in Eq. (4) for each projection plane, to the 2D centerlines in those projection planes, corrected 2D centerlines are generated [38].

III. PROPOSED 3D CORONARY ARTERIAL TREE RECONSTRUCTION
As already described in Section I, one of the main objectives of this paper is to solve the challenge of falsely reconstructed vessel segments. This subsection theoretically derives the existence of this specific challenge and then develops a novel 3D centerline reconstruction algorithm that can solve this problem, as well as several other X-ray imaging artifacts.

A. Theoretical Proof of the Existence of Falsely Reconstructed Vessel Segments
Since the objective of this subsection is to prove the existence of this reconstruction challenge, the generation of an incorrect vessel segment as well as its true counter-part using a simple vessel-like structure would suffice. The problem is most visible during reconstruction of curved, C-shaped vessels such as the RCA. So, in order to mimic a C-shaped vessel in 3D space, a parabola is first drawn in the 2D plane (z = 0) and then an arbitrary rotation is applied. Although the following derivation has been presented on a simple parabolic curve, it can be shown for other curves, such as higher-order polynomials. Let, the 2D curve be: y 2 = ax. So, the C-shaped vessel in R 3 is represented as To contain the curve in a finite-window, y is assumed to lie in [−k, k], where k is any positive real number. An arbitrary is applied to the curve, which generates the final 3D object: Let us now assume the X-ray source points be S 1 and S 2 , the mid-points of the projection planes be M 1 and M 2 , and the normals to the projection planes be N 1 and N 2 . For line passing through S 1 and A Since, P 1 lies on the projection plane with mid-point M 1 and normal N 1 , Similarly, for the projection of A on plane (M 2 , N 2 ), Hence, the equations of projected curves are given by: Now, lines passing through source points and points on projected curves are given by: In case of pointwise correspondence, The derivation of (10) is provided in the Appendix file [36].
A second solution will exist only if a Q 3 −y 1 Q 2 Q 2 −ay 1 Q 1 ≤ k for |y 1 | ≤ k. Hence, in case the 3D curve is quadratic in the initial X −Y plane, there will always exist a mathematically accurate but anatomically infeasible solution. For higher even-order polynomial curves, at least one mathematically accurate but anatomically infeasible solution will always exist (assuming it lies within finite window). However, these anatomically infeasible solutions may or may not exist for higher odd-order polynomial curves, since all other solutions can be imaginary.
The above proof of existence of falsely reconstructed vessel segments considers a simple example without any motion artifact. However, if motion artifacts (rigid or non-rigid) exist between projected image planes, the perfect pointwise correspondence between vessel segments, as formulated in (10), will sometimes not exist. The 3D projection lines, connecting the vessel centerlines on projection planes with the X-ray source, will often not intersect in that case and hence, create ambiguity in reconstructing the corresponding 3D point on the vessel tree. In the next subsection, a novel procedure is developed that can solve the problem of falsely reconstructed vessel segments during 3D centerlines reconstruction, even in the presence of significant motion artifacts.

B. Novel Point-Cloud Based Approach for 3D Centerline Reconstruction
Although the proposed motion correction algorithms are able to significantly reduce the effects of rigid and non-rigid motions from the 2D vessel centerlines, significant distortions still remain in the 2D vessel centerlines, which create ambiguity in identifying point-to-point correspondences as well as in generating intersections between projection surfaces. Our approach first attempts to identify all possible point correspondences from the intersections between projection surfaces and produces an initial dense point-cloud in the 3D space. The point-cloud of 3D points is a data set containing the nearest orthogonal points of the projection lines (both intersecting and non-intersecting) generated from the points on 2D vessel centerlines in projection planes to the corresponding X-ray sources. Since the point correspondences are not established at the initial stage of the algorithm, all possible nearest orthogonal points are considered at this stage. Hence, if there are, respectively, M and N number of points on 2D vessel centerlines in two projection planes, the initial point-cloud will consider M.N number of points.
From this initial point-cloud, redundant and irrelevant representative points are then removed. Those 3D points are identified where the distance between projection lines connecting X-ray sources to the corresponding points on 2D centerlines, that is, − −− → F j c j i , j = 1, · · · , s, i = 1, · · · , n j , is larger than a fixed value ; i.e. the 3D point ( j, i ; j i ) is removed if We estimate the value of the threshold as the average maximum diameter of the vessel being imaged.
To consider the issue of foreshortening of vessels, the proposed method allows many-to-one point correspondences between centerlines on projected planes. However, considering all many-to-one correspondences from each plane unnecessarily crowds the point-cloud. Hence, for each point along the centerline from every projection plane, at most k nearest point correspondences (with distance less than ) are included. This assumption significantly reduces the number of points in the point-cloud and provides a more compact outer hull for the generation of 3D centerline. The final refinement of the point cloud takes into account the relative location of potentially corresponding points on each projection curve. If these relative locations differ by more than a predefined threshold, the corresponding 3D orthogonal point is removed from the point cloud.

Algorithm 1 Generating 3D Centerline From Point-Cloud
Input : 2D centerline of vessel of interest from each projection plane Output: 3D centerline of the vessel 1 For every point along the 2D centerline in each projection plane, find the point of intersection between projection lines with every point along centerline in other planes, or the nearest orthogonal point in case of non-intersecting projection lines; 2 Retain points with orthogonal distance between projection lines less than the average maximum diameter of the vessel; 3 Compress the point-cloud by retaining some of the nearest points (in terms of orthogonal distance between projection lines) for each point along 2D centerline in each of the projection planes and discarding the rest; 4 Initialize b to a small integer; 5 do 6 Draw a B-spline curve with b number of breaks through the point-cloud; 7 If any point in the cloud lies at more than 3 times average distance from the fitted spline curve at that point, consider it as an outlier and discard it from the point cloud; 8 b ← b + 5; 9 while there are outliers; 10 Return the generated 3D curve as the 3D centerline of the vessel In this way, we avoid false reconstructions generated from the intersection of proximal with distal locations. Through the optimized point-cloud, a cubic B-spline curve is fitted with b 0 number of breakpoints (knots). The least squares method is used to fit the spline to the noisy data, where the smoothing effect is controlled by the selection of uniform breaks. After the first B-spline curve has been fitted to the initial set of points in the 3D point-cloud, the outliers are removed from the 3D point-cloud. In our model, a 3D point is considered an 'outlier' if the orthogonal distance between the point and the spline curve is greater than a fixed distance, which we set as d times the average of the Euclidean distances of the points on the 3D spline curve to the points in the point-cloud that are orthogonal to the curve at that point. In the following step, another B-spline curve is drawn through the reduced point-cloud set with the number of breakpoints b increased by a small number z, thus reducing the oversmoothing of the 3D centerline. The outlier points are again removed from the point-cloud set and the process is repeated until no outliers are found. The pseudo-code of the proposed 3D reconstruction technique is presented in Algorithm 1.
The generated 3D centerlines provide a visually satisfactory representation of the 3D skeleton of the vessel tree. However, the 3D skeleton generated using Algorithm 1 can still contain falsely reconstructed vessel segments, even after the point-cloud refinement steps described in this Section.
Regions with falsely reconstructed vessel segments can be identified because the generated 3D B-spline will fall between true and false nearest orthogonal points. As a result,

Algorithm 2 Generating Optimum 3D Centerline From Point-Cloud by Minimizing Reprojection Error
Input : 2D centerline of vessel of interest in each projection plane Output: 3D centerline of the vessel 1 do 2 Run Algorithm 1 to generate the 3D centerline of the vessel; 3 Back-project the 3D centerline on the original projection planes; 4 Compare with 2D centerlines: if the back-projection of a point perfectly matches with 2D centerline at all projection planes, assign 1, otherwise assign 0; Construct the point-cloud set combining both correctly and incorrectly reconstructed segments; 9 while back-projected points do not coincide with the 2D centerlines on projection planes and the maximum number of iterations has not reached; 10 Return the generated 3D centerline of the vessel the backprojection of this curve will fail to accurately match the centerlines in 2D projection planes. Thus, we backproject the generated 3D centerlines on the original 2D projection planes and mark the curve segments, where the distance between them is larger than a threshold, as incorrectly reconstructed. The correctly reconstructed segments can now be used to find one-to-one correspondences between projection points, which are used to further refine the point-cloud at these regions. This in turn allows the fitting of a B-spline with an increased number of break points. The process is iterated, gradually reducing the incorrectly reconstructed regions until the back-projections of generated 3D centerlines accurately match with the corresponding 2D centerlines at all projection planes, thus minimizing the average reprojection error. The pseudo-code of the proposed optimum 3D centerlines reconstruction is presented in Algorithm 2.
The algorithm described up to this point reconstructs the geometry of the centerlines only. For the complete 3D CA tree reconstruction, we use the vessel surface reconstruction method introduced by Galassi et al [4]. This uses the non-uniform rational basis splines (NURBS) method for the generation of luminal cross-sections along centerlines, parameterizing a luminal contour over boundary points Q i of a vessel in different projections, as follows: where p is the degree of the curve and B i, p (·) is a rational basis function of Q i of pth degree. The final 3D lumen surface is generated by applying the NURBS surface formulation: The cross-sections of diseased coronary arteries vary from roughly circular to various arbitrary shapes, depending on the extent and orientation of atherosclerotic plaques. Hence, we have used the NURBS model and parameterized the luminal contours and lumen surface to produce a flexible and versatile model, capable of representing from simple curves (e.g. circles, ellipses) to more complex free-form ones, depending on the boundary points in the lumen. In our implementation, the initial smooth B-spline curve is generated using 20 control points (b 0 = 20). The set distance for outlier removal is defined as 3 times the average of the Euclidean distances of the points on the 3D spline curve to the points in the point-cloud that are orthogonal to the curve at that point (d = 3). Algorithm 1 produces the final smooth 3D centerline from point-cloud by iteratively increasing the number of control points b by 5 (z = 5), while simultaneously removing the outliers from the point-cloud. Although this procedure usually converges after 10 − 15 iterations, the iterative increase of control points makes the algorithm susceptible to the overfitting of the B-splines. Hence, in order to produce anatomically plausible representation of 3D centerlines, we restrict the maximum number of control points to 50 in our implementation of Algorithm 1. Algorithm 2 for generating the optimal 3D centerlines by minimizing the reprojection errors does not take more than 3 iterations to converge.
The proposed approach is executed in Ubuntu 16.04 LTS 64-bit OS on an Intel(R) Core(TM) i5-7200U CPU @ 2.50GHz × 4 and 16 GB RAM. In our non-optimized Matlab implementation of the complete 3D CA tree reconstruction, the initial 2D preprocessing step (consisting of tubular structures identification, vessel segmentation, and 2D centerlines extraction) requires 45−60 seconds of runtime. The next step, motion correction (both rigid and non-rigid) takes approximately 10 − 15 seconds. The proposed point cloud-based 3D centerlines reconstruction procedure takes 2.5 − 3 minutes to run, where approximately 90% of the run-time is used to generate the point-cloud. The final steps of reconstructing 3D luminal cross-sections and 3D lumen surface require relatively longer time (about 10 − 12 minutes) depending on the length and number of coronary vessels. If, instead of using the NURBS fitting, the luminal cross-sections and lumen surface are fitted with circular contours, the run-time would be brought down to 1 − 1.5 minutes. An example of the 3D reconstructed LAD, LCx, and RCA, along with the complete reconstructed 3D CA tree (of Patient 2), is qualitatively depicted in Fig. 5.

IV. EXPERIMENTAL ANALYSIS
The performance of the proposed 3D CA tree reconstruction algorithm has been qualitatively and quantitatively evaluated on 45 patients enrolled in clinical studies (including 5 patients  with OCT imaging). For each patient, 2−3 X-ray angiographic projections are captured for imaging RCA and 3−6 projections are captured to image LCA. Since the left anterior descending (LAD) and left circumflex (LCx) arteries are usually not visible optimally on the same angiographic projections, different sets of LCA projections are used for reconstructing LAD and LCx. Among the total of 45 patients used for evaluation of 3D centerlines reconstruction, the LAD artery is visible from at least two projections for 40 patients. In case of the LCx artery, at least two X-ray projections are available from 37 patients. Similarly, in 32 out of 45 patients, the RCA is captured in at least two angiographic projection planes. Since our objective is to reconstruct the complete CA trees, the 3D reconstruction procedure of LAD, LCx, and RCA vessels also incorporates the attached sub-branches. Hence, the reconstructed 3D tree of LAD artery includes the diagonal and septal vessels, while the reconstructed 3D tree of LCx artery includes the marginal vessels. Similarly, the 3D reconstructed RCA tree includes the posterior descending artery, posterior left ventricular branch, acute marginal arteries, etc.
Since the objective of the proposed approach is to accurately reconstruct the geometry of the 3D coronary arterial tree, its performance has been evaluated in two separate ways. The 3D centerlines reconstruction is evaluated by backprojecting the 3D centerlines on the original projection planes, as well as on any additional projection plane not used for the 3D reconstruction. The reconstruction performance of 3D lumen volume is evaluated by comparing the reconstructed luminal cross-sections with the same derived from OCT.

A. Evaluating the Performance of Rigid Motion Correction
The qualitative performance of the proposed rigid-motion correction on two LCA image projections is visible in Fig. 2. The performance of the proposed rigid-motion correction method is quantitatively evaluated by comparing the 2D locations of landmark points (in this case, bifurcation points) on each projection plane with the back-projected locations of rigid motion-corrected 3D landmark points. For quantitative evaluation, the average reprojection errors, i.e. the average rms errors of the landmark points from all X-ray projections, are measured and the results are presented in Fig. 6 for LAD and LCx arteries and in Fig. 7 for RCA. The average reprojection error for LAD artery in 40 patients is 0.498 mm, while the same for LCx artery in 37 patients is 0.336 mm. Since the effect of non-rigid motion artifacts is relatively higher in the RCA, the average reprojection error for RCA in 32 patients is 0.517 mm. In total, the average reprojection error, combining the LAD, LCx, and RCA trees in 45 patients, is measured as 0.448 ± 0.396 mm.
The effectiveness of the proposed rigid motion correction approach is evaluated for 3D centerlines reconstruction. The existing 3D centerlines reconstruction method proposed in [4] by Galassi et al. is used for both rigid motion corrected and not rigid motion corrected X-ray angiographic projections and the comparative performance analysis is depicted using violin plots in Fig. 9, which demonstrate the kernel density of the data distribution allowing the visualization of multimodal distributions of reprojection errors. In all cases, the tails of the violins are trimmed to the range of the data distribution and scaled to maximum width of 1.
From the violin plots presented in Fig. 9, it is clearly observed that the use of proposed rigid motion correction significantly improves the accuracy of 3D centerlines reconstruction, in terms of reprojection errors, for all arterial trees (LAD, LCx, and RCA). The improvements in the reconstruction performance are statistically significant with respect to parametric paired-t test, where 0.01 is the desired level of significance. The performance improvement for all arterial trees of each of the patients is individually presented as violin plots in the Appendix file [36]. After the proposed rigid motion correction, the 2D centerlines are warped to correct the remaining non-rigid deformation in each plane, qualitatively depicted in Fig. 3 for Patient 2.

B. Performance Evaluation of 3D Centerlines Reconstruction
The performance of the proposed 3D centerlines reconstruction approach is evaluated by comparing the reprojected 3D centerlines on each of the angiographic projection planes with the original 2D centerlines on those planes. The performance is quantitatively measured using the reprojection errors and presented as box plots in Fig. 10. In our analysis over 45 patients, the proposed approach produces average reprojection error of 0.081 mm for 3D centerlines reconstruction of 40 LAD arteries and average error of 0.085 mm for reconstruction of 37 LCx arteries. The average reprojection error for the 3D reconstruction of 32 right CA trees is 0.113 mm. Overall, the average reprojection error for 3D reconstruction of 109 vascular trees from a total of 45 patients is 0.092 ± 0.055 mm. The centerlines reconstruction performance on all 45 patients is compared with the existing 3D CA reconstruction algorithm by Galassi et al [4], after motion correction using our proposed rigid motion correction algorithm. The overall result on 109 arteries, presented in Fig. 10, clearly demonstrates the substantial reduction in the reprojection errors. The improve- Fig. 9.
Comparison between reprojection errors, measured after reprojecting the reconstructed 3D centerlines as in Galassi et al., with and without the rigid motion correction. ments are also statistically significant with respect to paired-t test. For visual comparison, the upper limit of Y-axis in Fig. 10 is fixed to 10. The reconstruction performance on all arterial trees of each of the patients is individually presented as violin plots in the Appendix file [36].
The qualitative performance of the proposed 3D centerlines reconstruction algorithm is presented in Fig. 8. The first column presents the reconstructed 3D skeleton at one 3D view, while the next two columns depict the reprojection of the 3D centerlines on two original projection planes. The images show an accurate reprojection of the reconstructed centerlines. The quality of the reconstructed vasculatures has also been visually assessed and confirmed by two expert cardiologists. It should be noted that, in order to remove the non-rigid motion artifacts, the proposed approach warps the 2D vessel centerlines in each plane and slightly deviates from the vessel locations. Hence, if significant amount of non-rigid motion artifacts exist between projection planes, the reprojections of generated 3D centerlines on 2D projection planes will slightly deviate from the vessel centerlines, as visible in the third row of Fig. 8. However, the reprojection errors in this case will not be affected, since they are computed based on motion-corrected 2D centerlines.

C. Performance Evaluation Over Projection Planes Not Used for 3D Reconstruction
Since the proposed 3D reconstruction algorithm can effectively generate the 3D CA tree based on only two angiographic projections, any additional image projection, if available, can be applied to qualitatively, as well as quantitatively, validate the performance of the proposed technique. Measurements of the reprojection errors are obtained by back-projecting the reconstructed 3D centerlines on the additional projection planes and the results are quantitatively depicted as violin plots in Fig. 11. The qualitative results are presented in column 4 of Fig. 8. From the total of 45 patients, more than two angiographic projections were available in only 16 cases (4 LAD, 7 LCx, and 5 RCA). For objective qualitative evaluation of the reconstruction performance, we present a sample of 4 cases for visualization in Fig. 8. We have selected these 4 cases so that their reprojection accuracy on the third plane is representative of the whole 16-case set (0.907 ± 0.354 mm vs 0.910 ± 0.352 mm.) From the results presented in column 4, it can be observed that, even without using any information from that projection plane, the proposed 3D reconstruction algorithm is able to quite accurately produce the centerlines on that image plane. In this regard, it should be noted that the additional projection planes were not corrected for either rigid or non-rigid motion artifacts, which explains minor deviations in the reconstructed centerlines. The application of those steps should be able to significantly reduce the reprojection errors in each case. The remaining 3D reconstruction results are qualitatively depicted in the Appendix file [36].

D. Validation of the 3D Reconstructed Lumen Volume Over Optical Coherence Tomography
In order to evaluate the accuracy of the reconstructed 3D lumen volumes, we compared our approach with optical coherence tomography (OCT), an invasive medical imaging modality providing high-resolution endoluminal visualization of the lesion anatomy [39]. OCT applies back-reflected infrared light to perform in situ micron scale tomographic imaging of the vessel anatomy and internal microstructure of plaques [40]. For 5 patients, out of our total set of 45 patients enrolled in clinical studies, the OCT catheter (2.7F) was automatically pulled back at speed of 20mm/sec and rate of 100 frames per second, with manual injection of 20ml of contrast solution. The lumen cross-sectional areas were manually measured from every frame of the OCT video of each patient. From the corresponding X-ray angiographic projections of each patient, the starting position of the OCT catheter is visually located at the distal part of the vessel in one projection and the end position is identified based on the length of the pull-back of OCT catheter on the 3D reconstructed vessel centerline. Since the OCT frames were generated by pull-back with uniform speed, the length between the estimated starting and end positions on our 3D reconstructed vessel is uniformly partitioned and the cross-sectional areas from the 3D reconstructed lumen (Eqn. (12)) are computed. The quantitative results are presented in Figs. 12 and 13 for both OCT derived luminal areas and luminal cross-sections derived by the proposed 3D reconstruction.
From the results presented in Figs. 12 and 13 for 5 patients, it is observed that the reconstruction from the proposed algorithm accurately matches the OCT derived luminal crosssections. For statistical analysis of the comparison between OCT and our method, the χ 2 -test is applied to test the null hypothesis: the luminal cross-sectional areas derived by our proposed reconstruction is in agreement with the same derived from OCT in each frame. The calculated p-values are substantially above the usual significance threshold (they are higher than 0.85 in all cases), inferring our method is not significantly different from OCT.

V. DISCUSSION AND CONCLUSION
Several 3D CA tree reconstruction methods have been developed in the past two decades. However, many of the methods have not provided any quantitative evaluation of their reconstruction performance. Many of the existing methods followed breath-holding and no patient movement protocol during image acquisition. A few reconstruction methods that have quantitatively demonstrated their performance on clinical images suffer from sub-optimal reconstruction in terms of average reprojection error. Also, to the best of our knowledge none of them has explicitly addressed the problem of falsely reconstructed vessel segments.
In this regard, the objective of our work is to develop a new 3D CA tree reconstruction algorithm, maintaining the standard clinical image acquisition procedure, without the need for any additional constraint to reduce motion between acquisitions. The proposed approach first aims to remove the effects of both rigid and non-rigid motion, so that the point correspondences can be identified for 3D centerline generation. Next, it introduces and theoretically proves the problem of falsely reconstructed vessel segments in 3D CA reconstruction. A novel point-cloud based approach is then proposed for reconstruction of 3D centerlines that takes into consideration the problem of false vessel reconstruction, as well as foreshortening. The proposed algorithm minimizes the average reprojection errors to produce the optimally reconstructed 3D vessel skeleton. Experimental analysis over 45 patients has generated the average reprojection error of 0.092 ± 0.055 mm for 3D centerline reconstruction. The performance has been further evaluated by reprojecting the 3D centerlines on projection planes that were not involved in the 3D reconstruction procedure and it has produced average reprojection error of 0.910 ± 0.352 mm in 16 cases available from our set of 45 patients, even without discarding the motion artifacts.
One of the main advantages of our proposed reconstruction method is that it can be applied when only two image projections are available. Most state-of-the-art 3D reconstruction methods, such as [24] and [25], only work when at least three image projections are available, which is often infeasible in clinical practice (mostly for RCA). Although the method of [26] is applicable on two projections, it is suggested to be used on multiple angiographic projections. The second advantage of our proposed approach is that it can take care of falsely reconstructed vessel segments, as well as non-rigid deformation in the vessels, by identifying optimal point correspondences iteratively from a point-cloud of all possible correspondences. In contrast, both [24] and [25] use epipolar constraints for identifying point correspondences. So, even though they apply reprojection errors for correspondence identification using additional projection planes, they generate a significant number of incorrect vessel centerlines. The method of [26] first produces large number of false reconstructions using their approach and then tries to remove them based on a significance score. However, approaches like this are limited in that, when motion is present, vessel branches cannot be classified into true or false: a single branch might contain parts from the true reconstruction and parts from the false one.
Our proposed approach partitions the reconstructed centerlines in Algorithm 2 into accurate and incorrectly reconstructed segments, which creates a natural ordering in the 2D centerlines on projected planes. This, in turn, solves the problem of falsely reconstructed vessel segments generated from intersecting projection surfaces. In reality, a 3D vessel segment is falsely reconstructed, most commonly, in regions with small radius of curvature. This is because the projection surfaces generated from the 2D projections of that segment are likely to intersect at multiple 3D locations. For simplicity, if we divide a curved vessel segment (often C-shaped) into two parts: from the start of the segment to the 'vertex' and from 'vertex' to the end of the segment, the false solutions are generated when a projection line from the first part intersects with a projection line from the second part and vice versa. Since the proposed reconstruction technique partitions the vessel into accurately and incorrectly reconstructed segments, accurate correspondences between vessel segments are established. So, the projection lines generated from the first part of C-shaped vessel segment in one projection are only matched against the same from other projections and produce the accurately reconstructed 3D vessel segment. In addition, the proposed 3D reconstruction technique assumes manyto-one point correspondences between 2D centerlines from all projection planes in order to produce the point-cloud set. Hence, it automatically considers the problem of vessel overlapping and foreshortening in the model. In this regard, it should be mentioned that the proposed algorithm only uses the manually identified landmarks (bifurcation points) for rigid and non-rigid motion correction. The proposed 3D reconstruction algorithm is completely automated and makes no use of these landmarks. In case no motion artifact is present between image projections (e.g images generated from biplane system), the landmark identification step is not required.
The third advantage of our proposed reconstruction method is its performance compared to the existing literature. The results presented in Table I, shows clearly lower reprojection errors on original projection planes. The experimental analysis on a larger (more than 100) dataset, compared to the state-of-the-art methods, also proves the robustness of our method. Note that the proposed method does not depend on identifying point correspondences between projections, but produces point correspondences from the generated 3D centerlines. The proposed algorithm develops a pipeline where the 3D centerlines generation and point correspondences  I  COMPARATIVE 3D RECONSTRUCTION PERFORMANCE OF PROPOSED  ALGORITHM WITH STATE-OF-THE-ART TECHNIQUES identification are simultaneously optimized in an iterative fashion.
Since the proposed approach requires only two angiographic projections for the reconstruction, it generates the luminal contours by NURBS fitting over only 4 boundary points (2N boundary points from N projections). Hence, although the proposed approach can generate a quite accurate representation of the luminal contours, which are also statistically similar to the same derived from the OCT, the availability of any additional projection plane will provide a better representation of lumen contours and in turn, improve the 3D reconstruction of lumen surface. For the 3D centerlines reconstruction, the availability of any additional projection plane should not make much of a difference, as proven in Section IV-C. The proposed non-rigid motion correction algorithm deforms the 2D centerlines based only on rigid motion corrected landmarks and hence, it is capable of producing undesirable warping along the centerlines. However, the aim of non-rigid motion correction step is to deform the centerline, guided by the rigid motion corrected landmarks, for improved estimation of the point-cloud. Since the initial point-cloud considers all points on the 2D centerlines and finds optimal 3D centerlines by minimizing reprojection errors, any incorrect warping is automatically rectified by discarding non-correspondences from the point-cloud, as evidenced by its improved performance.
The development of a completely automated 3D CA tree reconstruction has larger clinical significance. A completely automated method can generate the 3D view of the coronary vasculatures during cardiac interventions and hence, reduce the likelihood of subjective interpretation of 3D CA tree from 2D projections. The generated 3D model can also be applied for quantification of coronary stenosis and hence, can be useful for measuring its severity [4]. The proposed 3D reconstruction method has been developed to solve the reconstruction of motion induced CA trees, but its applicability can be wide, including the reconstruction of cerebral vessels and branching non-cardiovascular structures, such as the biliary tree.