Abstract:
A method is disclosed to accurately determine spatial parameters (for example, the diameter, cross-sectional area and segment length) of bodily structures like human arteries. This is done by determining the three-dimensional position of these structures from two monoplane X-ray angiographic images obtained at arbitrary but different orientations with respect to the structure of interest. The images may be taken at a different moment, so a monoplane X-ray system can be used. Movement of the structure of interest in between the two X-ray exposures is allowed. The three-dimensional information obtained this way is subsequently used to correct the measurements of above mentioned parameters for commonly encountered error sources.

Description:
BACKGROUND  
         [0001]    The present invention relates generally to image analysis, and more particularly to the analysis of images of bodily structures, where those images are used for spatial measurements of those structures.  
           [0002]    The main arteries, like the aorta, the coronary arteries, the carotid arteries, the renal arteries and the femoral arteries provide the body with blood. Artery stenosis is an abnormal condition characterized by the constriction or narrowing of these vital arteries by a substance known as plaque, that prevents proper circulation of the blood. Plaque is a localized area of arteriosclerosis. Arteriosclerosis is a pathological condition that deposits lipids (fatty substances) and a proliferation of fibrous connective tissue on the inner walls of the arteries. Also anatomical (inherited) restrictions of the arteries exist (for instance because of tortuous vessels) or enlargements of the arteries (aneurysms) that strongly influence the flow of blood or, as in the latter case, could even endanger the life of the patient by spontaneous rupture of the artery concerned.  
           [0003]    Standard treatments for these pathological conditions consist of opening up the stenosed arteries by the use of drugs, or when this proves ineffective, by mechanical means. For example, one may force the opening of the artery by inflation of a little balloon inside the stenosed area (“dotter” or “PTA” procedure) often followed by the placement of a metallic frame (“stent”) inside the artery to prevent it from re-occlusion. Stent-like structures can also be used to treat the dangerous local enlargements of arteries. Important for these treatments is to know the exact position, shape, and extent of the stenoses or aneurysms. This will strongly influence the choice of treatment.  
           [0004]    An angiography-based diagnosis is the current standard for determining the extend of stenosis. Angiography is a special X-ray procedure that takes pictures (“angiograms”) of blood vessels. This diagnostic technique makes use of a radiopaque contrast medium, which is a chemical substance that strongly absorbs passing X-rays. Angiography is the X-ray visualization of the internal anatomy of the heart and/or blood vessels after the introduction of a radiopaque contrast medium into the blood. The contrast medium may be injected into an artery or a vein or introduced in a peripheral artery through a catheter (hollow tube) inserted in the artery. The radiologist carefully threads it into the blood vessel and guides it to the area to be studied, under continuous X-ray vision. When the catheter reaches the site under investigation, X-ray contrast medium is injected through the catheter and makes the artery with all its irregularities and blockages clearly visible.  
           [0005]    The term “angiogram” refers to the radiographic image of a blood vessel produced by angiography. Angiograms have darkened areas that represent open channels in blood vessels caused by the contrast medium blocking the X-rays. Digitizing the resulting images makes it possible to apply image processing techniques to the images. One of these techniques is to perform semi-automatic quantitative measurements of the vascular system such as vessel segment length, diameter, cross-sectional area, and the amount of narrowing of a vessel.  
           [0006]    Nowadays most examinations of vessel morphology are done using angiographic images acquired with a monoplane X-ray system. The term monoplane refers to the fact that such a system can only acquire an X-ray image from one direction at a time. In contrast, biplane X-ray systems can acquire images from two directions simultaneously.  
           [0007]    Performing quantitative measurements based upon a monoplane angiographic image has two important shortfalls that arise from the near absence of information about the three-dimensional position of the vascular structure. The image is in fact a projection that projects all three dimensional structures onto a two-dimensional plane. The two shortfalls are out-of-plane calibration errors and “foreshortening,” each of which will now be described.  
           [0008]    When using a monoplane X-ray image, one must “calibrate” the image in order to make absolute measurements. That is, the relationship between pixel size and real-world size needs to be determined. This can be done by including an object of known size in the recorded image. Problems arise when the object used for calibration is not in the same plane (parallel to the input screen of the image intensifier) as the structure under investigation. This will result in the calibration object being magnified differently than the structure, and hence in an incorrect measurement. In other words, every part of the vessel that is located in a different plane than the calibration object is magnified differently. This error source is called “out-of-plane calibration.” 
           [0009]    Generally, because of its intrinsic shape, the structure under investigation typically will not lie exactly in the image plane. If a structure, e.g., a vessel, has a directional component normal to the image plane, the length of a segment of that vessel, when projected on the image plane, will not equal its real three-dimensional length. This phenomenon, leading to errors in length measurements, is called “foreshortening.” 
           [0010]    To overcome these problems one needs a more accurate, three-dimensional (3-D) representation of the position and shape of the structure under investigation.  
           [0011]    Several methods have been developed that derive three-dimensional information from two digital, two-dimensional images. Stereoscopic digital angiography has been used in the calculation of three-dimensional position and orientation information of vessels (L. E. Fencil et al., Investigative Radiology, December 1987). However, stereoscopic determination of three-dimensional vessel position becomes less accurate if the main direction of the vessel is perpendicular to the direction of the stereoscopic shift. Thus, the reliability of this method in determining three-dimensional vessel structure depends on the orientation of the vessels themselves. This is clearly undesirable.  
           [0012]    In U.S. Pat. No. 4,630,203, Szirtes describes a technique for the three-dimensional localization of linear contours appearing in two stereoscopic images. However, this method also suffers from the limitation that the contour must not lie in the direction of the stereoscopic shift. In addition, a separate calibration step is required in this method to determine the three-dimensional locations of the X-ray sources relative to the imaging plane.  
           [0013]    Several workers have developed methods to derive three-dimensional structure from two radiographic images that are obtained in exactly orthogonal directions. The prerequisite that the images must be obtained in exactly orthogonal directions is a clear drawback of this method. This may be difficult to achieve in practice. In addition, determination of the positions of vessel segments perpendicular to one of the imaging planes is difficult or impossible with these methods.  
           [0014]    To address these problems, a method has been developed that allows calculation of three-dimensional vascular structure from two images obtained at arbitrary orientations (see S. A. MacKay et al., Computers and Biomedical Research, Vol. 15, p. 455, 1982). This method requires a calibration step on an object of known dimensions in the same X-ray system configuration as is used to image the patient. This calibration can be done before or after imaging the patient. This method is also referred to as the “calibrated epipolar geometry” method.  
           [0015]    In U.S. Pat. No. 4,875,165, Fencil et al. describe a method that requires a pair of biplane images as input and no calibration step, based upon the theoretical work of H. C. Longuet-Higgings, Nature, Vol. 293, p. 133, 1981. K. R. Hoffmann, U.S. Pat. No. 5,859,922, describes a method using several biplane images with the use of a calibration object. Both methods, however, have the drawback that a biplane X-ray system (that is, an X-ray system capable of making two recordings from different orientations simultaneously) is needed.  
           [0016]    Some other methods are aimed at reconstruction of the complete 3-D vascular tree from two or more images, but often require biplane X-ray systems. See for instance U.S. Pat. No. 6,047,080 to Chen et al. and U.S. Pat. No. 6,169,917 to Masotti et al.  
         SUMMARY  
         [0017]    Various embodiments of the present invention address problems encountered with some other parameter determination systems. For example, with many prior multi-image systems, no movement of the patient table is allowed between the collection of the two images. However, if a radiologist performing an angiographic examination has to reposition the patient table height or horizontal position between recording of the images in order to get the vessel under investigation in optimal view, measurements using other techniques can be adversely affected. The present method, therefore, allows such movements without compromising the accuracy of the results. Further, in some embodiments of the present invention, the calibration step can be done simultaneously with the acquisition of the image containing the vessel.  
           [0018]    A method according to one embodiment of the present invention allows two independently acquired monoplane angiographic images as input for calculation of the three-dimensional midline position of the three-dimensional vascular structures. This means that the two images can be obtained using a monoplane X-ray system (though a biplane system may also be used) at arbitrary orientations. Furthermore, the images may also be obtained with an intervening movement of the patient table.  
       
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0019]    [0019]FIG. 1 is a schematic diagram of a testing room having various components of an apparatus for practicing the present invention.  
         [0020]    [0020]FIG. 2 is a flowchart showing GUI-operator interaction in one method of collecting data for use according to the present invention.  
         [0021]    [0021]FIG. 3 is a flowchart showing one method of analyzing data according to the present invention.  
         [0022]    [0022]FIG. 4 is a diagram showing spatial features referred to in calculations herein.  
         [0023]    [0023]FIG. 5 is a diagram of the relationship among selected distances discussed herein.  
     
    
     DESCRIPTION  
       [0024]    For the purpose of promoting an understanding of the principles of the present invention, reference will now be made to the embodiment illustrated in the drawings and specific language will be used to describe the same. It will, nevertheless, be understood that no limitation of the scope of the invention is thereby intended; any alterations and further modifications of the described or illustrated embodiments, and any further applications of the principles of the invention as illustrated therein are contemplated as would normally occur to one skilled in the art to which the invention relates.  
         [0025]    Generally, the image capture and analysis system and method illustrated in FIGS.  1 - 5  provides for the capture of two digital images, each of a reference object and a bodily structure, such as an artery to be stenosed. In the example presented, the two images are processed to determine the size, shape, and position of the structure.  
         [0026]    An example arrangement for capturing images for use in one embodiment of the present invention will now be discussed with reference to FIG. 1. System  100  includes a C-arm  110  that supports X-ray tube  120  and image intensifier  130 . Patient table  140  supports the patient whose vessel  160  is under investigation in this example. Calibration object  150  is placed so that it appears in the captured images, and will be discussed in further detail below. A normal cath-lab imaging procedure can be followed, as will be understood by those skilled in the art, with the exceptions discussed in the following paragraph, to assure that the required input for the three-dimensional reconstruction is available.  
         [0027]    First, calibration object  150  is placed in the imaging scene. For example, calibration object  150  may be placed beneath the patient, on or under patient table  140 , in a pocket of the patient&#39;s gown, or other location. Then two images are taken of the vessel that is being investigated. The spatial angle between the two normal vectors for the imaging is preferably at least 30°, and the angulation angle γ, rotation angle θ, and source-to-intensifier distance SID should be known. Again, the vessel  160  under investigation and at least a portion of calibration object  150  must both be visible in both images. Between the recording of the two images, the customary patient table movement is allowed. The only requisite is that the vessel position relative to the calibration object does not change. In some embodiments, additional images may be captured, then a pair of images be selected for further analysis.  
         [0028]    Before the three-dimensional reconstruction can be used in practice, a few parameters should be determined. In most embodiments, this only needs to be done once and will in general take place at the installation of the software and hardware. In this example, the following parameters will be determined at installation: geometric image distortion, pixel size at the surface of the image intensifier, and systematic error in the X-ray system geometry.  
         [0029]    The geometric image distortion will be determined for a number of different orientations of C-arm  110 . Geometric image distortion is a distortion resulting from imperfect imaging by the X-ray imaging chain (image intensifier, video camera digitizer) and may be a non-negligible contribution to image errors. One typical distortion factor inherent to the use of image intensifiers is the so-called “pincushion” distortion. Geometric distortion is also dependent on the orientation of the X-ray system with respect to the earth because of the influence of the earth&#39;s magnetic field and gravity. In this embodiment, therefore, the geometric distortion is determined for different orientations of C-arm  110  by making a recording of an image grid of known dimensions (not shown, but well known in the art). Correction parameters are automatically derived, as is within the ability of those skilled in the art. See, for example, Gronenschild, E., J. Janssen, et al. “The accuracy and reproducibility of a global method to correct for geometric image distortion in the x-ray imaging chain.” Med. Phys. 24(12), pp. 1875-88 (1997).  
         [0030]    The pixel size at the surface of the image intensifier will be determined, as is within the ability of those skilled in the art.  
         [0031]    The systematic error in the X-ray system geometry as reported by this system needs to be determined to make correction possible.  
       Validation Data  
       [0032]    To determine the accuracy under different conditions of the method described herein, an internal validation study was performed. This consisted of three parts, which tested length measurement on simulated images, length measurement on phantom images, and diameter measurement on phantom images, respectively. During each test, table movement between the acquisition of the two required images was allowed, and used.  
       EXAMPLE 1  
       [0033]    The goal of the first test was to determine a minimum angle between the reference and analysis image to achieve a predefined maximum error of 5%, and to determine the theoretical accuracy of the three-dimensional reconstruction. The method described below was applied to find a particular length in simulated images, and the minimum angle between reference and analysis image to achieve the desired accuracy of 5% was found to be 30°. In fact, at this difference in projection angle, the accuracy is found to be better than 1.5%.  
       EXAMPLE 2  
       [0034]    The goal of the second test was to determine the accuracy of length measurements based upon X-ray images of a phantom vessel. The method described below was applied to find a particular length in phantom images. Phantom images were acquired with an angle difference of at least 30°, and the accuracy of the length measurement was found to be better than 2%.  
       EXAMPLE 3  
       [0035]    The goal of the third test was to determine the accuracy of the diameter measurements based upon X-ray images of a phantom vessel. The absolute accuracy of the diameter was better than −0.1 mm, and the precision was better than 0.2 mm based on phantom images that were acquired with an angle difference of at least 30°. This experiment used contour detection based on the CAAS II quantitative analysis software (Gronenschild, E., J. Janssen, et al. “CAAS. II: A second generation system for off-line and on-line quantitative coronary angiography.” Cathet Cardiovasc Diagn 33(1): 61-75 (1994)).  
       Image Analysis  
       [0036]    After the images are captured, an operator in this example embodiment performs an analysis using procedure  200 , which will now be discussed with reference to FIG. 2 and continuing reference to objects shown in FIG. 1, beginning at START point  201 . The operator first selects ( 205 ) a reference image. This image should represent both the part of vessel  160  that is being investigated and at least a part of calibration object  150 . The reference image may be captured with or without contrast liquid; in the latter case there is preferably a guide wire or catheter visible in the part of the vessel  160  under investigation.  
         [0037]    The reference image is then automatically corrected (step  210 ) for geometric distortion, as was mentioned above and is understood by those skilled in the art. The operator then indicates (step  215 ) the midline of the vessel in the reference image by indicating some control points. This step is preferably performed using a graphical user interface to a computer program that displays the reference image, accepts “click” actions by a pointing device to indicate control points, and accepts “drag” actions by a pointing device to move previously indicated control points. When the operator indicates that control point selection is complete, the computer program fits a curve through the indicated points, preferably by using a mathematical parametric curve definition known as splines. Preferably six control points are indicated to define the midline of the vessel.  
         [0038]    The operator performs (step  220 ) a calibration of the pixel size on the sphere in the reference image. This calibration is preferably performed automatically (using edge- and curve-detection techniques known to those in skilled in the art) after the operator indicates a point near the center of the sphere. This calibration method is highly reproducible.  
         [0039]    The operator then selects (step  225 ) a second image, to which this disclosure refers as the analysis image. The analysis image should be acquired using contrast liquid in the part of the vessel  160  that is being investigated and with calibration object  150  visible. The analysis image may be a subtracted image, as will be understood by those skilled in the art. The analysis image is automatically corrected (step  230 ) for geometric distortion.  
         [0040]    The operator then performs (step  235 ) a calibration on the sphere in the analysis image. Like the calibration in step  220 , this is preferably performed automatically after the operator indicates a point near the center of the sphere. If the analysis image is a subtracted image, the calibration should be performed without subtraction.  
         [0041]    The computer program then displays (step  240 ) two lines in the analysis image that indicate the corresponding position of the start and end points of the vascular segment as the operator indicated (step  215 ) in the reference image.  
         [0042]    The operator indicates (step  245 ) the vascular segment of interest by indicating the start and end points of the approximate centerline. When the operator finishes, the computer program automatically detects the borders of the vessel. This detection is preferably based on the CAAS II contour detection algorithm (see Gronenschild, Janssen, et al., supra). When the contour detection is finished, a number of physical parameters are calculated, such as vessel diameter, vessel segment length, cross sectional area, and the like. To this point these parameters are still based on the two 2-dimensional images only, and are likely to include out-of-plane and/or foreshortening errors. Results will be corrected for these errors in step  260 , below.  
         [0043]    The midline of the vessel  160  is computed (step  250 ) based on the borders detected in step  245 . A three-dimensional representation of the midline of vessel  160  is then computed (step  255 ), as will be discussed in conjunction with method  300 , below. The geometric results of step  245  are then corrected (step  260 ) according to the three-dimensional spatial position of the vessel, and method  200  ends at END point  299 .  
       Reconstruction Algorithm  
       [0044]    The three-dimensional reconstruction technique applied in step  255 , above, will now be discussed with reference to FIG. 3, and continuing reference to objects and step in FIGS. 1 and 2. The method  300  is based on three assumptions. First, it is assumed that different orientations of the image recording plane are used for the reference and analysis (2-D) images. The angle between the perspective vectors (the normal vectors to the image recording plane as the image is captured) for the two images is preferably between 30° and 150°. Second, at least part of calibration object  150 —enough to determine its center and size in the image plane—and the vessel  160  under investigation must be visible in both images. Finally, the distance SID between the X-ray source and the image intensifier, the rotation angle θ, the angulation angle γ, and the pixel size at the image intensifier should be known, as discussed above in connection with the initialization process.  
         [0045]    At the beginning of method  300  (START point  301 ), the two-dimensional midline from the reference image (indicated in step  215 ) is transformed (step  310 ) into the coordinate system for the three-dimensional world using a rotation matrix R(−γ, −θ), scaling matrix S(p x , p y ), and translation matrix T(−x calibration , −y calibration ). The transformation matrix from reference image to world coordinates is thus:  
           M   I→W     —     ref   =R (−γ,−θ)· S ( p   x   ,p   y )· T (− x   calibration   ,−y   calibration )  [Equation 1] 
         [0046]    which converts the two-dimensional coordinates of the reference image into the three-dimensional world coordinate system (see FIG. 1).  
         [0047]    Likewise, the two-dimensional midline from the analysis image is then transformed (step  320 ) into the three-dimensional world coordinate system. This transformation matrix is defined as  
           M   I→W     —     ana   =R (−γ,−θ)· S ( p   x   ,p   y )· T (− x   calibration   ,−y   calibration )  [Equation 2] 
         [0048]    which converts the two-dimensional coordinates of the analysis image into the three-dimensional world coordinate system.  
         [0049]    For each point in the analysis world contour-wire W a , starting (step  330 ) with the first point V 1   W     a   , the corresponding point in three-dimensional space will now be determined with reference to the points, lines, and planes illustrated in FIG. 4.  
         [0050]    Assume that the contour-wire W a  of the analysis image consists of n elements, and the contour-wire W r  of the reference image consists of m elements. The contour-wires W a  and W r  can now be computed in world coordinates using Equation 1 and Equation 2. For the reference image, the position of each point V j   W     r    can be calculated using  
           V   j   W     r     =M   I→W     —     ref   ·V   j   I     r     ;j={ 1, . . . , m}   [Equation 3] 
         [0051]    and for the analysis image, each point V i   W     a    can be calculated using  
           V   i   W     a     =M   I→W     —     ana   ·V   i   I     a     ;i={ 1, . . . , n}   [Equation 4] 
         [0052]    The real three-dimensional position of the midline wire W is the intersection of the normal vector at W a  with the normal vector at W r  for each segment of wire W.  
         [0053]    The plane Plane 1  perpendicular to the analysis world is defined (step  340 ) by  
         Plane 1   =P   ana     j     +μ·R   reference [2] T   +λ·R   analysed [2] T   [Equation 5] 
         [0054]    ∀μ, λ∈         . Because the X-ray system can be considered a perspective viewing system, the normal vectors should be corrected. The corrected normal is the normalized vector from X-ray source (P 1 ) to the center of the image intensifier (P 2 ). Thus,  
                     P   1     3      D       =            (         0           0             -   SID           )                   P   2     3      D       =            (             (         P   2     2      D          x     -     Po   x       )     ·     p   ii                   (         P   2     2      D          y     -     Po   y       )     ·     p   ii               0         )                     n   -&gt;     corr     =                  (       R     ref   r       ·     R     ref   a         )       -   1       ·     P   2     3      D         -         (       R     ref   r       ·     R     ref   a         )       -   1       ·     P   1     3      D                           n   _     corr     =                n   -&gt;     corr              n   -&gt;     corr                          [     Equation                 6     ]                               
 
         [0055]    The origin of the system can, of course, be shifted out of the origin of the world view, which is the center of calibration object  150 . This is done by correction of the position vector of Equation 5. P ana     j    should be corrected with the correction (offset) vector. The correction vector is the vector from the center of the image to the center of calibration object  150 , rotated over the rotation vector R. That is:  
               V   corr     =              (             (         P   sphere2     2      D          x     -     Po   x       )     ·     p   sphere                   (         P   sphere2     2      D          y     -     Po   y       )     ·     p   sphere               0         )     ·   R             [     Equation                 7     ]                               
 
         [0056]    Equation 5 can then be rewritten as:  
         Plane 1   =V   corr     1     ·P   ana     j     +μ·R   reference [2] T   +λ·R   analysis [2] T   [Equation 5.1] 
         [0057]    Now find (step  350 ) the intersection point between Plane 1  and the three-dimensional world reference midline:  
           P   ref     j   ={( l   refcontour ∩Plane 1   j )}  [Equation 8] 
         [0058]    The plane perpendicular to the analysis world is defined (step  360 ) by:  
         Plane 2   =V   corr     2     ·P   ref     j   +α·( R   reference [2] T   {circle over (x)}R   analysed [2] T )+β·R reference [2] T   [Equation 9] 
         [0059]    ∀α, β∈         .  
         [0060]    Define (step  370 ) the line l perpendicular to the reference plane at the position found in step  360  for ∀α∈         :  
           l=P   ana     j     +α·R   analysed [2] T   [Equation 10] 
         [0061]    The intersection of Plane 2  from step  360  with line l from step  370  is calculated ( 380 ) as the three-dimensional coordinate P 3D :  
           P   3D ={( l ∩Plane 2   j )}  [Equation 11] 
         [0062]    Now that the three-dimensional reconstructed midline of the vessel is available, it is possible to correct the length and diameter for the out-of-plane calibration and foreshortening errors.  
         [0063]    If there are more points to handle in the analysis image (decision step  390 ), the next point in the analysis world contour-wire is retrieved (step  395 ), and process  300  continues at step  340 . When there are no more points (negative result at decision step  390 ), method  300  ends at END point  399 .  
         [0064]    The projection of the three-dimensional midline resulting from the three-dimensional reconstruction process is oriented in the world Cartesian coordinate system. Before performing correction based on the three-dimensional midline of the vessel  160  under investigation, it is preferably positioned in the same view as the acquired image. The three-dimensional model, therefore, is preferably rotated to the same orientation as the projection of the analysis image. Spatial parameters of the vessel  160  are then calculated.  
         [0065]    The length L of the real-world contour-wire (representing the center line of the vessel  160 ) can be calculated by summing the lengths of the segments of the three-dimensional midline over the desired segment. For example, when the points on the contour-wire are denoted as P k  for k=1, . . . m, and the desired segment is from P a  to P b  for some integers a and b, where 1≦a&lt;b≦m, then the length is:  
             L   =       ∑     k   =   a       b   -   1                     P   _     k     -       P   _       k   +   1                        [     Equation                 12     ]                               
 
         [0066]    In addition, one can calculate a diameter of the vessel  160 . After repositioning the three-dimensional midline into the projection of the analysis image, the z-component indicates the out-of-plane distance as shown in FIG. 5. SID represents the source-to-intensifier distance, SOD is the source-to-calibration object distance, P o  is the actual profile of the object, and P ii  is the projected profile of the object.  
         [0067]    Now the pixel size at each position of the midline is computed. Based on the detected contour we can compute the diameter at each position along the vessel. Summarized, the calculation is:  
               SOD   =     SID   ·     (       p   calibration       p   ii       )              
            SOD   i     =     SOD   +            P     3        D   i         ·              0           0           1                              
            p   i     =       p   ii     ·       SOD   i     SID              
            D   i     =              Edge   left     -     Edge   right            ·     p   i                 [     Equation                 13     ]                               
 
         [0068]    Both length and diameter parameters are then available without out-of-plane or foreshortening errors.  
         [0069]    All publications, prior applications, and other documents cited herein are hereby incorporated by reference in their entirety as if each had been individually incorporated by reference and fully set forth.  
         [0070]    While the invention has been illustrated and described in detail in the drawings and foregoing description, the same is to be considered as illustrative and not restrictive in character, it being understood that only the preferred embodiments have been shown and described and that all changes and modifications that would occur to one skilled in the relevant art are desired to be protected.