Abstract:
An image reconstruction method for cone beam x-ray attenuation data acquired over a super-short-scan, short-scan or full-scan includes backproejcting over three adjacent segments of the arc traversed by the x-ray source. Each backprojection consists of a weighted combination of ID Hilbert filtering of the modified cone-beam data along both horizontal and non-horizontal directions.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS  
       [0001]     This application is based on U.S. Provisional Patent Application Ser. No. 60/630,919 filed on Nov. 24, 2004 and entitled “CONE-BEAM FILTERED BACKPROJECTION IMAGE RECONSTRUCTION METHOD FOR SHORT TRAJECTORIES.” 
     
    
     STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH  
       [0002]     This invention was made with government support under Grant No. CA109992 and EB001683 awarded by the National Institute of Health. The United States Government has certain rights in this invention 
     
    
     BACKGROUND OF THE INVENTION  
       [0003]     The present invention relates to computed tomography (CT) imaging apparatus; and more particularly, to a method for reconstructing images from divergent beams of acquired image data.  
         [0004]     In a current computed tomography system, an x-ray source projects a fan-shaped beam which is collimated to lie within an X-Y plane of a Cartesian coordinate system, termed the “imaging plane.” The x-ray beam passes through the object being imaged, such as a medical patient, and impinges upon an array of radiation detectors. The intensity of the transmitted radiation is dependent upon the attenuation of the x-ray beam by the object and each detector produces a separate electrical signal that is a measurement of the beam attenuation. The attenuation measurements from all the detectors are acquired separately to produce the transmission profile.  
         [0005]     The source and detector array in a conventional CT system are rotated on a gantry within the imaging plane and around the object so that the angle at which the x-ray beam intersects the object constantly changes. A group of x-ray attenuation measurements from the detector array at a given angle is referred to as a “view” and a “scan” of the object comprises a set of views made at different angular orientations during one revolution of the x-ray source and detector. In a 2D scan, data is processed to construct an image that corresponds to a two dimensional slice taken through the object. The prevailing method for reconstructing an image from 2D data is referred to in the art as the filtered backprojection technique. This process converts the attenuation measurements from a scan into integers called “CT numbers” or “Hounsfield units”, which are used to control the brightness of a corresponding pixel on a display.  
         [0006]     The term “generation” is used in CT to describe successively commercially available types of CT systems utilizing different modes of scanning motion and x-ray detection. More specifically, each generation is characterized by a particular geometry of scanning motion, scanning time, shape of the x-ray beam, and detector system.  
         [0007]     As shown in  FIG. 1 , the first generation utilized a single pencil x-ray beam and a single scintillation crystal-photomultiplier tube detector for each tomographic slice. After a single linear motion or traversal of the x-ray tube and detector, during which time 160 separate x-ray attenuation or detector readings are typically taken, the x-ray tube and detector are rotated through 1° and another linear scan is performed to acquire another view. This is repeated typically to acquire 180 views.  
         [0008]     A second generation of devices developed to shorten the scanning times by gathering data more quickly is shown in  FIG. 2 . In these units a modified fan beam in which anywhere from three to 52 individual collimated x-ray beams and an equal number of detectors are used. Individual beams resemble the single beam of a first generation scanner. However, a collection of from three to 52 of these beams contiguous to one another allows multiple adjacent cores of tissue to be examined simultaneously. The configuration of these contiguous cores of tissue resembles a fan, with the thickness of the fan material determined by the collimation of the beam and in turn determining the slice thickness. Because of the angular difference of each beam relative to the others, several different angular views through the body slice are being examined simultaneously. Superimposed on this is a linear translation or scan of the x-ray tube and detectors through the body slice. Thus, at the end of a single translational scan, during which time 160 readings may be made by each detector, the total number of readings obtained is equal to the number of detectors times 160. The increment of angular rotation between views can be significantly larger than with a first generation unit, up to as much as 36°. Thus, the number of distinct rotations of the scanning apparatus can be significantly reduced, with a coincidental reduction in scanning time. By gathering more data per translation, fewer translations are needed.  
         [0009]     To obtain even faster scanning times it is necessary to eliminate the complex translational-rotational motion of the first two generations. As shown in  FIG. 3 , third generation scanners therefore use a much wider fan beam. In fact, the angle of the beam may be wide enough to encompass most or all of an entire patient section without the need for a linear translation of the x-ray tube and detectors. As in the first two generations, the detectors, now in the form of a large array, are rigidly aligned relative to the x-ray beam, and there are no translational motions at all. The tube and detector array are synchronously rotated about the patient through an angle of 180-360°. Thus, there is only one type of motion, allowing a much faster scanning time to be achieved. After one rotation, a single tomographic section is obtained.  
         [0010]     Fourth generation scanners feature a wide fan beam similar to the third generation CT system as shown in  FIG. 4 . As before, the x-ray tube rotates through 360° without having to make any translational motion. However, unlike in the other scanners, the detectors are not aligned rigidly relative to the x-ray beam. In this system only the x-ray tube rotates. A large ring of detectors are fixed in an outer circle in the scanning plane. The necessity of rotating only the tube, but not the detectors, allows faster scan time.  
         [0011]     Divergent fan-beam scanning modes have the potential to allow fast data acquisition. But image reconstruction from divergent-beam projections poses a challenge. In particular, the projection-slice theorem is not directly applicable to the divergent-beam projections since the shift-invariance in a single view of projections is lost in the divergent-beam case. One way to bypass this problem is to explicitly rebin the measured divergent-beam projections into parallel beam projections. This is the basic method currently used in solving the problem of fan-beam image reconstruction. After the rebinning process, one can take the advantages of the fast Fourier transforms (FFT) to efficiently reconstruct images.  
         [0012]     Recently, the divergent cone-beam reconstruction problem in x-ray CT has attracted increased attention due to the rapid development of multi-row detectors. In the cone-beam case, it is much more complicated to rebin cone-beam projections into parallel-beam projections. The huge cone-beam data set also poses a big challenge to the potential data storage during the rebinning process. The main stream of the developments in cone-beam reconstruction has been focused on the development of approximate or exact reconstruction methods. For circular-based source trajectories, methods disclosed by L. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical Cone Beam Algorithm,” J. Opt. Soc. Am. A 1, 612-619(1984); G. Wang, T. H. Lin, P. Cheng, and D. M. Shinozaki, “A general cone-beam reconstruction algorithm,” IEEE Trans. Med. Imaging 12, 486-496 (1993); generate acceptable image quality up to moderate cone angles (up to 10° or so). Exact reconstruction algorithms have also been proposed and further developed for both helical source trajectory and more general source trajectories. Most recently, a mathematically exact and shift-invariant FBP reconstruction formula was proposed for the helicauspiral source trajectory A. Katsevich, “Theoretically exact filtered backprojection-type inversion algorithm for spiral CT,” SIAM (Soc. Ind. Appl. Math.) J. Appl. Math. 62, 2012-2026 (2002).  
         [0013]     There are two types of analytic cone-beam reconstruction algorithms: exact and approximate. The first and most practical approximate cone-beam reconstruction algorithm was proposed by Feldkamp, Davis and Kress (FDK) for a circular x-ray source trajectory. It is a heuristic extension of standard fan-beam reconstruction to the cone-beam case by introducing a “cosine” factor. Extensions to this approximate algorithm have been developed. Due to its one-dimensional shift invariant filtering kernel it has been applied to various medical imaging systems. The original FDK algorithm is applied to a complete circular scanning path. Recently, the fan-beam super-short-scan FBP reconstruction formulae has been extended to the cone-beam case to obtain super-short-scan FDK-type algorithms.  
         [0014]     There is an alternative method to obtain an approximate cone-beam reconstruction algorithm; namely, apply a mathematically exact algorithm to an incomplete source trajectory. In the case of a single arc used to reconstruct a 3D volume the data is incomplete according to Tuy&#39;s data sufficiency condition. The problem has been addressed in the Grangeat framework using a shift variant filtered backprojection reconstruction formula. When a rebinning scheme is introduced, a short scan image reconstruction algorithm in the Grageat framework results.  
         [0015]     Recently, a general procedure to generate shift-invariant cone-beam reconstruction algorithms has been developed for a general source trajectory provided that the data sufficiency condition is fulfilled. There is a need for such an image reconstruction method where the data sufficiency condition has not been met.  
       SUMMARY OF THE INVENTION  
       [0016]     The present invention is a new cone-beam FBP reconstruction method in which a shift-invariant FBP algorithm is applied to an arc scanning path. The scanning path may be a full circle (full-scan), an arc satisfying the short scan condition (short-scan), viz, the angular range of the scanning path is π+fan angle, or an arc with a shorter scanning path (super-short-scan), viz., the angular range of the scanning path is less than the above short scan condition. The resulting method is not mathematically exact; however, the shift-invariant feature is preserved. The method includes backprojections from three adjacent segments of the arc defined by T 1 ({right arrow over (x)}), T 2 ({right arrow over (x)}) and T 3 ({right arrow over (x)}). Each backprojection step consists of a weighted combination of ID Hilbert filtering of the modified cone-beam data along both horizontal and non-horizontal directions. The redundant projections were equally weighted to achieve the optimal noise variance. 
     
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0017]      FIGS. 1-4  are schematic representations of different CT system geometries;  
         [0018]      FIG. 5  is a pictorial representation of a 3D, or volume, CT system;  
         [0019]      FIG. 6  is a pictorial view of a CT system which employs the present invention; and  
         [0020]      FIG. 7  is a block diagram of the CT system of  FIG. 6 ;  
         [0021]      FIG. 8  is a pictorial representation of the coordinate system used to derive the equations for the present invention;  
         [0022]      FIG. 9  is another pictorial representation of parameters used to derive the equations for the present invention;  
         [0023]      FIGS. 10   a ,  10   b ,  10   c  and  10   d  are pictorial representations of parameters used to derive equations for the present invention;  
         [0024]      FIGS. 11   a ,  11   b ,  11   c  and  11   d  are pictorial representations of parameters used to derive equations for the present invention;  
         [0025]      FIGS. 12   a  and  12   b  are pictorial representations of the x-ray cone beam and detector array used to acquire data according to the present invention;  
         [0026]      FIG. 13  is a flow chart which illustrates the processing steps in a preferred embodiment of the present invention; and  
         [0027]      FIG. 14  is a side elevation view of a C-arm x-ray system which employs the present invention. 
     
    
     GENERAL DESCRIPTION OF THE INVENTION  
       [0028]     In the following discussion a third-generation cone-beam source-detector configuration is assumed where the detector is constantly opposed to the source. The coordinate system and the parameters used herein are demonstrated graphically in  FIG. 8 . Letter O labels the iso-center of rotation. The coordinate system x-y-z located at O is the laboratory Cartesian coordinate system, and the z axis is the axis about which the detector and source rotate. The detector plane is perpendicular to the iso-ray originating from the source point {right arrow over (y)}(t). The distance along the iso-ray between the source and detector is D. In the laboratory coordinate system, the point {right arrow over (x)} within the object is written as {right arrow over (x)}=(x,y,z). The density of a 3D object to be reconstructed is denoted as f({right arrow over (x)}). The cone-beam projection of the object from the source point {right arrow over (y)}(t) in the direction {circumflex over (r)} is written as: 
 
 g ( {circumflex over (r)},t )=∫ o   ∞   ds f[{right arrow over (y)} ( t )+ s{circumflex over (r)}]   (1) 
 
         [0029]     A unit vector {circumflex over (β)} will be defined in the direction from a source point {right arrow over (y)}(t) to an object point {right arrow over (x)}.  
                 β   ^     =           x   →     -       y   →     ⁡     (   t   )                  x   →     -       y   →     ⁡     (   t   )                =     (       cos   ⁢           ⁢   β   ⁢           ⁢   sin   ⁢           ⁢   α     ,     sin   ⁢           ⁢   β   ⁢           ⁢   sin   ⁢           ⁢   α     ,     cos   ⁢           ⁢   α       )         ,           (   2   )             
 
 where β and α are the longitudinal and azimuthal angles respectively. We next introduce a plane which contains the vector {circumflex over (β)}, and which is characterized by its normal vector {circumflex over (k)} φ . In a local coordinate system this normal vector may be locally parameterized by an angle φ. This plane will be denoted as Π({right arrow over (x)},t;φ) and may be written as: 
 
{circumflex over ( r )}(γ,φ)={circumflex over (β)} cos γ+{circumflex over (β)}×{circumflex over (k)} φ  sin γ, γε[−π,π],  (3) 
 
 where γ is the angle measured from {circumflex over (β)} to {circumflex over (r)} as shown in  FIG. 9 . 
 
 The shift-invariant FBP cone-beam reconstruction formula for a general source trajectory that fulfills the data sufficiency condition (H. K. Tuy, “An inverse formula for cone-beam reconstruction,”  SIAM J. Appl. Math  43, pp. 546-552, 1983.) is given by:  
               f   ⁡     (     x   →     )       =         -     1     4   ⁢     π   2           ⁢       ∑   m     ⁢     ∫       ⅆ   γ     ⁢         c   m     ⁡     (       x   →     ,   t     )                x   →     -       y   →     ⁡     (   t   )                ⁢       ∫     -   π     π     ⁢           ⁢       ⅆ   γ     ⁢     1     sin   ⁢           ⁢   γ       ⁢     ∂     ∂   q       ⁢     g   ⁡     [         r   ^     ⁡     (     γ   ,     ϕ   m       )       ,       y   →     ⁡     (   q   )         ]                   ⁢     ❘     q   =   t       .             (   4   )             
 
 with φ m  ε[0,π), [A. Katsevich, “A general scheme for constructing inversion algorithms for cone-beam CT,”  Int. J Math. Sci.  21, pp. 1305-1321, 2003 and G. H. Chen, “An alternative derivation of Katsevich&#39;s cone-beam reconstruction formula,”  Med. Phys.  30, pp. 3217-3226, 2003.] and where c m ({right arrow over (x)},t) is termed the structure factor and is explicitly dependent on the source trajectory and the selected weighting function. Thus, one may summarize the steps necessary to reconstruct a single point as: 1) Extract the measured projection data from the desired planes Π({right arrow over (x)},t;φ m ) which will be referred to as the critical planes and will be defined below, 2) Differentiate the projection data with respect to the source parameter, 3) Perform a one-dimensional filtering operation, 4) Sum the contributions over the weighted critical planes, 5) Perform a distance weighted backprojection. Thus, in deriving a reconstruction algorithm for any given trajectory the crucial step is to find and weight these critical planes. 
 
         [0030]     The structure factor is calculated as follows:  
                   c   m     ⁡     (       x   →     ,   t     )       =         sgn   ⁡     [         k   ^     ϕ     ·       y   →     ⁡     (   t   )         ]       ⁢     ω   ⁡     [       x   →     ,         k   ^     ϕ     ;   t       ]         ⁢     ❘     ϕ   =         ϕ   m     -     ∈         ϕ   =         ϕ   m     +     ∈             ,           (   5   )             
 
 where at φ m  either the function sgn[{circumflex over (k)} φ ·{right arrow over (y)}(t)] or the weighting function ω[{right arrow over (x)},{circumflex over (k)} φ ;t] as a function of φ is discontinuous. When a discontinuity occurs in either of these functions the plane Π({right arrow over (x)},t;φ) will be referred to as a critical plane. It is important to note that to reconstruct the density of a given image point, for each view angle, only the cone-beam projection data lying in these critical planes is necessary. 
 
         [0031]     Since the structure factor is highly dependent upon the selected weighting function we will also review here the conditions imposed on the weighting function. The weighting function properly accounts for data redundancy, i.e., for a given image point there may be multiple intersections between the plane Π({right arrow over (x)}, t;φ) and the source trajectory (i.e., {circumflex over (k)} φ ·[{right arrow over (x)}−{right arrow over (y)}(t)]=0). Obviously, the weighting function must be properly normalized,  
                 ∑   i     ⁢     ω   ⁡     [       x   →     ,         k   ^     ϕ     ;   t       ]         =   1.           (   6   )             
 
         [0032]     To obtain the above reconstruction formula, a necessary condition on the weighting function must be imposed  
                 ∂     ∂   ϕ       ⁢     ω   ⁡     (       x   →     ,       k   ^     ;   t       )         =   0.           (   7   )             
 
         [0033]     Namely, that the weighting function must be piecewise constant with respect to the angle of rotation, φ, of the plane Π({right arrow over (x)},t;φ) about {circumflex over (β)}. There are many possible weighting functions which satisfy both Eq. (6) and Eq. (7). In this work we have selected an equal weighting function which is desirable in medical CT since equal weighting is optimal for dose utilization. The equal weighting function which satisfies the conditions expressed in Eq. (6) and Eq. (7) is given as  
                 ω   ⁡     [       x   →     ,         k   ^     ϕ     ;   t       ]       =     1   N       ,           (   8   )             
 
 where N is the number of intersection points between plane Π({right arrow over (x)},t;φ) and the source trajectory. 
 
         [0034]     Up to this point, we have described the exact cone-beam FBP type framework for a complete source trajectory. In the following discussion, this is extended to develop an image reconstruction method for a single arc source trajectory.  
         [0035]     An arc source trajectory with radius R can be parameterized by the view angle variable t. thus, the position of any given view angle on the source trajectory will be given by 
 
{right arrow over ( y )}( t )= R (cos  t , sin  t,  0),  tε[t   i   ,t   f ],  (9) 
 
 where t i  and t f  are the parameters for the two ending points of the arc. The source-to-detector distance has been assumed to be twice the gantry radius, D=2R. 
 
         [0036]     In the geometrical setup introduced above and shown in  FIG. 8 , the detector plane is always perpendicular to the iso-ray. The analysis presented here assumes a flat equi-spaced detector, however a similar analysis may be conducted for a curved equi-angular detector. For convenience we define the following rotating system: 
 
{circumflex over (ω)}=(−cos  t , −sin  t,  0), û=(sin  t , −cos  t,  0), {circumflex over (v)}=(0,0,1). 
 
 The origin of the detector coordinate system is specified by the piercing point of the iso-ray onto the detector plane. The native detector coordinates u, v which are the signed distance along the directions û, {circumflex over (v)} respectively, will be introduced,  
               u   =     D   ⁢         r   ^     ·     u   ^           r   ^     ·     ω   ^             ,     v   =     D   ⁢           r   ^     ·     v   ^           r   ^     ·     ω   ^         .                 (   10   )             
 
         [0037]     The derivative  
       ∂     ∂   q         
 
 can be re-expressed using the chain rule of differentiation:  
             ∂     ∂   q       ⁢     g   ⁡     [       r   ^     ,       y   →     ⁡     (   q   )         ]         ⁢     ❘     q   =   t         =       (         -         D   2     +     u   2       D       ⁢     ∂     ∂   u         -       uv   D     ⁢     ∂     ∂   v         +     ∂     ∂   t         )     ⁢       g   ⁡     (     u   ,   v   ,   t     )       .           
 
         [0038]     By specifying the unit vector {circumflex over (β)} as one specific {circumflex over (r)} in this native geometry, i.e., substituting the previous expression for {circumflex over (β)} (Eq. (2)) into Eq. (10) one may obtain the cone-beam projection of the image point onto the detector, as given below:  
                   u   β     ⁡     (   t   )       =     D   ⁢         x   ⁢           ⁢   sin   ⁢           ⁢   t     -     y   ⁢           ⁢   cos   ⁢           ⁢   t           D   2     -     x   ⁢           ⁢   cos   ⁢           ⁢   t     -     y   ⁢           ⁢   sin   ⁢           ⁢   t             ,         v   β     ⁡     (   t   )       =     D   ⁢       z       D   2     -     x   ⁢           ⁢   cos   ⁢           ⁢   t     -     y   ⁢           ⁢   sin   ⁢           ⁢   t         .                 (   11   )             
 
 The detector position of the cone-beam projection of each image point from each view angle as given by u β (t) and v β (t) will later be used directly in performing the backprojection step. In addition, the angle γ ( FIG. 9 ) may also be expressed in terms of the detector coordinates u and v. 
 
         [0039]     To do so, we note that Eq. (3) dictates that: 
 
cos γ={circumflex over (β)}·{circumflex over (r)}.  (12) 
 
 Note that, in the local coordinate system, if the angle φ is fixed, the vector r is confined to the “critical fan” within the plane Π({right arrow over (x)}, t;φ) shown in  FIG. 9 . In other words, the point (u,v) where the ray vector {circumflex over (r)} is incident on the detector traverses the filtering line. In the detector plane, the slope of the filtering line is denoted as κ({right arrow over (x)},t):  
               κ   ⁡     (       x   →     ,   t     )       =         v   -       v   β     ⁡     (   t   )           u   -       u   β     ⁡     (   t   )           .             (   13   )             
 
 Note that the {right arrow over (x)} and t dependence in κ originates from the x and t dependence in u β (t) and v β (t). Therefore, one can see that the same filtering line may be used for the reconstruction of any image point lying in the critical plane Π({right arrow over (x)},t;φ). 
 
         [0040]     Finally we may re-express the filtering kernel in terms of the detector coordinates introduced above. Thus, using Eqs. (11)-(13), we obtain:  
                   d   ⁢           ⁢   γ       sin   ⁢           ⁢   γ       =       1          sin   ⁢           ⁢   α   ⁢           ⁢     cos   ⁡     (     β   -   t     )                ⁢     1         D   2     +     u   2     +     v   2           ⁢     du     u   -     u   β             ,           (   14   )             
 
 where v satisfies Eq. (13), (i.e. the point lies on the filtering line). 
 
         [0041]     Finally, substituting Eqs. (11) and (14) into Eq. (4), the following reconstruction formula is obtained:  
                 f   ⁡     (     x   →     )       =       -     1     4   ⁢     π   2           ⁢       ∑   m     ⁢     ∫       ⅆ   t     ⁢         c   m     ⁡     (       x   →     ,   t     )         L   ⁡     (       x   →     ,   t     )         ⁢       q   m     ⁡     [         u   β     ⁡     (   t   )       ,         v   β     ⁡     (   t   )       ;   t       ]                 ,     
     ⁢   where           (   15   )                   L   ⁡     (       x   →     ,   t     )       =       D   2     =       x   ⁢           ⁢   cos   ⁢           ⁢   t     -     y   ⁢           ⁢   sin   ⁢           ⁢   t           ,     
     ⁢         q   m     ⁡     [       u   _     ,       v   _     ;   t       ]       =       ∫     -   ∞       +   ∞       ⁢       ∫     -   ∞       +   ∞       ⁢           ⁢       ⅆ   u     ⁢           ⁢     ⅆ   v     ⁢           ⁢     δ   ⁡     [     v   -     v   _     -       κ   m     ⁡     (     u   -     u   _       )         ]       ⁢       h   H     ⁡     (     u   -     u   _       )       ⁢       g   _     ⁡     (     u   ,   v   ,   t     )               ,     
     ⁢         g   _     ⁡     (     u   ,   v   ,   t     )       =       D         D   2     +     u   2     +     v   2           ⁢     (             D   2     +     u   2       D     ⁢     ∂     ∂   u         +       uv   D     ⁢     ∂     ∂   v         -     ∂     ∂   t         )     ⁢     g   ⁡     (     u   ,   v   ,   t     )           ,           (   16   )             
 
 where  
           h   H     ⁡     (   x   )       =     1   x         
 
 is the kernel of Hilbert filter, and κ m  is the slope of the filtering line in the detector plane. 
 
         [0042]     In applying the Katsevich-type cone beam reconstruction formula described above a key task is to calculate the structure factor for the specific arc geometry. Namely, for each given reconstruction point {right arrow over (x)} and source point t one must find and weight the critical planes. From Eq. 7 the weighting function is piecewise constant with respect to φ. In order to calculate the structure factor one must find the transition points where the weighting function changes as the plane Π({right arrow over (x)},t;φ) rotates about β. In order to facilitate the analysis, one divides the source trajectory into several different regions based on the number of intersection points between the plane Π({right arrow over (x)},t;φ) (the plane that is perpendicular to the scanning plane) and the source trajectory as shown in  FIGS. 10   a - 10   c . For each image point the source trajectory may be segmented into three regions T 1 , (tε[t i , t a ]), T 2 , (tε(t a , t b )) and T 3 , (tε[t b , t f ]). The values t a  and t b  are easily determined analytically by projecting the image point to the x-y plane, and then respectively drawing lines from t f  and t i  through the projected image point as shown in  FIG. 10   d . The intersection between these lines and the source trajectory define t a  and t b , and are used to segment the source trajectory into T 1 , T 2 , and T 3  as given below in Eq. 20.  
         [0043]     The three sets T 1 ({right arrow over (x)}), T 2 ({right arrow over (x)}) and T 3 ({right arrow over (x)}) are defined as:  
                 T   1     =       (     x   →     )     ⁢     {     t   |       t   i     ≤   t   ≤       t   a     ⁡     (     x   →     )           }         ,           (   17   )                   T   2     =       (     x   →     )     ⁢     {     t   |         t   a     ⁡     (     x   →     )       &lt;   t   &lt;       t   b     ⁡     (     x   →     )           }         ,           (   18   )                   T   3     =       (     x   →     )     ⁢     {     t   |         t   b     ⁡     (     x   →     )       ≤   t   ≤     t   f         }         ,     
     ⁢   where           (   19   )                   t   n     =       [       t   m     +   π   -     2   ⁢           ⁢       tan     -   1       ⁡     (           -   x     ⁢           ⁢   sin   ⁢           ⁢     t   m       +     y   ⁢           ⁢   cos   ⁢           ⁢     t   m           R   -     x   ⁢           ⁢   cos   ⁢           ⁢     t   m       -     y   ⁢           ⁢   sin   ⁢           ⁢     t   m           )           ]     ⁢   mod   ⁢           ⁢   2   ⁢   π       ,           (   20   )             
 
 where n=a, b corresponding to m=i, f respectively and R is the gantry radius. 
 
         [0044]     It is important to know the manner in which the weighting function changes for a given image point and a given source point as the plane Π({right arrow over (x)},t;φ) rotates about {circumflex over (β)}. As this plane rotates about {circumflex over (β)} there are three critical angles where the weighting function abruptly changes. These three examples are illustrated for a sample image point {right arrow over (x)} and a sample source point tεT 2  in  FIG. 11 . The first transition occurs when the plane Π({right arrow over (x)},t;φ) contains {right arrow over (y)}(t i ) and thus this critical angle will be referred to as φ. In this case there is a transition from N=1 (before the plane Π({right arrow over (x)},t;φ) reaches {right arrow over (y)}(t i )) to N=2 (after the plane Π({right arrow over (x)},t;φ) passes {right arrow over (y)}(t i ), therefore the weighting function has a transition from 1 to ½ as shown in  FIG. 11   a . The second transition occurs at the critical angle when the plane Π({right arrow over (x)},t;φ) is tangential to the source trajectory and thus there is an abrupt jump from N=2 to N=1 and back to N=1 as shown in  FIG. 11   b . When this transition occurs the associated filtering line is horizontal and thus this critical angle will be referred to as φ h . The final case shown in  FIG. 11   c  is similar to the first case. Namely, the critical angle occurs when the plane Π({right arrow over (x)},t;φ) contains {right arrow over (y)}(t f ) as there is an abrupt transition from N=2 to N=1. It is clear that the discontinuities will be similar for source points in T 1  and T 3 ; namely, due to transitions as the plane Π({right arrow over (x)},t;φ) intersects either of the end points or becomes tangential to the source trajectory. To summarize, for all view angles there will be an abrupt jump in the weighting function when φ=φ m , m=h, i, f.  
         [0045]     In order to calculate the structure factor we need to consider sgn[{circumflex over (k)} φ ·{right arrow over (y)}(t)]. The argument of this function is calculated to be:  
                         k   ^     ϕ     ·       y   →     ⁡     (   t   )         =     r   ⁡     [       cos   ⁢           ⁢     ϕcos   ⁡     (     β   -   t     )         -     sin   ⁢           ⁢   ϕcos   ⁢           ⁢     αsin   ⁡     (     β   -   t     )           ]                   =     r   ⁢           cos   2     ⁡     (     β   -   t     )       +       cos   2     ⁢   α   ⁢           ⁢       sin   2     ⁡     (     β   -   t     )             ⁢       sin   ⁡     (     ϕ   -     ϕ   n       )       .                     (   21   )             
 
 Then we have:  
               sgn   ⁡     [         k   ^     ϕ     ·       y   →     ⁡     (   t   )         ]       =     {             -   1     ,               if   ⁢           ⁢   ϕ     ∈     [     0   ,     ϕ   h       )       ,               1   ,             if   ⁢           ⁢   ϕ     ∈       (     0   ,   π     ]     .                       (   22   )             
 
 Using the definition of the structure factor and the results in the preceding subsections, for z≠one may obtain:  
                 c   h     =   1     ,     
     ⁢       c   i     =     {                 -     1   2       ,               if   ⁢           ⁢   t     ∈       T   3     ⁡     (     x   →     )         ,               1   2               if   ⁢           ⁢   t     ∈         T   1     ⁡     (     x   →     )       ⋃       T   2     ⁡     (     x   →     )           ,           ⁢     
     ⁢     c   f       =     {                 -     1   2       ,               if   ⁢           ⁢   t     ∈       T   1     ⁡     (     x   →     )         ,               1   2               if   ⁢           ⁢   t     ∈         T   2     ⁡     (     x   →     )       ⋃       T   3     ⁡     (     x   →     )           ,           ⁢     
     ⁢   and   ⁢           ⁢   for   ⁢           ⁢   z     =     0   ⁢     :                         (   23   )               c   =     {           2   ,               if   ⁢           ⁢   t     ∈       T   2     ⁡     (     x   →     )         ,             1           if   ⁢           ⁢   t     ∈         T   1     ⁡     (     x   →     )       ⋃         T   3     ⁡     (     x   →     )       .                         (   24   )             
 
         [0046]     After substituting ( 23 ) into (24), we obtain the image reconstruction formula for an arc scanning path:  
                       -   4     ⁢     π   2     ⁢     f   ⁡     (     x   →     )         =       ⁢         ∫     t   i       t   f       ⁢           ⁢       ⅆ   t     ⁢     1     L   ⁡     (       x   →     ,   t     )         ⁢       q   h     ⁡     [         u   β     ⁡     (   t   )       ,         v   β     ⁡     (   t   )       ;   t       ]           +                     ⁢         1   2     ⁢     (       ∫     t   i       t   b       ⁢     -     ∫     t   b       t   f           )     ⁢           ⁢     ⅆ   t     ⁢     1     L   ⁡     (       x   →     ,   t     )         ⁢       q   f     ⁡     [         u   β     ⁡     (   t   )       ,         v   β     ⁡     (   t   )       ;   t       ]         -                     ⁢         1   2     ⁢     (       ∫     t   i       t   a       ⁢     -     ∫     t   a       t   f           )     ⁢           ⁢     ⅆ   t     ⁢     1     L   ⁡     (       x   →     ,   t     )         ⁢       q   f     ⁡     [         u   β     ⁡     (   t   )       ,         v   β     ⁡     (   t   )       ;   t       ]         ,                   (   25   )             
 
 where t a , t b  are shown in  FIG. 10   d  and filter functions q h , q i  and q f  were defined in Eq. 16. The expression for the slopes κ used in the filtering process are given below,  
                 κ   h     =   0     ,         κ   m     ⁡     (       u   _     ,       v   _     ;   t       )       =       v   _         u   _     -     D   ⁢           ⁢     cot   ⁡     (       t   -     t   m       2     )               ,           ⁢     m   =   i     ,     f   .             (   26   )             
 
         [0047]     To implement the preferred embodiment of the invention we rewrite the reconstruction formula (Eq. 25) into a more convenient form:  
                   -   4     ⁢     π   2     ⁢   f   ⁢     (     ⁢     x   -&gt;     ⁢     )       =         ∫     t   i       t   a       ⁢       ⅆ   t     ⁢     1       L   ⁢     (     ⁢     x   -&gt;       ,     t   ⁢     )           ⁢       Q   1     ⁡     [         u   _     =       u   β     ⁡     (   t   )         ,         v   _     b     =       v   β     ⁡     (   t   )         ,   t     ]           +       ∫     t   a       t   b       ⁢       ⅆ   t     ⁢     1       L   ⁢     (     ⁢     x   -&gt;       ,     t   ⁢     )           ⁢       Q   2     ⁡     [         u   _     =       u   β     ⁡     (   t   )         ,         v   _     b     =       v   β     ⁡     (   t   )         ,   t     ]           +       ∫     t   b       t   f       ⁢       ⅆ   t     ⁢     1       L   ⁢     (     ⁢     x   -&gt;       ,     t   ⁢     )           ⁢       Q   3     ⁡     [         u   _     =       u   β     ⁡     (   t   )         ,           v   _     b     =       v   β     ⁡     (   t   )         ;   t       ]               ,     
     ⁢   where           (   27   )                     q   h     ⁡     (       u   _     ,     v   _     ,   t     )       =       ∫     -   ∞       +   ∞       ⁢       ⅆ   u     ⁢           ⁢       h   H     ⁡     (     u   -     u   _       )       ⁢       g   _     ⁡     [     u   ,     v   _     ,   t     ]             ,     
     ⁢         q   m     ⁡     (       u   _     ,     v   m     ,   t     )       =       ∫     -   ∞       +   ∞       ⁢       ⅆ   u     ⁢           ⁢       h   H     ⁡     (     u   -     u   _       )       ⁢       g   _     ⁡     [     u   ,     v   m     ,   t     ]             ,     
     ⁢         Q   1     ⁡     (       u   _     ,     v   _     ,   t     )       =         q   h     ⁡     (       u   _     ,     v   _     ,   t     )       +       1   2     ⁢       q   i     ⁡     (       u   _     ,     v   i     ,   t     )         -       1   2     ⁢       q   f     ⁡     (     u   ,     v   f     ,   t     )             ,     
     ⁢         Q   2     ⁡     (       u   _     ,     v   _     ,   t     )       =         q   h     ⁡     (       u   _     ,     v   _     ,   t     )       +       1   2     ⁢       q   i     ⁡     (       u   _     ,     v   i     ,   t     )         -       1   2     ⁢       q   f     ⁡     (     u   ,     v   f     ,   t     )             ,     
     ⁢         Q   3     ⁡     (       u   _     ,     v   _     ,   t     )       =         q   h     ⁡     (       u   _     ,     v   _     ,   t     )       +       1   2     ⁢       q   i     ⁡     (       u   _     ,     v   i     ,   t     )         -       1   2     ⁢       q   f     ⁡     (     u   ,     v   f     ,   t     )             ,     
     ⁢     where   ⁢           ⁢     v   m     ⁢           ⁢   is   ⁢           ⁢   given   ⁢           ⁢   by             (   28   )                   v   m     =         v   _         u   _     =     D   ⁢           ⁢     cot   ⁡     (       t   -     t   m       2     )             ⁡     [     u   -     D   ⁢           ⁢     cot   ⁡     (       t   -     t   m       2     )           ]         ,     m   =   i     ,     f   .             (   29   )             
 
         [0048]     In the implementation of this reconstruction formula it is important to note that the {right arrow over (x)} dependence of Q m  (m=1,2, or 3) in Eq. 25 (or Eq. 27) may be expressed directly in terms of detector coordinates u β (t) and v β (t). Thus, the filtering operation is not explicitly voxel dependent, and this reconstruction formula may be implemented in a manner similar to the standard FDK method. For the full-scan case in which the x-ray source and detector array revolve completely around the subject, the filtering direction is identical to the FDK method and in the case of a short and super-short scan there are three separate filtering directions, linear combinations of which makeup three separate sets of filtration data to be backprojected.  
         [0049]     The reconstruction steps are summarized below: 
 
 Step 1: Modify the cone-beam data by computing the weighted derivative of the measured cone-beam data,  
           g   _     ⁡     (     u   ,   v   ,   t     )       =       D         D   2     +     u   2     +     v   2           ⁢     (             D   2     +     u   2       D     ⁢     ∂     ∂   u         +         u   ⁢           ⁢   v     D     ⁢     ∂     ∂   v         -     ∂     ∂   t         )     ⁢       g   ⁡     (     u   ,   v   ,   t     )       .           
 
 Step 2: Compute the ID Hilbert filter of the modified cone-beam data to obtain Q 1 , Q 2 , Q 3  according to Eq. 28. 
 
 Step 3: backproject Q 1 , Q 2 , Q 3  according to Eq.27 to obtain f({right arrow over (x)}) for each point {right arrow over (x)}. 
 
       DESCRIPTION OF THE PREFERRED EMBODIMENT  
       [0050]     With initial reference to  FIGS. 6 and 7 , a computed tomography (CT) imaging system  10  includes a gantry  12  representative of a “third generation” CT scanner. Gantry  12  has an x-ray source  13  that projects a cone beam of x-rays  14  toward a detector array  16  on the opposite side of the gantry. The detector array  16  is formed by a number of detector elements  18  which together sense the projected x-rays that pass through a medical patient  15 . Each detector element  18  produces an electrical signal that represents the intensity of an impinging x-ray beam and hence the attenuation of the beam as it passes through the patient. During a scan to acquire x-ray projection data, the gantry  12  and the components mounted thereon rotate about a center of rotation  19  located within the patient  15 . The path y(t) of the x-ray source  13  is thus a circular arc.  
         [0051]     The rotation of the gantry and the operation of the x-ray source  13  are governed by a control mechanism  20  of the CT system. The control mechanism  20  includes an x-ray controller  22  that provides power and timing signals to the x-ray source  13  and a gantry motor controller  23  that controls the rotational speed and position of the gantry  12 . A data acquisition system (DAS)  24  in the control mechanism  20  samples analog data from detector elements  18  and converts the data to digital signals for subsequent processing. An image reconstructor  25 , receives sampled and digitized x-ray data from the DAS  24  and performs high speed image reconstruction according to the method of the present invention. The reconstructed image is applied as an input to a computer  26  which stores the image in a mass storage device  29 .  
         [0052]     The computer  26  also receives commands and scanning parameters from an operator via console  30  that has a keyboard. An associated cathode ray tube display  32  allows the operator to observe the reconstructed image and other data from the computer  26 . The operator supplied commands and parameters are used by the computer  26  to provide control signals and information to the DAS  24 , the x-ray controller  22  and the gantry motor controller  23 . In addition, computer  26  operates a table motor controller  34  which controls a motorized table  36  to position the patient  15  in the gantry  12 .  
         [0053]     The above-described CT system is designed to rapidly perform a full-scan in which the x-ray source revolves in a circular arc completely around the subject being imaged. The present invention may be employed in such a CT system, but because the invention may also be employed in short-scans and super-short-scans, it is particularly useful when used with a C-arm x-ray system depicted in  FIG. 14 .  
         [0054]     Referring particularly to  FIG. 14 , the C-arm x-ray system includes a C-arm  60  to which the two-dimensional detector  16  and the x-ray source  13  are mounted. Here, again the patient  15  is positioned on a table  62  that extends between the source  13  and the detector  16 . The C-arm is rotationally mounted to a base  64 , and cone beam data is acquired by rotating the x-ray source  13  and detector  16  around a defined axis  66 . C-arm scanners are particularly useful in image-guided interventions and they are characterized by scans in which the x-ray source and detector are rotated over a very short arc.  
         [0055]     As shown best in  FIG. 12   a , in the preferred embodiment of the present invention the detector array  16  is a flat array of detector elements  18 , having N r  (e.g. 1024) elements  18  disposed along the in-plane (x,y) direction, and N z  (e.g. 1024) elements  18  disposed along the z axis. The x-ray beam emanates from the x-ray source  13  and fans out as it passes through the patient  15  and intercepts the detection array  16 . Each acquired view is a N r  by N z  array of attenuation measurements as seen when the gantry is oriented in one of its positions during the scan. As shown in  FIG. 12B , the object of the present invention is to reconstruct a set of 2D image slices  35  from the 3D array of acquired data produced by the x-ray cone beam during the scan.  
         [0056]     It can be seen that because the cone beam diverges as it passes through the patient  15 , the accurate reconstruction of the parallel image slices  35  is not possible with a straight forward fan beam filtering and backprojection process. The present invention enables an accurate reconstruction of the image slices  35  from this acquired cone beam data. In addition, the present invention enables the accurate reconstruction of image slices  35  even when the circular path of the x-ray source  13  is less than a complete circle.  
         [0057]     This reconstruction method is implemented in the image reconstructor  25 . Referring particularly to  FIG. 13 , the cone beam projection data is received from the DAS  24  as a two-dimensional array of values which are preprocessed in the standard manner at process block  40 . Such preprocessing includes correcting for known errors and offsets and calculating the minus log of the data to convert it to x-ray attenuation values. The preprocessed cone beam attenuation profiles g(u,v,t) are then modified by calculating the weighted derivative {overscore (g)}(u,v,t) as indicated at process block  42 . This is performed as set forth above in Eq. (16). The ID Hilbert filter of the weighted derivative of the cone beam data {overscore (g)}(u,v,t) is then computed as set forth above in Eq. (28) to produce filtered cone beam data sets Q 1 ({overscore (u)},{overscore (v)},t), Q 2 ({overscore (u)},{overscore (v)},t) and Q 3 ({overscore (u)},{overscore (v)},t) as indicated at process blocks  44 ,  46  and  48 . As indicated by process blocks  50 ,  52  and  54  these filtered cone beam data sets Q 1 ({overscore (u)},{overscore (v)},t), Q 2 ({overscore (u)},{overscore (v)},t) and Q 3 ({overscore (u)},{overscore (v)},t) are each backprojected as set forth above in Eq. (27) and the backprojected data is summed as indicated at process block  56  to produce a 3D image f(x). The desired 2D image slices  58  are produced from the 3D image and output to computer  26 .  
         [0058]     The present invention is valid when the x-ray source arc covers an angular range of 180°+fan angle, namely a short scan mode. The backprojection segments T 1 , T 2  T 3  is correspondingly changed for each image voxel. The present invention is also valid even if the source arc is shorter than 180°+fan angle. This is a super-short scan mode.  
         [0059]     A number of variations are possible from the preferred embodiment described above. The differentiations in preweighted data may be computed using a fast Fourier transform (FFT) method. Also, the Hilbert convolution process may be implemented using an FFT method. Parameters such as source to detector distance, source to isocenter distance, the source to image point distance may be measured in a geometrical calibration process and incorporate into the method. The convolution process may be implemented using either a horizontal detector coordinate or vertical coordinate. The convolution procedure may be numerically implemented by combining the convolution procedure in terms of horizontal coordinate and in terms of vertical coordinate. The present invention may be applied to a curved detector and circular source trajectory.  
         [0060]     It should be apparent to those skilled in the art that the convolution procedure may be numerically implemented by combining the convolution procedure in terms of horizontal coordinate and in terms of vertical coordinate. The present invention may also be applied to a curved detector and circular source trajectory. The present invention may be utilized to reconstruct an image for other clinical applications such as radiation therapy where an x-ray source and a flat-panel imager slowly rotates around the patient. While the invention is most advantageously applied to cone beam data sets, the present invention may also be utilized to reconstruct fan-beam CT images. This is the case of Nz=1.