Abstract:
In one embodiment, the present invention is a method for reconstructing a scanned CT image that includes steps of: acquiring projection data of an object utilizing a flat-panel detector; and filtering and backprojecting the projection data, without radially interpolating, to reconstruct a three-dimensional image of the object. This embodiment provides direct fan-parallel, three-dimensional reconstruction of computed tomographic images without loss of resolution associated with radial interpolation, while retaining imaging quality and backprojection speed gains associated with parallel reconstruction.

Description:
BACKGROUND OF THE INVENTION 
     This invention relates to methods and apparatus for volumetric computed tomograph (CT) reconstruction and more specifically to methods and apparatus for fan parallel cone beam volumetric CT reconstruction utilizing a flat panel detector. 
     In at least one known computed tomography (CT) imaging system configuration, an x-ray source projects a fan-shaped beam which is collimated to lie within an X-Y plane of a Cartesian coordinate system and generally referred to as the “imaging plane”. The x-ray beam passes through the object being imaged, such as a patient. The beam, after being attenuated by the object, impinges upon an array of radiation detectors. The intensity of the attenuated beam radiation received at the detector array is dependent upon the attenuation of the x-ray beam by the object. Each detector element of the array produces a separate electrical signal that is a measurement of the beam attenuation at the detector location. The attenuation measurements from all the detectors are acquired separately to produce a transmission profile. 
     In known third generation CT systems, the x-ray source and the detector array are rotated with a gantry within the imaging plane and around the object to be imaged so that the angle at which the x-ray beam intersects the object constantly changes. A group of x-ray attenuation measurements, i.e., projection data, from the detector array at one gantry angle is referred to as a “view”. A “scan” of the object comprises a set of views made at different gantry angles, or view angles, during one revolution of the x-ray source and detector. In an axial scan, the projection data is processed to construct an image that corresponds to a two dimensional slice taken through the object. One method for reconstructing an image from a set of projection data is referred to in the art as the filtered back projection 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 cathode ray tube display. 
     Three-dimensional reconstruction of scanned CT images is possible. A direct fan parallel reconstruction method is used for data acquired using a flat panel detector. The use of parallel reconstruction advantageously provides an amount of aliasing cancellation. However, in at least one known method, a radial interpolation step is required when rebinning from fan beam to parallel data. This interpolation step produces a loss of resolution. 
     It would therefore be desirable to provide methods and apparatus for direct fan parallel, three dimensional reconstruction of computed tomographic images that avoid the loss of resolution associated with radial interpolation and that retain the image quality and backprojection speed gains associated with parallel reconstruction. 
     BRIEF SUMMARY OF THE INVENTION 
     There is therefore provided, in one embodiment, a method for reconstructing a scanned CT image that includes steps of: acquiring projection data of an object utilizing a flat-panel detector; and filtering and backprojecting the projection data, without radially interpolating, to reconstruct a three-dimensional image of the object. 
     The above described embodiment provides direct fan-parallel, three-dimensional reconstruction of computed tomographic images without loss of resolution associated with radial interpolation, while retaining imaging quality and backprojection speed gains associated with parallel reconstruction. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a pictorial view of a CT imaging system. 
     FIG. 2 is a block schematic diagram of the system illustrated in FIG.  1 . 
     FIG. 3 is a representation of fan beam and fan parallel geometries for one fan formed by a row of a planar two-dimensional detector orthogonal to the system z-axis and a fan-vertex V. 
     FIG. 4 is a representation of geometry for a fan-beam projection defined by readings of a row of detector cells. 
     FIG. 5 is a geometric representation of a rebinning process of one embodiment of the present invention in which rebinning is to a fan-parallel geometry for a flat panel detector and a circular source trajectory. 
     FIG. 6 is a different representation of the geometry of the rebinning process illustrated in FIG.  5 . 
     FIG. 7 is a geometric representation of a rebinning process of one embodiment of the present invention, where rebinning is to a fan parallel geometry for a flat parallel detector and a helical source trajectory. 
     FIG. 8 is a drawing that shows a representation of an orthogonal projection of a scanner isocenter onto a plane defined by an xy line of detector cells and a cone vertex defined by an x-ray source and that defines a tilted plane defined by the vertex and a an axis x′ parallel to and offset from an x axis. 
     FIG. 9 is a geometric representation of a voxel M to be reconstructed. 
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     Referring to FIGS. 1 and 2, a computed tomograph (CT) imaging system  10  is shown as including a gantry  12  representative of a “third generation” CT scanner. Gantry  12  has an x-ray source  14  that projects a beam of x-rays  16  toward a detector array  18  on the opposite side of gantry  12 . Detector array  18  is formed by detector elements or cells  20  which together sense the projected x-rays that pass through an object  22 , for example a medical patient. Detector array  18  may be fabricated in a single slice or multi-slice configuration, and also, in various shapes, such as an arc of a circle, an arc of a cylinder, or flat panels. Each detector element or cell  20  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 patient  22 . During a scan to acquire x-ray projection data, gantry  12  and the components mounted thereon rotate about a center of rotation  24 . 
     Rotation of gantry  12  and the operation of x-ray source  14  are governed by a control mechanism  26  of CT system  10 . Control mechanism  26  includes an x-ray controller  28  that provides power and timing signals to x-ray source  14  and a gantry motor controller  30  that controls the rotational speed and position of gantry  12 . A data acquisition system (DAS)  32  in control mechanism  26  samples analog data from detector elements  20  and converts the data to digital signals for subsequent processing. An image reconstructor  34  receives sampled and digitized x-ray data from DAS  32  and performs high speed image reconstruction. The reconstructed image is applied as an input to a computer  36  which stores the image in a mass storage device  38 . 
     Computer  36  also receives commands and scanning parameters from an operator via console  40  that has a keyboard. An associated cathode ray tube display  42  allows the operator to observe the reconstructed image and other data from computer  36 . The operator supplied commands and parameters are used by computer  36  to provide control signals and information to DAS  32 , x-ray controller  28  and gantry motor controller  30 . In addition, computer  36  operates a table motor controller  44  which controls a motorized table  46  to position patient  22  in gantry  12 . Particularly, table  46  moves portions of patient  22  through gantry opening  48 . 
     In one embodiment of the present invention, a fan-parallel weighted-convolution reconstruction kernel decomposition is utilized. Let f(r,φ) and  f (x,y) denote respectively polar and cartesian representations of an x-ray attenuation coefficient distribution to be reconstructed, and p(s,θ) a polar representation of a Radon transform of  f :                p        (     s   ,   θ     )       =       ∫     -   ∞     ∞            ∫     -   ∞     ∞              f   _          (     x   ,   y     )            δ        (       x                 cos                 θ     +     y                 sin                 θ     -   s     )               x             y                   (   1   )                                
     Under conditions satisfied by any patient or scanned object, and referring to FIG. 3, the Radon inversion formula is written as:                f        (     r   ,   φ     )       =       1     4        π   2                ∫   0     2      π              lim     ɛ   →   0              ∫     -   ∞     ∞              F   ɛ          (   t   )            p        (     s   ,   θ     )               s             θ                       (   2   )                                
     with: 
     
       
         t=s−r cos(θ−φ)  (3) 
       
     
     and:                  F   ɛ          (   t   )       =       (           1     ɛ   2                  t        ≤   ɛ                 -   1       t   2                  t        ≥   ɛ           )     .             (   4   )                                
      where FIG. 3 illustrates fan-beam and fan-parallel geometries for one fan formed by a horizontal row of a planar two-dimensional detector and fan vertex V. Two source positions (from vertices) are shown at two angles, β and β′. 
     In the following, the Radon space parameterization (s,θ) is changed into the fan-beam parameterization (β, η(ν)), and the integral is then expressed as a weighted convolution. Reference is made to FIG. 3 discussed above, and to FIG. 4, which is a representation of geometry for one fan beam projection defined by the readings of a row of detector cells in plane x′Oy′. S and D are, respectively, distances from fan-vertex V and the detector plane to O; E=S+D, and s is a Radon distance (in the plane of FIG. 4) associated with a ray detected at detector cell position υ. Introducing β, an angle from the y axis to fan-vertex V, S, a distance from vertex V to O and P, a point to be reconstructed, projection data p(s,θ) is parameterized via angle η from VO to VP (see FIG.  3 ), with η(0)=0 and η limited to [η min ,η max ] so that an entire object cross-section (bounded in an image plane) is covered. Parameters (β,ν) are independent variables, with β being a source (vertex) angle and ν a detector cell coordinate on a straight line defined by an intersection of an imaging plane and detector  18  plane. Changing variables (s,θ)→(β,ν) via the parameterization η(ν) produces:                [         θ           S         ]     =       [           β   +     η        (   ν   )                   S                   sin        [     η        (   ν   )       ]               ]     .             (   5   )                                
     Parameter limits are defined by β min , β max , ν min , ν max  so that: 
     
       
         β max −β min =2π; η(ν min )=η min ; η(ν max )=η max .  (6) 
       
     
     For a flat-panel detector  18 , and considering all data to be filtered to be acquired along a single “horizontal” detector row, fan-beam parameterization is given by:                  η        (   ν   )       =       arctan        [     v     S   +   D       ]       =     arctan        [     ν   E     ]           ,           (   7   )                                
     where E=S+D is a sum of distances from source V to isocenter O, and D from O to detector  18  (as measured along a ray from V to detector  18  through O). 
     The Jacobian of the transformation is written as:              J   =         S   ×   E   ×     cos        [     arctan        (     ν   E     )       ]             E   2     +     v   2         =       S   E            1       (     [     1   +       (     v   E     )     2       ]     )       3   /   2         .                 (   8   )                                
     Variable t, the argument of F is written as: 
     
       
         t=s−r cos(θ−φ) 
       
     
     so that: 
     
       
         t=S×sin[η(ν)]−r cos(θ−φ).  (9) 
       
     
     In fan-parallel reconstruction, as in parallel reconstruction, angle θ is constant for one fan-parallel projection, unlike the case in fan-beam reconstruction. Defining a fan-angle for a ray passing through point f(r,φ) to be reconstructed: 
     
       
         {tilde over (η)}=η({tilde over (ν)}) with: S×sin[{tilde over (η)}]=r cos(θ−φ),  (10) 
       
     
     argument t is written as:              t   =       S   ×     {       sin        [     η        (   ν   )       ]       -     sin        [     η        (     ν   ~     )       ]         }       =       {       sin   [       η        (   ν   )       -     sin        [     η        (     ν   ~     )       ]             ν   -     ν   ~         }     ×   S   ×     (     ν   -     ν   ~       )                 (   11   )                                
     The following relation is used below:                  Lim     ɛ   →   0              ∫     -   ∞     ∞              F   ɛ          [     λ                   b        (   λ   )         ]            a        (   λ   )               λ           =       Lim     ɛ   →   0              ∫     -   ∞     ∞              F   ɛ          (   λ   )              a        (   λ   )           b   2          (   λ   )                   λ     .                   (   12   )                                
     Considering F, with t as above (11), then, under the summation sign:                  F   ɛ          (   t   )       =         F   ɛ          (     ν   -     ν   ~       )       ×     1     S   2       ×         {       ν   -     ν   ~           sin        [     η        (   ν   )       ]       -     sin        [     η        (     ν   ~     )       ]           }     2     .               (   13   )                                
     Let K be written as:                K        (     ν   ,     ν   ~       )       =     {         sin        [     η        (   ν   )       ]       -     sin        [     η        (     ν   ~     )       ]           ν   -     ν   ~         }             (   14   )                                
     If K can be written as a weighted convolution kernel:                K        (     ν   ,     ν   ~       )       =       {         sin        [     η        (   ν   )       ]       -     sin        [     η        (     ν   ~     )       ]           ν   -     ν   ~         }     =       E        (   ν   )            H        (     ν   -     ν   ~       )            G        (     ν   ~     )                   (   15   )                                
     where H is a convolution kernel, and E and G are pre- and post-convolution weights (E(ν) as used here is not to be confused with the constant E introduced in equation (7)), then:                  F   ɛ          (   t   )       =         F   ɛ          (     ν   -     ν   ~       )           S   2     ×       [       E        (   ν   )            H        (     ν   -     ν   ~       )            G        (     ν   ~     )         ]     2                 (   16   )                                
     is also in the form of a weighted-convolution (shift-invariant filtering). Note that Jacobian J (equation (8)) is a function of η, and therefore contributes a pre-convolution weight only. The reconstruction equation is thus written in the form of a weighted-convolution backprojection:                f        (     r   ,   φ     )       =       1     4        π   2          S   2                ∫   0     2      π              1       G        (     ν   ~     )       2            {       Lim     ɛ   →   0              ∫     ν   inf       ν   sup                  f   ɛ          (     ν   -     ν   ~       )           H   2          (     ν   -     ν   ~       )         ×       J        (   ν   )           E        (   ν   )       2       ×     p        (     β   ,   ν     )               ν           }             β                   (   17   )                                
     In equation (17), the          J        (   ν   )           E        (   ν   )       2                            
     term corresponds to a pre-filtering weighting; the bracketed {} expression corresponds to a filtering operation; the term        1       G        (     ν   ~     )       2                            
     corresponds to a post-filtering weighting; the integral over dβ corresponds to a backprojection operation; and the term        1     4        π   2          S   2                              
     corresponds to a normalization factor. Also: 
     r, φ, and θ are polar coordinates in a polar coordinate system; 
     f(r,φ) is a polar representation of an x-ray attenuation coefficient distribution to be reconstructed; 
     ν is a detector cell coordinate on a straight line defined by an intersection of an imaging plane and a plane defined by the flat panel detector; 
     β is an angle from a y-axis to a radiation source; 
     p(β,ν) is a representation of a radon transform of  f , which itself is a Cartesian representation of an x-ray attenuation distribution to be reconstructed; 
     S is a distance from a radiation source to an isocenter; 
     ν sup  is a maximum detector cell coordinate on flat panel  18 ; 
     ν inf  is a minimum detector cell coordinate; 
     H(ν−{tilde over (ν)}) and F ε (ν−{tilde over (ν)}) are convolution kernels; 
     η({tilde over (ν)}) is a fan angle for a ray passing through a voxel being reconstructed, and {tilde over (ν)} is a cell index for the ray; 
     E(ν) is a preconvolution weight that is a function of detector cell position ν; 
     G({tilde over (ν)}) is a post convolution weighting; and 
     J(ν) is a Jacobian function. 
     Standard techniques for adapting integrals to algorithms operating on discrete data are known and need not be described in detail here. In embodiments of the present invention, such techniques are used to approximate integrals such as those in equation (17). Thus, filtering and backprojection comprises a discrete reconstruction of an image as an approximation to the continuous weighted-convolution backprojection of equation (17). 
     In one embodiment of the invention, after acquiring projection data of an object utilizing a flat panel detector, image reconstruction proceeds as follows, according to equation (17), looping over projections corresponding to parameter β for each fixed image point. For a given projection, apply a pre-convolution weighting to the projection data using function E and Jacobian J; noting that pre-convolution weights are independent of an image point described by (r,φ)). Next, perform a convolution using kernels F and H, and apply a post-convolution weighting using function G; independent of the image point. Finally, perform a backprojection, and, in one embodiment, apply a normalization factor. Note that a pixel-dependent weighted fan-beam backprojection is not required for fan-parallel backprojection. 
     Necessary condition for reconstruction 
     A necessary condition for the fan-parallel kernel K to be decomposed as:                K        (     ν   ,     ν   ~       )       =       {         sin        [     η        (   ν   )       ]       -     sin        [     η        (   ν   )       ]           ν   -     ν   ~         }     =       E        (   ν   )            H        (     ν   -     ν   ~       )            G        (     ν   ~     )                   (   18   )                                
     can be expressed in term of the fan-beam parameterization η as:                K        (     ν   ,     ν   ~       )       =       {         sin        [     η        (   ν   )       ]       -     sin        [     η        (   ν   )       ]           ν   -     ν   ~         }     =         {           η   ′          (   ν   )            cos        [     η        (   ν   )       ]              η   ′          (     ν   ~     )            cos        [     η        (     ν   ~     )       ]               η   ′          (   0   )              η   ′          (     ν   -     ν   ~       )            cos        [     η        (     ν   -     ν   ~       )       ]           }       1   /   2                sin        [     η        (     ν   -     ν   ~       )       ]         ν   -     ν   ~         .                 (   19   )                                
     Equivalently, the following necessary condition is written for η:                  sin        [     η        (   ν   )       ]       -     sin        [     η        (     ν   ~     )       ]         =         {           η   ′          (   ν   )            cos        [     η        (   ν   )       ]              η   ′          (     ν   ~     )            cos        [     η        (     ν   ~     )       ]               η   ′          (   0   )              η   ′          (     ν   -     ν   ~       )            cos        [     η        (     ν   -     ν   ~       )       ]           }       1   /   2              sin        [     η        (     ν   -     ν   ~       )       ]       .               (   20   )                                
     In equation (20), η(ν) is a fan angle as a function of detector cell  20  coordinate ν. 
     From (4) and (5) it is written that: 
     
       
         ν=E×tan[arc sin(α)],  (21) 
       
     
     where α is a coordinate describing equidistant increments, is a solution. As is shown below, an equidistant parameterization written as:                  η        (   ν   )       =       arctan        [     ν     S   +   D       ]       =     arctan        [     ν   E     ]           ,           (   22   )                                
     associated with flat-panels with equidistant detector cells, is not a solution to equation (20). Yet an approximate solution provided by the necessary condition (and henceforth referred to as a “natural approximation”) is a very close one that leads to excellent imaging results. 
     Fan-parallel CT reconstruction from flat-panel detector 
     Substituting the parameterization (22) into the necessary condition (20) leads to an approximate relation written as:                  K        (     ν   ,     ν   ~       )       =       {         sin        [     η        (   ν   )       ]       -     sin        [     η        (   ν   )       ]           ν   -     ν   ~         }     =         ν       1   +       (     ν   E     )     2           -       ν   ~         1   +       (       ν   ~     E     )     2               ν   -     ν   ~             ,           (   23   )                 K        (     ν   ,     ν   ~       )       ≈         {           E   2          [       E   2     +       (     ν   -     ν   ~       )     2       ]           (       E   2     +     ν   2       )          (       E   2     +     ν   2       )         ×       cos                 arctan        (     ν   E     )        cos                   arctan        (       ν   ~     E     )           cos                   arctan        (       ν   -     ν   ~       E     )             }       1   /   2              sin                 arctan        (       ν   -     ν   ~       E     )         ν   -     ν   ~                   (   24   )                                
     which, using trigonometric relations, is rewritten as:                K        (     ν   ,     ν   ~       )       =           ν       1   +       (     ν   E     )     2           -       ν   ~         1   +       (       ν   ~     E     )     2               ν   -     ν   ~         ≈       1   E              [     1   +       (       ν   -     ν   ~       E     )     2       ]       1   /   4               [     1   +       (     ν   E     )     2       ]       3   /   4            [     1   +       (       ν   ~     E     )     2       ]         3   /   4                     (   25   )                                
     Accordingly, by substitution into (17), to the following image reconstruction algorithm                f        (     r   ,   φ     )       ≈       1     4        π   2              E   S            ∫   0     2      π              Lim     ɛ   →   0              ∫     ν   min       ν   max                [     1   +       (       ν   ~     E     )     2       ]       3   /   2       ×         F   ɛ          (     ν   -     ν   ~       )           1   +       (       ν   -     ν   ~       E     )     2           ×     p        (     β   ,   ν     )               ν             β                       (   26   )                                
     Recalling equation (4) for the Jacobian, an expression for pre-convolution weights is identically constant and equal to 1.0, and a filter expression is written as:                    F   ɛ          (     ν   -     ν   ~       )           [     1   +       (       ν   -     ν   ~       E     )     2       ]       1   /   2         =         F   ɛ          (     ν   -     ν   ~       )              cos        [     η        (     ν   -     ν   ~       )       ]       .               (   27   )                                
     Post-convolution weights are written:                  [     1   +       (       ν   ~     E     )     2       ]       3   /   2       =     1       cos   3          (     η   ~     )                 (   28   )                                
     For comparison purposes, a 2D fan-beam reconstruction formula for straight-line of detector cells is:          f        (     r   ,   φ     )       =       1     4        π   2                ∫   0     2      π              Lim     ɛ   →   0              ∫     ν   min       ν   max                (     S   L     )     2          1       cos   2          (     η   ~     )         ×       F   ɛ          (     ν   -     ν   ~       )       ×     p        (     β   ,   ν     )            cos        [     η        (   ν   )       ]               ν             β                                      
     The pre- and post-convolution weight expressions and filter expression of this embodiment differ from those for Feldkamp reconstruction (as well as fan-beam reconstruction from a straight line of detectors) and are unique to direct fan-parallel image reconstruction from a straight line of detector. For example, the 1/L 2  of fan-beam reconstruction is absent from the equation for fan-parallel reconstruction. 
     Cone-beam fan-parallel reconstruction for a circular source trajectory 
     In one embodiment, rebining is done on a row-by-row basis, to define an associated fan-parallel geometry. In addition, the dβ infinitesimal is changed to dβ′ to account for a fan average cone-angle and fan geometry is adapted to account for tilting of the fan beam away from a central gantry plane, i.e., an equivalent source rotation in a tilted plane. See FIGS. 5 and 6, which are graphical illustration of the process of rebinning to a fan-parallel geometry for a flat panel detector and a circular source trajectory. See also FIG. 8, in which O′ is an orthogonal projection of a scanner  10  isocenter onto a plane defined by an xy line of detector cells  20  and cone vertex V represents x-ray source  14 . A tilted fan is defined by V and an axis x′ parallel to and offset from z. O is on a central (or gantry  12 ) plane. The central plane is orthogonal to z and passes through V. 
     Cone-beam fan-parallel reconstruction for a helical source trajectory 
     In one embodiment employing a helical source trajectory, and referring to FIG. 7, rays are rebined to fan-parallel by interpolating a projection along z, so that the rays intersect an associated fan-parallel x′=t axis at a constant elevation. A weight is included to account for individual ray cone-angles, and an average cone-angle of the fan-parallel projection (equal to the cone-angle for the ray at γ=0) is used to calculate dβ′ in the tilted plane. Backprojection then proceeds on a voxel basis, as in the FDK algorithm. 
     From the preceding description of various embodiments of the present invention, it is evident that fan-parallel, three-dimensional reconstruction of computed tomographic images is achieved without loss of resolution associated with radial interpolation and without the loss of aliasing cancellation due to the 1/L 2  factor in fan-beam reconstruction. Referring to FIG. 9, in fan beam reconstruction, the two conjugate rays contributing to M (from source  14  positions V and V′) have unequal weights (except on x′); therefore limiting aliasing cancellation effectiveness away from isocenter  24  (as L≠L′,1/L 2 ≠1/L′ 2 ). This same consideration explains why noise in reconstructed images increases away from isocenter  24  in fan-beam reconstruction (assuming a constant noise level in backprojection data). Both of these effects (aliasing cancellation and noise) are improved with fan-parallel backprojection, due to the absence of the 1/L 2  factor in reconstruction equation (26). 
     Although particular embodiments of the invention have been described and illustrated in detail, it is to be clearly understood that the same is intended by way of illustration and example only and is not to be taken by way of limitation. In addition, the CT system described herein is a “third generation” system in which both the x-ray source and detector rotate with the gantry. Many other CT systems including “fourth generation” systems wherein the detector is a full-ring stationary detector and only the x-ray source rotates with the gantry, may be used if individual detector elements are corrected to provide substantially uniform responses to a given x-ray beam. Accordingly, the spirit and scope of the invention are to be limited only by the terms of the appended claims and legal equivalents.