Method and apparatus for three-dimensional reconstruction of coronary vessels from angiographic images and analytical techniques applied thereto

[0002] This patent application is a continuation in part of U.S. Pat. No. 6,047,080 issued Apr. 4, 2000 and filed Jun. 16, 1996 as application Ser. No. 08/665,836.

[0003] A portion of the disclosure of this patent document contains material which is subject to copyright protection. The copyright owner has no objection to the facsimile reproduction by anyone of the patent document or the patent disclosure, as it appears in the Patent and Trademark Office patent file or records, but otherwise reserves all copyright rights whatsoever.

[0004] 1. Field of the Invention

[0005] The present invention relates generally to a method for reconstructing images of vascular structures and more specifically to an improved method for three-dimensional (3-D) reconstruction of vascular structures from two two-dimensional biplane projection images and methods and structures for quantitative analysis of such a reconstructed structure.

[0006] 2. Discussion of Related Art

[0007] Several investigators have reported computer assisted methods for estimation of the 3-D coronary arteries from biplane projection data. These known methods are based on the known or standard X-ray geometry of the projections, placement of landmarks, known vessel shape, and on iterative identification of matching structures in two or more views. Such methods are described in a publication entitled “3-D digital subtraction angiography”, IEEE Trans. Med. Imag.,vol. MI-1, pp. 152-158, 1982 by H. C. Kim, B. G. Min, T. S. Lee, et. al. and in a publication entitled “Methods for evaluating cardiac wall motion in 3-D using bifurcation points of the coronary arterial tree”, Invest. Radiology,January-February pp. 47-56, 1983 by M. J. Potel, J. M. Rubin, and S. A. Mackay, et al. Because the computation was designed for predefined views only, it is not suitable to solve the reconstruction problem on the basis of two projection images acquired at arbitrary and unknown relative orientations.

[0008] Another known method is based on motion and multiple views acquired in a single-plane imaging system. Such a method is described in a publication entitled “Reconstructing the 3-d medial axes of coronary arteries in single-view cineangiograms”, IEEE Trans. MI, vol.13, no. 1, pp. 48-60, 1994 by T. V. Nguyen and J. Sklansky uses motion transformations of the heart model. However, the motion transformations of the heart model consist only of rotation and scaling. By incorporation of the center-referenced method, initial depth coordinates, and center coordinates, a 3-D skeleton of the coronary arteries was obtained. However, the real heart motion during the contraction involves five specific movements: translation, rotation, wringing, accordion-like motion, and movement toward the center of the ventricular chamber. Therefore, the model employed is not general enough to portray the true motion of the heart, especially toward the end-systole.

[0009] Knowledge-based or rule-based systems have been proposed for 3-D reconstruction of coronary arteries by use of a vascular network model. One such knowledge-based system is described in a publication entitled “An expert system for the labeling and 3-D reconstruction of the coronary arteries from two projections”, International Journal of Imaging,Vol. 5, No. 2-3, pp. 145-154, 1990 by Smets, Vandewerf, Suctens, and Oosterlinck. Because the rules or knowledge base were organized for certain specific conditions, it does not generalize the 3-D reconstruction process to arbitrary projection data. In other knowledge-based systems, the 3-D coronary arteries were reconstructed from a set of X-ray perspective projections by use of an algorithm from computed tomography. Due to the motion of the heart and only a limited number of projections (four or six), accurate reconstruction and quantitative measurement are not easily achieved.

[0010] Closed-form solutions of the 3-D reconstruction problem using a linear approach was a significant development and is described in, for example, a publication entitled “Determining 3-d motion and structure of a rigid body using the spherical projection”, CVGIP,vol. 21, pp. 21-32, 1983 by B. L. Yen and T. S. Huang. Unfortunately, actual data is always corrupted by noise or errors and the linear approach based techniques may not be sufficiently accurate when using noisy data. Hence, optimal estimation has been explicitly investigated. Additionally, U.S. Pat. No. 4,875,165 entitled Method for Determination of 3-D Structures in Biplane Angiography issued in the name of Fencil et al. also has significant drawbacks.

[0011] Use of a two-step method is known for producing an optimal estimation for a 3-D structure based on maximum-likelihood and minimum-variance estimation. In these techniques, for example, two publications entitled “Optimal motion and structure estimation”, IEEE Trans. on PAMI, Vol. 15, no. 9, September 1993, pp. 864-884, and “Structure from motion using the reconstruction and projection technique”, Proc. IEEE Workshop Computer Vision, November 1987, pp. 345-348, image error was employed in the objective function for a non-constricted minimization process. Preliminary estimates computed by a linear algorithm were used as initial estimates for the process of optimal estimation. However, if the initial solution from the linear approach is not sufficient, (e.g., with more than 2 pixels=0.6 mm error in the input 2-D image data), the iterative minimization process at the second step may become trapped in a local minimum due to a lack of prior information concerning the variables to be optimized.

[0012] Quantitative coronary analysis (“QCA”) of a reconstruction of an arterial tree was known and developed in the 1970's to quantify vessel geometry and the effects of drugs on the regression and progression of coronary artery disease. In the mid-1980's, digital systems were introduced into the catheterization laboratory to support the angiographer during the interventional procedures. With the advent digital angiographic technology, on-line QCA has been widely used predominantly for the selection of the optimal sizes of the interventional devices, and for the assessment of the efficacy of the procedure. However, current QCA techniques are performed on the 2-D projection views in which the vessel overlap and foreshortening are subjectively minimized in a “trial-and-error” manner by the interventionist. FIGS. 56(a)-(d) demonstrate the problems resulting from use of such 2-D projections in QCA techniques. FIG. 56cis an RCA image as known in the art used for QCA and resulting in a 30% narrowing measurement as shown in FIG. 56d.With 2-D projection views there is no way to know or estimate how much error occurs in the QCA process due to foreshortening with respect to the stentotic segment. In FIGS. 56aand 56b,the same vessel segment between two bifurcation points as marked by two dots at its proximal and distal ends depicts 77% and 52% foreshortening, respectively.

[0013] It is evident from the above discussion that a need exists for improved reconstruction of 3-D images from 2-D image data and that a further need exists for improved QCA techniques utilizing such 3-D reconstruction to provide needed analysis in the intervention process.

[0014] The disadvantages of present methods known for visual reconstruction of vascular structures are substantially overcome with the present invention by providing a novel method for three-dimensional reconstruction of vessels using two-dimensional angiograms and further by providing improved QCA techniques utilizing such 3-D reconstructions of vascular structures.

[0015] This patent application teaches the methods and structures of the invention with reference to vascular structures (i.e., coronary arteries or coronary arterial trees). Those skilled in the art will readily recognize that the improved reconstruction methods and the quantitative analysis methods and structures are equally applicable to any vascular structure. References herein to “coronary” or “arterial” applications of the invention are intended merely as exemplary of common vascular structures where such reconstruction and analysis techniques are beneficially applied. Nothing in this patent application should be construed as limiting the reconstruction and analysis features to coronary applications. In particular, reference to “quantitative coronary analysis” (or “QCA”) should be understood to refer broadly to quantitative analysis of any vascular structures. Likewise, reference to an “artery” or “arteries” should be broadly understood to refer to any vessel or vessels. “Artery tree”, “arterial tree”, “coronary arterial tree” or “coronary artery tree” should all be understood to broadly refer to any vascular tree structure. Coronary applications of the present invention are therefore merely one common exemplary application of the features, methods and structures of the present invention.

[0016] In the present inventive method, a novel optimization technique minimizes the image point errors and directional vector errors in both imaging systems subject to the constraints derived from the individual intrinsic imaging parameters of the employed imaging system. Given five or more corresponding object points in both views (projection image), a constrained nonlinear optimization algorithm is applied to obtain an optimal estimate (transformation) of the biplane imaging geometry in the form of R and {right arrow over (t)} which characterize the position and orientation of one imaging system (imaging portion of the imaging system) relative to the other. The initial solution is estimated on the basis of the individual intrinsic imaging parameters.

[0017] Angiograms of eight hundred patients were analyzed in which two cases are selected for discussion hereinafter. The biplane imaging geometry was first determined without a calibration object, and the 3-D coronary arterial trees were reconstructed, including both left and right coronary artery systems. Various two-dimensional (2-D) projection images of the reconstructed 3-D coronary arterial tree were generated and compared to other viewing angles obtained in the actual patient study. Similarity between the real and reconstructed arterial structures was excellent. The accuracy of this method was evaluated by using a computer-simulated coronary arterial tree. Root-mean-square (RMS) errors in the 3-D position and the 3-D configuration of vessel centerlines and in the angles defining the R matrix and {right arrow over (t)} vector were 0.9-5.5 mm, 0.7-1.0 mm, and less than 1.5 and 2.0 degrees, respectively, when using 2-D vessel centerlines with RMS normally distributed errors varying from 0.4-4.2 pixels (0.25-1.26 mm).

[0018] More specifically, the method for three-dimensional reconstruction of a target object from two-dimensional images involves a target object having a plurality of branch-like vessels. The method includes the steps of: a) placing the target object in a position to be scanned by an imaging system, the imaging system having a plurality of imaging portions; b) acquiring a plurality of projection images of the target object, each imaging portion providing a projection image of the target object, each imaging portion disposed at an unknown orientation relative to the other imaging portions; c) identifying two-dimensional vessel centerlines for a predetermined number of the vessels in each of the projection images; d) creating a vessel hierarchy data structure for each projection image, the data structure including the identified two-dimensional vessel centerlines defined by a plurality of data points in the vessel hierarchy data structure; e) calculating a predetermined number of bifurcation points for each projection image by traversing the corresponding vessel hierarchy data structure, the bifurcation points defined by intersections of the two-dimensional vessel centerlines; f) determining a transformation in the form of a rotation matrix and a translation vector utilizing the bifurcation points corresponding to each of the projections images, the rotation matrix, and the translation vector representing imaging parameters corresponding to the relative orientations of the imaging portions of the imaging system; g) utilizing the data points and the transformation to establish a correspondence between the two-dimensional vessel centerlines corresponding to each of the projection images such that each data point corresponding to one projection image is linked to a data point corresponding to the other projection images, the linked data points representing an identical location in the vessel of the target object after the projections; h) calculating three-dimensional vessel centerlines utilizing the two-dimensional vessel centerlines and the correspondence between the data points of the two-dimensional vessel centerlines; and i) reconstructing a three-dimensional visual representation of the target object based on the three-dimensional vessel centerlines and diameters of each vessel estimated along the three-dimensional centerline of each vessel; and j) determining the optimal view of the vessel segments with minimal vessel foreshortening.

[0019] Further features of the present invention utilize such 3-D reconstructions to provide improved QCA techniques for analysis of coronary artery images.

[0020] The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee.

[0021] FIG. 1 is a flowchart of a high level view of the reconstruction and analysis procedures of the present invention.

[0022] FIG. 2 is a flowchart providing additional details of the collection of gantry system parameters and image acquisition.

[0023] FIG. 3 is a diagram of a typical bi-plane imaging system.

[0024] FIG. 4 is a diagram of a typical single plane imaging system.

[0025] FIG. 5 is a diagram of intrinsic imaging parameters of a typical single plane imaging system with 7″ by 7″ field-of-view and 94 cm SID characterized by mathematical model.

[0026] FIG. 6 depicts a mathematical model of two projections, showing a 3-D object point and its projection in two views.

[0027] FIG. 7 is exemplary of estimated pixel size on 7″ field-of-view.

[0028] FIG. 8 is exemplary of the determination of the image centroid by calculation of the intersection of two lines.

[0029] FIG. 9 is a flowchart providing additional details of the identification of 2-D coronary arterial tree and feature extraction.

[0030] FIG. 10 depicts a manual indication of initial vessel centerline by selection of a sequence of points.

[0031] FIG. 11 depicts a manual indication of pseudo-branches by selection of a sequence of points originating from a desired landmark.

[0032] FIG. 12 depicts resultant ridge-point map where each dot corresponds to a local minima grey level.

[0033] FIG. 13 depicts final results of deformation process after an initial curve is aligned with the actual vessel centerlines.

[0034] FIG. 14 depicts a bifurcation point calculated by finding the intersection between a side branch and the precedent vessel.

[0035] FIGS. 15a-cdepict a vessel hierarchical digraph constructed from 2-D coronary tree identified in two images.

[0036] FIG. 16 is a flowchart providing additional details of the determination of the transformation defining spatial relationship of the employed two views.

[0037] FIG. 17 depicts the actual and initially calculated bi-plane geometry.

[0038] FIG. 18 depicts the calculation of bi-plane transformation based on bifurcation points and directional vectors.

[0039] FIG. 19 depicts cone-shaped bounding regions used to constrain the calculated 3-D points.

[0040] FIG. 20 is a flowchart providing additional details of the establishment of vessel centerline correspondence for each pair of arteries in the coronary trees.

[0041] FIG. 21 depicts the selection of corresponding pairs of arteries using vessel hierarchical digraphs.

[0042] FIG. 22 depicts epipolar plane geometry constraints.

[0043] FIG. 23 depicts epipolar line geometry corresponding to FIG. 22.

[0044] FIGS. 24a-bdepict typical examples of finding corresponding points using epipolar lines.

[0045] FIG. 25 depicts a refinement process for establishment of vessel centerline correspondences.

[0046] FIGS. 26a-ddepict comparison of previous methods for establishing correspondence with the methods of the present invention.

[0047] FIG. 27 is a flowchart providing additional details of calculation of 3-D coronary arterial skeleton in terms of 3-D vessel centerlines and cross-sectional diameters.

[0048] FIG. 28 is a flowchart providing additional details of rendering the reconstructed 3-D coronary arterial tree with gantry angulation.

[0049] FIG. 29 is an exemplary skeleton image of 3-D coronary arteries.

[0050] FIG. 30 is an exemplary surface based modeling image of coronary arteries.

[0051] FIGS. 31a-fdepict 3-D reconstruction of RCA images viewing at four different angles of sample acquired images.

[0052] FIGS. 32a-fdepict 3-D reconstruction of LCA images viewing at four different angles of sample acquired images.

[0053] FIG. 33 is a typical example of the vessel skeleton structure in terms of 3-D curves.

[0054] FIGS. 34a-ddepict schematic diagrams of curvature and torsion.

[0055] FIGS. 35a-bdepict schematic diagrams of 3-D and 2-D take-off angles.

[0056] FIG. 36 depicts a schematic diagram of insertion angle of a bypass graft.

[0057] FIG. 37 is a flowchart providing additional details of QCA processing.

[0058] FIGS. 38a-bdepict a selected region of interest on a bypass graft and its QCA result.

[0059] FIGS. 39a-bdepict a selected region of interest on a bypass graft and its regional QCA result.

[0060] FIGS. 40a-cdepict 3-D and 2-D take-off angles of LCA and an example of insertion angle of a bypass graft.

[0061] FIGS. 41a-cdepict the percent foreshortening.

[0062] FIGS. 42a-bdepict typical vessel foreshortening on LCA segments.

[0063] FIGS. 43a-ddepict accumulated percent foreshortening on the entire LCA reconstructed from an image pair.

[0064] FIGS. 44a-ddepict a selected single segment of a region of interest and its maps.

[0065] FIGS. 45a-ddepict a selected bifurcation region of a region of interest and its maps.

[0066] FIGS. 46a-ddepict a selected bypass graft segment of a region of interest and its three regional maps.

[0067] FIG. 47 is a flowchart providing additional detail of the processing to assess length of a selected arterial segment.

[0068] FIG. 48 is a flowchart providing additional detail of the processing to assess tortuosity of a selected arterial segment.

[0069] FIG. 49 is a flowchart providing additional detail of the processing to assess take-off angles of a selected arterial segment.

[0070] FIG. 50 is a flowchart providing additional detail of the processing to assess insertion angles of a selected arterial segment.

[0071] FIG. 51 is a flowchart providing additional detail of the processing to verify connectivity of a selected arterial segment.

[0072] FIG. 52 is a flowchart providing additional detail of the processing to optimal view of a selected arterial segment to minimize foreshortening.

[0073] FIG. 53 is a flowchart providing additional detail of the processing to determine accumulated percent foreshortening of a selected arterial segment.

[0074] FIG. 54 is a flowchart providing additional detail of the processing to optimal view of a selected arterial segment to minimize overlap.

[0075] FIG. 55 is a flowchart providing additional detail of the processing to optimal view of a selected arterial segment to minimize both foreshortening and overlap.

[0076] FIGS. 56a-ddepict exemplary vessel foreshortening and existing 2-D QCA process.

[0077] While the invention is susceptible to various modifications and alternative forms, a specific embodiment thereof has been shown by way of example in the drawings and will herein be described in detail. It should be understood, however, that it is not intended to limit the invention to the particular form disclosed, but on the contrary, the invention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention as defined by the appended claims.

[0078] Overview

[0079] A mono-plane or bi-plane angiographic system can be employed for cardiac angiographic acquisition. The pair of angiograms are acquired based on any two selected gantry orientation in terms of left/right anterior oblique (RAO/LAO) angle, caudal/cranial (CAUD/CRAN) angle. The angulation is defined with respect to the iso-center to which the gantry arm rotates. The gantry system information including focal spot to image intensifier distance (SID), field of view (FOV, e.g., 5, 7, or 9 inch mode), and gantry orientation as in FIG. 1, and the acquired images are automatically recorded. The employed 3-D reconstruction protocol involves the following steps:

[0080] 1. Routine cardiac catheterization is initiaeeed and a standard coronary artenogram is completed in two views; one injection in the biplane system, two injections in the single plane system.

[0081] 2. A pair of image sequences, consisting of 8 images (including from end-diastole to end-systole frames) at each view, are chosen and transferred to a computer workstation.

[0082] 3. The 2-D coronary arterial trees are identified in the selected pair of angiograms (e.g., images corresponding to end-diastole) and the vessel features are extracted including vessel centerlines, bifurcation points, and hierarchy.

[0083] 4. The processing of 3-D reconstruction is completed by use of the information of identified vessel features yielding a 3-D patient-specific coronary arterial tree model.

[0084] 5. The computer-simulated projection from reconstructed 3-D coronary model along with the corresponding gantry angles are displayed on a computer graphics workstation.

[0085] 6. The process of 3-D quantitative coronary analysis (QCA) process including optimal-view strategy map (OSM) and assessment of coronary morphology relative to any selected arterial segment(s) is completed. The results are used as guidance for minimizing the vessel foreshortening and overlap, and for facilitating selection of other intra-coronary devices, i.e., stents.

[0086] The proposed method of reconstruction consists of seven major steps shown in the flowchart figures and descriptive text to follow. The major steps are:

[0087] (1) acquisition of two standard angiogram sequences by use of a single-plane or biplane imaging system,

[0088] (2) identification of 2-D coronary arterial trees and feature extractions including bifurcation points, vessel centerlines, and construction of vessel hierarchy in the two images,

[0089] (3) determination of transformation in terms of a rotation matrix R and a translation vector {right arrow over (t)} based on the extracted vessel bifurcation points,

[0090] (4) establishment of vessel centerline correspondence,

[0091] (5) calculation of the skeleton of 3-D coronary arterial trees,

[0092] (6) rendering of 3-D coronary arterial tree with associated gantry angulation, and

[0093] (7) calculation of optimal view(s) and 3-D QCA associated with the selected arterial segment(s) of interest.

[0094] A. Image Acquisition

[0095] Details of the methods for acquiring images of the arterial tree are shown in the flowchart of FIG. 2.

[0096] A pair of angiographic views are acquired by a typical single-plane or biplane imaging system as shown in FIG. 3 and FIG. 4 with intrinsic imaging parameters as illustrated in FIG. 5. The spatial relationship of the acquired two views can be characterized by use of a mathematical model as shown in FIG. 6. The imaging system parameters including pixel size and the imaging center coordinates are calculated as shown in FIG. 7 and FIG. 8, respectively. Once the information are calculated, they will be utilized in al the subsequent 3-D reconstruction process. After routine cardiac catheterization is initiated, a standard coronary arteriogram is completed in two standard views based on one injections in a bi-plane imaging system. With the acquired two angiographic sequences, a pair of images were chosen corresponding to the same or close acquisition time instance (e.g., end-diastole) in the respective cardiac cycles. They were transferred to the computer workstation through the fiber optic link for analysis and reconstruction. The transfer rate is 2 seconds per image with the size of 512×512 pixels and one byte storage per pixel.

[0097] B. Segmentation and Feature Extraction

[0098] Details of the methods for defining segments and extracting features from the arterial tree are shown in the flowchart of FIG. 19.

[0099] A semi-automatic system based on the technique of deformation model and segmentation technique was employed for identification of the 2-D coronary arterial tree in the given angiograms. The required user interaction in the identification process involves only the indication of several points inside the lumen of each artery, near the projection of vessel centerline in the angiogram as shown in FIG. 10. The artery with stenotic lesion or tortuous shape can be selected as landmarks to improve the reconstruction results as shown in FIG. 11. Afterwards, a spline-curve is formed based on the selected points, which serves as the initial centerline of the vessel. Meanwhile, an operator was automatically applied to the image to identify those pixels which are local minimum on intensity (i.e., the intensity of coronary arterial structures is darker than that of background). These identified pixels indicate the location of vessel centerline as shown in FIG. 12. By use of the deformation process, the identified pixels serve as the external forces to act on the initial model curve such that it will be gradually deformed and finally resides on the actual centerline of vessel as shown in FIG. 13.

[0100] The identified centerlines and the branching relationships are used for construction of the vessel hierarchy in each angiogram by a labeling process according to the anatomy and connectivity among the primary and secondary coronary arteries as shown in FIG. 14. The labeling process on the coronary tree is performed automatically by application of the algorithm of the breadth-first search to traverse identified vessel centerlines and associated bifurcations. From each vessel of the coronary tree that is currently visited, this approach tries to search as broadly as possible by next visiting all of the vessel centerlines that are adjacent to it, and finally results in a vessel hierarchically directed graph (digraph) containing a set of nodes with corresponding depth levels and two types of arcs (descendant and sibling arcs) defining the coronary anatomy. In the hierarchical digraph, each node represents a vessel in the coronary arterial tree, and the descendent arc connects the nodes between the vessel and its branch. If a vessel has multiple branches, the nodes associated with these branches are connected as a linked list by sibling arcs according to the branching locations relative to the proceeding vessel.

[0101] With the constructed hierarchical digraph, the remaining vessel features including bifurcation points and diameters as shown in FIGS. 15a-care then calculated for all of the coronary arterial segments, and the results are saved into the nodes associated with the branching vessels in the hierarchical digraph.

[0102] The deformation model is controlled under the influence or energy E_defof the internal forces characterized by smoothness constraint and external forces derived from image features. $E_{\mathrm{def}}=\sum _{i=1}^n\{\sum _{j=1}^{m_1}{E}_{\mathrm{int}}(l_i\left(v\left(s_j\right)\right)+E_{\mathrm{ext}}\left(l_i\left(v\left(s_j\right)\right)\right\}$

[0103] where n denotes the number of coronary arteries in the coronary arterial tree, and m_i, i=1, . . . ,n denotes the number of centerline points v(s_j)=(x(s_j), y(s_j)) on each artery l_i(v),The internal curve energy can be written as

[0104] E_int(l_i(v(s_j))=(α(s_j)|v_s(s_j)|²+β(s_j)|v_ss(s_j)|²)/2,

[0105] where α(s) and the first order term v_s(s) make the model curve like a membrane and β(s) and the second order term v_ss(s) make the it act like a thin plate. The external energy is derived from the ridge points of image I denoting the minimal (or maximal) gray level in a local window (e.g., k pixels by k pixels) region as

[0106] E_ext(l_i(v(s_j))=∇_k*(G_k(v(s_j))⊕I))

[0107] where ⊕ denotes the convolution operator, G_k(v(s)) is a Gaussian based filter centered at v(s) with k by k pixels of window, ∇_kdenotes the ridge point filter defined within a k by k pixels of window.

[0108] C. Calculation of Transformation Characterizing the Two Views

[0109] Details of the methods for transforming the views of the arterial tree are shown in the flowchart of FIG. 16.

[0110] By use of a mono-plane or biplane imaging system for image acquisitions, the locations and dimensions of imaging plane and x-ray source can be defined by use of a mathematical model as shown in FIG. 16. The spatial relationship between the employed two views can be characterized by a transformation in the forms of a rotation matrix R and a translation vector {right arrow over (t)} on the basis of the respective locations of x-ray source (or focal spot) served as the origin in 3-D coordinate space.

[0111] In the first view, let (u_i, v_i) denote the image coordinates of the ith object point, located at position (x_i, y_i, z_i). We have u_i=Dx_i/z_i, v_i=Dy_i/z_i, where D is the perpendicular distance between the x-ray focal spot and the image plane when the view was acquired. Let (ξ_i, η_i) denote the scaled image coordinates, defined as ξ_i=u_i/D=x_i/z_i, η_i=v_i/D=y_i/z_i. The second projection view of the employed bi-plane imaging system can be described in terms of a second pair of image and object coordinate systems u′v′ and x′y′z′ defined in an analogous manner. Scaled image coordinates (ξ_i, η′_i′) in the second view for the ith object point at position (x_i′, y_i′, z_i′) are given by ξ_i′=u_i′/D′=x_i′/z_i′, η_i′=v_i′/D′=y_i′/z_i′. The spatial relationship between the two views can be characterized by $\left[\begin{array}{c}{x}_i^{\prime }\\ y_i^{\prime }\\ z_i^{\prime }\end{array}\right]=R·\left\{\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right]-\overrightarrow{t}\right\}=\lfloor \begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\rfloor ·\left[\begin{array}{c}{x}_i-t_x\\ y_i-t_y\\ z_i-t_z\end{array}\right].$

[0112] In our previously proposed technique, the transformation was calculated based on the identified bifurcation points of the 2-D coronary arterial trees in the two views and can be only effective to the images acquired from the biplane systems or a single-plane imaging system without table panning or movement during the acquisition.

[0113] In this invention, a new 3-D reconstruction technique is developed by incorporation of the corresponding bifurcation points and the vessel directional vectors of bifurcations. The proposed method can also be employed to obtain accurate results of 3-D reconstruction without the need of ECG-triggered acquisition even when there is table panning or moving during the injection. The required prior information (i.e., the intrinsic parameters of single-plane imaging system) for 3-D reconstruction include: (1) the distance between each focal spot and its image plane, SID (focal-spot to imaging-plane distance), (2) the field of view in terms of the pixel size, p_size(e.g., 0.3 mm/pixel), (3) the distance {overscore (ff′)} between the two focal spots or the known 3-D distance between two points in the projection images, and (4) iso-center distance (with respect to which the rotary motion of the gantry arm rotates) relative to the focal spot. Items (1) and (2) are prerequisite and (3) and (4) are optional.

[0114] In the single-plane imaging system utilized, the gantry angulation associated with different view is always recorded and accompany with the acquired image file. However, the transformation directly calculated from the recorded gantry angles may not accurately define the spatial relationship between the two views due to (i) images non-simultaneously acquired from two different cardiac cycles and (ii) table panning or movement during the acquisitions. Although such a directly calculated transformation can not truly characterize the two views, the resultant estimate is generally good enough to serve as an initial guess to the non-linear minimization process for obtaining an optimal solution. In the biplane imaging system utilized, the gantry angulations associated with different views may or may not be recorded with the acquired image file. If the gantry angulations of the employed two views are unknown, orthogonal angles will be used as the initial solutions of transformation as shown in FIG. 17 by use of the equations herein below.

[0115] If the gantry angulations are known in advance, however, the transformation directly calculated from the recorded gantry angles may not accurately define the spatial relationship between the two views due to (i) patient's respiration, (ii) non-synchronous x-ray source firing on the x-ray tubes, and (iii) mis-alignment of the respective iso-centers corresponding to the individual x-ray gantry instruments of the bi-plane system. Although such a directly calculated transformation can not truly characterize the two views, the resultant estimate is generally good enough to serve as an initial guess to the non-linear minimization process for obtaining an optimal solution. With the initial solution calculated from the pre-defined orthogonal angels or recorded gantry angles, the global minimum resulting from the non-linear optimization can be guaranteed since the initial estimate is close to the actual solution.

[0116] On the basis of the initial solution, the final accurate transformation is calculated from an objective function by minimization of (1) image point errors in term of the square of Euclidean distance between the 2-D input data and the projection of calculated 3-D data points and (2) directional vector errors between the 2-D vessel directional vectors and the projection of calculated 3-D vessel directional vectors based on equations below and as shown in FIG. 18 subject the constraints shown below and in FIG. 19.

[0117] The initial estimates of the z and z′ axes of the two imaging systems are taken to be orthogonal in this situation, and the unit translation vector are assumed to lie on the x-z plane. Two different initial solutions ({overscore (R₁)}, {right arrow over (t)}_{{overscore (u1)}}) and ({overscore (R₂)}, {right arrow over (t)}_{{overscore (u2)}}) are employed as follows: ${\stackrel{\_}{R}}_1=\left[\begin{array}{ccc}0& 0& -1\\ 0& 1& 0\\ 1& 0& 0\end{array}\right],{\overrightarrow{t}}_{{\stackrel{\_}{u}}_1}=\lfloor \begin{array}{c}-D^{\prime }\\ {\alpha }^{\prime }·t_d\\ 0\\ D\\ \alpha ·t_d\end{array}\rfloor ,\mathrm{and}$ ${\stackrel{\_}{R}}_2=\left[\begin{array}{ccc}0& 0& 1\\ 0& 1& 0\\ -1& 0& 0\end{array}\right],{\overrightarrow{t}}_{{\stackrel{\_}{u}}_2}=\left[\begin{array}{c}{D}^{\prime }\\ {\alpha }^{\prime }·t_d\\ 0\\ D\\ \alpha ·t_d\end{array}\right],$

[0118] where D and D′ denote the SID distances, α and α′ are the MF factors, and t_drepresents the magnitude of {right arrow over (t)}. If the magnitude of {right arrow over (t)} is not available, an approximate measure is taken as: $t_d=\sqrt{{\left(\frac{D}{\alpha }\right)}^2+{\left(\frac{D^{\prime }}{{\alpha }^{\prime }}\right)}^2}.$

[0119] Let α₁,β₁) and (α₂,β₂) denote the gantry orientations of the 1st and 2nd views, where α_xand β_xare the LAO/RAO, and CRAN/CAUD angles, respectively. Let (x₁, y₁, z₁) and (x₂, Y₂, z₂) denote the respective position vectors defined in the 1st and 2nd views. The spatial relationship between the two views can be characterized based on the recorded gantry information as $\left[\begin{array}{c}{x}_2\\ y_2\\ z_2\end{array}\right]=R·\left\{\left[\begin{array}{c}{x}_1\\ y_1\\ z_1\end{array}\right]-\overrightarrow{t}\right\}=\lfloor \begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\rfloor ·\left[\begin{array}{c}{x}_1-t_x\\ y_1-t_y\\ z_1-t_z\end{array}\right].$

[0120] where

[0121] R=R_x(−α₁)·R_y(β₁)·R_y(−β₂)·R_x(α₁),

[0122] and $\lfloor \begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\rfloor =(R_x\left(-{\alpha }_1\right)·R_y\left({\beta }_1\right)·R_y\left(-{\beta }_2\right)·R_x\left({\alpha }_2\right)·\left[\begin{array}{c}0\\ 0\\ -{\mathrm{Iso}}_2\end{array}\right])-\left[\begin{array}{c}0\\ 0\\ -{\mathrm{Iso}}_1\end{array}\right].$

[0123] The matrices R_x(θ) and R_y(φ) denote the rotation with respect to x-axis and y-axis with θ and φ angles, respectively, and Iso₁and Iso₂denote the two iso-center distances associated with the 1st and 2nd views.

[0124] The scaled 3-D points ({overscore (x)}′_i, {overscore (y)}′_i, {overscore (z)}′_i) defined in the x′y′z′ coordinate system are initialized as ${\stackrel{\_}{x}}_i^{\prime }=\frac{u_i^{\prime }}{{\alpha }^{\prime }·t_d},{\stackrel{\_}{y}}_i^{\prime }=\frac{v_i^{\prime }}{{\alpha }^{\prime }·t_d},{\stackrel{\_}{z}}_i^{\prime }=\frac{D^{\prime }}{{\alpha }^{\prime }·t_d},i=1,\dots ,n,$

[0125] where (u′_i, v′_i) denotes the 2-D input coordinates of the input points. In other words, the initial scaled 3D points are taken to lie in a plane that is parallel to the image plane in the x′y′z′ system and passes through the iso-center.

[0126] Let {right arrow over (μ)}_i={μ_i, μ_i} and {right arrow over (μ)}_i′={μ′_i, μ′_i} denote the respective identified 2-D vessel directional unit vectors of bifurcations in each views. Let v={{right arrow over (v)}₁, {right arrow over (v)}₂, . . . , and {right arrow over (v)}_n} and v′={{right arrow over (v)}₁′, {right arrow over (v)}₂′, . . . , {right arrow over (v)}_n′} denote the calculated 3-D vessel directional unit vectors of bifurcations in two views. The initial estimate of the 3-D directional vectors defined in the two views are characterized as

[0127] {right arrow over (μ)}_i=[μ_i, μ_i, 0]^t, {right arrow over (μ)}_i′=[μ′_i,μ′_i,0]^t

[0128] With the known estimates of the scaled image coordinates (ξ_i, η_i) and (ξ_i, η_i′), transformation in terms of R and {right arrow over (t)}, the object coordinates (x_i,y_i,z_i) can be calculated as follows: $\lfloor \begin{array}{ccc}{r}_{11}-r_{31}{\xi }_i^{\prime }& r_{12}-r_{32}{\xi }_i^{\prime }& r_{13}-r_{33}{\xi }_i^{\prime }\\ r_{21}-r_{31}{\eta }_i^{\prime }& r_{22}-r_{32}{\eta }_i^{\prime }& r_{23}-r_{33}{\eta }_i^{\prime }\\ 1& 0& -{\xi }_i\\ 0& 1& -{\eta }_i\end{array}\rfloor \left[\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right]=\left[\begin{array}{c}\overrightarrow{a}·\overrightarrow{t}\\ \overrightarrow{b}·\overrightarrow{t}\\ 0\\ 0\end{array}\right],$

[0129] where $\overrightarrow{a}=\lfloor \begin{array}{c}{r}_{11}-r_{31}{\xi }_i^{\prime }\\ r_{12}-r_{32}{\xi }_i^{\prime }\\ r_{13}-r_{33}{\xi }_i^{\prime }\end{array}\rfloor ,\overrightarrow{b}=\lfloor \begin{array}{c}{r}_{21}-r_{31}{\eta }_i^{\prime }\\ r_{22}-r_{32}{\eta }_i^{\prime }\\ r_{23}-r_{33}{\eta }_i^{\prime }\end{array}\rfloor .$

[0130] The scaled coordinates of the initial 3D points ({right arrow over (x)}_i, {right arrow over (y)}_i, {right arrow over (z)}_i) defined in the x′y′z! coordinate system are calculated by ®

[0131] ([x_i,y_i,z_it]^t−{right arrow over (t)})/|{right arrow over (t)}|.

[0132] Let {right arrow over (μ)}_i={μ_i, μ_i} and {right arrow over (μ)}′_i={μ′_i, μ′_i} denote the respective identified 2-D vessel directional unit vectors of bifurcations in each view. Let v={{right arrow over (v)}₁, {right arrow over (v)}₂, . . . , {right arrow over (v)}_n} and v′={{right arrow over (v)}′₁, {right arrow over (v)}′₂, . . . , {right arrow over (v)}′_n} denote the calculated 3-D vessel directional unit vectors of bifurcations in two views. The initial estimate of the 3-D directional vectors defined in the two views are characterized as

[0133] {right arrow over (μ)}_i=[μ_i, μ_i,0]^t,{right arrow over (μ)}′_i=[μ′_i,μ′_i,0]^t

[0134] When there occurs a table panning or movement during the image acquisitions, a translation vector {right arrow over (r)}=[r_x,r_y,r_z]^tthat defines the initial solution of such a movement can be estimated if the transformation in terms of a rotation matrix and a translation vector is known as $\underset{\overrightarrow{r}}{\min }ℱ\left(\overrightarrow{r}\right)=\sum _{i=1}^n{\left({\xi }_i^{\prime }-\frac{x_i^{\prime }}{z_i^{\prime }}\right)}^2+{\left({\eta }_i^{\prime }-\frac{y_i^{\prime }}{z_i^{\prime }}\right)}^2+{\left({\xi }_i-\frac{{\overrightarrow{c}}_1·{\overrightarrow{p}}_i^{\prime }+t_x+r_x}{{\overrightarrow{c}}_3·{\overrightarrow{p}}_i^{\prime }+t_z+r_z}\right)}^2+{\left({\eta }_i-\frac{{\overrightarrow{c}}_2·{\overrightarrow{p}}_i^{\prime }+t_y+r_y}{{\overrightarrow{c}}_3·{\overrightarrow{p}}_i^{\prime }+t_z+r_z}\right)}^2,$

[0135] where {right arrow over (c)}_kdenotes the respective kth column vectors of matrix R. The initial estimates of the 3-D bifurcation points {right arrow over (p)}′_i, i=1, . . . , n can be calculated based on the formula as shown in the above equations.

[0136] Given the set of 2-D bifurcation points and directional vectors extracted from the pair of images, “optimal” estimates of the transformation and 3-D bifurcation points are obtained by finding the minimum of the following equation: $\underset{p,p^{\prime },v,v^{\prime }}{\min }{ℱ}_1\left(P,P^{\prime },v,v^{\prime }\right)=\sum _{i=1}^n{\left({\xi }_i-\frac{x_i}{z_i}\right)}^2+{\left({\eta }_i-\frac{y_i}{z_i}\right)}^2+{\left({\xi }_i^{\prime }-\frac{x_i^{\prime }}{z_i^{\prime }}\right)}^2+{\left({\eta }_i^{\prime }-\frac{y_i^{\prime }}{z_i^{\prime }}\right)}^2+{\Vert {\overrightarrow{\mu }}_i-{\left[{\overrightarrow{v}}_i\right]}_{\mathrm{xy}}\Vert }^2+{\Vert {\overrightarrow{\mu }}_i^{\prime }-{\left[{\overrightarrow{v}}_i^{\prime }\right]}_{x^{\prime }{y}^{\prime }}\Vert }^2,$

[0137] where n denotes the number of pairs of corresponding points extracted from the two images and P and P′ denote the sets of 3-D object position vectors (i.e., bifurcation points), {right arrow over (p)}_i=(x_i,y_i,z_i) and {right arrow over (p)}′_i=(x′_i,y′_i, z′_i), respectively, {right arrow over (μ)}_iand {right arrow over (μ)}′_i, denote the respective 2-D vessel directional vectors of bifurcations in each views, and v={{right arrow over (v)}₁, {right arrow over (v)}₂, . . . , {right arrow over (v)}_n} and v′={{right arrow over (v)}′₁, {right arrow over (v)}′₂, . . . , {right arrow over (v)}′_n} denote the calculated 3-D vessel directional vectors of bifurcations in two views, and [{right arrow over (v)}_i]_xyand [{right arrow over (v)}′_i]_x′y′, i=1, . . . ,n, denote the calculated 2-D projection vectors from {right arrow over (v)}_iand {right arrow over (v)}_ionto the image planes, respectively.

[0138] Since the relationship between the two views can be characterized by a rotation matrix R and a translation vector {right arrow over (t)}=[t_x, t_y, t_z]^t, the above equation can be expressed as $\underset{R,\overrightarrow{t},p^{\prime },v^{\prime }}{\min }{ℱ}_2\left(R,\overrightarrow{t},p^{\prime },v^{\prime }\right)=\sum _{i=1}^n{\left({\xi }_i^{\prime }-\frac{x_i^{\prime }}{z_i^{\prime }}\right)}^2+{\left({\eta }_i^{\prime }-\frac{y_i^{\prime }}{z_i^{\prime }}\right)}^2+{\left({\xi }_i-\frac{{\overrightarrow{c}}_1·{\overrightarrow{p}}_i^{\prime }+t_x}{{\overrightarrow{c}}_3·{\overrightarrow{p}}_i^{\prime }+t_z}\right)}^2+{\left({\eta }_i-\frac{{\overrightarrow{c}}_2·{\overrightarrow{p}}_i^{\prime }+t_y}{{\overrightarrow{c}}_3·{\overrightarrow{p}}_i^{\prime }+t_z}\right)}^2+{\Vert {\overrightarrow{\mu }}_i-{\left[R^{-1}·{\overrightarrow{v}}_i^{\prime }\right]}_{\mathrm{xy}}\Vert }^2+{\Vert {\overrightarrow{\mu }}_i^{\prime }-{\left[{\overrightarrow{v}}_i^{\prime }\right]}_{x^{\prime }{y}^{\prime }}\Vert }^2,$

[0139] where {right arrow over (c)}_kdenotes the respective kth column vectors of matrix R and R¹is the inverse matrix of R. When there occurs a table panning or movement during the image acquisitions, however, the above equation can not truly characterize the spatial relationship of the two views. Instead, a translation vector {right arrow over (r)}=[r_x, r_y, r_z]^tneeds to be incorporated to characterize such a situation, and the above equation can be rewritten as $\underset{R,\overrightarrow{t},p^{\prime },v^{\prime },\overrightarrow{r}}{\min }{ℱ}_3\left(R,\overrightarrow{t},p^{\prime },v^{\prime }\overrightarrow{r}\right)=\sum _{i=1}^n{\left({\xi }_i^{\prime }-\frac{x_i^{\prime }}{z_i^{\prime }}\right)}^2+{\left({\eta }_i^{\prime }-\frac{y_i^{\prime }}{z_i^{\prime }}\right)}^2+{\left({\xi }_i-\frac{{\overrightarrow{c}}_1·{\overrightarrow{p}}_i^{\prime }+t_x+r_x}{{\overrightarrow{c}}_3·{\overrightarrow{p}}_i^{\prime }+t_z+r_z}\right)}^2+{\left({\eta }_i-\frac{{\overrightarrow{c}}_2·{\overrightarrow{p}}_i^{\prime }+t_y+r_y}{{\overrightarrow{c}}_3·{\overrightarrow{p}}_i^{\prime }+t_z+r_z}\right)}^2+{\Vert {\overrightarrow{\mu }}_i-{\left[R^{-1}·{\overrightarrow{v}}_i^{\prime }\right]}_{\mathrm{xy}}\Vert }^2+{\Vert {\overrightarrow{\mu }}_i^{\prime }-{\left[{\overrightarrow{v}}_i^{\prime }\right]}_{x^{\prime }{y}^{\prime }}\Vert }^2,$

[0140] It is well known that any 3D rigid motion can be decomposed uniquely into a translation and a rotation by an angle θ around a unit vector {right arrow over (v)}_upassing through the origin of the coordinate system. In the proposed method, a quaternion denoted by {tilde over (q)}=(s,{right arrow over (w)})=(s, w₁, w₂, w₃) of norm equal to 1 is employed to represent the rotation transformation as

[0141] {right arrow over (w)}=sin(θ/2){right arrow over (v)}_u, s=cos(θ/2).

[0142] For any quaternion {tilde over (q)}=(s,{right arrow over (w)})=(s, w₁, w₂, W₃) of norm 1, there exists a rotation matrix defined as: $R\equiv 2\lfloor \begin{array}{ccc}{s}^2+{\left(w_1\right)}^2-\frac{1}{2}& w_1{w}_2-{\mathrm{sw}}_3& {\mathrm{sw}}_2+w_1{w}_3\\ w_1{w}_2+{\mathrm{sw}}_3& s^2+{\left(w_2\right)}^2-\frac{1}{2}& w_2{w}_3-{\mathrm{sw}}_1\\ w_1{w}_3-{\mathrm{sw}}_2& {\mathrm{sw}}_1+w_2{w}_3& s^2+{\left(w_3\right)}^2-\frac{1}{2}\end{array}\rfloor$

[0143] With this quaternion representation, the above equations can be rewritten as $\underset{\tilde{q},\overrightarrow{t},{\hat{p}}^{\prime }}{\min }{ℱ}_4\left(\tilde{q},{\overrightarrow{t}}_u{\hat{p}}^{\prime }\right)=\sum _{i=1}^n{\left({\xi }_i^{\prime }-\frac{{\hat{x}}_i^{\prime }}{{\hat{z}}_i^{\prime }}\right)}^2+{\left({\eta }_i^{\prime }-\frac{{\hat{y}}_i^{\prime }}{{\hat{z}}_i^{\prime }}\right)}^2+{\left\{{\xi }_i-\frac{\begin{array}{c}2\left(s^2+w_1^2-1/2\right){\hat{x}}_i^{\prime }+2\left(w_1{w}_2+{\mathrm{sw}}_3\right){\hat{y}}_i^{\prime }+\\ 2\left(w_1{w}_3-{\mathrm{sw}}_2\right){\hat{z}}_i^{\prime }+t_{u_x}+{\stackrel{\_}{r}}_x\end{array}}{\begin{array}{c}2\left({\mathrm{sw}}_2+w_1{w}_3\right){\hat{x}}_i^{\prime }+2\left(w_2{w}_3-{\mathrm{sw}}_1\right){\hat{y}}_i^{\prime }+\\ 2\left(s^2+w_3^2-\frac{1}{2}\right){\hat{z}}_i^{\prime }+t_{u_z}+{\stackrel{\_}{r}}_z\end{array}}\right\}}^2+{\left\{{\eta }_i-\frac{\begin{array}{c}2\left(w_1{w}_2-{\mathrm{sw}}_3\right){\hat{x}}_i^{\prime }+2\left(s^2+w_2^2-1/2\right){\hat{y}}_i^{\prime }+\\ 2\left({\mathrm{sw}}_1+w_2{w}_3\right){\hat{z}}_i^{\prime }+t_{u_y}+{\stackrel{\_}{r}}_y\end{array}}{\begin{array}{c}2\left({\mathrm{sw}}_2+w_1{w}_3\right){\hat{x}}_i^{\prime }+2\left(w_2{w}_3-{\mathrm{sw}}_1\right){\hat{y}}_i^{\prime }+\\ 2\left(s^2+w_3^2-\frac{1}{2}\right){\hat{z}}_i^{\prime }+t_{u_z}+{\stackrel{\_}{r}}_z\end{array}}\right\}}^2+{\Vert {\overrightarrow{\mu }}_i-{\left[R^{-1}·{\overrightarrow{v}}_i^{\prime }\right]}_{\mathrm{xy}}\Vert }^2+{\Vert {\overrightarrow{\mu }}_i^{\prime }-{\left[{\overrightarrow{v}}_i^{\prime }\right]}_{x^{\prime }{y}^{\prime }}\Vert }^2,$

[0144] where ({circumflex over (x)}′_i, ŷ′_i, {circumflex over (z)}′_i) denotes the scaled 3D coordinates, within a scale factor of the magnitude of the translation vector |{right arrow over (t)}| and where [tu_x,tu_y,tu_z]^tdenotes the unit translation vector that corresponds to {right arrow over (t)}, and {right arrow over (r)}/|{right arrow over (t)}|=[{right arrow over (r)}_x, {right arrow over (r)}_y, {right arrow over (r)}_z]. $\begin{array}{cc}{C}_1\text{:}& s^2+w_1^2+w_2^2+w_3^2=1\\ C_2\text{:}& t_{u_x}^2+t_{u_y}^2+t_{u_z}^2=1\\ C_3\text{:}& 0<{\hat{z}}_i^{\prime },i=1,\dots ,n\\ C_4\text{:}& 0<2\left({\mathrm{sw}}_2+w_1{w}_3\right){\hat{x}}_i^{\prime }+2\left(w_2{w}_3-{\mathrm{sw}}_1\right){\hat{y}}_1^{\prime }+2\left(s^2+w_3^2-\frac{1}{2}\right){\hat{z}}_i^{\prime }+t_{u_z}+{\stackrel{\_}{r}}_z,\\ & i=1,\dots ,n.\end{array}$

[0145] The constraint C₁characterizes the quaternion norm; the constraint C₂ensures a unit translation vector; and the constraints C₃and C₄force the scaled coordinates {circumflex over (z)}′_iand {circumflex over (z)}_ito be in front of the respective focal spots.

[0146] If the iso-center distances of the biplane imaging systems or the MF factors are available, the constraints C₃and C₄above can be modified as: $\frac{D^{\prime }/{\alpha }^{\prime }-{\delta }_h}{\left|\overrightarrow{t}\right|}<{\hat{z}}_i^{\prime }<\frac{D^{\prime }/{\alpha }^{\prime }+{\delta }_h}{\left|\overrightarrow{t}\right|},\dots ,n$ $\frac{D/\alpha -{\delta }_h}{\left|\overrightarrow{t}\right|}<{\hat{z}}_i=2\left[\left({\mathrm{sw}}_2+w_1{w}_3\right){\hat{x}}_i^{\prime }+\left(w_2{w}_3-{\mathrm{sw}}_1\right){\hat{y}}_i^{\prime }+\left(s^2+w_3^2-\frac{1}{2}\right){\hat{z}}_i^{\prime }\right]+t_{u_z}+{\stackrel{\_}{r}}_z\dots ,n,$

[0147] where α and α′ denote the MF factors and δ_h(≈12.5±2.0 cm at end-diastole) denotes the maximal length of the heart in a long-axis view. For each 3D object point, the ray connecting that point and the focal spot intersects the image plane near the associated 2D image point even when the input data are corrupted by noise. In addition to the constraints imposed on the z and z′ coordinates, two other constraints are incorporated to confine the x, x′, y, and y′ coordinates of each calculated 3D point as follows: $C_5\text{:}{\left(\frac{{\hat{x}}_i}{{\hat{z}}_i}-{\xi }_i\right)}^2+{\left(\frac{{\hat{y}}_i}{{\hat{z}}_i}-{\eta }_i\right)}^2\le {\left(\frac{{\delta }_c}{p_{\mathrm{size}}}\right)}^2,i=1,\dots ,n,C_6\text{:}{\left(\frac{{\hat{x}}_i^{\prime }}{{\hat{z}}_i^{\prime }}-{\xi }_i^{\prime }\right)}^2+{\left(\frac{{\hat{y}}_i^{\prime }}{{\hat{z}}_i^{\prime }}-{\eta }_i^{\prime }\right)}^2\le {\left(\frac{{\delta }_c}{p_{\mathrm{size}}}\right)}^2,i=1,\dots ,n,$

[0148] where δ_cdefines the radius of a circular disk (e.g., 20 pixels) centered at (ξ_i, η_i) or (ξ′_i, η′_i). The constraints C₃to C₆form a cone-shape bounding region for each 3D point with the focal spot as its apex and the a circular disk (of radius δ_c) centered at (ξ_i, η_i) as its base to confine the searching space of the calculated 3D object point which is shown in FIG. 19.

[0149] D. Establishment of Vessel Centerline Correspondence

[0150] Details of the methods for establishing the vessel centerline correspondence of the arterial tree are shown in the flowchart of FIG. 20.

[0151] After the transformation (R, {right arrow over (t)}) that defines the two views is obtained, this information will be used for establishing the point correspondences on 2-D vessel centerlines in the pair of images and calculating 3-D morphologic structures of coronary arterial tree. The calculated transformation in conjunction with the epi-polar constraints were employed as the framework for establishing the point correspondences on the vessel centerlines based on the two identified 2-D coronary arterial trees and the vessel hierarchy as shown in FIG. 21.

[0152] If the information on relative positions of two cameras is known (i.e., the known locations and orientations of two focal spots F_aand F_b), the correspondences of image points can be solved by use of “epi-polar constraints.” Let P denote a 3-D point, and let P_aand P_b, the projections of P, denote the pair of corresponding points in two images. Such constraints state that these points, F_a, F_b, P_a, P_b, and P, are all on one plane called the epi-polar plane as shown in FIG. 22. The epi-polar plane and each image plane intersect along a straight line called the epi-polar line (see FIG. 23). Therefore, the location of point P_bin the second image must lie on the epi-polar line resulting from the intersection between the second image plane and the plane defined by point P_aand the two focal spots. A typical example of local correspondence establishment by use of epi-polar constraints where points A, B, C, and D at the first view are defined by finding the intersections between the 2-D curve associated with left anterior descending artery (LAD) artery and respective epi-polar lines defined by a, b, c, and d at the second view as shown in FIG. 24.

[0153] If the epi-polar line is “tangential” relative to the 2-D vessel centerline curve, the accuracy of calculated intersection point becomes sensitive to any small errors in epi-polar line derived based on the calculated transformation in terms of R and {right arrow over (t)}. To overcome the problem, a refinement process as shown in FIG. 25 will be employed to calculate the global correspondences using optimal parametric arguments t₁, t₂, . . . t_n−1, t_nbased on the equations below. $\underset{\hat{t},p^{\prime },v^{\prime }}{\min }{ℱ}_{\mathrm{epi}}\left(\hat{t},P^{\prime },v^{\prime }\right)=\sum _{i=1}^n{\left({\xi }_i^{\prime }-\frac{x_i^{\prime }}{z_i^{\prime }}\right)}^2+{\left({\eta }_i^{\prime }-\frac{y_i^{\prime }}{z_i^{\prime }}\right)}^2+{\left(f_x\left(t_i\right)-\frac{{\overrightarrow{c}}_1·{\overrightarrow{p}}_i^{\prime }+t_x}{{\overrightarrow{c}}_3·{\overrightarrow{p}}_i^{\prime }+t_z}\right)}^2+{\left(f_y\left(t_i\right)-\frac{{\overrightarrow{c}}_2·{\overrightarrow{p}}_i^{\prime }+t_y}{{\overrightarrow{c}}_3·{\overrightarrow{p}}_i^{\prime }+t_z}\right)}^2+{\Vert f^{\prime }\left(t_i\right)-{\left[R^{-1}·{\overrightarrow{v}}_i^{\prime }\right]}_{\mathrm{xy}}\Vert }^2+{\Vert {\overrightarrow{\mu }}_i^{\prime }-{\left[{\overrightarrow{v}}_i^{\prime }\right]}_{x^{\prime }{y}^{\prime }}\Vert }^2,$

[0154] where {circumflex over (t)}={t₁, t₂, . . . , t_n}, P′={{right arrow over (p)}′_i=(x′₁,y′_i,z′_i) | i=1, . . . , n}, and v′={{right arrow over (v)}₁, {right arrow over (v)}₂, . . . , {right arrow over (v)}_n} denote the set of parametric variables, the set of 3-D object position vectors, and the set of 3-D vessel directional vectors to be optimized, respectively, n denotes the number of points on the 2-D vessel centerline, (f_x(t_i), f_y(t_i)) and (ξ_i, η′_i) are the respective spline-curve modeled and extracted 2-D vessel centerline points in two views, f′ (t_i)=(f′_x(t_i), f′_y(t_i)) and {right arrow over (μ)}′_idenote the respective spline-curve modeled and extracted 2-D directional vectors of vessel centerlines in each view, [R⁻¹{right arrow over (v)}_i]_xyand [{right arrow over (v)}_i]_x′y′ denote the projection of the calculated 3-D directional vector on the respective image planes, {right arrow over (c)}_k, k=1,2,3, denotes the respective kth column vectors of matrix R, and R⁻¹is the inverse matrix of R. The 3-D position vectors P′≡{{right arrow over (p)}_i=(x′_i,y′_i,z′_i) | i=1, . . . , n} on vessel centerline are computed by use of equations described above where the 2-D pair correspondence (ξ,η_i) are substituted by the 2-D spline-curve modeling function f(t_i)=(f_x(t_i),f_y(t_i)). Based on the equations, the vessel centerline correspondences are globally established by incorporating the entire vessel shape in terms of directions and locations that will yield more accurate results than those obtained by only use of epi-polar constraints with local vessel segment points; especially when epi-polar line and vessel segment are tangential. FIGS. 26(c) and (d) show the results of 3-D LCA coronary reconstructed from a pair of angiograms as shown in FIGS. 26(a) and (b) based on the simple epi-polar technique and the refinement process, respectively. The reconstruction of left main artery apparently illustrates inaccurate results based on the simple epi-polar process which were corrected after employing the refinement process.

[0155] E. Calculation of 3-D Skeleton of Coronary Arterial Tree

[0156] Details of the methods for calculation of the 3-D skeleton of the arterial tree are shown in the flowchart of FIG. 27.

[0157] With the point correspondences on 2-D vessel centerlines (ξ_j, η_j) and (ξ′_j, η′_j) and the transformation (R, {right arrow over (t)}), the 3-D vessel centerline points of coronary arteries (x_j, y_j, z_j) can then be calculated as shown in equations below. For each artery, the cross-sectional diameters at the two end-points can be calculated as $d_i={\stackrel{\_}{d}}_i·\left(\sqrt{x_i^2+y_i^2+z_i^2}/\sqrt{{\xi }_i^2+{\eta }_i^2+D^2}\right),$

[0158] where {right arrow over (d)} is the extracted 2-D vessel diameter centered at (ξ_i, η_i), and D is the distance between the focal spot and the image intensifier. The remaining cross-sectional diameters on a artery are linearly interpolated based on the calculated diameters at two end-points as shown in FIG. 29.

[0159] In the first view, let (u_i, v_i) denote the image coordinates of the ith object point, located at position {right arrow over (p)}_i=(x_i, y_i, z_i) relative to the focal spot of the first view. Then $u_i=\frac{{\mathrm{Dx}}_i}{z_i},v_i=\frac{{\mathrm{Dy}}_i}{z_i},$

[0160] where D is the SID (focal-spot to imaging-plane distance) in the first imaging system. Let (ξ_i, η_i) denote the scaled image coordinates, defined as ${\xi }_i=\frac{u_i}{D}=\frac{x_i}{z_i},{\eta }_i=\frac{v_i}{D}=\frac{y_i}{z_i}.$

[0161] The second projection view of the biplane imaging system can be described in terms of the SID D′, and a second pair of image and object coordinate systems u′v′ and x′y′z′ can be defined in an analogous manner. Scaled image coordinates (ξ′_i, η′_i) in the second view for the ith object point at position {right arrow over (p)}_i=(x′_i, y′_i, z′_i) are then given by ${\xi }_i^{\prime }=\frac{u_i^{\prime }}{D^{\prime }}=\frac{x_i^{\prime }}{z_i^{\prime }},{\eta }_i^{\prime }=\frac{v_i^{\prime }}{D^{\prime }}=\frac{y_i^{\prime }}{z_i^{\prime }}.$

[0162] With the known estimates of the scaled image coordinates (ξ_i, η_i) and (ξ′_i, η′_i), transformation in terms of R and {right arrow over (t)}, the object coordinates (x_i, y_i, z_i) can be calculated as follows: $\lfloor \begin{array}{ccc}{r}_{11}-r_{31}{\xi }_i^{\prime }& r_{12}-r_{32}{\xi }_i^{\prime }& r_{13}-r_{33}{\xi }_i^{\prime }\\ r_{21}-r_{31}{\eta }_i^{\prime }& r_{22}-r_{32}{\eta }_i^{\prime }& r_{23}-r_{33}{\eta }_i^{\prime }\\ 1& 0& -{\xi }_i\\ 0& 1& -{\eta }_i\end{array}\rfloor \left[\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right]=\left[\begin{array}{c}\overrightarrow{a}·\overrightarrow{t}\\ \overrightarrow{b}·\overrightarrow{t}\\ 0\\ 0\end{array}\right],$

[0163] where $\overrightarrow{a}=\lfloor \begin{array}{c}{r}_{11}-r_{31}{\xi }_i^{\prime }\\ r_{12}-r_{32}{\xi }_i^{\prime }\\ r_{13}-r_{33}{\xi }_i^{\prime }\end{array}\rfloor ,\overrightarrow{b}=\lfloor \begin{array}{c}{r}_{21}-r_{31}{\eta }_i^{\prime }\\ r_{22}-r_{32}{\eta }_i^{\prime }\\ r_{23}-r_{33}{\eta }_i^{\prime }\end{array}\rfloor .$

[0164] F. Rendering of 3-D Coronary Arterial Tree with Associated Gantry Angulation

[0165] Details of the methods for rendering the 3-D arterial tree with associated gantry angulation are shown in the flowchart of FIG. 28.

[0166] To obtain the geometrical continuity of the vessel centerline, a spline-based smoothing process is applied to the originally calculated data sets in accord with the equations below. The surface-based modeling technique is employed to recover the complete morphology of the 3-D arterial lumen. For every two consecutive cross sections centered at (x_i, y_i, z_i) and (x_i+1, y_i+1, z_i+1), a sequence of polygonal patches are formed with the corresponding cross sectional diameters as shown in FIG. 30.

[0167] When an arbitrary computer-generated image is produced, the gantry information defining the current projection is calculated in the form of left anterior oblique/right anterior oblique (LAO/RAO) and caudal/cranial (CAUD/CRAN). These gantry angles are defined in a spatial coordinate system with the iso-center as the origin. The LAO/RAO angulation is defined on the y-z plane, while the CAUD/CRAN angulation is defined on the x-z plane. With the spatial coordinates, the position of focal spot (x_f, y_f, z_f) can be formulated by use of two rotations R_x(α) and R_y(−β) as presented in equations to follow.

[0168] The reconstructed 3-D coronary arterial tree can be rotated to any selected viewing angle yielding multiple computer-generated projections to improve the quality and utility of the images subsequently acquired as shown in FIGS. 31a-fand 32a-f.

[0169] Due to the poor imaging resolution in the depth or z-axis, the 3-D shape of coronary vessel tree become joggling. To overcome the drawbacks, the technique of Bezier curve for data smoothing is adopted. The Bezier curve can be defined as a parametric function as follows: $P\left(u\right)=\sum _{k=0}^n{p}_k{B}_{k,n}\left(u\right),p_k=\left(x_k,y_k,z_k\right).$

[0170] where p_kdenote the n+1 input control points. Each B_k,n(u) is a polynomial function of degree one less than the number of control points used and is defined as  B_k,n(u)=C(n, k)u^k(1−u)^n−k

[0171] and the C(n,k) represent the binomial coefficients $C\left(n,k\right)=\frac{n!}{k!\left(n-k\right)!}$

[0172] By applying the Bezier equation, a new representation of 3-D coronary vessels can be obtained.

[0173] When an arbitrary computer-generated image is produced, the gantry information defining the current projection is calculated in the form of left anterior oblique/right anterior oblique (LAO/RAO)and caudal/cranial (CAUD/CRAN). These gantry angles are defined in a spatial coordinate system with the iso-center as the origin. The LAO/RAO angulation is defined on the y-z plane, while the CAUD/CRAN angulation is defined on the x-z plane. With the spatial coordinates, the position of focal spot (x_f, y_f, z_f) can be formulated by use of two rotations R_x(α) and R_y(−β) as, $\lfloor \begin{array}{c}{x}_f\\ y_f\\ z_f\end{array}\rfloor =R_y\left(-\beta \right)R_x\left(\alpha \right)·\lfloor \begin{array}{c}{x}_n\\ y_n\\ z_n\end{array}\rfloor =\left[\begin{array}{ccc}\cos \left(-\beta \right)& \sin \left(\alpha \right)\sin \left(-\beta \right)& \cos \left(\alpha \right)\sin \left(-\beta \right)\\ 0& \cos \left(\alpha \right)& -\sin \left(\alpha \right)\\ -\sin \left(-\beta \right)& \sin \left(\alpha \right)\cos \left(-\beta \right)& \cos \left(\alpha \right)\cos \left(-\beta \right)\end{array}\right]·\left[\begin{array}{c}{x}_n\\ y_n\\ z_n\end{array}\right],$

[0174] where (x_n, Y_n, z_n)=(0, 0, D_fc) denotes the neutral position of focal spot (or anterior-posterior (AP) view with 0° LAO and 0° CAUD, D_fcis the distance between iso-center and focal spot which is provided by the manufacture), R_xand R_ydenotes the rigid rotations with respect to x-axis and y-axis, and α(−α) and β(−β) denote the LAO (RAO) and CAUD (CRAN) angles, respectively.

[0175] G. 3-D Coronary Quantitative Analysis: Optimal Views and Morphology Assessment

[0176] FIG. 37 provides a flowchart describing the high level operations included in the improved QCA of the present invention. Specifically, the method of FIG. 37 presumes the availability of 3-D reconstruction image data of a coronary arterial tree. The methods presented above are useful for generation of such 3-D reconstructed data from 2-D coronary images. Those skilled in the art will recognize a variety of other sources for such 3-D image data as well as alternate techniques for reconstructing such data from 2-D images.

[0177] Element 3700 is first operable to determine whether the user of the QCA method desire to assess the geometry of the coronary artery reconstruction or whether the user wishes to determine optimal view of the 3-D image for particular parameters. If the user is assessing the geometry, element 3702 is operable to determine what aspect of the arterial image geometry is to be assessed. Element 3702 is operable to evaluate take-off angles of each arterial branch in the arterial tree image presently under study. Element 3706 is operable to evaluate the insertion angle a bypass graft selected in the arterial tree image presently under study. Element 3708 is operable to evaluate the length of a selected arterial segment from the arterial tree image presently under study. And element 3710 is operable to evaluate the tortuosity of a selected arterial pathway from the arterial tree image presently under study. Following completion of the selected evaluative task processing is completed or may continue by looping back to element 3700 to await a next selection from the arteriographer.

[0178] If the user wishes to determine optimal view for particular parameters as determined by element 3700, element 3712 is next operable to determine which view the user desires for purposes of optimizing. Element 3714 is operable to present the user with a graphical map useful for viewing and optimizing foreshortening of the coronary arterial tree image presently under study. Element 3716 is operable to present the user with a graphical map useful for viewing and optimizing overlap of the arterial tree image presently under study. Element 3718 is operable to present the user with a graphical map useful for optimizing the combination offoreshortening and overlap in the arterial tree image presently under study. Element 3720 is operable to present to the user information indicative of the accumulated percent foreshortening of each arterial segment divided by bifurcation points in the arterial tree image presently under study. Following completion of the selected view optimization, processing is completed or may continue by looping back to element 3700 to await a next selection from the arteriographer.

[0179] Each of the above QCA tasks is further described in detail herein below with reference to additional figures as required for a clear understanding.

[0180] G.1. Assessment of the Coronary Geometry

[0181] The curvilinear nature of vessels is important in the pathophysiology of coronary artery disease and the performance of interventions. Therefore such a geometric property can not be quantified from 2-D projection data. The skeleton of a 3-D vessel can be mathematically defined as a curve function p(s)=(x(s), y(s), z(s)) connecting all the 3-D centerline points as shown in FIG. 33, where 0≦s≦1 is the parametric variable.

[0182] G.1.1. Evaluation of Length of Arterial Segment

[0183] As noted above, element 3708 is operable to evaluate the length of a selected artery segment of the arterial tree image presently under study. FIG. 47 is a flowchart describing additional details of the operation of element 3708 to evaluate the length of a selected artery segment in the arterial tree image presently under study. In general, the 3-D vessel centerline is characterized by a parametric curve function in terms of an ordered sequence of 3-D data points as described previously. The length of any selected artery segment on the vessel centerline can be calculated by the summation of the inter-distance of every pair of consecutive vessel centerline points.

[0184] More specifically, element 4700 is operable to accept user input defining two points on the displayed arterial tree. The user provides the two points by clicking a pointer device at each of the desired points. The two selected points define a segment of the displayed arterial tree for which the user wishes to evaluate the length. Element 4702 then determines whether a path exists in the arterial tree from the start of the displayed tree through bot selected points. Calculation of the segment length from the 3 -D data requires a path exists from the start of the image tree through both selected points. Details of this element are provided further herein below. Element 4704 then tests whether the required path is found. If not, processing loops back to element 4700 to permit the user to select different points for a desired segment. If the required path through the elected segment exists, processing continues with element 4706 to calculate the length between the two identified points along the path identified by element 4702.

[0185] FIG. 51 is a flowchart providing additional details of the processing of element 4702 (and 4802 discussed below) to determine the connectivity of a selected arterial pathway defined by two selected end points. In particular, element 5100 uses the vessel hierarchy data structure to start at the distal end point of the selected pathway. Element 5102 then traverses the vessel centerline point sequence from the present point in the upstream direction. Element 5104 then determines whether the selected proximal end point of the selected pathway is encountered in the next point traversed by element 5102. If the proximal end point is encountered, then the selected pathway is connected and element 5106 is operable to return a “CONNECTED” status to the invoking procedure. If the proximal end point is not yet encountered, element 5108 is operable to determine whether the beginning of the coronary tree is encountered in this next point traversed by element 5102. If so, element 5110 is operable to return a “NOT CONNECTED” status to the invoking procedure because the beginning of the tree is encountered without traversing through the proximal end point. If the beginning of the tree is not yet encountered, processing loops back to element 5102 to traverse a next point along the selected path and repeat the tests of elements 5104 and 5108.

[0186] G. 1.2. Evaluation of the tortuosity of artery

[0187] Element 3710 is operable to evaluate the tortuosity of a selected arterial segment pathway. FIG. 48 is a flowchart providing additional details of the processing to evaluate tortuosity of a selected arterial segment.

[0188] G. 1.2.1. Evaluation of the tortuosity of an artery by the F-S theory

[0189] Generally, by use of such spline-based curve modeling technique, we are able to apply the theory of differential geometry such as Frenet-Serret apparatus (or F-S theory) or its variation to study the geometrical nature of the 3-D coronary artery shape at any time point in the cardiac cycle. A technique based on the Frenet-Serret (F-S) theory of differential geometry has been developed to study the geometrical nature or tortuosity of the 3-D coronary artery shape at any time point in the cardiac cycle. The F-S theory consists of five components: three vector fields along the given curve (the tangent T(s), the normal N(s), and the bi-normal B(s) vectors) and two scalar valued functions (the curvature κ(s) and the torsion τ(s) where s denotes the parametric variable defining the location of point on the curve function α(s). The curvature κ(s₀) measures the rate of change of the angle defined by the two neighboring tangents T(s₀) and T(s₀+δ_s) associated with the two respective points α(s₀) and α(s₀+δ_s). In other word, it characterizes the degree of bending pertinent to a local curve segment in the neighborhood between s₀and s₀+δ_s(i.e., how rapidly the curve pulls away from the plane Ω_nperpendicular to the tangent vector at T(s₀)) as shown in FIG. 34a. Similarly, the torsion at τ(s₀) measures a rate of change (or a twisting) at a point α(s₀) how its trajectory twists out of the plane Ω_nperpendicular to the normal vector B(s₀)) as shown in FIG. 34b.

[0190] In FIG. 34c, a circle defining on the x-y plane of radius r>0 is illustrated, where the vector fields T(s) and N(s) vary along the curve and the binormal field B(s) is constant. Note that the curvature κ(s) (for all s) of a circle is equal to the reciprocal of its radius (i.e., 1/r). Since the binormal field B is a constant, the torsion at every point is equal to 0. The example illustrated in FIG. 34dshows a circular helix in which has constant curvature and torsion at every point in the curve. Given the curvature and torsion functions κ(s) and τ(s), they completely describe the curve α(s) geometrically.

[0191] The calculation of curvature and torsion at every vessel centerline point s₀is characterized by the following equations:

[0192] T(s₀)=α⁽¹⁾(s₀)/|α⁽¹⁾(s₀)|

[0193] B(s₀)=α⁽¹⁾(s₀)×α⁽²⁾(s₀)/|α⁽¹⁾(s₀)×α⁽²⁾(s₀)|

[0194] N(s₀)=B(s₀)×T(s₀)

[0195] κ(s₀)=|α⁽¹⁾(s₀)×α⁽²⁾(s₀)|/|α⁽¹⁾(s₀)|³

[0196] τ(s₀)=[α⁽¹⁾(s₀)×α⁽²⁾(s₀)]·α⁽³⁾(s₀)/|α⁽¹⁾(s_)×α⁽²⁾(s₀)|²,

[0197] where α⁽¹⁾(s) denotes the ith derivative with respect to s. In FIG. 38a, the length of the selected pathway is 10.0 cm. The narrowest diameter of the diseased arterial segment is 3.4 mm. The measurement of tortuosity associated with this pathway was illustrated in FIG. 38b.

[0198] 1.2.2. Evaluation of the Tortuosity of Artery by Point-to-Point Method

[0199] A new method, which is a variation of F-S theory, is employed to characterize the vessel tortuosity in order to facilitate the study of vessel straightening due to various intra-coronary devices. Three types of vessel tortuosity are defined: (1) point-to-point (or local) tortuosity, (2) point-to-line (or relative) tortuosity, and (3) tortuosity spectrum. Given an arterial segment, the point-to-point (P-to-P) tortuosity at each vessel centerline point is estimated by calculation of the differential of tangent vectors (i.e., the angle in degree between the two tangent vectors) associated with two consecutive arterial points on the 3-D centerline curve. The point-to-line (P-to-L) tortuosity at each vessel centerline point is defined by the differential of the tangent vector and the directional vector of the line connecting the two end points of the given arterial segment (i.e., the angle in degree between the tangent vector of a vessel centerline point and the directional vector of line). Afterwards, the mean values of the P-to-P tortuosity and P-to-L tortuosity are calculated. To evaluate the “global variation” of vessel geometry, the tortuosity spectrum is derived by calculation of the normalized integration of P-to-L tortuosity along the arterial segment.

[0200] FIG. 39ashows the reconstructed 3-D bypass graft. In FIG. 39bthe P-to-P (bottom) and P-to-L (the 2nd graph from the bottom) tortuosity measurements associated with the selected 3 arterial segments defined between the proximal and distal locations in FIG. 39aare illustrated, respectively. The mean values of regional P-to-P and P-to-Line are (7.6,6.8, 1.9) and (37.5,27.5,9.0) degrees, respectively. The total mean values of P-to-P and P-to-Line are 5.4 and 23.3 degrees, respectively. The tortuosity spectrum in terms of the normalized integration of P-to-L measurements and the diameters along the arterial segments are shown at the 2nd graph from the top and the top graphs in FIG. 39b.

[0201] Referring again to FIG. 48, elements 4800 . . . 4804 are operable in a manner analogous to elements 4700 . . . 4704 of FIG. 47. Specifically, the user is prompted to enter two points on the arterial tree image presented to define a segment of interest. When a valid segment is so selected, element 4806 is then operable to determine whether the user wishes to evaluate the tortuosity in accordance with the point-to-point method or the F-S method both described above. If the user selects the point-to-point (p-to-p) method, elements 4808 and 4810 are operable to calculate the point-to-point regional curvature and point-to-line regional curvature, respectively, of each point on the vessel centerline as described above. If the F-S method is selected by the user, elements 4812 and 4814 are operable to calculate the curvature and torsion, respectively, of each vessel centerline point along the selected segment pathway as described above.

[0202] G. 1.3. Evaluation of the Take-off Angle

[0203] Element 3704 is operable to evaluate the take-off angles of a selected arterial segment pathway. FIG. 49 is a flowchart providing additional details of the processing to evaluate take-off angles of a selected arterial segment. The 3-D take-off angle of a branch vessel is estimated by calculating the angle between the two tangent vectors originating from the bifurcation; i.e., one toward the direction of precedent vessel and the other toward the direction of branching vessel as shown in FIG. 35a. The 2-D take-off angle of a branch vessel is obtained by first computing the branch vessel and preceding vessel centerlines projected onto the imaging plane based on the gantry orientation (α,β) followed by calculating the angle of the 2-D tangent vectors as shown in FIG. 35b. When an arbitrary computer-generated image is produced, the gantry information defining the current projection is calculated in the form of left anterior oblique/right anterior oblique (LAO/RAO) and caudal/cranial (CAUD/CRAN). Mathematically, these gantry angles are defined in a spatial coordinate system with the iso-center as the origin. The LAO/RAO angulation is defined on the y-z plane, while the CAUD/CRAN angulation is defined on the x-z plane. With the gantry viewing angle (α,β), the 2-D projection coordinates (u_(α,β),ν_(α,β)of a vessel centerline point can be calculated by use of the following equations: $\lfloor \begin{array}{c}{u}_{\left(\alpha ,\beta \right)}\\ v_{\left(\alpha ,\beta \right)}\end{array}\rfloor =\lfloor \begin{array}{c}{x}_{\left(\alpha ,\beta \right)}/z\left(\alpha ,\beta \right)\\ x_{\left(\alpha ,\beta \right)}/z\left(\alpha ,\beta \right)\end{array}\rfloor \left[\begin{array}{c}{x}_{\left(\alpha ,\beta \right)}\\ y_{\left(\alpha ,\beta \right)}\\ z_{\left(\alpha ,\beta \right)}\end{array}\right]=R_y\left(-\beta \right)R_x\left(\alpha \right)·\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{ccc}\cos \left(-\beta \right)& \sin \left(\alpha \right)\sin \left(-\beta \right)& \cos \left(\alpha \right)\sin \left(-\beta \right)\\ 0& \cos \left(\alpha \right)& -\sin \left(\alpha \right)\\ -\sin \left(-\beta \right)& \sin \left(\alpha \right)\cos \left(-\beta \right)& \cos \left(\alpha \right)\cos \left(-\beta \right)\end{array}\right]·\left[\begin{array}{c}x\\ y\\ z\end{array}\right],$

[0204] where R_xand R_ydenotes the rigid rotations with respect to x-axis, (x_(α,β), γ_(α,β), z_(α,β)) denotes the 3-D coordinates of a vessel centerline point defined at the current gantry coordinate system, (x,y,z) denotes the 3-D coordinates of a vessel centerline point in the AP viewing angle (0 degree LAO and 0 degree CAUD) coordinate system, and α (−α) and β (−β) denote a gantry position at LAO (RAO) and CAUD (CRAN) angles, respectively. In FIGS. 40(a) and (b), the 3-D and 2-D take-off angles of individual coronary arteries are shown based on two different gantry viewing angles.

[0205] The above equation can be employed to model the majority of existing imaging systems having only two degrees of freedom. Another variable R_z(γ), where R_zdenotes the rigid rotation with respect to z-axis by γ degrees, can be readily added to emulate the imaging system if a rotational C-arm base is employed. With the new third degree of freedom in term of the new angle variable, the above equation can be re-written as follows. $\lfloor \begin{array}{c}{u}_{\left(\alpha ,\beta ,\gamma \right)}\\ v_{\left(\alpha ,\beta ,\gamma \right)}\end{array}\rfloor =\lfloor \begin{array}{c}{x}_{\left(\alpha ,\beta ,\gamma \right)}/z\left(\alpha ,\beta ,\gamma \right)\\ x_{\left(\alpha ,\beta ,\gamma \right)}/z\left(\alpha ,\beta ,\gamma \right)\end{array}\rfloor \left[\begin{array}{c}{x}_{\left(\alpha ,\beta ,\gamma \right)}\\ y_{\left(\alpha ,\beta ,\gamma \right)}\\ z_{\left(\alpha ,\beta ,\gamma \right)}\end{array}\right]=R_y\left(-\beta \right)R_x\left(\alpha \right)R_z\left(\gamma \right)·\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=\left[\begin{array}{ccc}C\left(-\beta \right)C\left(\gamma \right)+S\left(\alpha \right)S\left(-\beta \right)S\left(\gamma \right)& -S\left(\gamma \right)C\left(-\beta \right)+S\left(\alpha \right)S\left(-\beta \right)C\left(\gamma \right)& C\left(\alpha \right)S\left(-\beta \right)\\ C\left(\alpha \right)S\left(\gamma \right)& C\left(\alpha \right)S\left(\gamma \right)& -S\left(\alpha \right)\\ -S\left(-\beta \right)C\left(\gamma \right)+S\left(\alpha \right)C\left(-\beta \right)S\left(\gamma \right)& S\left(-\beta \right)S\left(\gamma \right)+S\left(\alpha \right)C\left(-\beta \right)C\left(\gamma \right)& C\left(\alpha \right)C\left(-\beta \right)\end{array}\right]·\left[\begin{array}{c}x\\ y\\ z\end{array}\right],$

[0206] where C(x) and S(x) denote the cos and sin functions respectively.

[0207] FIG. 49 is a flowchart providing additional details of the processing of element 3704 of FIG. 37 comprising a method for computing take-off angles in accordance with a 2-D or 3-D model. Element 4900 is first operable to determine which model the user has requested. If a 3-D model is used, elements 4902 and 4904 compute the two required tangent vectors from the bifurcation point downstream and from the bifurcation point upstream, respectively. If a 2-D model is used, elements 4906 and 4908 calculate the required two tangent vectors as originated from the bifurcation point downstream in the direction of the branch and from the bifurcation point upstream in the direction of the preceding vessel, respectively. In both models, element 4910 then computes the take-off angel as the angle between the two tangent vectors as determined in accordance with the requested model.

[0208] G. 1.4. Evaluation of the Insertion Angle

[0209] Element 3706 is operable to evaluate the insertion angle of a selected arterial segment pathway. FIG. 50 is a flowchart providing additional details of the processing to evaluate the insertion angle of a selected arterial segment. The insertion angles (upstream and downstream angles) of bypass graft at the end-to-side anastomosis are estimated by calculation of the angles between the 3-D tangent vectors at the native vessels and the tangent vector at the graft as shown in FIG. 36. Element 5000 of FIG. 50 is operable to receive operator input selecting two end points on the arterial tree having the anastomosis located between them. Element 5002 then receives operator input to select one point on the native artery connecting to the anastomosis. Element 5004 then calculates the insertion angle as the upstream and downstream angles defined by the three selected points.

[0210] G.2. Determination of Optimal View Strategy

[0211] When an arbitrary computer-generated image is produced, the gantry information defining the current projection is calculated in the form of left anterior oblique/right anterior oblique (LAO/RAO)and caudal/cranial (CAUD/CRAN). These gantry angles are defined in a spatial coordinate system with the iso-center as the origin. The LAO/RAO angulation is defined on the y-z plane, while the CAUD/CRAN angulation is defined on the x-z plane. With the spatial coordinates, the position of focal spot (x_f, y_f, z_f) can be formulated by use of two rotations R_x(α) and R_y(−β) as, $\lfloor \begin{array}{c}{x}_f\\ y_f\\ z_f\end{array}\rfloor =R_y\left(-\beta \right)R_x\left(\alpha \right)·\lfloor \begin{array}{c}{x}_n\\ y_n\\ z_n\end{array}\rfloor =\left[\begin{array}{ccc}\cos \left(-\beta \right)& \sin \left(\alpha \right)\sin \left(-\beta \right)& \cos \left(\alpha \right)\sin \left(-\beta \right)\\ 0& \cos \left(\alpha \right)& -\sin \left(\alpha \right)\\ -\sin \left(-\beta \right)& \sin \left(\alpha \right)\cos \left(-\beta \right)& \cos \left(\alpha \right)\cos \left(-\beta \right)\end{array}\right]·\left[\begin{array}{c}{x}_n\\ y_n\\ z_n\end{array}\right],$

[0212] where (x_n, y_n, z_n)=(0, 0, D_fc) denotes the neutral position of focal spot or anterior-posterior (AP) view with 0° LAO and 0° CAUD, D_fcis the distance between iso-center and focal spot which is provided by the manufacture), R_xand R_ydenotes the rigid rotations with respect to x-axis andy-axis, and α (−α) and β (−β) denote the LAO (RAO) and CAUD (CRAN) angles, respectively.

[0213] The above equation can be employed to model the majority of existing imaging systems having only two degrees of freedom. Another variable R_z(γ), where R_zdenotes the rigid rotation with respect to z-axis by γ degrees, can be readily added to emulate the imaging system if a rotational C-arm base is employed. With the new third degree of freedom in term of the new angle variable, the above equation can be re-written as follows. $\lfloor \begin{array}{c}{x}_f\\ y_f\\ z_f\end{array}\rfloor =R_y\left(-\beta \right)R_x\left(\alpha \right)R_z\left(\gamma \right)·\lfloor \begin{array}{c}{x}_n\\ y_n\\ z_n\end{array}\rfloor =\left[\begin{array}{ccc}C\left(-\beta \right)C\left(\gamma \right)+S\left(\alpha \right)S\left(-\beta \right)S\left(\gamma \right)& -S\left(\gamma \right)C\left(-\beta \right)+S\left(\alpha \right)S\left(-\beta \right)C\left(\gamma \right)& C\left(\alpha \right)S\left(-\beta \right)\\ C\left(\alpha \right)S\left(\gamma \right)& C\left(\alpha \right)S\left(\gamma \right)& -S\left(\alpha \right)\\ -S\left(-\beta \right)C\left(\gamma \right)+S\left(\alpha \right)C\left(-\beta \right)S\left(\gamma \right)& S\left(-\beta \right)S\left(\gamma \right)+S\left(\alpha \right)C\left(-\beta \right)C\left(\gamma \right)& C\left(\alpha \right)C\left(-\beta \right)\end{array}\right]·\left[\begin{array}{c}{x}_n\\ y_n\\ z_n\end{array}\right],$

[0214] where C(x) and S(x) denote the cos and sin functions respectively.

[0215] Given the 3-D character of the coronary artery tree, it is expected that any projection will foreshorten a variety of segments. The visual perception provided by such a projection shows a variety of arterial segments with a degree of foreshortening that may or may not be appreciated. In FIG. 41a, the graphical example shows that the projected object has 0% degree of foreshortening since the tangential direction of such a line segment is perpendicular to the normal vector of image plane. Note that the projections of the two line segments have various magnifications due to the different locations relative to the image intensifier. If the two line segments are rotated by an arbitrary angle, e.g., 45 degrees relative to their original orientation in FIG. 41a, there will be approximately 30% degree of foreshortening as shown in FIG. 41b. In FIG. 41c, a synthetic example of vessel foreshortening for a typical segment of a curvilinear arterial centerline is illustrated.

[0216] Due to the problem of vessel overlap and foreshortening, multiple projections are necessary to adequately evaluate the coronary arterial tree with arteriography. The elimination or at least minimization of foreshortening and overlap is a prerequisite for an accurate quantitative coronary analysis (QCA) such as determination of intra coronary lengths in a 2-D display. Hence, the patient might receive substantial radiation and contrast material during diagnostic and interventional procedures. This traditional trial-and-error method provides views in which overlapping and foreshortening are minimized, but only in terms of the subjective experience-based judgment of the angiographer. In FIG. 42aand 42b, the same vessel segment between two bifurcation points as marked by two dots at its two distal ends depicts 77% and 52% foreshortening in the respective angiograms.

[0217] The reconstructed 3-D coronary arterial tree can be rotated to any selected viewing angle yielding multiple computer-generated projections to determine for each patient which standard views are useful and which are of no clinical value due to excessive overlap. Therefore, the 3-D computer assistance provides means to improve the quality and utility of the images subsequently acquired.

[0218] G.2. 1. Calculation of Percent Foreshortening Map

[0219] Element 3714 is operable to present foreshortening information for an arterial tree to permit a user to optimize the view for foreshortening. FIG. 52 is a flowchart providing additional details of the processing to optimize the view of an arterial tree 3-D projection for reducing foreshortening. Let p_i, i=0, 1, . . . , m denote the points on the 3-D vessel centerline curve {haeck over (L)}. Let {right arrow over (l)}_j=[l_j, l_j, l_j]^t, j=1, 2, . . . , m denote the vector formed by every two consecutive vessel centerline points p_j−1and p_j. The percent foreshortening of a vessel segment {haeck over (L)} based on any specific gantry orientation α and β angles can be obtained as follows: ${ℱ}_{\mathrm{fores}}\left(\alpha ,\beta ,\stackrel{⋓}{L}\right)=\frac{\sum _{j=1}^m\left|{\overrightarrow{l}}_j\right|\cos \left({\theta }_j\right)}{\sum _{j=1}^m\left|{\overrightarrow{l}}_j\right|},0°\le {\theta }_j\le 180°,-120°<\alpha <120°,-60°<\beta <60°.$

[0220] where θ is the angle between the directional vector {right arrow over (l)}_jand projection vector {right arrow over (z)}_p=[−cos(α)sin(β), −sin(α),cos(α)cos(β)]^t, and the constraints associated with the angle parameters α and β pertinent to the achievable gantry angles of imaging system. If there are k multiple arterial segments to be evaluated (e.g., a bifurcation region), the aforementioned equation is extended as follows: ${ℱ}_{\mathrm{fores}}\left(\alpha ,\beta ,\stackrel{⋓}{L}\right)=\frac{\sum _{i=1}^k\sum _{j=1}^{m_i}\left|{\overrightarrow{l}}_j^i\right|\cos \left({\theta }_j\right)}{\sum _{i=1}^k\sum _{j=1}^{m_i}\left|{\overrightarrow{l}}_j^i\right|},0°\le {\theta }_j^i\le 180°,-120°<\alpha <120°,-60°<\beta <60°.$

[0221] where {haeck over (L)}={L₁, L₂, . . . , L_k} denote the set of arterial segments and each arterial segment of interest L_i={{right arrow over (l)}₁ⁱ, {right arrow over (l)}₂ⁱ, . . . , {right arrow over (l)}ⁱ_m} consists of (m_i+1) centerline points, i=1, 2, . . . , k.

[0222] The above equations can be employed to model the majority of existing imaging systems having only two degrees of freedom. Another variable γ, denoting the rigid rotation with respect to z-axis by γ degrees, can be readily added to emulate the imaging system if a rotational C-arm base is employed. With the new third degree of freedom in term of the new angle variable, the above equation can be re-written as follows. ${𝔉}_{\mathrm{fores}}\left(\alpha ,\beta ,\gamma ,\breve{L}\right)=\frac{\sum _{j=1}^m\left|{\stackrel{->}{1}}_j\right|\cos \left({\theta }_j\right)}{\sum _{j=1}^m\left|{\stackrel{->}{1}}_j\right|},0°\le {\theta }_j\le 180°,-120°<\alpha <120°,-60°<\beta <60°,-90°<\gamma <90°.{𝔉}_{\mathrm{fores}}\left(\alpha ,\beta ,\gamma ,\breve{L}\right)=\frac{\sum _{i=1}^k\sum _{j=1}^{m_1}\left|{\stackrel{->}{1}}_j^i\right|\cos \left({\theta }_j^i\right)}{\sum _{i=1}^k\sum _{j=1}^{m_1}\left|{\stackrel{->}{1}}_j^i\right|},0°\le {\theta }_j^i\le 180°,-120°<\alpha <120°,-60°<\beta <60°,-90°<\gamma <90°.$