Estimating the 3D skeletons and transverse areas of coronary arteries from biplane angiograms.

-A method for estimating the three-dimensional (3-D) skel- etons and transverse areas of the lumens of coronary arteries from digital X-ray angiogram is described. The method is based on the use of a 3-D generalized cylinder (GC). The GC consists of a series of 3-D elliptical disks transverse to and centered on a 3-D skeleton (medial axis) of the coronary arteries. The estimates of the transverse areas are based on a nonlinear least-squares-error estimation technique de- scribed by Marquardt. This method exploits densitometric profiles, boundary estimates, and the orientation of the arterial skeleton in 3- space. Our method includes an automatic artery tracking procedure. This procedure applies an adaptive window to the densitometric profile data that is used in the parameter estimation. The tracking procedure in conjunction with Marquardt’s technique yields the areas of sections of lumens of coronary arteries intercepted by planes containing the target focal point and the video scan lines. By applying the method to two views, as in biplane angiography, the 3-D skeletons and the transverse areas of these lumens are estimated. This procedure requires little human intervention. Preliminary experimental tests of our procedure on angiogram of in vivo human coronaries and on synthetic images yield encouraging results. The angiogram were obtained by intra-arterial injection of a contrast agent.


I. INTRODUCTION LTHOUGH the incidence of coronary artery disease
A is decreasing in the United States, as many as 1.5 million Americans had a heart attack in 1985 and about 550 000 of those died. To counteract this epidemic, there have been significant improvements in medical and surgical therapy. The technique of percutaneous transluminal coronary angioplasty is developing as an alternative to coronary bypass surgery in selected patients. In addition, it may soon be feasible to use laser energy systems to ablate atherosclerotic lesions. As these modalities become more sophisticated, it becomes increasingly important to apply accurate quantitative analysis to coronary stenoses in order to assess the results of these clinical interventions.
In this paper, we present a method for estimating the three-dimensional (3-D) skeletons (medial axes) and transverse areas of the lumens of coronary arteries from digitally acquired X-ray angiograms obtained by intra-ar-Manuscript received January 14, 1986; revised May 17, 1988. This work was supported in part by the Focused Research Program on Image Engineering at the University of California, Irvine.
K. Kitamura and J. Sklansky  terial injection of a contrast agent. This method also provides estimates of the boundaries of the images of the arterial lumens in the angiograms. Consequently, this method is a step toward three-dimensional geometric modeling of coronary arteries. In clinical applications our method will enhance the visibility of stenotic lesions and will facilitate the assessment of atheroma mass before and after balloon or laser angioplasty .
Our method has several benefits: 1) it introduces an automatic tracking procedure that suppresses both random noise and line-like structures in the nonarterial background; 2) it provides a robust estimate of the 3-D skeleton (medial axis) of the coronary arteries; and 3) it provides estimates of the cross sections of arteries in planes perpendicular to the local 3-D orientation of the skeleton.
Our method is based on several assumptions, which seem to be supported by experiments. Some of these experiments are Freudenberg's [38]. The remaining experiments are ours, reported in Section V of this paper. The assumptions are: 1) the central projections produced by the X-rays may be approximated by parallel-ray projections perpendicular to the image planes; 2) the extended focal spot may be approximated by a point source of Xrays; 3) the blurring of the image may be approximately represented by a Gaussian point spread function; 4) the attenuation coefficients for soft tissue, bone, and iodineenhanced blood are approximately constant for all wavelengths of the incident X-rays; and 5) the cross sections of the arteries, including stenotic portions, are approximately circular or elliptic.
Several earlier methods for estimating the boundaries of X-ray images of the lumens of coronary arteries have been reported. (In the remainder of this paper we usually use "artery" in place of "arterial lumen" for simplicity of expression.) Selzer [l], [2] developed a method of detecting these boundaries based on the first-and secondorder derivatives of the densitometric profile perpendicular to the medial axis ("skeleton") [3]- [5] of the artery image. Similar boundary detection methods were described by Spears [6], Fukui [7], and Kooijman [8]. Because differentiation amplifies the deleterious effects of noise, the profile is smoothed by polynomial curve fitting before the derivatives are estimated. Unfortunately, when the artery is narrowed, as at a stenosis, or when the densitometric profile of the lumen is nonunimodal, polynomial curve fitting of the lumen boundary may not be sufficiently accurate.
0278-0062/88/0900-0173$01 .OO @ 1988 IEEE Parametric modeling offers an approach to improved estimates of the boundaries [9], [ 101. In this approach one first forms a parametric model of the profile across the artery based on a priori information about its shape. Then each parameter involved in the model is trained recursively based on the observed profile. Shmueli [9] applied a circular cross section model of the artery to digital subtraction angiography images. Pappas [lo] used a more general model with an elliptic cross section and accounted for the background of soft tissue and bones. Both of these parametric methods make more effective use of the profile data than earlier methods. In particular, they exploit the global characteristics of the profile data involved, whereas earlier methods make use of only local changes in the profile data. As a result, the estimated boundaries are more accurate even for the stenosis portion in a noisy image.
Stansfield [ 113 proposed a rule-based expert system (a form of artificial intelligence) to identify coronary arteries in digital subtracted angiograms. This system was demonstrated to be effective for recognizing arteries in a noisy image.
We used the above body of techniques as a point of departure for estimating the shapes and sizes of coronary arterial cross sections from one or two images. A major impediment toward developing an effective estimation technique is our understanding of the likely shapes of the lumen of the artery.
A great majority of sclerotic arteries have elliptic cross sections [12], [13], [37], [38]. With an elliptic cross section the apparent diameter of an artery in a single image depends on the X-ray projection angle. Consequently, under the assumption of an elliptic cross section, an estimate of the diameter of an artery at a point P on its medial axis must take into account the shape of its cross section and the orientation of its medial axis at P in three-dimensional space. From a clinical standpoint, estimates of the transverse areas are usually more significant than estimates of the boundaries because the transverse areas are more directly related to the ability of the arteries to transport blood in sufficient quantities and speeds.
High-speed computerized tomography may avoid many of the problems inherent in monoplane conventional angiography because it provides cross-sectional scans of coronary arteries. Improved quantitative analyses of coronary arteries of dogs have been presented based on this approach [14], [15]. However, this approach is computationally very expensive because it processes large amounts of projection data and requires complex computations such as a filtered projection [16].
Recently, important advances in coronary angiography have been brought about with the help of medical imaging techniques such as digital subtraction angiography (DSA) [ 171-[ 191 and biplane angiography [20]- [22]. In these advances, more effective use of angulated views may enhance our ability to reconstruct coronary arteries in 3space. An idea of three-dimensional reconstruction of coronary arteries from a few angiographic views was proposed by Bresler and Macovski [23]. They introduced a statistical volume model of coronary arteries, assuming, as we do, an elliptical cross section. They used a form of Kalman filtering to suppress the background noise. They tested their method on a software phantom and simulated noise. However, the method was not examined on a real artery image. Recently, Parker [24] reported the use of tracking the leading edge of the bolus and medial-axis angle-corrected densitometry to obtain the entire vascular bed assuming circular cross sections. For sclerotic arteries, however, the restriction to circular cross sections is usually not correct [37], [38].
Our method of skeleton and transverse area estimation [25] is in part an extension of a procedure used by N. Kehtarnavaz in cardiac positron emission tomography [28]. In that procedure the Marquardt nonlinear leastsquare-error estimator [26], [27] plays a key role, and showed good accuracy and convergence rate.
Most earlier methods require substantial human intervention for artery tracking. Automatic tracking of the artery is difficult because the artery usually curves a great deal, becomes narrower at stenoses, and furcates into smaller branches. Our method includes an automatic tracking procedure which applies an adaptive window to the profile data to be used in the parameter estimation. This procedure is designed to suppress subtraction artifacts and edge-like structured noise in the background of the artery images, and contributes to the accuracy of the skeleton and area estimates.
We tested our method on real coronary angiograms as well as on synthetic images with encouraging results.
The remainder of this paper is organized as follows. The 3-D model of coronary arteries is presented in Section 11. In Section 111, our method for estimating the Xray intercepted areas and boundaries of coronary arteries from a single view is described. In Section IV, our method for estimating the 3-D skeletons and transverse areas of coronary arteries from two views is presented. In Section V the results of preliminary experiments to test our method are presented. Our conclusions are presented in Section VI.
11. THREE-DIMENSIONAL MODEL OF CORONARY ARTERIES There are several methods for quantitative representation of 3-D structures. We suggest that coronary arteries may be represented efficiently by a generalized cylinder (GC) [29], [30] because the volume of an artery is naturally described as the swept volume of a 2-D disk (not necessarily circular) moved along a 3-D curve. A GC is a solid having a 3-D curve as a medial axis. At every point on the medial axis, there is a closed and bounded cross section perpendicular to it. The GC model was introduced by Binford [31] as a way of representing 3-D shapes in a computer program.
Brown [12], [13], [37] and Freudenberg [38] have observed that cross sections of coronary arteries are usually elliptic or circular: healthy arteries have circular cross sections, sclerotic arteries usually have elliptic or circular cross sections. Accordingly, our GC model of coronary arteries is restricted to elliptic cross sections, viewing the circle as a form of ellipse. For the sake of geometrical and mathematical simplicity, we approximate the central projection of X-ray imaging by parallel-ray projection orthogonal to the image plane. This is a good approximation in coronary angiograms because the distance between the X-ray source and the object (human heart) is much longer than that between the image plane and the object, and also much longer than the object size. In subsequent refinements of our technique, the parallel-ray projections will be replaced by central projections.
With the use of this GC model, the lumens of coronary arteries can be represented by a 3-D skeleton and a series of elliptical disks centered on the skeleton, as shown in Fig. 1. In this figure the direction of the X-ray beam is assumed to be horizontal. Fig. 1 illustrates that the plane of the cross section does not necessarily coincide with the X-ray beam. In such cases the cross section does not correspond to a gray-level profile along any straight line in the X-ray image of the artery. This fact makes it difficult to use the profile data in the artery image to estimate the parameters of the elliptic cross section.
To overcome this difficulty, we approximate the GC model by a series of truncated elliptic cylinders and/or elliptic cones, as shown in Fig. 2(a). Note that the cross section of an elliptic cylinder or an elliptic cone formed by an arbitrary plane is elliptic. Consequently, we assume that the cross sections of the coronary artery formed by the planes parallel to the X-ray beams are elliptic. Suppose we cut the component cylinders and cones by horizontal planes parallel to the X-ray beams. Then each section intercepted by each of these planes is elliptic, provided the plane passes through just one cylinder or cone and provided the axis of the cylinder or cone is not horizontal. This is shown in Fig. 2(a). In Fig. 2(a), each hatched elliptic cross section is produced by a horizontal plane parallel to the X-ray projection. These intercepted sections correspond to the densitometric profiles along horizontal video scan lines in the X-ray image of the artery. Under the assumption that these intercepted sections are elliptic, we can estimate the parameters of these ellipses based on the corresponding horizontal profile lines in the artery image.
Based on the above observations we model the coronary arteries by a 3-D skeleton and a series of elliptic sections intercepted by the X-ray beams, as shown in Fig. 2 The most convenient plane of these intercepted ellipses is the plane determined by a video scan line (which in this paper is assumed to be horizontal). This orientation of the plane, however, will not be feasible when the image of the local portion of the skeleton containing the intercepted ellipse is also horizontal. In such cases, we must compute profile data along lines perpendicular or nearly perpendicular to the image of the local portion of the skeleton. The details of this technique of estimating intercepted areas are described in Section 111. In Section 11, a 3-D model of a coronary artery is formed based on a 3-D GC. The model assumes elliptic cross sections perpendicular to an artery image plane. Under this assumption, we devised a method which can estimate the X-ray intercepted areas, 2-D skeletons, and boundaries of a coronary artery simultaneously from a single view of the artery image. In this section, the method is described in detail.
The method consists of an application of Marquardt's nonlinear least-square-error estimation and an adaptivewindow tracker of the image of the artery.
Experimental results for this method are presented in Section V.

A. Densitometric ProJle Model of Coronary Arteries
In X-ray images, an incident X-ray beam of intensity IO passing through a length z of material results in a transmitted beam of intensity (1) where U is an attenuation coefficient determined by the material 1321, and the z-axis is parallel to the X-ray beam. When forming angiographic X-ray images of coronary arteries, there are three major types of materials along the X-ray paths: soft tissue, bones, and the iodine-enhanced blood in the arteries. In this case, (1) becomes where U,, U,, and U, are attenuation coefficients for soft tissue, bone, and iodine-enhanced blood, respectively. We assume here that these coefficients are constant within Because of the variations of thickness from one point to the next, we see a meaningful image. The projected image p (x, y ) is given by In our GC model, the thickness L , ( x ) of an artery at x is given by the following equation because the shape of the arterial cross section is elliptic: In This approximation is justified because only a limited region of the profile data along the x-axis is selected for the estimation', as discussed later in this section.
Finally, we assume that the blurring of the X-ray angiogram may be accounted for approximately by convolution by a Gaussian point spread function (PSF). Earlier research indicates that this assumption is correct primarily in portions dominated by low spatial frequencies; in these portions the blurring is caused primarily by scatter. The blurring produced by scatter is well approximated by a Gaussian PSF [34]. Scatter is dependent on the imaged object; hence, this PSF vqries with the spatial coordinates of the image in a manner determined primarily by the imaged object.
In portions of the image dominated by high spatial frequencies, for example at sharp boundaries of high-contrast objects, the blurring is caused primarily by the nonzero extent of the focal spot. The focal spot distribution is usually nonisotropic and may be distinctly non-Gaussian [34]. Since our estimation technique minimizes the mean square error over a window that includes an image of the arterial lumen and some of its neighboring background, and since the boundary of the image of the arterial lumen consumes only a small portion of this window, we believe the assumption of a Gaussian PSF to represent the blurring is reasonable. Both our experiments and those of Pappas [lo] seem to support this belief.
Let g(x) denote this PSF, and let S denote its standard deviation. Then g(x) may be expressed as follows: By convolving the projected image p ( x ) with g ( x ) , we obtain the following densitometric profile estimate f ( x ) in the X-ray intercepted plane: where the asterisk * denotes convolution, v(x) denotes the projected image of the contrast-enhanced artery defined by (9), the artery model, associated with parameters 0, y, and C , , w ( x ) denotes the projected image of the soft tissue and bones defined by (111, the background model, associated with parameters Qo and Q,, and g(x) denotes the PSF defined by (12), the blurring model, associated with parameter S .
Using the densitometric profile estimate f ( x ) , each profile parameter among { 0, y, C , , S , Qo, and Q, } in f ( x ) is estimated based on the comparison to the corresponding observed densitometric profilef (x) in the artery image. The estimation of C , for successive profile lines in the artery image gives us the 2-D skeleton of the artery image (i.e., an estimate of the image of the 3-D skeleton).
The estimation of C, -y and C , + y for successive profile lines in the artery image gives us the boundary of the image of the arterial lumen. The product Py is the maxi- proportional to the maximum projected thickness [cry in Fig. 3(a)] of an X-ray intercepted ellipse of the artery.
The product by2 gives us an estimate of the area of each ellipse.within a scale factor. Thus by2 = 2kUcAB.
Equation (14) shows that Py2 is proportional to the area of the ellipse, because the area of the ellipse is TAB. Thus Py2 = (2kUC/7r) * (area of the ellipse).
Below we describe a training technique for estimating the profile parameters 0 , y, C, , S , Qo, and Q,. In this technique it is assumed that one may calibrate the area.
We suggest two approaches to doing this. Method I : Include an object of known size in the image-often a catheter or a pacemaker may provide this reference image. Method 2: Compute the ratio Z defined by AB Z = cos2 e + p2 sin2 e in a series of successive planes along a relatively large, healthy, clearly imaged artery. Find the average of these ratios as well as the average of the corresponding values of P. Call these two averages ( Z ) and ( P ), respectively. By (lo), This gives us a means of estimating 2kUc. From the discussion following (14), we see that the area can be estimated by the ratio a p y 2 / 2 k U c .

B. Preliminary Estimate of ProBle Parameters
To estimate the parameters P, y, C,, S , Qo, and Q , in f ( x ) , we have used a parameter training technique. There are many such techniques available for this purpose [35]. From them, we chose the Marquardt nonlinear leastsquares estimation technique [ 2 6 ] , [27]. The Marquardt technique minimizes the distance between f ( x ) a n d f ( x ) in the least-squares sense. In Referring to (15) and (16), we define the error vector e as The parameter vector a is defined as The N x 6 Jacobian matrix J is defined as Parameter training by the Marquardt technique is carried out by the recursive vector equation where I is a 6 X 6 identity matrix, K is the iteration number in the training process, and X is a constant that is chosen empirically so as to ensure convergence of the iterative process. (We refer to this as an "iterative parametric training" estimation process.) In our experiments on real coronary angiograms, we found that ak may be kept inside the region of convergence throughout the training process by setting X at the convenient value of 1 . This is because: 1) the initial values of the other parameters may be chosen within the region of convergence by a careful matching to the profile data, and 2) because the final estimated parameters in the present plane are used as the initial parameters for the iterative estimation process in the next plane. Setting X = 1 did not cause any convergence problems in any of our experiments, and achieved a substantial reduc- We observe that E is proportional to the convergence rate I e 12/ I f l2 since f is constant during the iterative estimation process.
The results of using our version of the Marquardt technique for parameter estimation in an experiment are shown in Table I and Fig. 4. The data in this experiment are obtained from a real artery image. Table I  The time consumed in computing these parameters is relatively long: about 1 min for each scan line for the case where the number of pixels in each profile is 20 and the number of iterations is limited to 20. Since all the programs were written in Fortran (a language whose execution is relatively slow), and the computer used in the experiments was old (the Perkin-Elmer 7/32 and the PDP-11/73), there is strong evidence that the computation time will speed up substantially by the use of a more powerful programming language and more powerful computers.
The estimate of the profile parameter obtained by the above procedure is effective on angiograms in which the structured background noise has been suppressed (e.g., by subtraction). However, subtracted angiograms may not be available or may include substantial noise as a result of inexact subtraction. This noise may corrupt the computed transverse areas, typically by producing false oscillations in the curve of area versus distance along the arterial medial axis. To suppress this noise we devised a technique for producing noise-suppressing estimates of the profile diameters. This technique is described below.

C. Noise-Suppressing Estimate of Projle Parameters
In this section we derive a noise-suppressing estimate of the profile parameters. This estimate is based on the interaction among the noise and profile parameters described below.
In an examination of the effect of structured noise on the estimation of the profile parameters, we observed that the standard deviation S and the parameter 0 are sensitive to the noise. The parameter y is relatively insensitive to the noise. Based on these observations, we conjecture that the estimate of 0 is highly correlated to that of S , and that y is largely independent of changes in 0 and S.
From our experiments investigating the behavior of S, 0, and y in response to structured noise, we deduced the following mechanism of interaction among these parameters. Because S, the standard deviation of our Gaussian PSF, represents the blurring caused by the noise, S changes with the noise. (Since the noise is caused primarily by scatter, S is spatially varying and image-dependent.) With a change in S, 0 also changes to optimize the shape of the reconstructed densitometric profile. This is explained schematically in Fig. 5 . In Fig. 5 , we assume data profile (a) with the height of h to be reconstructed by the parameters 0, y, and S. If the standard deviation S is estimated as a larger value S' (b) because of the effect of the structured noise, then 0 is also estimated as a larger value 0' to maintain the same profile height h. Because both 0 and S determine the height of the reconstructed profile and our estimator minimizes the height difference between the data profile and the reconstructed profile at each point, the two parameters 0 and S are closely cor- The profile reconstructed by the parameters p', y, and S'. related. On the other hand, y primarily determines the width of the profile rather than the height. Thus, the effect on y of the change in S is relatively small.
Based on the above analysis, we devised the following noise-suppressing estimate of the profile parameter. The estimate is obtained in two steps.
Step I : For each profile, estimate 0, y, C,, S, Q,, and Q, by the preliminary estimation method described in Section 111-B. Then fix y and C,.
Step 2: Replace S by a running average of S in a neighborhood alcng the local skeleton. Denote this running average by S. Then estimate P, Q,, and Q , , keeping 3, y, and C, fixed. In Step 1, we let all the parameters change as in our preliminary estimation method. Then we fix the values of y and C,, which may be accurately estimated at this point.

In
Step 2, to eliminate the effect of S on 0, the value of S is first averaged in the neighborhood.
is estimated again. T o take the average of S in the neighborhood is reasonable, because the standard deviation of the scatter produced at location P is affected primarily by the portion of the object imaged within a neighborhood of P .

D. Tracking Procedure
Note that the observed densitometric profile f (xi) illustrated in Fig. 4 contains only one local maximum. Unfortunately, the parameter training technique given by (20) does not yield accurate results when the observed densitometric profile has two or more local maxima. Such profiles are often observed in the following three types of segments of an artery image: a junction where the artery furcates into smaller branches (a similar profile may be produced when an edge-like artifact or a bone trabecula is near the image of the artery); a stenosis, i.e., a relatively narrow section of the artery; and a section where the artery is oriented nearly parallel to the cross-sectional plane. Examples of such segments of arteries and their profiles are sketched in Fig. 6. Suppose we are scanning the artery at plane L. Then only the profile data f (x) between A and B specified in Fig. 6 should be used for the parameter estimation in each case. If our parameter training technique is applied to the full profiles shown in Fig.  6, the results will not be accurate.
To overcome this difficulty in the use of the parameter training technique, we have developed an ''adaptive window" that selects the appropriate region of the profile data for estimating the parameters. In addition, we have developed a procedure of automatic tracking along the artery. This tracking procedure facilitates the use of the adaptive window and the parameter estimator in successive scanning planes.
In our automatic tracking procedure, the selection of the appropriate region of the profile data in the next plane is based on the estimates of parameters in the present plane and the shape of the profile data in the next plane. We explain the details of the procedure by referring to Fig.   7. Assume that we are tracking the artery upward and have just finished the parameter estimation in plane L -1. Supposef(x) represents the profile data in plane L. Now we have to determine the region of the profile data f ( x ) to be used for parameter estimation in plane L. At this point, a pair of boundary points ( x u , xb) of the artery in plane L -1 have been obtained from the estimates of the parameters ( x u = c, -y, xb = C , + y ) . The algorithm to determine the adaptive window for the profile dataf(x) in plane L is as follows: Adaptive Window Algorithm: Step 1: Let [xu, xb] be the initial window o n f ( x ) .
Step 2: Find the point x, in the window [xu, xb] such to Step 2. Else go to Step 5 .
Step 5: Starting from x, and moving leftward on the xaxis, find the first point xI satisfyingf(xl) 5 f ( x r + 1 ) Step 6: If f ( x l ) < TH (a fixed threshold value), then go to Step 9. Else go to Step 7.
Step 8: Moving the left from x l , replace x I by the next point x1 satisfying f ( x l ) 5 f ( x l + 1 ) and f ( x r ) < f ( x l -1). Then go back to Step 6.
Step 9: Starting from x, and moving rightward on the x-axis, find the first point x, satisfyingf(x,) f ( x r -1 ) Step 10: Iff(x,) < TH, then go to Step 13. Else go to Step 11.
Step 12: Moving to the right from x,, replace x, by the next point x, satisfying f ( x , ) f ( x , -1 ) andf(x,) < f ( x r + 1). Then go back to Step 10.
Step 13: Set the adaptive window o n f ( x ) to [xl, x,]. The algorithm is illustrated in Fig. 7. In this figure x, is the maximum point in the initial window [ x u , xb], and x, is not on the border of this window. To the left of x,, xg is the first local minimum point andf(x,) 2 T H . Since the next left point xf is a local maximum point, the local minimum point xg is ignored. x, is the next local minimum point. Sincef(x,) < TH, x, becomes the left endpoint of the adaptive window. To the right of x,, xh is the first local minimum point andf(xh) 2 TH. However, since the next right point xi is not a local maximum point, x h becomes the right endpoint of the adaptive window. As a result, we obtain the adaptive window [x,, xh] o n f ( x ) . With this algorithm, the adaptive window selects the appropriate region of the profile data even for the difficult cases shown in Fig. 6.
The rationale justifying our algorithm for the adaptive window is as follows. Consider the angiogram as a superposition of foreground and background. The foreground is the image of the iodinated coronary lumens; the background is the image of the noniodinated tissue-primarily soft tissue (muscle and fat) and bone. The threshold TH establishes a criterion for separating foreground from background except when the background has unusual contrast (as may be caused by line-like subtraction artifacts, e.g., by inadequately subtracted bone trabeculae or by a furcation in the artery [ Fig. 6(a) local maximum, illustrated at point C in Fig. 6(a), obscures the true boundary of the arterial lumen; hence we use the preceding local minimum as an endpoint of the window in the presence of this strong local maximum.
Thus, the first minimum away from x, that is above TH is not made an endpoint of the window unless a strong local maximum is nearby. A weak local maximum one pixel beyond a local minimum above TH is ignored because that weak local maximum most likely represents random noise. In the absence of strong local maxima the first local minimum away from x, is made an endpoint of the window. Based on the adaptive window, automatic tracking along an artery is carried out as follows. Assume that we have an artery image as shown in Fig. 8. For the other branches of the artery, like 0, 0, and @ in Fig. 8, the only necessary manual inputs are the initial values of the profile parameters for the first planes (C, E, and G ). Each branch is automatically tracked up- ward just as for branch 0. In these branches, the tracking can stop automatically when it meets branch 0, which has been already tracked. For simplicity the automatic tracking procedure above is discussed for a special case where the basic arterial tree is vertical (or nearly vertical), and the artery image can be tracked upward by horizontal scanning planes. But horizontal scanning planes are not suitable for arterial segments that are horizontal or nearly so. For such cases our procedure computes a suitable scanning direction. This computation requires the manual marking of several points on the medial axis of the artery image. These points are illustrated in Fig. 9 by the vertices of the piecewise linear approximation of the medial axis. In this tracking procedure, we implemented eight discrete scanning directions (a-@ shown in Fig. 9). Among these eight directions, the procedure selects the direction closest to the perpendicular to each linear segment specified by two successive manually marked points on the medial axis. With the use of this tracking procedure, the problem encountered in Fig. 6(c) is avoided.
Because the scanning direction at each linear segment may be calculated from the manually marked points, it is easy to obtain the densitometric profile data. For each linear segment, the xand y-axes of the artery image are rotated according to the calculated scanning direction at the segment. After the rotation of the coordinate axes, the densitometric profile data are obtained for each scanning plane of the segment.

IV. ESTIMATING THE 3-D SKELETONS AND THE TRANSVERSE AREAS OF CORONARY ARTERIES FROM
Two VIEWS In Section I11 we described our parameter estimation method for a single angiogram. This method yields the Xray intercepted areas, the image of the arterial skeleton, and the boundaries of the image of the arteries. In this Fig. 10. The geometry for reconstructing the 3-D skeleton of an artery from two X-ray views.
section we extend this procedure to compute the 3-D skeleton and transverse areas of the arteries from biplane image data. This extension is important because the 3-D skeleton and the transverse areas are medically more meaningful than the 2-D image of the skeleton and the Xray intercepted areas.
The extended procedure requires two simultaneous images (Image 1 and Image 2), produced by two parallelray projections (Projection 1 and Projection 2), as shown in Fig. 10. In this figure the coordinate axes for Projection 2 are denoted by (x', y , z' ). The angle between Projection 1 and Projection 2 is denoted by 6. The directions of the X-ray beams for both Projection 1 and Projection 2 are assumed to be horizontal.
For each of these two images (Image 1 and Image 2), the parameter estimation method described in Section I11 is applied. This yields the 2-D skeleton for each of the artery images. From the obtained 2-D skeletons, the 3-D skeleton of the artery is reconstructed as follows. In Fig.   10, C, and C,. denote the xand x'-coordinates of the 2-D skeletons in the horizontal plane L , respectively. Because the geometrical relations between the two images are known, the position of the 3-D skeleton in plane L may be obtained by backprojecting C, and Cx, along the X-ray projections, as shown in Fig. 10 [36]. By repeating this procedure at successive horizontal planes, the 3-D skeleton may be reconstructed.
In carrying out the above procedure we encountered the practical problem of identifying corresponding points in the two images of the arterial skeleton. To compute corresponding points in pairs of images, we first identified manually those pairs of corresponding points that are easily identified by the human viewer. We refer to such points as landmarks. We marked two types of landmarks on the two images: junction points and stenotic points. A junction point is a point where the blood flow from two vessels joins to form a larger flow in a larger vessel. Using the reconstructed 3-D skeleton and the computed X-ray intercepted areas of the arteries, we obtain the transverse areas as follows. Let A , be a cross section of the GC model of an artery and a', is its normal 3-D vector.
Let AF represent the X-ray intercepted section of the artery at the corresponding position (the same y-coordinate) along the 3-D skeleton, and ZF is its normal 3-D vector. The vector ZF is known because we know the direction of A F . The vector a'R can be obtained for any point along the 3-D skeleton by calculating the local direction of the reconstructed 3-D skeleton. We estimate the angle 4 between a'F and a'R from (22) a'F a'R cos q!l = -.
In (22), ZF -ZR denotes the inner product of the vectors a F and ZR. I ZF I and I Z R I are the norms of these vectors.
The areas of AF and AR are related by the following approximation: where [ A R ] and [ A F ] stand for the areas of AR and AF, respectively. Finally, from the estimated area of AF using the method discussed in Section I11 and the obtained angle q!l, we estimate the area of AR from Experiments testing the above procedure are described in the next section.
V. EXPERIMENTS Our method for quantitative analysis of coronary arteries described in the preceding sections was tested in several experiments. These experiments are reported in this section. The experiments were performed on the following three types of images: a) subtracted artery images; b) biplane nonsubtracted artery images; and c) synthetic artery images superimposed on background noise in real angiograms. The results are consistent, and lend support to the validity of our method.

A. Subtracted Artery Images
We tested our preliminary estimation method (Section 111-B) on a subtracted X-ray image of a patient's left coronary artery at LAO 60" obtained in the Cardiac Catheterization Laboratory of the University of California Irvine Medical Center. This image is one of the series of cinematic frames of an in vivo heart shortly after the injection of iodinated contrast material into the coronary arterial network. These frames were acquired by an X-ray imaging system consisting of a Siemens Gigantos Optimatic X-ray generator with a tungsten target operated at 64 mA and 96 kVp, a Siemens photointensifierlvideo camera that produces 30 frames/s, and a Fischer Model DA-100 computer system that digitizes each frame into a 512 X 512 X 8-bit array, processes the image, and stores the processed image on 1 /4-in tape in a Scotch DC 600 HC tape cartridge. The width of the focal spot was 0.6 mm. The magnification of the X-ray image projected onto the input screen was approximately 1.3. Fig. 11 shows the subtraction-enhanced image that we tested. This image was obtained by subtraction of a noniodine image from an iodine-enhanced image, followed by the addition of a constant gray level (105 in this case) to each pixel in order to eliminate negative-valued pixels. The tracking order and direction are shown in Fig. 8, which illustrates the structure of the arterial network. In this experiment, the artery is tracked upward by horizontal scanning planes without manual marking of points on the medial axis, because in the image at LAO 60" the left coronary artery is oriented within 45 O of the vertical. The number of scanning planes between A and B is 235. The value of TH (the threshold in the adaptive window algorithm) used here is 130, which is about half of the full range of gray level (from 0 to 255 ), and obtained empirically. Fig. 12(a) shows the estimated boundaries (C, -y, C , + y). The method worked automatically from A to B , from C to D, from E to F , and from G to H, respectively, as discussed in Section 111. The method also worked at a stenosis and at a large junction ( S and T i n Fig. 8). The branch @ the estimated boundary points included a few outlying points. This could happen at a scanning plane where the local direction of the artery is almost horizontal-as at position U in Fig. 8-because the profile shape in such a plane is complicated. These outlying boundary points were deleted manually, and the remaining gaps and points were filled and smoothed by linear interpolation. Fig. 12(b) shows the estimated 2-D skeleton produced by the values of C,. Fig. 12(c) shows the estimate of /'3y centered on the 2-D skeleton. Finally, Fig. 12(d) shows the estimate of Py2 centered on the 2-D skeleton. In Fig. 12(c) and (d), the length of each horizontal black line shows the magnitude of /'3y and or2, respectively. In Fig. 12(d), the computed intercepted area of the artery is quite small at the stenosis S . These results encourage us to believe that our procedure for finding intercepted areas is a good way of finding stenoses in coronary angiograms.

B. Biplane Nonsubtracted Artery Images
Next we tested our method on real biplane images. Four iodine-enhanced biplane images of a right coronary artery are used: two from LAO 30" (labeled as Image L-1 and Image L-2 ), the other two from R A 0 30" (Image R-1 and Image R-2). These four images are taken at the same phase of the heartbeat: at ED (end of diastole). These four images are shown in Fig. 13 (LAO 30") and Fig. 14 ( R A 0 30" ). These images were obtained from the Cardiac Catheterization Laboratory at the University of California Irvine Medical Center using the system described earlier in this section. Because no background images taken before iodine injection were available to us, these images are nonsubtracted. Fig. 15 shows perspective views of the estimated 3-D skeletons from two views using the method described in Section IV. As shown in Fig. 15, the two skeletons are reconstructed from two different combinations of the images: (a) the 3-D skeleton estimated from Image L-1 and Image R-1; and (b) the 3-D skeleton estimated from Image L-2 and Image R-1. Here the skeletons of the artery are reconstructed between the two points 0 and 0' shown in Fig. 13, using the automatic tracking procedure with manually marked points on the medial axis. The estimated skeletons in Fig. 15 are consistent with each other. Based on these reconstructed 3-D skeletons, the transverse areas of the 3-D GC model are estimated between the two points A and B shown in Fig.  13. The results are shown in Fig. 16. Fig. 16(a) shows the estimated transverse areas from the X-ray intercepted areas of Image L-1 and the 3-D skeleton shown in Fig.  15(a). Fig. 16(b) shows the estimated transverse areas from the X-ray intercepted areas of Image L-2 and the 3-D skeleton shown in Fig. 15(b). Fig. 16(c) shows the estimated transverse areas from the X-ray intercepted areas of Image R-1 and the 3-D skeleton shown in Fig.  15(a). Fig. 16(d) shows the estimated transverse areas from the X-ray intercepted areas of Image R-2 and the 3-D skeleton shown in Fig. 15(b). These results are consistent among the four images. This consistency argues for the validity of our estimation method. These estimated transverse areas decrease to a minimum at the stenosis labeled S in Fig. 13(a).

C. Synthetic Artery Images
Next we tested our method on a synthetic artery image superimposed on a portion of the background of a real artery image. Here a software phantom of an elliptic cyl- --  Fig. 13(a) (Image L-1 ) is superimposed on the image of the phantom. The standard deviation of this background was estimated to be 6.89 within a 50 X 100-pixel window placed so that the phantom is contained within the window. This window is marked by dashed lines in Fig. 13(a). In contrast to our lack of knowledge of the true 3-D structure and dimensions of the real artery in Figs. 13 and 14, this phantom image provides a ground truth to evaluate our In order to analyze the effect of the structured noise on the parameter estimation quantitatively, our simplified method for the profile parameter estimation was tested on the synthetic images generated by the procedure described above. Fig. 17 shows the estimates of y and Py of the synthetic images for various values of r: (a) r = 0.0; (b) r = 0.5; and (c) r = 1.0. These three values of r correspond to standard deviations of background noise of CJ = 0, 1.57, and 3.14, respectively. The corresponding signal-to-noise ratios are approximately by/. = 16/y = 00, 10.2, and 5.1, respectively. In Fig. 17, the length of each horizontal black line shows the magnitude of the estimate of y or Py at each cross section. As shown in Fig. 17, the estimate of p is susceptible to the structured noise. When there is no such structured noise [ Fig. 17(a)], both and y are accurately estimated with no errors. However, as the noise parame- ter r increases [ Fig. 17(b) and (c)], the error of the estimate of increases as well. The estimate of y is much more stable.
To evaluate the effectiveness of our noise-suppressing estimate quantitatively, it was applied to the synthetic images. Fig. 18 shows the experimental result. Compared to Fig. 17, the effectiveness of the method is evident. Fig.   18 indicates that the estimate of p is approximately as stable as the estimate of y.
In summary, these experiments support the validity of our method for estimating the 3-D skeletons and transverse areas of coronary arteries from biplane angiograms.

VI. CONCLUSIONS
In this paper we describe a computer-based method that helps cardiologists to visualize, detect, identify, and measure coronary lesions from digitized single-plane and biplane angiograms. In particular: We showed how to obtain estimates of the arterial boundaries and the areas intercepted by planar X-ray beams from a single view.
We showed how to obtain estimates of the three-dimensional skeletons of coronary arteries and the transverse areas locally perpendicular to the skeletons from two views.
Our method for estimating the transverse areas suppresses much of the structured noise in nonsubtracted images.
Our experimental results, using both coronary angiograms and synthetic images, support the validity of our estimation methods.
Consequently, we believe this paper is a step forward in the automation of the detection and measurement of coronary lesions and in quantitative geometric modeling of coronary arteries from just one or two X-ray image planes.
The effectiveness of the technique presented here is limited by 1) the sampling imposed by the video link and A/D converter, 2) the complexity of the arterial network structure and its image, e.g., the existence of overlapping images of vessels, 3) the deviation of the lumen cross section from an ellipse, 4) the contrast of the arterial lumen against its background, 5) the parallel-ray approximation, 6) the point source approximation, 7) the Gaussian point spread function approximating the blurring operation, and 8) the polyenergetic dependency of attenuation of X-rays by human tissue.