3D reconstruction of coronary arteries from 2D angiographic projections using non-uniform rational basis splines (NURBS) for accurate modelling of coronary stenoses

Objective Assessment of coronary stenosis severity is crucial in clinical practice. This study proposes a novel method to generate 3D models of stenotic coronary arteries, directly from 2D coronary images, and suitable for immediate assessment of the stenosis severity. Methods From multiple 2D X-ray coronary arteriogram projections, 2D vessels were extracted. A 3D centreline was reconstructed as intersection of surfaces from corresponding branches. Next, 3D luminal contours were generated in a two-step process: first, a Non-Uniform Rational B-Spline (NURBS) circular contour was designed and, second, its control points were adjusted to interpolate computed 3D boundary points. Finally, a 3D surface was generated as an interpolation across the control points of the contours and used in the analysis of the severity of a lesion. To evaluate the method, we compared 3D reconstructed lesions with Optical Coherence Tomography (OCT), an invasive imaging modality that enables high-resolution endoluminal visualization of lesion anatomy. Results Validation was performed on routine clinical data. Analysis of paired cross-sectional area discrepancies indicated that the proposed method more closely represented OCT contours than conventional approaches in luminal surface reconstruction, with overall root-mean-square errors ranging from 0.213mm2 to 1.013mm2, and maximum error of 1.837mm2. Comparison of volume reduction due to a lesion with corresponding FFR measurement suggests that the method may help in estimating the physiological significance of a lesion. Conclusion The algorithm accurately reconstructed 3D models of lesioned arteries and enabled quantitative assessment of stenoses. The proposed method has the potential to allow immediate analysis of the stenoses in clinical practice, thereby providing incremental diagnostic and prognostic information to guide treatments in real time and without the need for invasive techniques.


Results
Validation was performed on routine clinical data. Analysis of paired cross-sectional area discrepancies indicated that the proposed method more closely represented OCT contours than conventional approaches in luminal surface reconstruction, with overall root-meansquare errors ranging from 0.213mm 2 to 1.013mm 2 , and maximum error of 1.837mm 2 . Comparison of volume reduction due to a lesion with corresponding FFR measurement suggests that the method may help in estimating the physiological significance of a lesion. PLOS

Introduction
The approach presented in this paper to reconstruct 3D vessel centrelines has similar characteristics to the work presented by Cardenes et al. [10]. Similarly, it considers a centreline segment as a curve joining two nodes (bifurcation points) and it labels segments using the connectivity matrix of the coronary tree image. Additionally, in our method, the ends (nodes) of corresponding branches (segments) are employed as landmark points for estimating the optimal geometry parameters. Importantly, no assumption based on index correspondences is presumed, which can be inaccurate for perspective projections of a tortuous structure, and 3D centrelines of paired segments are reconstructed as intersection of two surfaces.
The reconstruction of the 3D vessel lumen is a crucial step for the assessment of stenosis severity. A common approach is to approximate vessels with tubular objects that are reconstructed by sweeping a circular cross-section along a spine (3D centreline) [4,[10][11][12]. Diameter values are computed on a view as the Euclidean distance between two border points of the vessel centred at a centreline point and then corrected for geometric magnification. More accurate approaches account for the fact that, due to the rotational movements of the X-ray source, the projections used for the calculation are in general not perpendicular to the normal of the cross sectional planes of the centreline [5,8]. Some investigators use diameter information from two views to fit an ellipse to the cross section [13]. Other investigators use the diameter vectors from two views directly as the axes of an elliptic cross-section [14]. However, because of the lack of triple orthogonality, the ellipse can be ambiguous, and ambiguity increases with a decreasing angle between the two diameters. In addition, the non-orthogonality of contours to the vessel centreline usually yields irregularities in the 3D surface (i.e. contours intersecting each other).
Generally, the circular cross-section approximation simplifies and speeds up the reconstruction process, and it may be sufficiently accurate in the case of healthy vessels. However, when coronary arteries are diseased and may contain a significant extent of atherosclerotic plaque, their cross-section can change from roughly circular to various arbitrary shapes, depending on the extent and relative orientation of the plaque [15]. To address this issue, we propose a strategy that aims to model a luminal template cross-section to an extent that is consistent with the stenosis profile on each view, using Non-uniform Rational Basis Splines (NURBS) [16].
In order to provide validation for the approach, we compared lesion assessment using this novel 3D method with Optical Coherence Tomography (OCT), an invasive imaging modality that enables high-resolution endoluminal visualization of coronary anatomy [17]. Additionally, we investigated the correlation between the volume reduction due to a lesion and the Fractional Flow Reserve (FFR) measurement, an invasive technique that uses vessel physiology to gauge the functional significance of a stenotic lesion [18].

Materials and methods Algorithm
The algorithm comprises two main steps: (1) 2D pre-processing of the projections to identify the stenotic vessel (in terms of centreline and profiles), and (2) 3D processing to obtain a 3D model of the vessel and stenotic segment (Fig 1). The main contribution of our work resides in the second step.
2D projections processing: Image enhancement and vessel extraction. Images were firstly pre-processed using a classic Hessian-based multiscale filter, which detects tubular structures while removing background anatomy, such as bones and muscle tissues [19]. The filter was calculated for different scales in order to enhance vessels with different diameters (seven linearly increasing scales ranging between 0.5mm and 4mm, which approximately corresponded to the radius range of projected coronary arteries). The resulting vesselness image, i.e. the maximum response for each pixel, provided a measure for the probability of each pixel to belong to a vessel (Fig 2).
2D vessel boundaries were extracted using a fast marching level set method on the vesselness image, by interpreting vessel boundaries as the propagating front final position [20,21]. The front expanded outward from a seed point towards an end point selected by the user (by simple mouse clicking), and with speed inversely proportional to the gradient magnitude of the vesselness image i.e. within a vessel, the gradient is low, and the front propagates fast, while at the boundaries the gradient is high and the front slows down. A constrained map was defined such that the front propagation was constrained to values within the range 0.1-1 (a low boundary >0 was required to avoid losing regions that show low vessel-like structures due to insufficient contrast). The exact position of the selected initialization points did not influence the behaviour of the fast marching process, as long as they were located within a vessel. Selection of one start and one end points worked robustly for a targeted arterial segment. The user could select multiple initialization points, for detection of multiple vessels. The subsequent analysis proceeded automatically i.e. no additional manual interaction. To compute the vessel centreline a subvoxel precise skeletonization method was applied [22]. The centreline was computed as the minimum-cost path from an end to point to a start point by performing a gradient descent on the fast marching distance map (cost function).
Given the skeleton of a coronary tree in the form of a binary image, then nodes (i.e. branching points) and segments (i.e. subsets of the skeleton that are bifurcation free) were detected and labelled on each projection by exploiting the connectivity matrix of the binary image. To this end, connected pixels were found using an eight-connectivity seed-fill algorithm [23]. The seed-fill algorithm started from a seed point (a foreground pixel) and then iteratively searched its neighbours to detect connected components. For different projections, labels could propagate differently. To ensure correspondence of labels, the initialization points previously selected for vessels extraction (i.e. fast marching seed points) were employed. These points were indicated so that they corresponded across projections. An initialization point was matched to a node based on the Euclidean distance. Only nodes having a matching initialization point were retained. Therefore, node labels across projections were ensured to correspond, and so were branches. Finally, the 2D centrelines were represented as cubic B-splines. A cubic B-spline representation allowed to combine smoothing with the ability to evaluate derivatives, and hence vessel orientations, at any location, thus enabling the steps that follow [24] [25]. 3D reconstruction required definition of a 3D global reference system in which a 3D vessel will be generated. As the geometrical relationship between a radiographic plane (local reference system) and C-arm gantry can be easily described, we defined as global reference system the C-arm acquisition system (an exhaustive description of the C-arm system geometry can be found in [26]). The transformation from the local to the global reference system comprises a translation t followed by a rotation R.
Precisely, an initial estimate of the imaging geometry was obtained from the available image acquisition settings recorded by the X-ray system and stored in the file header of a radiographic image. Then, because small perturbation errors could affect the imaging geometry, i.e. possible table translation during the image acquisition or errors in the parameters, an optimization of the system geometry was necessary [27]. Corresponding anatomical landmarks, i.e. corresponding nodes extracted in previous step, were employed to refine the initial estimate [6]. With this approach, four corresponding landmarks were required (if no initial estimate is used and two projection images are at arbitrary and unknown orientation, determining the geometry requires eight or more corresponding landmark points [28], which can be unrealistic in many clinical cases). A 3D landmark point was reconstructed as the nearest point to corresponding non-intersecting projection lines [23]. Then, the obtained 3D locations were projected back onto a plane. The transformation, in the form of a rotation matrix R and a translation vector t that minimized the sum of distance errors of newly projected points from the input set of image coordinates was estimated using least-squares. Refinement of the geometry parameters was performed iteratively until convergence, such that at each iteration a new set of 3D landmark points was calculated using the newly computed transformation. For the cases in our study, we found that after 25 iterations, the accuracy improved relatively slowly with increasing number of iterations (after optimization, the maximum Root-Mean-Square (RMS) Euclidean distance between the set of image landmarks and corresponding projected reconstructed landmark points was 1.41mm, and the average RMS Euclidean distance was 0.85mm). Resulting transformation was applied to relative centrelines and profiles prior to 3D reconstruction.
The number of planes acquired in a routine ICA can occasionally vary from a minimum of two planes up to four planes, depending mostly on the lesion location and its visibility. It is important to make the method adaptable to a case in which more than two planes are acquired so to exploit all available information. If more than two projection images were available, 3D landmarks were reconstructed for all pairs of two projection images. Planes were sorted based on increasing RMS Euclidean distance between the set of image landmarks and corresponding projected 3D reconstructed landmarks, after optimization. According to the smallest RMS Euclidean distance, one plane was referred to as reference plane, and the pair of planes P n;m with the smallest RMS distance was referred to as reference pair and employed to generate a 3D centreline C. The imaging geometry of the additional planes was optimized with respect to the reference plane.
(a) 3D centreline reconstruction. A 3D vessel centreline was reconstructed as intersection of surfaces defined by corresponding arterial segments. Notice that from step 1) the segments correspondences between two projection planes P n;m were known. Given paired transformed centrelines C n,m and their respective focal points f n,m , the surfaces S n,m representing projectional beams for C n,m were formed by connecting C n,m to the respective focal point (Fig 3). The problem of finding the 3D centreline C from the projected curves C n,m was thus reduced to finding the intersection curve of the two surfaces S n,m : Precisely, S n,m were triangulated surfaces and their intersection was calculated by applying the fast mutual triangle intersection test [29]. The basis of this test is that if two triangles (elements of the surfaces S n,m ) intersect, then they overlap along the line of intersection of their respective planes. Briefly, the test computes the signed distance of the three vertices of a triangle from the plane containing the other triangle. If all the distances have the same sign, then the two triangles do not intersect. Otherwise, they may intersect, and the problem is reduced to an overlap test of two segments lying on the intersection line. The algorithm computes the scalar intervals for which the intersection line lies inside each triangle; if the intervals overlap, the intersection line segment is computed (i.e. two unique points). Robustness problems may rise when the triangles are nearly co-planar or an edge is nearly co-planar to the other triangle (especially when an edge is close to an edge of the other triangle). To handle these cases, a tolerance value was defined for the distances between the vertices of the two triangles by experimentally observing the highest 3D Euclidean distance between corresponding landmarks achieved i.e. if a point is relatively close to the other triangle's plane, it is considered as being on the plane. All the intersections (line segments or single point) for the two triangulated surfaces S n,m gave an ordered set of points that, connected, formed a 3D centreline C.
If the dataset included more than two projections, the pair of planes P n;m with the smallest RMS distance was employed to generate the reference 3D centreline C. To incorporate the information in vessel diameter from an additional image, the intersection of the reference surface with the additional surface was computed. A point (x k ,y k ,z k ) 0 on the newly computed intersection curve was associated to a point (x k ,y k ,z k ) on the reference centreline C as the closest point to the projectional line f n , (x k ,y k ,z k ) (Fig 4). Hence, the diameter vector information at (x k ,y k ,z k ) 0 was stored at the point (x k ,y k ,z k ) for next 3D luminal cross-section reconstruction. (b) 3D luminal cross-section reconstruction. In 3D reconstruction from X-ray angiograms, the accuracy of the luminal contour will depend on the number N of available projection planes, their mutual orientation and their orientation with respect to the arterial segment of interest. Due to the limited number N of projection planes, a conic function could provide an acceptable approximation for luminal contours. A visual inspection of various OCT datasets indicated that healthy arterial segments have nearly circular luminal cross-sections. However, when a build-up of plaque occurred, the circular contour deformed in various and unpredictable ways: in some cases, the plaque distribution was uniform along a cross-section, and the lumen still appeared as roughly circular, in other cases the distribution was less uniform, yielding various shapes (e.g. semi-circular, oval, triangular). Samples of OCT cross-sections from diseased arterial segments are reported in the next section. Based on such observations, design of a luminal cross-section consisted of two stages: first, a cross-section was designed using the circularity approximation to handle the limited number of projections, and then, local adjustment of the initial cross-section was carried out in accordance with 2D projected information (2D boundaries), as follows.
The computed 3D centreline C was the spine along which the lumen cross-sections were generated. Cross-sectional planes were defined in 3D space at equidistant points k along a 3D reconstructed centreline curve C. Given a 3D point (x k ,y k ,z k ) on a centreline curve C, a crosssectional plane S k was described by the tangential vectorn S k of C at (x k ,y k ,z k ) (Fig 5A). 2D image point (u k,n , v k,n ) was the projection point of (x k ,y k ,z k ) on a projection plane P n , and d k,n was the diameter vector orthogonal to the 2D vessel centreline direction at (u k,n , v k,n ). The vector d k,n was then scaled by the magnification factor, computed as the ratio of the distances d(f n , (x k ,y k ,z k ))/d(f n , (u k,n , v k,n )), with f n the focal point. Generally, the diameter vector d k,n was not  orthogonal to the 3D centreline curve C at (x k ,y k ,z k ); thus, the scaled diameter d M k;n was projected on the cross-sectional plane S k at (x k ,y k ,z k ) as ( Fig 5B) [23]: At this stage, a luminal contour was designed as a circle, with diameter given by the average length of reconstructed diameter vectors d 0 k . We used Non-uniform Rational Basis Splines (NURBS) to parametrize a luminal contour. Parameterization using NURBS allowed us to obtain a flexible and versatile modelling of a cross-section, capable of representing simple curves (e.g. circles) and more complex free-form ones. A luminal contour was represented as: where Q i were the control points, n+1 was the number of control points, p was the degree of the curve, and B i,p (u) was a rational basis function of i th control points and p th -degree, defined on the knot vector. For a comprehensive treatment of NURBS, refer to Piegl and Tiller [16]. Initially, a circular luminal contour γ k (u) was uniformly sampled along its length; the order of a curve γ k (u) was three and n+1 control points were uniformly distributed along its length (n+1 = 17 control points were found to be suitable for all our cases with 2N = 4,. . ., 8 reconstructed 3D boundary points). The shape of the curve was then adjusted as follows. From Eq 2, a set of two 3D coordinate points describing a 3D luminal cross-section at (x k ,y k ,z k ) was obtained per plane. Given N projection planes, Q k = 2N points describing a 3D luminal crosssection at (x k ,y k ,z k ) were obtained. 3D points Q k were labelled as interpolatory control points for the luminal contour γ k (u). For each Q k , the Euclidean distances from the control points of γ k (u) were computed. Each point Q k replaced its closest control point. To achieve interpolation at a point Q k , the multiplicity of its knot value was increased (more precisely, a knot value was repeatedly inserted until its multiplicity was equal to the degree of the curve) [16]. The remaining control points of γ k (u) were labelled as non-interpolatory. In order to avoid unrealistic spikes, a neighbouring control point of an interpolatory point was removed if the difference between their radial distances was at least 25% of the contour mean radii. Otherwise, γ k (u) was locally adjusted by spacing apart the knots of non-interpolatory control points. Fig 6 shows two examples of initial contour (blue curve), computed interpolatory points Q k (green squares), and resulting adjusted contour γ k (u) (green curve).
(c) 3D Surface reconstruction. A technique called lofting (or skinning) can be used to generate the surface passing through a set of luminal contours {γ k (u)}, k = 0,. . .,K [16]. Lofting was performed as an interpolation across the control points in each luminal contour, yielding the new control points P j of the lofted vessel surface (Fig 7).
In order to preserve continuity in the luminal surface, some precautions were taken prior to interpolation across contours. Control points were sorted according to their angles on the plane and in the same direction (clockwise), so that the luminal surface had no flips or twists. In addition, the set of first control points of all luminal contours formed a continuous curve; without this precaution, regions of torsion in the luminal surface would appear in correspondence of discontinuities along the curve.
Importantly, lofting required that the curves were compatible, that is they had common degree and number of control points and were defined over the same knot vector [30]. To ensure compatibility, degrees were elevated to a maximum degree, then knot vectors were merged, and finally knots were inserted so that contours had the same knot vector and, hence, the same number of control points [31]. This process did not change the shape of luminal contours.
Then, the parameters fv k g and the knot vector V were computed [16], and they were used to interpolate degree-q curves through the control points of Eq 4: so that Eq 5 interpolated Q i at certain v values. These newly obtained control points were the control points of the lofted surface (luminal vessel surface) defined over the knot vectors U and V. By applying the lofting technique, 3D luminal vessel surfaces were reconstructed as in Eq 5.

Dataset
Study population. The study population comprised patients enrolled in clinical studies at the catheterization laboratory, of the John Radcliffe University Hospital, Oxford, UK. Ethics approval (Oxford Regional Ethics Committee) was obtained and all patients provided informed written consent for all study procedures. OCT cohort was enrolled from the Plaque Imaging and Biomarker Study-Ethics Reference 08/H0603/41. Patients who were scheduled for non-emergency coronary angiography with the view to possible percutaneous coronary intervention were recruited. The angiographic images were acquired prior to inserting the intracoronary imaging catheter. OCT of the diseased coronary artery segment was performed prior to the deployment of any stent. The FFR cohort was part of the Oxford Acute Myocardial Infarction Study (OxAMI)-Ethics Reference 10/H0408/24. Patients were recruited in 2010/ 2011 for the Biomarkers Study (OCT cohort) while the OxAMI study (FFR cohort) is still ongoing (patients were selected in 2013/2014.) In this study, results are obtained on twelve patients. Six patients underwent ICA and OCT; six patients underwent ICA and FFR.
Data acquisition. X-ray arteriographic images were recorded with single-flat-panel mobile C-arm (Axiom Artis, Siemens; contrast agent Niopam, Bracco), at 15 frames/sec, with the detector at angles selected by the interventional cardiologist. Each frame represents a projection of the coronary tree onto the single-plane system. The frames were related to the cardiac cycle by electrocardiographic (ECG) gating, which enabled selection of views at the same cardiac phase. Image acquisition settings, including focal point to image plane distance, field of view, image plane orientation, image pixel spacing, were automatically stored with each image file in DICOM format.
Serial OCT images were acquired using the Ilumien system (St. Jude Medical, Minneapolis, MN, USA). Optical Coherence Tomography uses back-reflected infrared light to perform insitu micron scale tomographic imaging of the vessel anatomy and internal microstructure of plaques [32,33]. With manual injection of 20ml of contrast solution, the OCT catheter (2.7 F) was automatically pulled back at speed of 20mm/sec and rate of 100 frames per second spanning approximately 50mm (a high frame-rate is desirable for many reasons, including higher longitudinal resolution and cardiac motion-free intracoronary imaging [34,35]).
Fractional Flow Reserve is defined as the ratio between the distal pressure P d (measured close to the distal stenosis) and the proximal pressure P a (equivalent to mean aortic pressure) of the stenosis [36]. A cut-off value of 0.8 indicates ischemia. Calibrated and equalised pressure wire (Certus, St Jude Medical, St Paul, MN, USA) was advanced to cross the lesion. Hyperaemia was induced using an intravenous infusion of adenosine at a rate of 140 μg/kg/min. Measurements of absolute mean aortic pressure, P a , and mean distal pressure, P d , were recorded when reaching steady-state hyperaemia (under the assumption that during hyperaemic condition, myocardial resistances are minimal and remain constant).

Experiments
In order to investigate the accuracy of the proposed method, a series of experiments were performed on routine clinical data.
Assessment of 2D vessel detection. Accurate vessel extraction is crucial to accurate 3D vessel reconstruction. Obtained vessel boundaries were compared with expert manual tracing, in terms of reduction in diameter, R, due to a lesion (stenosis severity): where D o is the minimum occluded diameter, D p and D d are the proximal and distal reference diameters, respectively (D d ; D p is their average value). 2D expert manual tracing was performed off-line on each radiographic image using apposite 2D-QCA software (McKesson Cardiology). The profiles of a segment of interest were traced to measure the above quantities of interest (calibration was performed automatically by means of the size of the size of the guiding catheter).

Assessment of 3D reconstruction.
OCT images were analysed off-line by one author (MA) who was completely blinded to the 3D reconstruction method. Analysis was performed using specialist software Ilumien Analysis Package (St. Jude Medical). Manual measurements of lumen mean, minimum, and maximum diameters and cross-sectional areas were done on every individual frame of the lesion (calibration was performed automatically by means of the size of the OCT catheter).
The comparison of the OCT analysis and corresponding 3D reconstructed segment analysis demanded matching of the two modalities. First, corresponding segments were visually identified on OCT frames and 2D angiograms. Then, a common anatomical landmark was identified on both imaging modalities (e.g. a side branch or the ostium for the main epicardial vessels); from the landmark, a segment length of interest was established based on OCT data. Cross-sectional areas were computed along this length, for both OCT and 3D reconstructed data (with sampling rate as OCT).
Comparison with OCT derived analysis involved the proposed method and conventional approaches, namely i) circular cross-section approximation (with diameter d k,n ) and ii) ellipse fitting of the computed 3D profile points Q k . A least-squares method was employed that minimized the Euclidean distances of the computed Q k = 2N points to the ellipse [37]. Comparison was carried out by computing, for each approach, the difference in every pair of corresponding cross-sections between OCT and the approach, followed by a RMS of all the discrepancies derived from every individual pair of cross-sections.
For a 3D reconstructed segment, luminal cross-sectional areas were calculated using Stokes Theorem for computing the area of a planar polygon [38]. Parametric points {p i = (u i ,v i )}, i = 0, . . ., n on a luminal cross-sectional curve were computed following a uniform distribution, and the area A enclosed by the curve was calculated as: Comparison with FFR measurements. Progressively increasing plaque burden causes increasing severity of luminal stenosis and can lead to hemodynamically restriction of flow across a lesion [39]. ICA-derived luminal measurements might be employed for the prediction of hemodynamically significant stenosis as determined by FFR [40,41]. Correlation between the FFR measurement (i.e. derived from pressure gradient across a stenotic lesion) and the volume reduction of that lesion, as calculated with the proposed 3D reconstruction method, was investigated. [18].
For comparison, the matching segment was visually identified on radiographic images, acquired while performing the FFR procedure, as between the proximal location of the pressure sensor and its distal location. The same segment was identified on corresponding diagnostic 2D angiograms.
The reduction in volume was computed as the ratio of the real volume (V R ) to the interpolated volume (V I ), as if there were no lesion between the proximal and the distal portions of a segment of interest. As before, a set of luminal cross-sectional areas along a segment was computed by applying Eq 7. For the interpolated lumen, cross-sections were assumed circular and diameters were estimated through linear interpolation between proximal and distal diameters.

Results
The algorithm was implemented in MATLAB (The MathWorks Inc., Natick, MA, USA). 3D reconstruction and analysis were performed off-line after the completion of angiographic study, using an Intel Core 2 i7-480MQ CPU @ 2.70 GHz, with a computational time of 1-2 min. In this section, coronary artery segments for a patient are abbreviated as: LCx, Left Circumflex; LAD, Left Anterior descending; Marg, Marginal; IntM, Intermediate; RCA, Right Coronary Artery.
The results of the 2D detection are given in Fig 8. A total of 24 stenotic segments from 12 patients were included to compare reduction in diameter by the fully manual and the semiautomatic segmentation method. The diseased segments had various degrees of severity: the average reduction in diameter was of 53.2% as assessed by the manual segmentation, the minimum reduction was 28.5%, and the maximum reduction was 77.4%. There was an excellent correlation between the manual segmentation and the implemented method, with a coefficient of determination r 2 = 0.969 (Fig 8, left). The Bland-Altman plot (Fig 8, right) also depicts a high agreement between the manual and semi-automatic driven values with a mean difference close to zero and with only one value outside the 1.96 Ã SD cut-offs (SD = standard deviation, value corresponding to 55.0% reduction as assessed by manual segmentation against 49.6% as assessed by the implemented method).
A total of 6 diseased segments from 6 OCT patients (one segment per patient) were examined (LAD = 2, RCA = 1, Marg = 1, IntM = 1, LCx = 1). The number of analysed cross-sections varied between 18 and 134 per patient, this depending on the acquired OCT data (i.e. quality of the OCT run, lack of anatomical landmarks, lesion length). The computed volume of the lumen based on the OCT analysis was compared to the computed volume of the lumen based on the proposed 3D reconstruction method, for corresponding segments. The results show an excellent correlation between the two methods with a coefficient of determination r 2 = 0.986 (Fig 9,  left). Bland-Altman plot (Fig 9, right) confirmed the high agreement between the two methods with all the values being within the 1.96 Ã SD cut-offs. Corresponding quantitative data are presented in Table 1. Mean diameter, minimum diameter, maximum diameter are also reported, with a maximum distance error of 0.79mm and a mean error of 0.29mm. Overall, Patient 1 (LAD) shows the highest percentage error in volume (18.10%) and in maximum diameter (24.18%), while Patient 3 (Marg) shows the highest percentage error in minimum diameter (28.59%). For these patients, we observed a lower contrast in correspondence of the lesion segment (possible thrombus), which can result in slight underestimation of the diameter. Analysis of paired cross-sectional areas involved the proposed method, circular approximation and ellipse fitting. Two cases of diseased luminal cross-sections, acquired with OCT and as reconstructed with the three approaches, are shown in Fig 10. The reconstructed luminal cross-sections were found to be a closer representation of the true OCT contours both visually and quantitatively. Maximum and minimum discrepancies, standard deviation and RMS discrepancy of cross-sectional areas are reported in Table 2. The proposed method showed the lowest RMS discrepancy, with overall values ranging from 0.213mm 2 to 1.013mm 2 , and maximum error of 1.837mm 2 . Maximum errors with circular and elliptical fitting were, respectively, 2.322mm 2 and 5.078mm 2 .
A total of 6 diseased segments from 6 FFR patients were analysed (LAD = 3, RCA = 2, Marg = 0, IntM = 0, LCx = 1). Fig 11 illustrates two 3D reconstructed lesions, their  interpolated volume and proximal and distal pressure measurement locations. The measured FFR values were compared to segment volume reduction, V R /V I , calculated based on the 3D reconstructed lesion. The results are illustrated in Fig 12. An FFR measurement of 1.0 indicates normal blood flow, an FFR measurement lower than 0.8 is the clinically validated cut-off that is used in clinical practice to indicate a functionally significant lesion that may cause ischemia. FFR measurements varied between 0.64 and 0.89, and 3/6 lesions were evaluated as physiologically significant. The results showed a good correlation between FFR and estimated volume reduction, V R /V I , for 5/6 patients. However, volume reduction for Patient 6 seemed to indicate a not severely obstructive stenosis, while the FFR measurement equal to 0.64 indicates a clinically significant stenosis.

Discussion and conclusions
This paper proposes a novel method to generate in an effective manner a 3D computational patient-specific model of lesioned vessel, directly from 2D projections acquired while performing invasive coronary X-ray angiography, and appropriate for immediate geometric and quantitative analysis. The technical novelty of our work resides in the 3D reconstruction strategy and the cross-section NURBS parameterization. The proposed method i) computes a 3D centreline as intersection of surfaces defined by corresponding branches, overcoming the issues related to point-to-point matching based reconstruction; ii) it parametrizes a 3D luminal cross-sectional contour by using NURBS, and adjusts a model to interpolate computed contour points of the luminal surface, thus avoiding the tubular approximation; iii) it generates a 3D model by lofting of the obtained cross-sections, thus it allows us to generate any complex shape from any set of cross-sections. Importantly, the method is not restricted to a specific number of 2D projections and additional information can be integrated sequentially. The design of a 3D cross-section consisted of two main steps: first, a template contour was built on a cross-sectional plane with the circularity approximation and, second, an alteration of the contour was performed accordingly to the computed 3D luminal contour points. In the first step, the use of 3D cross-sectional planes circumvents the problems with degenerated ellipses and non-parallel neighbouring cross-sections [14]. Importantly, the use of a template contour handles the limited number of projections in ICA, often as low as two, i.e. lack of information on the vessel boundary. The choice of a circular template contour derived from the observation on OCT data that healthy cross-sections are roughly circular. When computing the diameter vectors, i.e. 3D lumen boundary points, we accounted for the fact that the projections are in general not perpendicular to the normal of the cross sectional planes of the centreline [5,8]. In the second step, parameterization of the luminal contour by means of NURBS functions permitted a smooth representation of the contour, and control over changes on its form. Initial control points of the template circular contour were adjusted or deleted so that the contour curve interpolated the computed points representing 3D lumen boundary points. The influence on the shape of the initial control points depends on their location with respect to the interpolatory control points. Generally, for a higher number of acquired projections, i.e. more information on the luminal surface is provided, their influence will decrease, while for two projections only, their influence will be higher.
Validation of the 3D reconstruction method involved comparison with independent clinical measurements that are used routinely in the clinical assessment of coronary artery stenosis severity. To the best of our knowledge, no previously proposed method was validated using OCT data. Importantly, clinical data included stenoses of various degree of severity and comparison was performed on a frame-to-frame basis. This comparison is of great clinical value, as it suggests that 3D renderings of coronary arteries from ICA could help assessment of coronary artery stenosis severity, in combination or without the need for expensive and potentially hazardous invasive techniques.
The measured lumen volumes based on the OCT analysis were similar to the computed lumen volumes based on the 3D reconstruction method, indicating a very good performance of the algorithm. For the six patients, the maximum difference in minimum diameter and maximum diameter values were about 29% and 24%, respectively. For these two cases, we observed a lower image contrast in correspondence of the lesion segment, possibly due to a thrombus, which could result in an underestimation of the diameter. Correlation of volumes showed very good correspondence between our proposed method and OCT; however, the relative errors in diameter between our method and OCT are likely to be underestimated as they are averaged along the length of the segment.
For a thorough comparison, paired luminal cross-sectional areas were analysed. When compared with circular approximation and elliptical fitting, the proposed method showed a closer representation of OCT luminal shapes and a higher accuracy in computed cross-sectional areas. The highest improvement of the proposed method over the circular and elliptical approaches was, respectively, about 40% (Patient 6), and about 70% (Patient 3). It was observed that for small angles between reconstructed diameter vectors, the ellipse fitting approach tends to underestimate cross-sectional areas and can result in a degenerate ellipse. It is important to underline that the discrepancy between the cross-sectional areas generally depends on the number of projections employed in the reconstruction, as well as on their mutual orientation and their orientation with respect to the arterial segment of interest. Fig 13 shows the RSM discrepancy between OCT and our approach (Patient 5), with varying number of employed angiographic planes in the reconstruction. The lowest RSM error of 0.213mm 2 was observed when all the four planes were employed in the 3D reconstruction; this error increased above 0.4mm 2 when subsets of three and, then, two planes only were employed. Similarly, we observed that the RSME increased to 2.534mm 2 for Patient 1, when two out of the three acquired planes were employed.
We also explored whether the proposed method may reflect an FFR measurement, a widely used physiological estimate of lesion severity. This is of potential clinical relevance; however, to extend this to a clinical validation would require a clinical study powered adequately and specifically designed for that purpose. We suggest that the level of validation provided in the current manuscript provides a rationale and justification for such a study. The comparison of  Table 2.
https://doi.org/10.1371/journal.pone.0190650.g013 volume reduction due to a plaque with the FFR measurements showed a good agreement in establishing the clinical significance of a stenosis, except for one patient. This was the patient with the lowest FFR value, who showed a highly irregular luminal surface. Such disagreement might indicate that the plaque morphology has implication on the fluid dynamics and that for the same degree of stenosis (volume reduction due to a lesion) hemodynamic parameters might dependent on the plaque shape, as claimed by recent works [42]. Overall, these results seems to be in line with a recent clinical study that investigates the relationship between 2D quantities measured on ICA images and FFR measurement, showing that the 2D minimum diameter correlates well with FFR value in intermediate coronary lesions [43].

Future work
The 3D model obtained in this paper is suitable for isogeometric analysis (IGA) of blood flow in the coronary arteries [44]. Also, the model intends to preserve the lumen morphology of a lesion, and recent work has shown the implications of it on the flow dynamics [42]. In the future, computational fluid dynamics approaches may supplement the 3D anatomical information obtained by ICA data by providing accurate estimation of FFR [44][45][46].
Moreover, NURBS properties (e.g. local support property) make the 3D model suitable for a co-registration problem. For instance, co-registration of the 3D coronary artery tree with corresponding 3D volume rendering of the epicardial surface (extracted from cardiac magnetic resonance images) may help relating location and status of a stenosis with area of myocardium at-risk [47].