Patent Publication Number: US-8116426-B2

Title: Computed tomography device and method using circular-pixel position-adaptive interpolation

Description:
BACKGROUND OF THE INVENTION 
     1. Field of the Invention 
     The present invention relates to helical and circular x-ray computed tomographic (CT) imaging, and in particular to CT imaging with adaptive interpolation. 
     2. Discussion of the Background 
     There are several methods of backprojection. A common method is pixel-driven backprojection illustrated in  FIG. 1 . A ray  11  is drawn from the source  10  through the center of the reconstruction pixel at (i,j) to the detector array  12 . Ray  11  intersects detector  12  at floating point channel index Ch i,j  on channel c k . The value assigned to the pixel is a linear interpolation of the values of the detector the ray intersects and the adjacent elements:
 
 V ( i,j )= w   k−1   ·D ( c   k−1 )+ w   k   ·D ( c   k )+ w   k+1   ·D ( c   k+1 )  (1)
 
where c k  is the channel of intersection and w k−1 +w k +w k+1 =1. For typical linear interpolation, the interpolation kernel width ΔCh is fixed at one channel and only two detector elements are used (either w k−1  or w k+1  is zero). A nearest neighbor implementation has w k+1 =w k−1 =0 and w k =1. Another pixel-driven method adjusts the weights such that interpolation kernel locations on the detector are always adjacent, as described in U.S. Pat. No. 6,724,856.
 
     Pixel-driven methods are easy to implement. However, the disadvantage is that actual projection of the pixel onto the detector array is not used to compute V(i,j). For large pixel sizes, the projection of the pixel onto the array will be larger than one channel, meaning that too few detector elements contribute to V(i,j), and the signal-to-noise ratio (SNR) will be poorer. For small pixel sizes the projection of the pixel onto the array is smaller than one channel, causing a loss of lateral resolution. 
     Another backprojection method is ray-driven backprojection. This technique is shown in  FIG. 2  and described in Barrett, Harrison and Swindell, William, “ Radiological Imaging: The Theory of Image Formation, Detection, and Processing ,” New York: Academic Press, 1981, pp. 421-422. In this case the backprojection ray  21  is drawn from the source  20  to the centers (indicated as c) of detector elements in the detector array  22 , and the contribution weight w ijc  of the detector element to pixels in the ray path is determined from the length of the ray within the pixel, as follows:
 
 V ( i,j )= V ( i,j )+ w   ijc   ·D ( c )  (2)
 
     Other ray driven methods perform linear interpolation between two pixels for each row or column intersected by the projection line, as described in Joseph P, “ An improved algorithm for reprojecting rays through pixel images ,” IEEE Trans. Med. Imaging 1, pp. 192-196 (1982). 
     A third technique, distance-driven backprojection, maps the boundaries of the pixels and the detector elements onto a common axis (referred to as centering), and the amount of overlap is used as the backprojection weight, as shown in  FIG. 3 . Rays  31  from source  30  strike array  32 . See U.S. Pat. No. 7,227,982 and De Man, Bruno, and Basu, Samit, “ Distance - driven projection and backprojection in three dimensions ,” Phys Med Biol 49, pp. 2463-2475 (2004). 
     The pixel  33  is assumed to be square with a half-width δ, and the edges of the square pixel are in the x-ray beam  31  from source  30  are projected to points x pc1  and x pc2  on the common centering x axis. Similarly, the boundaries of detectors c are projected to positions xcc on the x axis. In  FIG. 3 , channel boundaries c s  to c e+1  are projected to positions xcc cs  to xcc ce+1  on the centering axis. The interpolation kernel width is the overlap of x pc1  to x pc2  with xcc cs  to xcc ce+1 : 
                     V   ⁡     (     i   ,   j     )       =           f   iscc     ·     D   ⁡     [     c   s     ]           (       ycc   ks     -     ycc     ks   +   1         )       +         f   iecc     ·     D   ⁡     [     c   e     ]           (       ycc   ke     -     ycc     ke   +   1         )       +       ∑       c   ′     =       c   s     +   1           c   e     -   1       ⁢     D   ⁡     [     c   ′     ]                   (   3   )               
The advantage of this method is that the projection of the pixel onto the detector is used to determine the interpolation kernel width, yielding improved SNR at large pixel sizes and improved resolution at smaller pixel sizes. The disadvantage to this method is that both the channels and pixels must be projected onto a centering axis (typically the y or x axis), which makes processing more complex. Furthermore, the resulting interpolation kernel width is a function of view.
 
     SUMMARY OF THE INVENTION 
     The present invention, in one embodiment, is directed to a computed-tomography method including determining a channel position located on an x-ray detector using a linear path from an x-ray source through a backprojection pixel, and determining an interpolation kernel width using a distance from the source to the backprojection pixel. 
     In another embodiment, the present invention is directed to an apparatus having an x-ray source, an x-ray detector having channels and segments, a signal processing unit connected to receive data collected the detector, and configured to determine a channel position on the detector located on a linear path from the x-ray source through a backprojection pixel and determine an interpolation kernel width using a distance from the source to the backprojection pixel. 
     The present invention may also be directed to a computer-readable medium containing instructions, wherein the instructions, when executed by a computer, perform a method including determining a channel position located on a linear path from an x-ray source through a backprojection pixel, and determining an interpolation kernel width using a distance from the source to the backprojection pixel. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       A more complete appreciation of the invention and many of the attendant advantages thereof will be readily obtained as the same becomes better understood by reference to the following detailed description when considered in connection with the accompanying drawings, wherein: 
         FIG. 1  is a diagram illustrating conventional pixel-driven backprojection; 
         FIG. 2  is a diagram illustrating conventional ray-driven backprojection; 
         FIG. 3  is a diagram illustrating conventional distance-driven backprojection; 
         FIG. 4  is a diagram of a system according to the invention; 
         FIG. 5  is a diagram illustrating possible axial spatial orientations for adjacent views with resulting interpolated views; 
         FIG. 6  is a diagram illustrating backprojection according to the invention; 
         FIG. 7  is a diagram illustrating calculation of the backprojection value; 
         FIG. 8  is a graph illustrating backprojection according to the invention for different field of views; 
         FIG. 9  is a graph illustrating a comparison between backprojection according to the invention and distance-driven backprojection at different field of views; 
         FIGS. 10A-10E  are graphs illustrating a comparison between backprojection according to the invention and distance-driven backprojection at different field of views; and 
         FIG. 11  is a diagram illustrating a second embodiment of the invention. 
     
    
    
     DETAILED DESCRIPTION 
       FIG. 4  shows an x-ray computed tomographic imaging device according to the present invention. The projection data measurement system constituted by gantry  1  accommodates an x-ray source  3  that generates a cone-beam of x-ray flux approximately cone-shaped, and a two-dimensional array type x-ray detector  5  consisting of a plurality of detector elements arranged in two-dimensional fashion, i.e., a plurality of elements arranged in one dimension stacked in a plurality of rows. X-ray source  3  and two-dimensional array type x-ray detector  5  are installed on a rotating ring  2  in facing opposite sides of a subject, who is laid on a sliding sheet of a bed  6 . Two-dimensional array type x-ray detector  5  is mounted on rotating ring  2 . Each detector element will correspond with one channel. X-rays from x-ray source  3  are directed on to subject through an x-ray filter  4 . X-rays that have passed through the subject are detected as an electrical signal by two-dimensional array type x-ray detector  5 . 
     X-ray controller  8  supplies a trigger signal to high voltage generator  7 . High voltage generator  7  applies high voltage to x-ray source  3  with the timing with which the trigger signal is received. This causes x-rays to be emitted from x-ray source  3 . Gantry/bed controller  9  synchronously controls the revolution of rotating ring  2  of gantry  1  and the sliding of the sliding sheet of bed  6 . System controller  10  constitutes the control center of the entire system and controls x-ray controller  8  and gantry/bed controller  9  such that, as seen from the subject, x-ray source  3  executes so-called helical scanning, in which it moves along a helical path. Specifically, rotating ring  2  is continuously rotated with fixed angular speed while the sliding plate is displaced with fixed speed, and x-rays are emitted continuously or intermittently at fixed angular intervals from x-ray source  3 . The source may also be scanned circularly. 
     The output signal of two-dimensional array type x-ray detector  5  is amplified by a data collection unit  11  for each channel and converted to a digital signal, to produce projection data. The projection data output from data collection unit  11  is fed to processing unit  12 . Processing unit  12  performs various processing using the projection data. Unit  12  performs interpolation, backprojection and reconstruction, as described in more detail below. Unit  12  determines backprojection data reflecting the x-ray absorption in each voxel. In the helical scanning system using a cone-beam of x-rays, the imaging region (effective field of view) is of cylindrical shape of radius o) centered on the axis of revolution. Unit  12  defines a plurality of voxels (three-dimensional pixels) in this imaging region, and finds the backprojection data for each voxel. The three-dimensional image data or tomographic image data compiled by using this backprojection data is sent to display device  14 , where it is displayed visually as a three-dimensional image or tomographic image. 
     The basic geometry of the backprojection according to the invention is illustrated in  FIG. 5 . The backprojection according to the invention may by used with helical and circular scanning, and with cone-beam or fan beam x-ray CT. It incorporates the pixel-driven attribute of tracing a line  51  from the source  50  through the center of the pixel  53  (having radius P) to the detector  52  to determine the backprojection position Ch i,j , and adaptively determines the interpolation kernel width ΔCh as a function of the source to pixel distance L assuming the pixel has a circular shape. The backprojection according to the invention does not require projection of pixels and detector boundaries onto a centering axis, yet retains the advantage of distance-based backprojection; i.e., improved signal-to-noise ratio (SNR) at large pixel sizes and improved resolution at smaller pixel sizes. 
     For a 2-D implementation, the adaptive kernel is applied only in the channel direction, and bilinear interpolation in the segment direction. The backprojection according to the invention in the channel direction is shown in more detail in  FIG. 6 , with segments omitted for clarity. The pixel is assumed to have a circular shape with center at position x i,j , y i,j  and radius P: 
                   P   =       F   ⁢           ⁢   O   ⁢           ⁢   V       2   ·   MATRIX               (   4   )               
where FOV is the reconstruction field of view and MATRIX is the reconstructed image matrix size. The floating point channel position Ch i,j  (in units of channels) is determined in the same manner as for a pixel-driven method:
 
                     Ch     i   ,   j   ,   k       =       1     Δ   ⁢           ⁢   γ       ⁢       tan     -   1       ⁡     [           x     i   ,   j       ⁢   cos   ⁢           ⁢     β   k       +       y     i   ,   j       ⁢   sin   ⁢           ⁢     β   k           R   +       x     i   ,   j       ⁢   sin   ⁢           ⁢     β   k       -       y     i   ,   j       ⁢   cos   ⁢           ⁢     β   k           ]                 (   5   )               
where Δγ the angular channel width, β k  is the view angle at view k, and R is the source to isocenter distance, where the isocenter is the gantry rotation axis. The projection of the circular pixel determines the adaptive channel interpolation kernel width (in units of channels), which is given by:
 
                       Δ   ⁢           ⁢     Ch     i   ,   j   ,   k         =       2   ⁢           ⁢   P       Δ   ⁢           ⁢     γ   ·     L     i   ,   j   ,   k               ⁢     
     ⁢   where           (   6   )                   L     i   ,   j   ,   k       ⁡     (       x     i   ,   j       ,     y     i   ,   j       ,     β   k       )       =           (       R   ⁢           ⁢   sin   ⁢           ⁢     β   k       +     x     i   ,   j         )     2     +       (       R   ⁢           ⁢   cos   ⁢           ⁢     β   k       -     y     i   ,   j         )     2                 (   7   )               
Calculation of the channel interpolated backprojection value is illustrated in  FIG. 7  Here, c is an integer and represents the channel number that a backprojection ray intersects. In this example, it is assumed that the indexing of the channels is such that the first channel is located at c=0, and that a ray intersecting channel  0  exactly in the center will have floating point position 0.5 (i.e., possible range of positions for channel  0  is 0.0 to 0.9999 . . . ) Other indexing schemes may have the center of channel  0  at 0.0, with a possible range of positions from −0.5 to +0.49999 . . . . If a different indexing scheme is used, Eq. 5 will be accordingly altered by adding a shift amount.
 
                     PDI     CH     i   ,   j   ,   k         =     {           PD   ⁡     (   c   )               if   ⁢           ⁢     c   s       =       c   ⁢           ⁢   and   ⁢           ⁢     c   e       =   c                   1     Δ   ⁢           ⁢     Ch     i   ,   j   ,   k           ⁢     (         δ   s     ·     PD   ⁡     (     c   s     )         +         δ   e     ·   PD     ⁢     (     c   e     )       +       ∑       c   ′     =       c   s     +   1           c   e     -   1       ⁢     PD   ⁡     (     c   ′     )           )           otherwise                   (   8   )               
where PD( ) is raw projection data, PDI is interpolated value of projection data, c s , c, and c e  are integer starting, center, and ending channel indexes and
 
     
       
         
           
             
               
                 
                   λ 
                   = 
                   
                     
                       Δ 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         Ch 
                         
                           i 
                           , 
                           j 
                           , 
                           k 
                         
                       
                     
                     2 
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
             
               
                 
                   
                     f 
                     s 
                   
                   = 
                   
                     
                       Ch 
                       
                         i 
                         , 
                         j 
                       
                     
                     - 
                     λ 
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
             
               
                 
                   
                     f 
                     e 
                   
                   = 
                   
                     
                       Ch 
                       
                         i 
                         , 
                         j 
                       
                     
                     + 
                     λ 
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
             
               
                 
                   
                     c 
                     s 
                   
                   = 
                   
                     int 
                     ⁡ 
                     
                       ( 
                       
                         f 
                         s 
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   12 
                   ) 
                 
               
             
             
               
                 
                   
                     c 
                     e 
                   
                   = 
                   
                     int 
                     ⁡ 
                     
                       ( 
                       
                         f 
                         e 
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   13 
                   ) 
                 
               
             
             
               
                 
                   
                     δ 
                     s 
                   
                   = 
                   
                     
                       c 
                       
                         s 
                         + 
                         1 
                       
                     
                     - 
                     
                       f 
                       s 
                     
                   
                 
               
               
                 
                   ( 
                   14 
                   ) 
                 
               
             
             
               
                 
                   
                     δ 
                     e 
                   
                   = 
                   
                     
                       f 
                       e 
                     
                     - 
                     
                       c 
                       e 
                     
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
       
     
       FIG. 8  shows an example of the interpolation kernel ΔCh as a function of L for various FOV sizes. For standard pixel-driven backprojection, ΔCh is equal to 1.0. Smaller FOVs are typically used to visualize smaller structures, where a higher resolution is desirable. For small FOVs of 240 and less, ΔCh&lt;1, which provides a better resolution. For larger FOVs, resolution is usually not the issue. Here ΔCh&gt;1, which provides a better signal-to-noise ratio, since more channels are used in the interpolation as well as better aliasing artifact. 
       FIG. 9  shows a comparison of ΔCh between the present invention (solid lines) and distance driven backprojection (DDBPJ, dashed lines) at extremes of view angle β=0° and β=45° (which are extremes for DDBPJ for the centering process). These are shown as a function of angular pixel position φ at radii of various FOVs. Overall, ΔCh for the present invention is somewhat larger than DDBPJ for a given FOV and pixel position, but will have similar image quality. 
       FIGS. 8 and 9  show that the method according to the invention can provide similar performance in terms of images compared to the conventional methods, while the method according to the invention. 
     In the first embodiment, the method is implemented in the channel direction. In an alternate embodiment, the method can be implemented in the segment direction by assuming the voxel is circular in z and calculating the projection of the circle across the segments in the z direction, as shown in  FIG. 10 . Q is the pixel radius in z and is equal to the reconstruction slice spacing (though in other embodiments it can be calculated differently).  FIGS. 10A-E  show results for FOVs of 500, 400, 320, 240 and 180, respectively. The solid lines represent the present invention and the dashed lines show, for comparison, distance driven backprojection. 
     The floating point segment position is determined exactly the same as it is for the pixel-driven method: 
                     z   seg     =       D   ·     z   p       L             (   16   )               
where z seg  is the distance from the then center of the detector to the backprojection segment location, D is the source to detector distance and is z p  is the z distance of the reconstruction pixel from the source-detector plane at view k. For circular scanning this is independent of view and for helical it is a function of view. The floating point segment position Seg is given by
 
                   Seg   =       z   seg     w             (   17   )               
where w is the segment width at isocenter.
 
     The adaptive segment interpolation kernel width ΔSeg (in units of segments) given by: 
     
       
         
           
             
               
                 
                   
                     
                       Δ 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       Seg 
                     
                     = 
                     
                       
                         2 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           σ 
                           seg 
                         
                       
                       w 
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   where 
                 
               
               
                 
                   ( 
                   18 
                   ) 
                 
               
             
             
               
                 
                   
                     σ 
                     seg 
                   
                   = 
                   
                     
                       Q 
                       · 
                       cos 
                     
                     ⁢ 
                     
                         
                     
                     ⁢ 
                     θ 
                     ⁢ 
                     
                       
                         d 
                         seg 
                       
                       
                         d 
                         p 
                       
                     
                   
                 
               
               
                 
                   ( 
                   19 
                   ) 
                 
               
             
             
               
                 
                   
                     d 
                     seg 
                   
                   = 
                   
                     
                       
                         D 
                         2 
                       
                       + 
                       
                         z 
                         seg 
                         2 
                       
                     
                   
                 
               
               
                 
                   ( 
                   20 
                   ) 
                 
               
             
             
               
                 
                   
                     d 
                     p 
                   
                   = 
                   
                     
                       
                         L 
                         2 
                       
                       + 
                       
                         z 
                         p 
                         2 
                       
                     
                   
                 
               
               
                 
                   ( 
                   21 
                   ) 
                 
               
             
             
               
                 
                   θ 
                   = 
                   
                     
                       tan 
                       
                         - 
                         1 
                       
                     
                     ⁢ 
                     
                       
                         z 
                         p 
                       
                       L 
                     
                   
                 
               
               
                 
                   ( 
                   22 
                   ) 
                 
               
             
           
         
       
     
     Calculation of the segment backprojection value is calculated in the same manner as for the channels as in  FIG. 7  and Eq. 8, substituting segments for channels. 
     In a third embodiment, the method can be implemented for both the channel and segment directions at the same time. This is the 3-D case, where the circular pixel is now a spherical voxel with radius S. Equations (6-15) are used to calculate channel direction projection with P=S and equations (16-22) for segment direction, with Q=S. This embodiment can be used for full cone-beam detectors. 
     The invention may also be embodied in the form a computer-readable medium containing a stored program to cause a computer to carry out the various operations and functions described above. 
     Numerous other modifications and variations of the present invention are possible in light of the above teachings. This document and equations have been developed for a curved detector array. For example, a flat or other detector array shape can be implemented. It is therefore to be understood that within the scope of the appended claims, the invention may be practiced otherwise than as specifically described herein.