# Patent

## Title
Cardiac computed tomography methods and systems using fast exact/quasi-exact filtered back projection algorithms

## Abstract

The present invention provides systems, methods, and devices for improved computed tomography (CT). More specifically, the present invention includes methods for improved cone-beam computed tomography (CBCT) resolution using improved filtered back projection (FBP) algorithms, which can be used for cardiac tomography and across other tomographic modalities. Embodiments provide methods, systems, and devices for reconstructing an image from projection data provided by a computed tomography scanner using the algorithms disclosed herein to generate an image with improved temporal resolution.

## Background
CROSS-REFERENCE TO RELATED APPLICATIONS
This application claims priority to and the benefit of the filing date of U.S. Provisional Application No. 61/225,708, filed Oct. 28, 2009, which is incorporated by reference herein in its entirety.


STATEMENT OF GOVERNMENT INTEREST
This invention was made made with government support under contracts EB002667, EB004287 and EB007288 awarded by National Institutes of Health. The government has certain rights in the invention.


BACKGROUND OF THE INVENTION
1. Field of the Invention
The present invention provides systems, methods, and devices for improved computed tomography (CT). More specifically, the present invention includes methods for improved cone-beam computed tomography (CBCT) resolution using improved filtered back projection (FBP) algorithms, which can be used for cardiac tomography and across other tomographic modalities.
2. Description of Related Art
Cardiovascular diseases (CVDs) are pervasive (American Heart Association 2004). CVD is the number one killer in the western world. The cost of the health care for CVD is skyrocketing. In 2004, the estimated direct and indirect cost of CVD was $368.4 billion.
Coronary artery disease is a leading cause of death as a result of a myocardial infarct (heart attack) without any symptom. Tomographic equipment with high temporal resolution is needed in order to successfully perform a cardiac scan and understand the etiology and pathogenesis of CVD, such as high blood pressure, coronary artery diseases, congestive heart failure, stroke and congenital cardiovascular defects, as well as to develop effective prevention and treatment strategies. CT scanners are now considered instrumental for detecting early heart diseases and are a centerpiece of preventive cardiology programs.
Although there has been an explosive growth in the development of CT scanners for cardiac CT studies, the efforts are generally limited to regular heartbeats. When applying traditional CT algorithms for cardiac CT reconstruction, the cardiac images may be inaccurate or useless based on substantial motion blurring, especially seen in patients who have high and irregular heartbeats due to the fact that each projection sector covers a projection angular range of a substantial length. Within such an angular range, the heart moves appreciably, especially when it is not in a relative stationary phase. As a benchmark, a  ̃0.3 mm spatial resolution is routinely achieved in spiral CT of the temporal bone where the motion magnitude is much less than that of the heart (see M. Vannier and G. Wang, Spiral C T refines imaging of temporal bone disorders, Diagnostic imaging, vol. 15, p. 116-121, 1993 and G. Wang, et al., Design, analysis and simulation for development of the first clinical micro-CT scanner1, Academic Radiology, vol. 12, pp. 511-525, 2005, which is incorporated by reference herein in its entirety). Spatial resolution with cardiac CT is at best in the millimeter domain.
Over the last thirty years, computer tomography (CT) has gone from image reconstruction based on scanning in a slice-by-slice process to spiral scanning. From the mid-1980s to present day, spiral type scanning has become the preferred process for data collection in CT. Under spiral scanning, a table with the patient continuously moves through the gantry while the source in the gantry is continuously rotating about the table. At first, spiral scanning used a one-dimensional detector array, which received data in one dimension (a single row of detectors). Later, two-dimensional detectors, where multiple rows (two or more rows) of detectors sit next to one another, were introduced. In CT there have been significant problems for image reconstruction especially for two-dimensional detectors.
For three/four-dimensional (also known as volumetric/dynamic) image reconstruction from the data provided by a spiral scan with two-dimensional detectors, known groups of algorithms include: exact algorithms, quasi-exact algorithms, approximate algorithms, and iterative algorithms. While the best approximate algorithms are of Feldkamp-type, the state of the art of the exact algorithms is the recently developed Katsevich algorithm.
Under ideal circumstances, exact algorithms can provide a replication of a true object from data acquired from a spiral scan. However, exact algorithms can require a larger detector array, more memory, are more sensitive to noise, and run slower than approximate algorithms. Approximate algorithms can produce an image very efficiently using less computing power than exact algorithms. However, even under typical circumstances they produce an approximate image that may be similar to but still different from the exact image. In particular, approximate algorithms can create artifacts, which are false features, in an image. Under certain circumstances these artifacts can be quite severe.
To perform the long object reconstruction with longitudinally truncated data, the spiral cone-beam scanning mode and a generalized Feldkamp-type algorithm were proposed by Wang and others in 1991. However, the earlier image reconstruction algorithms for that purpose are either approximate or exact only using data from multiple spiral turns.
In 2002, an exact and efficient method was developed by Katsevich, which is a significant breakthrough in the area of spiral cone-beam CT. The Katsevich algorithm is in a filtered-backprojection (FBP) format using data from a PI-arc (scanning arc corresponding to the PI-line and less than one turn) based on the so-called PI-line and the Tam-Danielsson window. The principle is that any point inside the standard spiral or helical belongs to one and only one PI-line. Any point on the PI-line can be reconstructed from filtered data on the detector plane with the angular parameter from the PI-arc. In 2003, a slow FBP and a backprojected-filtration algorithm (BPF) were developed for helical cone-beam CT based on the Katsevich algorithm by exchanging the order of integrals. For important biomedical applications including application with movement present such as cardiac CT, generalization of the exact cone-beam reconstruction algorithms from the case of standard spirals to the case of nonstandard spirals and other scanning loci is desirable and useful. Although the current Katsevich-type algorithms are known for a standard spiral scan, there are no known fast algorithms, systems, devices and methods that can reconstruct an image exactly or quasi-exactly from data acquired in a CT scan with good temporal resolution.
Therefore, despite the impressive advancement of the CT technology, there are still unmet, critical and immediate needs such as those mentioned above for better image quality in many cardiac and other CT investigations wherein the motion magnitude is increased.
SUMMARY OF THE INVENTION
The numerous limitations inherent in the scanning systems described above provide great incentive for new, better systems and methods capable of accounting for one or more of these issues. If CTs are to be seen as an accurate, reliable therapeutic answer, then improved methods for reconstructing an image should be developed that can more accurately predict the image with improved temporal resolution and less artifacts.
The primary limitation to the above-mentioned system is its need to provide good temporal resolution and image reconstruction when movement is involved. However, as more complex applications for scanning are encountered, reconstruction of key subject areas such as the heart, lung, head and neck is cumbersome at best and may be inadequate to develop reliable diagnosis and therapies. Therefore, a more-advanced system that allows for the production of better object reconstruction would be ideal. This allows for the adaptation of exact and quasi-exact algorithms to provide better images.
Accordingly, embodiments of the invention provide methods, systems, and devices for reconstructing an image from projection data provided by a computed tomography scanner comprising: scanning an object in a cone-beam imaging geometry following a general triple helix path wherein projection data is generated; reconstructing the image, wherein the reconstructing comprises performing a filtered backprojection; using a fast exact or quasi-exact filtered back projection algorithm to generate the backprojected data; and using the backprojected data to generate an image with improved temporal resolution. Preferably, embodiments of the invention provide images with less than about 500 ms, e.g., about 100 ms temporal resolution or less, such as about 80 ms or less, or about 60 ms or less, or about 50 ms or less, or about 30 ms or less, or even about 10 ms or less, and so forth.
In the context of this disclosure, exact or quasi-exact means that the algorithm is theoretically exact for a good portion of voxels in the object or theoretically exact if a practically insignificant portion of data could be handled in a more complicated fashion. Said another way, quasi-exact means that the algorithm is derived from an exact three-dimensional reconstruction approach, in which deviations from exactness are introduced which are sufficiently small and lead to minor artifacts, but result in a numerically efficient algorithm. By way of example, these deviations may lead to inexact weighting of a small percentage of Radon planes at every voxel.
In preferred embodiments, the temporal resolution may be in the range of about 100 ms to less than about 10 ms.
The present invention includes a computed tomography (CT) imaging method comprising: scanning an object using multi-source helical cone-beam computed tomography (CBCT) to acquire projection data relating to the object; and reconstructing the scanned portion of the object into an image by performing a computationally efficient filtered backprojection (FBP) and theoretically exact/quasi-exact algorithm to generate image data. Such methods, for example, are included within the scope of the invention which are capable of achieving a temporal resolution in the image in the range of <500 ms.
In the context of this disclosure, what is meant by “computationally efficient”, we mean that the computational efficiency of a cone-beam reconstruction algorithm is in the same order of magnitude of that developed by Katsevich as reported in (A. Katsevich, “Theoretically exact filtered backprojection-type inversion algorithm for spiral CT,” SIAM Journal on Applied Mathematics, vol. 62, p. 2012, 2002), which is herein incorporated by reference in its entirety. Further, by “theoretically exact”, we mean that the reconstructed image would be exact if both the sampling interval and the measurement error approach zero.
Embodiments of the invention further include a computed tomography (CT) imaging system comprising: a multi-source helical cone-beam computed tomography (CBCT) scanner operably configured for scanning an object to acquire projection data relating to the object; a processing module operably configured for reconstructing the scanned portion of the object into an image by performing a filtered backprojection (FBP) with a fast exact or quasi-exact FBP algorithm to generate image data; and a processor for executing the processing module. Such systems can include software and hardware operably configured for performing the functions of the processing module.
Embodiments of the present invention provide for reconstructing the image using a fast exact or quasi-exact algorithm developed by defining the weight function, determining filtering directions, calculating the backprojection coefficients, and reconstructing the object with, for example:




  
    f
    ⁡
    
      (
      x
      )
    
  
  =
  
     
    
      
        -
        
          1
          
            4
            ⁢
            
              π
              2
            
          
        
      
      ⁢
      
         
        
          
            ∫
            I
          
          ⁢
          
            
              ∑
              m
            
            ⁢
            
              
                
                  
                    c
                    m
                  
                  ⁡
                  
                    (
                    
                      s
                      ,
                      x
                    
                    )
                  
                
                
                  
                  
                    x
                    -
                    
                      y
                      ⁡
                      
                        (
                        s
                        )
                      
                    
                  
                  
                
              
              ×
              
                 
                
                  
                
                ⁢
                
                   
                  
                    
                      ∫
                      0
                      
                        2
                        ⁢
                        π
                      
                    
                    ⁢
                    
                      
                        ∂
                        
                          ∂
                          q
                        
                      
                      ⁢
                      
                        
                      
                      ⁢
                      
                        
                          D
                          f
                        
                        (
                        
                          
                            y
                            ⁡
                            
                              (
                              q
                              )
                            
                          
                          ,
                          
                            
                          
                          ⁢
                          
                            
                               
                              
                                
                                  cos
                                  ⁢
                                  
                                    
                                  
                                  ⁢
                                  
                                    γβ
                                    ⁡
                                    
                                      (
                                      
                                        s
                                        ,
                                        x
                                      
                                      )
                                    
                                  
                                
                                +
                                
                                  sin
                                  ⁢
                                  
                                    
                                  
                                  ⁢
                                  
                                    
                                      γα
                                      ⊥
                                    
                                    ⁡
                                    
                                      (
                                      
                                        s
                                        ,
                                        x
                                        ,
                                        
                                          
                                        
                                        ⁢
                                        
                                          θ
                                          m
                                        
                                      
                                      )
                                    
                                  
                                
                              
                              )
                            
                            ⁢
                            
                              |
                              
                                q
                                =
                                s
                              
                            
                            ⁢
                            
                              
                                
                                  d
                                  y
                                
                                
                                  sin
                                  ⁢
                                  
                                    
                                  
                                  ⁢
                                  γ
                                
                              
                              ⁢
                              
                                
                                  d
                                  s
                                
                                .
                              
                            
                          
                        
                      
                    
                  
                
              
            
          
        
      
    
  




Such methods, systems, and devices can further be characterized in having the fast exact or quasi-exact algorithm implemented by differentiating each projection with respect to variable s; for each y1(s), i={1, 2, 3}, performing the Hilbert transform of derivative data along the given filtering directions on the corresponding detector plane; and backprojecting the filtered data on the inter-PI segments to reconstruct the object.
The features and advantages of the present invention will be apparent to those skilled in the art. While numerous changes may be made by those skilled in the art, such changes are within the spirit of the invention.



BRIEF DESCRIPTION OF THE DRAWINGS
These drawings illustrate certain aspects of some of the embodiments of the present invention, and should not be used to limit or define the invention.
FIGS. 1A-B are schematic diagrams showing geometry of triple-source helical CBCT. Three x-ray sources are rotated around the x3-axis along the helices y1(s), y2(s) and y3(s), respectively. The y1(s), y2(s) and y3(s) are on a cylinder of radius R. An object to be reconstructed is confined within a cylinder of radius r, where r<R. Parameter h denotes the pitch of each helix. The inter-helix distance along the x3-axis is h/3.
FIG. 2 is an illustration of the Zhao window bounded by solid lines Γ±1 and the Tam-Danielsson window bounded by dashed lines Γ±2. The detector plane is represented by the Cartesian coordinate system (u, v).
FIG. 3 is a schematic diagram of inter-PI arcs (thick solid curve-arcs).
FIG. 4 is a schematic diagram of the decomposition of the Zhao window into the regions G1, G2 and G3. Lit and Lib are the inflection lines at ŝu and ŝd, respectively.
FIG. 5A is a graphic representation allowing for the visualization of the domains delimited by the A-curves and T-curves on the surface of the unit sphere in spherical coordinates for x=(0.1, 0, 0).
FIGS. 5B, C, and D are graphic representations of the zoom-in versions of the areas bounded by the bottom left, bottom right, and top circles shown in FIG. 5A, respectively.
FIG. 6 is an illustration of the osculating plane πc,
FIGS. 7A and B are graphical representations of the close-up views of the diagram for x=(0, −0.15, 0) for 7A and x=(0, −0.3, 0) for 7B.
FIG. 8 is a graphical illustration of L-curves in the spherical coordinates θ1, θ2).
FIG. 9A is a graphic representation of the full diagram showing different regions split by A-, T- and L-curves for L-curves in x=(0.2, −0.3, 0).
FIGS. 9B, C, and D are graphic representations of the zoom-in versions of the regions bounded by the upper, bottom right, and bottom left circles shown in FIG. 9A.
FIG. 10 is a graphic representation of domains on the detector plane.
FIG. 11A is a graphical representation of the Bs-curve being tangent to a T-curve in D11 for the source on y1(s).
FIG. 11B is a graphical representation of the Bs-curve being tangent to a T-curve in D5 for the source on y1(s).
FIG. 11C is a graphical representation of the Bs-curve passing across the second T-curve for the source on y1(s).
FIG. 12A is a graphical representation of the determination of c0.
FIG. 12B is a graphical representation of the determination of c1.
FIGS. 13A-B are respectively graphical representations of filtering lines in the case of {circumflex over (x)}εG1∪G3 and {circumflex over (x)}εG2 for the first fast FBP algorithm.
FIG. 14 is a graphical representation showing that the required detector area is bounded by Γt, Γl, Lmax and Lmin for the first algorithm, and by Γt, Γl, L′max and L′min for the second algorithm.
FIG. 15 is a graphical representation of filtering lines for two fast FBP algorithms when {circumflex over (x)} is above where Lit, Lt2, and Lt1 are for the first and second algorithms, respectively.
FIGS. 16A and 16B are graphical representations illustrating the second fast FBP algorithm.
FIG. 17A is a reconstructed image of the Clock phantom with r=375 mm using the first fast FBP algorithm.
FIG. 17B is a reconstructed image of the Clock phantom with r=375 mm using the second fast FBP algorithm.
FIGS. 17C and D are images representing the differences between the reconstructed images in FIGS. 17A and 17B and the ground truth respectively in the display window [−0.5, 0.5].
FIG. 18 is a graphical representation showing projected inter-PI lines on the detector plane, where the thick curve segments denote the inter-PI arcs.
FIG. 19A is a graphical representation of a plot of EQUATION 20 with r0=0.495 R.
FIG. 19B is a graphical representation of a plot of Φs over a range of sε(0, 2π).
FIG. 19C is a graphical representation of a plot of ψ(s2s) over a range of 0≦s22−s0≦4.1773.
FIG. 20 is a graphical representation for possible locations of the “critical event” for Case 4.
FIG. 21 is a graphical representation illustrating regions G21 and G22.
FIG. 22A is a graphical representations of the relationship among the inter-PI line La, L-line, Ll and inflection line Lit (Lib) for {circumflex over (x)} in G2 and above ŝu.
FIG. 22B is a graphical representation of the relationship among the inter-PI line La, L-line, Ll and inflection line Lit (Lib) for {circumflex over (x)} in G2 and above ŝd.
FIG. 23 is a schematic representation of the angular transformation from a spherical coordinate system to a detector.



DETAILED DESCRIPTION OF VARIOUS EMBODIMENTS OF THE INVENTION
In accordance with embodiments of the present invention, a method of the present invention may comprise introducing two fast FBP algorithms for use with conventional cardiac CT technologies in order to obtain better reconstruction images. One of the many potential advantages of the methods of the present invention, only some of which are discussed herein, is that images with less blurring and improved temporal resolution may be obtained even when there is movement in the object being scanned. The current invention may provide benefits to various types of interior tomography including, but not limited to, cardiac, lung, head and neck tomography. In the medical field and in biomedical science, the methods disclosed herein may greatly reduce the production of unusable images and thereby potentially allow increased early detection of diseases, reduced amount of radioactive contrast used on the patients, and/or reduced costs associated with CTs. Better temporal resolution in the images may provide a cost savings by reducing the number of images needed to conclude a finding. This type of scanning may likewise provide more flexibility in designing experiments in small animals in order to better study these diseases and develop effective treatments.
Another potential advantage is that the two fast FBP proposed algorithms utilize the inter-PI line and inter-PI arcs, and have a shift-invariant filtering structure. Unlike our slow-FBP algorithm performing filtration spatial-variantly line by line, the proposed fast-FBP algorithms filter projection data spatial-invariantly view by view, representing a significant computational benefit. Since triple-source helical CBCT may triple temporal resolution, it seems a promising mode for cardiac CT and other CT applications, and our proposed algorithms may find applications in this context. The methods of the present invention allow for temporal resolution in the range of about 100 ms to less than about 10 ms.
Geometry of Triple-Source Helical CBCT.
In particular embodiments, the geometry of the triple-source helical. CBCT may be measured by allowing f(x) be an object function to be reconstructed. In embodiments where this function is smooth and vanishes outside the object cylinder EQUATION 1 may be applied as described below:

Ω={x=(x1,x2,x3)|x12+x22≦r2,x3min≦x3≦x3max},0<r<R,  (EQUATION 1)

where r is the radius of the object cylinder and R the radius of the scanning cylinder on which a scanning trajectory resides. In embodiments with the Cartesian coordinate system (x1, x2, x3), the triple-helix trajectories can be expressed as shown in EQUATION 2 below:




  
    
      
        {
        
          
            
              
                
                  
                    y
                    1
                  
                  ⁡
                  
                    (
                    s
                    )
                  
                
                =
                
                  (
                  
                    
                      R
                      ⁢
                      
                        
                      
                      ⁢
                      cos
                      ⁢
                      
                        
                      
                      ⁢
                      s
                    
                    ,
                    
                      R
                      ⁢
                      
                        
                      
                      ⁢
                      sin
                      ⁢
                      
                        
                      
                      ⁢
                      s
                    
                    ,
                    
                      
                        h
                        
                          2
                          ⁢
                          π
                        
                      
                      ⁢
                      s
                    
                  
                  )
                
              
            
          
          
            
              
                
                  
                    y
                    2
                  
                  ⁡
                  
                    (
                    s
                    )
                  
                
                =
                
                  (
                  
                    
                      R
                      ⁢
                      
                        
                      
                      ⁢
                      
                        cos
                        ⁡
                        
                          (
                          
                            s
                            +
                            
                              
                                2
                                3
                              
                              ⁢
                              π
                            
                          
                          )
                        
                      
                    
                    ,
                    
                      R
                      ⁢
                      
                        
                      
                      ⁢
                      
                        sin
                        ⁡
                        
                          (
                          
                            s
                            +
                            
                              
                                2
                                3
                              
                              ⁢
                              π
                            
                          
                          )
                        
                      
                    
                    ,
                    
                      
                        h
                        
                          2
                          ⁢
                          π
                        
                      
                      ⁢
                      s
                    
                  
                  )
                
              
            
          
          
            
              
                
                  
                    y
                    3
                  
                  ⁡
                  
                    (
                    s
                    )
                  
                
                =
                
                  (
                  
                    
                      R
                      ⁢
                      
                        
                      
                      ⁢
                      
                        cos
                        ⁡
                        
                          (
                          
                            s
                            +
                            
                              
                                4
                                3
                              
                              ⁢
                              π
                            
                          
                          )
                        
                      
                    
                    ,
                    
                      R
                      ⁢
                      
                        
                      
                      ⁢
                      
                        sin
                        ⁡
                        
                          (
                          
                            s
                            +
                            
                              
                                4
                                3
                              
                              ⁢
                              π
                            
                          
                          )
                        
                      
                    
                    ,
                    
                      
                        h
                        
                          2
                          ⁢
                          π
                        
                      
                      ⁢
                      s
                    
                  
                  )
                
              
            
          
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          2
        
        )
      
    
  




where h>0 is the pitch of each helix, and sε is the rotation angle. FIG. 1 illustrates the triple-source helical CBCT geometry.
Previously, the inter-helix PI-lines-were defined and extended the traditional Tam-Danielsson window to the Zhao window in the case of triple helices. The terms inter-helix PI-lines and inter-PI lines are the same and are used interchangeably throughout. Specifically, for each source position yj(s), jε {1, 2, 3}, the corresponding Zhao window is the region on the surface of the scanning cylinder bounded by the nearest helix turn of yjmod3+1(s) and the nearest helix turn of y(j+1)mod3+1(s), jε {1, 2, 3}. In FIG. 2, Γ±1 and Γ±2 denote the boundaries of the Zhao window and the Tam-Danielsson window on the detector plane, respectively. In certain embodiments, the algorithms described herein are designed for flat-panel detectors. However, in embodiments with detectors of other shapes, the arbitrary-shaped detector may be rebinned to a virtual flat-panel detector in a preprocessing step so that the algorithms of the current disclosure may be used.
The properties of the inter-PI lines and inter-PI arcs may be determined by recalling that an inter-PI line for yj(s) and yjmod3+1(s), jε{1, 2, 3}, is the line that (1) intersects yj(s) at one point and yjmod3+1(s) at another point; and (2) the absolute difference between the angular parameter values at the two intersection points is less than 2π. The existence and uniqueness of the inter-PI line is shown in Theorem 1 below.
Theorem 1 states that through any fixed xεΩ, there exists one and only one inter-PI line for any pair of the three helices defined by EQUATION 2. In the triple-helix case, there are three inter-PI lines for a fixed xεΩ and corresponding inter-helix PI-arcs whose end points may be along the corresponding helices and share the intersection points of the inter-PI lines. In some embodiments, the three inter-PI arcs represent the source trajectory arcs along which the sources illuminate the point x as shown in FIG. 3.
In 2003, Katsevich proposed a general scheme for constructing inversion algorithms for CBCT. It can be stated as follows in EQUATIONS 3-8:




  
    
      
        
          f
          ⁡
          
            (
            x
            )
          
        
        =
        
          
            
              -
              
                1
                
                  4
                  ⁢
                  
                    π
                    2
                  
                
              
            
            ⁢
            
              
                ∫
                I
              
              ⁢
              
                
                  ∑
                  m
                
                ⁢
                
                  
                    
                      
                        c
                        m
                      
                      ⁡
                      
                        (
                        
                          s
                          ,
                          x
                        
                        )
                      
                    
                    
                      
                      
                        x
                        -
                        
                          y
                          ⁡
                          
                            (
                            s
                            )
                          
                        
                      
                      
                    
                  
                  ×
                  
                    
                      ∫
                      0
                      
                        2
                        ⁢
                        π
                      
                    
                    ⁢
                    
                      
                        ∂
                        
                          ∂
                          q
                        
                      
                      ⁢
                      
                        
                          D
                          f
                        
                        ⁡
                        
                          (
                          
                            
                              y
                              ⁡
                              
                                (
                                q
                                )
                              
                            
                            ,
                            
                              
                                cos
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  γβ
                                  ⁡
                                  
                                    (
                                    
                                      s
                                      ,
                                      x
                                    
                                    )
                                  
                                
                              
                              +
                              
                                sin
                                ⁢
                                
                                  
                                
                                ⁢
                                
                                  
                                    γα
                                    ⊥
                                  
                                  ⁡
                                  
                                    (
                                    
                                      s
                                      ,
                                      x
                                      ,
                                      
                                        θ
                                        m
                                      
                                    
                                    )
                                  
                                
                              
                            
                          
                          )
                        
                      
                    
                  
                
              
            
          
          ⁢
          
            |
            
              q
              =
              s
            
          
          ⁢
          
            
              
                d
                y
              
              
                sin
                ⁢
                
                  
                
                ⁢
                γ
              
            
            ⁢
            
              d
              s
            
          
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          3
        
        )
      
    
  
  
    
      
        
          
        
        ⁢
        
          
            
              β
              ⁡
              
                (
                
                  s
                  ,
                  x
                
                )
              
            
            :=
            
              
                x
                -
                
                  y
                  ⁡
                  
                    (
                    s
                    )
                  
                
              
              
                
                
                  x
                  -
                  
                    y
                    ⁡
                    
                      (
                      s
                      )
                    
                  
                
                
              
            
          
          ,
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          4
        
        )
      
    
  
  
    
      
        
          
        
        ⁢
        
          
            
              
                α
                ⊥
              
              ⁡
              
                (
                
                  s
                  ,
                  x
                  ,
                  θ
                
                )
              
            
            :=
            
              
                β
                ⁡
                
                  (
                  
                    s
                    ,
                    x
                  
                  )
                
              
              ×
              
                α
                ⁡
                
                  (
                  
                    s
                    ,
                    x
                    ,
                    θ
                  
                  )
                
              
            
          
          ,
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          5
        
        )
      
    
  
  
    
      
        
          
        
        ⁢
        
          
            
              
                c
                m
              
              ⁡
              
                (
                
                  s
                  ,
                  x
                
                )
              
            
            :=
            
              
                lim
                
                  ɛ
                  ->
                  
                    0
                    +
                  
                
              
              ⁢
              
                (
                
                  
                    φ
                    ⁡
                    
                      (
                      
                        s
                        ,
                        x
                        ,
                        
                          
                            θ
                            m
                          
                          +
                          ɛ
                        
                      
                      )
                    
                  
                  -
                  
                    φ
                    ⁡
                    
                      (
                      
                        s
                        ,
                        x
                        ,
                        
                          
                            θ
                            m
                          
                          -
                          ɛ
                        
                      
                      )
                    
                  
                
                )
              
            
          
          ,
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          6
        
        )
      
    
  
  
    
      
        
          
        
        ⁢
        
          
            
              φ
              ⁡
              
                (
                
                  s
                  ,
                  x
                  ,
                  θ
                
                )
              
            
            :=
            
              
                sgn
                ⁡
                
                  (
                  
                    α
                    ·
                    
                      
                        y
                        .
                      
                      ⁡
                      
                        (
                        s
                        )
                      
                    
                  
                  )
                
              
              ⁢
              
                n
                ⁡
                
                  (
                  
                    s
                    ,
                    x
                    ,
                    α
                  
                  )
                
              
            
          
          ,
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          7
        
        )
      
    
  
  
    
      
        
          
        
        ⁢
        
          
            I
            :=
            
              
                
                  ⋃
                  
                    l
                    =
                    1
                  
                  
                    L
                    c
                  
                
                ⁢
                
                  [
                  
                    
                      a
                      i
                    
                    ,
                    
                      b
                      i
                    
                  
                  ]
                
              
              ->
              
                                3
              
            
          
          ,
          
            
              I
              ∋
              s
            
            ->
            
              
                y
                ⁡
                
                  (
                  s
                  )
                
              
              ∈
              
                                3
              
            
          
          ,
          
            

          
          ⁢
          
            
          
          ⁢
          
            
              
              
                
                  y
                  .
                
                ⁡
                
                  (
                  s
                  )
                
              
              
            
            ≠
            0.
          
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          8
        
        )
      
    
  




where Df(y, β) is the cone-beam transform of f, θ the polar angle in the plane perpendicular to β(s, x), α(s, x, θ) a unit vector perpendicular to β(s, x), θm a point where φ(s, x, θ) is discontinuous, n(s, x, α) a weight function, C a finite union of C∞ curves in 3, −∞<a1<b1<∞, and {dot over (y)}(s):=dy/ds.
The aforementioned general inversion formula can be applied to any trajectory that satisfies Tuy's condition, but only when the weight function n(s, x, α) is well designed can the inversion formula have a shift-invariant filtering structure. To derive fast exact FBP algorithms for triple-source helical CBCT, our general approach involves the following key concepts of and analyses on the inflection line, A-, T-, L-, and Bs-curves.
Inflection line. On the detector plane, the boundaries of the Zhao window may be expressed as EQUATION 9 below:




  
    
      
        
          
            u
            ⁡
            
              (
              s
              )
            
          
          =
          
            
              D
              ⁢
              
                
              
              ⁢
              sin
              ⁢
              
                
              
              ⁢
              s
            
            
              1
              -
              
                cos
                ⁢
                
                  
                
                ⁢
                s
              
            
          
        
        ,
        
          
            v
            ⁡
            
              (
              s
              )
            
          
          =
          
            
              Dh
              
                2
                ⁢
                π
                ⁢
                
                  
                
                ⁢
                R
              
            
            ⁢
            
              
                (
                
                  s
                  +
                  
                    Δ
                    ⁢
                    
                      
                    
                    ⁢
                    s
                  
                
                )
              
              
                (
                
                  1
                  -
                  
                    cos
                    ⁢
                    
                      
                    
                    ⁢
                    s
                  
                
                )
              
            
          
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          9
        
        )
      
    
  




where D is the distance between the detector and the source, s is the angular parameter relative to the corresponding source position, Δs=−2/3π and Δs=−4/3π are for the top and bottom boundaries respectively. Then, EQUATIONS 10-14 can be used.




  
    
      
        
          
        
        ⁢
        
          
            
              
                u
                .
              
              ⁡
              
                (
                s
                )
              
            
            =
            
              D
              
                
                  cos
                  ⁢
                  
                    
                  
                  ⁢
                  s
                
                -
                1
              
            
          
          ,
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          10
        
        )
      
    
  
  
    
      
        
          
        
        ⁢
        
          
            
              
                u
                 ̈
              
              ⁡
              
                (
                s
                )
              
            
            =
            
              
                D
                ⁢
                
                  
                
                ⁢
                sin
                ⁢
                
                  
                
                ⁢
                s
              
              
                
                  (
                  
                    
                      cos
                      ⁢
                      
                        
                      
                      ⁢
                      s
                    
                    -
                    1
                  
                  )
                
                2
              
            
          
          ,
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          11
        
        )
      
    
  
  
    
      
        
          
        
        ⁢
        
          
            
              
                v
                .
              
              ⁡
              
                (
                s
                )
              
            
            =
            
              
                Dh
                
                  2
                  ⁢
                  π
                  ⁢
                  
                    
                  
                  ⁢
                  R
                
              
              ⁢
              
                
                  [
                  
                    1
                    -
                    
                      cos
                      ⁢
                      
                        
                      
                      ⁢
                      s
                    
                    -
                    
                      
                        (
                        
                          s
                          +
                          
                            Δ
                            ⁢
                            
                              
                            
                            ⁢
                            s
                          
                        
                        )
                      
                      ⁢
                      sin
                      ⁢
                      
                        
                      
                      ⁢
                      s
                    
                  
                  ]
                
                
                  
                    (
                    
                      1
                      -
                      
                        cos
                        ⁢
                        
                          
                        
                        ⁢
                        s
                      
                    
                    )
                  
                  2
                
              
            
          
          ,
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          12
        
        )
      
    
  
  
    
      
        
          
        
        ⁢
        
          
            
              
                v
                 ̈
              
              ⁡
              
                (
                s
                )
              
            
            =
            
              
                Dh
                
                  2
                  ⁢
                  π
                  ⁢
                  
                    
                  
                  ⁢
                  R
                
              
              ⁢
              
                
                  
                    
                      
                        
                          
                            (
                            
                              
                                cos
                                ⁢
                                
                                  
                                
                                ⁢
                                s
                              
                              -
                              1
                            
                            )
                          
                          ⁢
                          
                            (
                            
                              s
                              +
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                s
                              
                            
                            )
                          
                          ⁢
                          cos
                          ⁢
                          
                            
                          
                          ⁢
                          s
                        
                        +
                      
                    
                  
                  
                    
                      
                        2
                        ⁢
                        sin
                        ⁢
                        
                          
                        
                        ⁢
                        
                          s
                          ⁡
                          
                            (
                            
                              
                                cos
                                ⁢
                                
                                  
                                
                                ⁢
                                s
                              
                              +
                              
                                s
                                ⁢
                                
                                  
                                
                                ⁢
                                sin
                                ⁢
                                
                                  
                                
                                ⁢
                                s
                              
                              +
                              
                                Δ
                                ⁢
                                
                                  
                                
                                ⁢
                                s
                                ⁢
                                
                                  
                                
                                ⁢
                                sin
                                ⁢
                                
                                  
                                
                                ⁢
                                s
                              
                              -
                              1
                            
                            )
                          
                        
                      
                    
                  
                
                
                  (
                  
                    
                      (
                      
                        1
                        -
                        
                          cos
                          ⁢
                          
                            
                          
                          ⁢
                          s
                        
                      
                      )
                    
                    3
                  
                  )
                
              
            
          
          ,
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          13
        
        )
      
    
  
  
    
      
        
          
            
              d
              2
            
            ⁢
            v
          
          
            d
            
              u
              2
            
          
        
        =
        
        ⁢
        
          
            
              
                
                  
                    v
                     ̈
                  
                  ⁡
                  
                    (
                    s
                    )
                  
                
                ⁢
                
                  
                    u
                    .
                  
                  ⁡
                  
                    (
                    s
                    )
                  
                
              
              -
              
                
                  
                    u
                     ̈
                  
                  ⁡
                  
                    (
                    s
                    )
                  
                
                ⁢
                
                  
                    v
                    .
                  
                  ⁡
                  
                    (
                    s
                    )
                  
                
              
            
            
              
                [
                
                  
                    u
                    .
                  
                  ⁡
                  
                    (
                    s
                    )
                  
                
                ]
              
              3
            
          
          =
          
            
              h
              
                2
                ⁢
                π
                ⁢
                
                  
                
                ⁢
                D
              
            
            ⁢
            
              
                (
                
                  s
                  +
                  
                    Δ
                    ⁢
                    
                      
                    
                    ⁢
                    s
                  
                  -
                  
                    sin
                    ⁢
                    
                      
                    
                    ⁢
                    s
                  
                
                )
              
              .
            
          
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          14
        
        )
      
    
  




The inflection point exists when




  
    
      
        d
        2
      
      ⁢
      v
    
    
      d
      
        u
        2
      
    
  
  =
  0.




Thus, we obtain su=2.6053 and sd=3.6779 when Δs=−2/3π and −4/3π. The slope of the tangent line at s can be computed as shown in EQUATION 15 below:





  
    
      
        
          
            
              
                
                  d
                  v
                
                
                  d
                  u
                
              
              =
              
              ⁢
              
                
                  
                    v
                    .
                  
                  ⁡
                  
                    (
                    s
                    )
                  
                
                
                  
                    u
                    .
                  
                  ⁡
                  
                    (
                    s
                    )
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                
                  h
                  ⁡
                  
                    [
                    
                      1
                      -
                      
                        cos
                        ⁢
                        
                          
                        
                        ⁢
                        s
                      
                      -
                      
                        
                          (
                          
                            s
                            +
                            
                              Δ
                              ⁢
                              
                                
                              
                              ⁢
                              s
                            
                          
                          )
                        
                        ⁢
                        sin
                        ⁢
                        
                          
                        
                        ⁢
                        s
                      
                    
                    ]
                  
                
                
                  2
                  ⁢
                  π
                  ⁢
                  
                    
                  
                  ⁢
                  
                    R
                    ⁡
                    
                      (
                      
                        
                          cos
                          ⁢
                          
                            
                          
                          ⁢
                          s
                        
                        -
                        1
                      
                      )
                    
                  
                
              
            
          
        
        
          
            
              =
              
              ⁢
              
                
                  h
                  
                    2
                    ⁢
                    π
                    ⁢
                    
                      
                    
                    ⁢
                    R
                  
                
                ⁢
                cos
                ⁢
                
                  
                
                ⁢
                
                  s
                  .
                
              
            
          
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          15
        
        )
      
    
  




Because cos su=cos sd=−0.8596, the slope is the same (−0.1368 h/R) at both inflection points. For practical medical applications, it is common that rFOV≦0.5 R, and a boundary limitation x12+x22≦r2 (r=0.495 R) may be included, which is shown as the vertical lines Γl and Γr in FIG. 4. Now, the inflection lines (the tangent lines at ŝd and ŝu, where ŝd and ŝu are the projection of y(j+1)mod3+1(sd) and yjmod3+1(su), jε{1, 2, 3} on the detector plane) and the boundary lines split the Zhao window into the following three regions: G1, G2 and G3. Only the points in G1 and G3 can have tangent lines with Γ±1.
A-Curve and T-Curve.
To construct an appropriate weight function, the understanding of how Radon planes intersect with the trajectories is important. The number of intersection points only changes when a Radon plane is tangent to the trajectory or contains one PI line/inter-PI line. Hence, if we find all such Radon planes, we can determine the distribution of the intersection points. Since each plane is uniquely determined by its normal vector, in the following sections we use unit vectors instead of the Radon planes. An A-curve consists of all unit vectors orthogonal to an inter-PI line. A T-curve consists of all unit vectors in EQUATION 16 as shown below:




  
    
      
        
          
            α
            ⁡
            
              (
              s
              )
            
          
          =
          
            ±
            
              
                
                  (
                  
                    x
                    -
                    
                      y
                      ⁡
                      
                        (
                        s
                        )
                      
                    
                  
                  )
                
                ×
                
                  
                    y
                    .
                  
                  ⁡
                  
                    (
                    s
                    )
                  
                
              
              
                
                
                  
                    (
                    
                      x
                      -
                      
                        y
                        ⁡
                        
                          (
                          s
                          )
                        
                      
                    
                    )
                  
                  ×
                  
                    
                      y
                      .
                    
                    ⁡
                    
                      (
                      s
                      )
                    
                  
                
                
              
            
          
        
        ,
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          16
        
        )
      
    
  




where s belongs to an inter-PI arc. Actually, the A-curve represents all Radon planes containing one inter-PI line, and the T-curve represents all Radon planes tangent to the trajectory. Since there are three inter-PI lines and three inter-PI arcs for a fixed x, there are accordingly three A-curves and three T-curves. The use of spherical coordinates (θ1, θ2) to describe these curves on the unit sphere is shown in EQUATION 17:

α=(cos θ1 sin θ2, sin θ1 sin θ2, cos θ1),−π≦θ1≦π,0≦θ2≦π.  (EQUATION 17)

With the identification (θ1, θ2)≅((θ1+π)mod2π, π−θ2), each a corresponds to a unique plane through x with the normal vector α.
As an example, the A-curves and T-curves of point x=(0.1, 0, 0) are illustrated in FIG. 5, where R=1, h=2π. T1, T2 and T3 stand for the T-curves corresponding to the inter-PI arcs ,  and  respectively. Similarly, A1, A2 and A3 are for the A-curves corresponding to the inter-PI lines s1ss2e, s2ss3e and ŝ{circumflex over (s1es3e)} respectively.
The A-curves and T-curves may divide the surface of the unit sphere into several connected domains, in each of which all the planes through x have the same number of intersection points (IPs) with the inter-PI arcs of x. Given an object point and one trajectory, the number of IPs only changes when a Radon plane is tangent to the trajectory or contains the endpoints of the trajectory. The A-curve represents all planes containing the endpoints of the trajectory, and the T-curve represents all planes tangent to the trajectory. If any Radon plane is chosen and rotated around one direction, the normal vector of this plane forms a curve on the unit sphere. Clearly, only when this curve intersects with the A-curve or T-curve does the number of IPs change. Thus, the A-curve and T-curve define the boundaries of different domains in which the number of IPs is constant. The distribution of IPs over the inter-PI arcs is listed in Table I. To determine the distribution of IPs, we first pick a vector α(θ1, θ2) in each domain, and then generate the plane through x and perpendicular to α(θ1, θ2), and compute numerically the number of IPs.







TABLE I







Distribution of IPs on Each of the Domains Delimited by the A-curves


and T-curves on the Surface of the Unit Sphere. The dash indicates


no IP on the corresponding inter-Pl arc.











   
   
   
















D1
1
—
—



D2
—
1
—



D3
—
—
1



D4
1
1
1



D5
2
1
—



D6
2
—
1



D7
1
2
—



D8
—
2
1



D9
1
—
2



D10
—
1
2



D11
3
1
1



D12
1
3
1



D13
1
1
3










By construction, a T-curve always starts from an A-curve and ends on another A-curve. It can be seen from FIG. 5 that a T-curve is may not be smooth at some point ac, but the limits of unit tangent vectors at ac− and ac+ are equal Such a point ac is defined herein as a “cusp”. The term cusp indicates that the two vectors determine the same plane, and ac is the normal vector to that plane. It has been proved in that the cusp is equivalent to the osculating plane Πc(x) which goes through y1(s1c(x)) iε{1, 2, 3}, is parallel to {dot over (y)}1(s1c(x)), ÿ1(s1c(x)) and contains x (see FIG. 6). On the detector plane, this corresponds to a point where the projected boundary has a point of zero curvature, i.e., the point of inflection.
The diagram plotted in spherical coordinates deforms smoothly as a function of x. The new diagram is equivalent to the old one in the embodiments where the distribution of IPs remain the same. An essential change could happen when three boundaries intersect each other at one point, which is defined herein as a “critical event”. The term “critical event” may happen in the following seven cases:
1. Three A-curves intersect at one point;
2. Three T-curves intersect at one point.
3. Two A-curves and one T-curve intersect at one point;
4. One A-curve and two T-curves intersect at one point;
5. A T-curve becomes tangent to an A-curve at a point of non-smoothness (i.e., cusp);
6. When the order of tangency (i.e., the zero derivatives of this order) at the beginning of the T-curve is increased, the T-curve re-emerges on the other side of the A-curve;
7. T-curve develops a smooth dent and becomes-tangent to an A-curve.
From Lemmas 2-3 described below in the Examples section, it is shown that Cases 1, 2, 4, 6, 7 do not occur for r<0.495 R and Case 5 does not occur for r<0.265 R. In some embodiments, Case 3 is possible. In embodiments where Case 3 takes place, the tangency of T-curve and A-curve will move across another A-curve, then one domain disappears. For example, when x=(0.1, 0, 0) gradually changes to x=(0, −0.15, 0), in FIG. 5D the tangency of T3 and A2 will move across A3, and domain D10 disappears (see FIG. 7A). In other embodiments, Case 5 is possible for r≧0.265 R. In these embodiments, a T-curve will only intersect A-curves at the endpoints. That is, the cusp of that T-curve and one domain disappear (see FIG. 7B).
L-Curve.
The L-curve may be used to split the domain D4 into sub-domains, making the weight function n continuous across all the A-curves. This is the key requirement, which may allow for the development of efficient reconstruction algorithms. Thus, the L-curve should not go across an A-curve. In embodiments where x is fixed and s is run over the three inter-PI arcs, {circumflex over (x)} forms a trajectory on the detector plane. Because at the endpoint of the inter-PI arc the line connecting yi(s) and x happens to be an inter-PI line, {circumflex over (x)} always starts from one endpoint of the inter-PI arc on the boundary of the Zhao window, and ends at the other one. Hence, whatever the trajectory of {circumflex over (x)} is, part of the trajectory is in G2. In other words, {circumflex over (x)} will run across one inflection line, then move in G2, and finally cross the other inflection line. Note that {circumflex over (x)} on the inflection line indicates a plane containing the inflection line, i.e., a cusp in one T-curve. From Lemma 3 described in the Example section below, the cusps always belong to the boundary of domain D4. Thus, they can be used as the endpoints of the L curve. A family of L-curves is formed as follows. Run s over the three inter-PI arcs of x. If {circumflex over (x)} is in G2 and above ŝu where Γ+1 intersects Lit (FIG. 4), find the plane through {circumflex over (x)} and ŝu. If {circumflex over (x)} is in G2 and below ŝd where Γ−1 intersects Lib, find the plane through {circumflex over (x)} and ŝd. If {circumflex over (x)} is in G2 and between ŝd and ŝu, find the plane through {circumflex over (x)} and parallel to the u-axis. A plot of all the normal unit vectors of these planes in the spherical coordinates (θ1, θ2) may then be constructed. This gives us three L-curves. The corresponding lines on the detector plane are called L-lines. FIG. 8 shows the L-curves on the diagram in spherical coordinates (θ1, θ2), where L1, L2 and L3 denote the L-curves corresponding to the inter-PI arcs ,  and  respectively. As is seen from the above construction, the L-curve always starts and ends on the cusps, and not defined for those parameter values when {circumflex over (x)} is not in G2.
For r≧0.265 R, one or more cusps will disappear if “critical event Case 5” occurs, then the L-curve may start from the intersection of T-curve and A-curve, and end at one A-curve. Also, the L-curve may start from one cusp and end at one A-curve or start from the intersection of T-curve and A-curve, and end at one A-curve. For example, see FIG. 9, L1 starts from the intersection of T3 and A2, and ends at the intersection of T2 and A2; L2 starts from one point on A2, and ends at the cusp of T1; L3 starts from the cusp of T1, and ends at one point on A2.
Whatever the endpoints of L-curves are, the L-curves intersects at one point, for example, in some embodiments, θ2=π or 0 in the spherical coordinates (which corresponds to the plane containing x and parallel to the x1-x2 plane). Then D4 is split into several sub-domains. If the endpoints of L-curves are cusps, by Lemma 5, each sub-domain contains only one A-curve. If not, small “line segments” on A-curves may appear and the sub-domains may contain more than one A-curve (see FIG. 9B).
Bs-curve. A Bs-curve may consist of all unit vectors perpendicular to x-yi(s), iε{0, 2, 3}. Each intersection of Bs- and A-curves corresponds to a plane containing an inter-PI line and yi(s). Each intersection of B, and T-curves corresponds to a plane tangent to an inter-PI arc and containing yi(s). For example, in certain embodiments, one may choose xεΩ with {circumflex over (x)}εG2, {circumflex over (x)} is above L0, where L0 is the projection of the helical tangent at the current position. If L(θ):=DP(s)∩Π(x, α(s, x, θ)) is denoted, where DP(s) is the detector plane corresponding to the source position s, and L(θ) is the projection of the plane through x with the normal vector α(s, x, θ) (FIG. 10), then as θ increases, α(s, x, θ)εβL(s, x) rotates clockwise on DP(s), and the following sequence of events takes place. First, Π(x, α(s, x, θ)) intersects ŝ1sŝ2c, and a pair of IPs is born. On the unit sphere, this is seen as an intersection of Bs and A1, after which Bs enters D5 (FIG. 11a). Second, Π(x, α(s, x, θ)) intersects ŝ3sŝ1c and another pair of IPs is born. On the unit sphere, this is seen as an intersection of Bs and A3, after which Bs enters D11. Third, a swap of two IPs takes place. On DP(s) this happens when θ=θ0, L(θ0) is parallel to the helical tangent. On the unit sphere, this means that Bs is tangent to T1 at α0=α(s, x, θ0). Fourth, Bs exits D11 by intersecting T1. On DP(s), this takes place when L(θ) is tangent to Γ+2. Finally, Π(x, α(s, x, θ)) intersects the L-line. This will not change the number of IPs but it will be useful for construction of the weight function. On the unit sphere, this is seen as an intersection of Bs and L1.
The jumps across an A-curve can only be of two types: from a 1-IP domain to a 3-IP domain and from a 3-IP domain to a 5-IP domain. Note that the Bs-curve is tangent from the inside to T1, which means a swap of two intersection points at α=α0 where sgn(α0·{dot over (y)}(s))=0 (see [15]). For a fixed s, if x is allowed to change slightly inside the Zhao window, the tangency point will move from D11 to D5 across A3 (or from D11 to D6 across A1) (FIGS. 12-13). If x projects into G1 or G3, the Bs-curve will pass not only through D5 (or D6) and D11 but also through D7 and D12 (FIG. 14). The similar results can be obtained if the source is on y2(s) or y3(s).
Two filtered-backprojection algorithms for triple-source helical cone-beam CT can be used to obtain images having higher temporal resolution. The first exemplary algorithm uses two families of filtering lines, which are parallel to the tangent of the scanning trajectory and the so-called L lines. The second algorithm uses two families of filtering lines tangent to the boundaries of the Zhao window and L lines, respectively, but it eliminates the filtering paths along the tangent of the scanning trajectory, thus reducing the detector size greatly. Additional information concerning these algorithms can be found in Lu, Yang, et al., “Fast Exact/Quasi-Exact FBP Algorithms for Triple-Source Helical Cone-Beam CT,” IEEE Transactions on Medical Imaging, Vol. 29, No. 3, March 2010, which is incorporated by references herein in its entirety.
First Fast FBP Algorithm.
In order to design an algorithm for triple-source helical CBCT useful in cardiac CTs and other CTs where movement exists, the weight function n(s, x, α) must be specified. The filtering directions by the discontinuities of φ(s, x, θ):=sgn(α·{dot over (y)}(s))n(s, x, α) must also be determined. Following the determination of the filtering directions, the backprojection coefficients can be calculated according to EQUATION 6. Once the filtering lines and the backprojection coefficients are determined, EQUATION 3 may be used to reconstruct the object.
In order to construct the weight function n(s, x, α) one should know the following. In certain embodiments, in order to have an efficient FBP structure, the weight function n(s, x, α) should be continuous across all A-curves. Thus, the weight function can be defined as shown in Table II. The values in Table II are the weights assigned to IPs. For example, in the D1 domain the Radon plane has only one IP on the inter-PI segment . Accordingly, a weight of 1 may be assigned to this IP and a dash used to indicate that there is no IP on the inter-PI segments  and . In the D11 domain the Radon plane has three IPs on , one IP on  and one IP on . Thus, a weight of −1 may be assigned to two IPs on  and a weight of 1 to all other IPs.







TABLE II







Weight Function for the First Fast FBP Algorithm











   
   
   
















D1
+1
—
—



D2
—
+1
—



D3
—
—
+1



D5
+1, 0
0
—



D6
+1, 0
—
0



D7
0
+1, 0
—



D8
—
+1, 0
0



D9
0
—
+1, 0



D10
—
0
+1, 0



D11
+1, 0, 0
0
0



D12
0
+1, 0, 0
0



D13
0
0
+1, 0, 0



D41
+1
0
0



D42
0
+1
0



D43
0
0
+1










To find the backprojection coefficients, a representative point in each area is selected in order to determine the discontinuities of φ(s, x, θ):=sgn(α·y(s))n(s, x, α) and extend the results by continuity to the entire area. A discontinuity of sgn occurs only when a Bs-curve is tangent to a T-curve from the inside. In other words, the Bs-curve does not go across the T-curve, but stays on one side in a neighborhood of the point of tangency. On the detector plane, L(θ) is parallel to the helical tangent for all {circumflex over (x)}εG1∩G2∩G3. Hence, this gives a family of filtering lines parallel to L0, where L0 is the projection of the helical tangent at y(s). The swap of two IPs changes the weight from n=0 to n=1. The backprojection coefficient is computed as c0=φ+−φ−=(+1)(0)−(−1)(1)=1 (FIG. 12A).
In certain-embodiments, by construction, the weight function n is continuous across all inter-PI lines. A discontinuity of n occurs only when a Bs-curve intersects a T-curve or an L-curve. Without loss of generality or wishing to be limited by theory, choosing y1(s0) on . For {circumflex over (x)}εG2, after the swap mentioned in the above paragraph the weight at the current position is zero. Hence, when the Bs-curve passes through a T-curve, i.e., from D11 to D4, n is continuous. Possible jumps of n may only occur when a Bs-curve passes through an L-curve, i.e., from D43 to D41 or from D42 to D41 in FIG. 8. On the detector plane, this occurs when L(θ) overlaps the L-line of {circumflex over (x)}. Then, the backprojection coefficients may be computed as c1=φ+−φ−=(+1)(1)−(+1)(0)=1 (FIG. 16.). For {circumflex over (x)}εG1∩G3, the Bs-curve will not enter D41. Instead, it passes through a second T-curve twice, i.e., from D43 to D12 and from D7 to D1. From Table II, the jumps of n may only occur in the latter intersection. On the detector plane, this happens when L(θ) overlaps the line tangent to Γ±1. Then, the backprojection coefficients may be computed as ct=φ+−φ−=(+1)(1)−(+1)(0)=1 (FIG. 15).
FIG. 13A and FIG. 13B summarize the filtering lines and the backprojection coefficients discussed above. In these figures, L0′ is the line parallel to L0 and L1 denotes the L-line. To implement the proposed algorithm, the filtering lines cannot be truncated. Thus, the detector size should be large enough to cover the area bounded by Γr, Γt, Lmax and Lmin, where Lmax and Lmin are the lines across the intersections of (1) Γ1 and Γ+2 and (2) Γ, and Γ−2 respectively, and parallel to L0 (FIG. 14). In certain embodiments, the required detector area can be determined by two factors: (1) the ratio of the pitch h and the scanning radius R and (2) the ratio of the object support radius r and the scanning radius R. If R is fixed, the required detector area grows as h or r increases.
Second Fast FBP Algorithm.
Again, the design of the second fast FBP algorithm starts with specifying new weights (Table III). By construction, n(s, x, α) is continuous across all inter-PI lines. More importantly, a swap of two IPs takes place when a Bs-curve becomes tangent to a T-curve, and n(s, x, α) changes from +1 to −1. The discontinuity of sgn(α·ŷ(s)) appears only when a Bs-curve is tangent to a T-curve from inside. Since both n(s, x, α) and sgn(α·{dot over (y)}(s)) are discontinuous at that point, the function φ(s, x, θ):=sgn(α·{dot over (y)}(s))n(s, x, α) is continuous. Thus, the filtering operation along the tangent of the scanning trajectory is eliminated.







TABLE III







Weight Function for the Second Fast FBP Algorithm











   
   
   
















D1
+1
—
—



D2
—
+1
—



D3
—
—
+1



D5
+1, −1
+1
—



D6
+1, −1
—
+1



D7
+1
+1, −1
—



D8
—
+1, −1
+1



D9
+1
—
+1, −1



D10
—
+1
+1, −1



D11
−1, +1, −1
+1
+1



D12
+1
−1, +1, −1
+1



D13
+1
+1
−1, +1, −1



D41
+1
0
0



D42
0
+1
0



D43
0
0
+1










A discontinuity of n can occur only when a Bs-curve intersects a T-curve or an L-curve. Follow the discussion in Section IV, jumps of n may occur when (1) a Bs-curve passes through a T-curve, i.e., from D11 to D4 or from D4 to D12 in FIGS. 11A and 11C, and (2) a Bs-curve passes through an L-curve, i.e., from D42 to D41 or from D43 to D41 in FIG. 8. On the detector plane, this gives two families of filtering lines: the lines tangent to Γ±2 or Γ±1 and the L-lines. Note that the filtering lines tangent to γ±1 are different from those for our first fast FBP algorithm (FIG. 15), because the discontinuity of n(s, x, α) occurs on the different side of the cusp. Then, the backprojection coefficients can be calculated a ct=φ+−φ−=(+1)(0)−(+1)(−1)=1 and c1=φ+−φ−=(+1)(1)−(+1)(0)=1.
The reconstruction formula for the second algorithm is the same as that for the first algorithm. The only difference lies in the selection of the filtering lines. For clarity, our second fast FBP algorithm is illustrated in FIGS. 16A-B. Because the filtering paths along the tangent of the scanning trajectory are eliminated, the required detector area is reduced by at least 30% (FIG. 14).
By Lemma 3 described below in the Examples section, there are two types of “line segments” according to different critical events. First, let us consider the “line segment” related to a critical event in Case 3. Recall that before entering D4 the Bs-curve will be tangent to a T-curve. For the first algorithm, at the tangency the weight n changes from 1 to 0, then it does not change whether the Bs-curve enters D4 across an A-curve or a T-curve. For the second algorithm, if the weight n changes from 1 to −1 at the tangency, then it will jump from −1 to 0 when the Bs-curve enters D4. If the Bs-curve enters D4 across an A-curve (i.e., the “line segment”; see FIG. 7), the FBP structure is ruined. Thus, the critical event in Case 3 will only affect the second algorithm, without damaging the FBP structure of the first algorithm. Then, let us consider the “line segment” related to a critical event in Case 5, which only occurs when r≧0.265 R. From the discussion in Section III.D such “a line segment” is the boundary of the one-IP and three-IP domains. Hence, for both the algorithms the weight n will jump from 1 to 0 if the Bs-curve enters D4 across the “line segment”, and the FBP structure is ruined. Consequently, the first algorithm is theoretically exact for r<0.265 R and not exact for 0.495 R≧r≧0.265 R, and the second algorithm is not exact for r≦0.495 R.
Since the algorithms are not always exact, it may be important to estimate what percentage of the Radon planes is incorrectly calculated. If the Radon planes with approximate weighting only have a small percentage, the algorithms can be considered quasi-exact, and we can still reconstruct with high image quality.
First, one must consider the incorrectly weighting planes caused by the critical event in Case 3. It appears in the area r<0.495 R. Let us fix x for r<0.495 R, denote the intersection of the Radon plane and the detector plane as L(θ), θε[0, 2π], run s over the three inter-PI arcs, and see what happens with {circumflex over (x)} and L(θ). Based on the discussion in Section III. E, for {circumflex over (x)} in G2, if the critical event occurs, the Bs-curve will first intersect a T-curve, and then go cross an A-curve. For example, in FIG. 12 the Bs-curve will first intersect T1, then enter D4 across A3. On the detector plane, this corresponds to that L(θ) intersects the tangent of Γ+2 before the inter-PI line ŝ1cŝ3s while L(θ) is rotated clockwise. Therefore, the Radon planes between the tangent of Γ+2 and ŝ1cŝ3s are not exactly weighted. Because the slope of ŝ1cŝ3s is positive and the slope of the tangent of Γ+2 is less than h/2πR, the percentage of the incorrectly weighted Radon planes is less than




  
    ρ
    =
    
      β
      π
    
  
  ,
  
    
      where
      ⁢
      
        
      
      ⁢
      α
    
    <
    β
    <
    
      2
      ⁢
      
        arctan
        (
        
          
            2
            
              3
            
          
          ⁢
          tan
          ⁢
          
            α
            2
          
        
        )
      
    
  
  ,
  
    α
    =
    
      arctan
      ⁢
      
        h
        
          2
          ⁢
          π
          ⁢
          
            
          
          ⁢
          R
        
      
    
  




(See Appendix). It is common that h/R<0.2 in practical applications, hence ρ<1.17%. For {circumflex over (x)} in G1 (G3), if the critical event occurs, the Bs-curve will first go across an A-curve, and then over a T-curve. On the detector plane, this corresponds to the case when L(θ) intersects the tangent to Γ+1(Γ−1) before the inter-PI line Ŝ3eŜ2s while L(θ) is rotated clockwise. Hence, the Radon planes between the tangent of Γ+1 (Γ−1) and Ŝ3eŜ2s are not exactly weighted. Because the slope of ŝ3eŝ2s is negative and the slope of the tangent of Γ+1 (Γ−1) is more than −0.35 h/R, the percentage of the incorrectly weighted Radon planes is less than ρ=2.57% for h/R=0.2.

Then, one must further consider the incorrectly weighting planes caused by the critical event in Case 5. Recall that the L-curves are used to split the domain D4 into sub-domains, making the weight function n continuous across all the A-curves, and the cusps are the starting and ending points of the L-curves. If the endpoints of the L-curves are not the cusps, there will be small fractions (or “line segments”) on the A-curves, making the weight function n discontinuous across them and ruining the FBP structure of our algorithms. It possibly occurs for r≧0.265 R. As discussed above, the Bs-curve will first enter a 1-IP domain from a 3-IP domain across the “line segment”, and then pass through an L-curve. On the detector plane, this corresponds to that L(θ) intersects the inter-PI line ŝ3eŝ2s before the L-line while L(θ) is rotated clockwise. Recall that if the cusp is not in D4, Ŝ2s is possibly to the left of Ŝu or Ŝ3e is to the right of Ŝd. Thus, the slope of Ŝ3eŜ2s is always more than −0.35 h/R.
Because the slope of the L-curve is never positive, the percentage of the incorrectly weighted Radon planes is less than ρ=2.57% for h/R=0.2. On the other hand, based on the discussion on Lemma 3, one or more cusps may possibly remain even when r≧0.265 R, which means that less “line segments” related to critical events will appear in Case 5, and in fact more Radon planes may be correctly weighted.
The implementation of these algorithms consists of one or more, and preferably all, of the following steps: Step 1) Differentiate each projection with respect to variable s; Step 2) For each yi(s), i={1, 2, 3}, perform the Hilbert transform of derivative data along the given filtering directions on the corresponding detector plane; Step 3) Backproject the filtered data on the inter-PI segments to reconstruct the object point. Differences between the algorithms described herein and some of the ones previously described include, but are not limited to, differences in triple-helices geometry the filtered data are backprojected on inter-PI segments and that there are two families of filtering lines for each algorithm so that each point on the detector plane will be filtered twice. Also, since the algorithms described herein allow shift-invariant filtration, all results are in Cartesian coordinates directly, and there is no coordinate transform necessary similar to what was used in the slow-FBP algorithm or BPF algorithm.
Previously published BPF algorithms for triple-source helical CBCT can indeed produce excellent image quality, FBP algorithms (either “slow” or “fast”) are computationally desirable for several reasons, such as being amendable for parallel processing. In particular, while the computational structures of our BPF algorithm and FBP algorithms are quite similar, the FBP algorithms avoid densely sampled intermediate reconstruction in the PI-line-based coordinate system, and more importantly they can reconstruct a region of interest (ROI) or volume of interest (VOI) much more efficiently than the BPF counterpart. Note that ROI/VOI reconstruction is very common in medical imaging. A related technology called “interior tomography” is being actively developed to target this type of problems. Then, an interesting possibility would be to develop tripe-source interior CBCT.
The inventive two fast exact/quasi-exact FBP algorithms for triple-source helical CBCT have their advantages and disadvantages. From the perspective of exact reconstruction, the first algorithm is more desirable than the second algorithm because it is not affected by critical events in Case 3. However, in terms of efficient data acquisition, it may require a larger detector area than the second algorithm. In the medical CT field, the rectangular detector shape is most popular, and the helical pitch may be varied case by case. Therefore, it is practically possible to have projection data for reconstruction using either or both of the two fast FBP algorithms.
The methods disclosed herein can be practiced on any CT system. An example of a CT system and apparatus capable of implementing the methods is provided is an electron beam CT. In that framework, a curvilinear tungsten material or target can be arranged along a non-standard curve to be traced by an electro-magnetically driven electron-beam for formation of an X-ray source and collection of cone-beam data.
An exemplary electron beam CT comprises a vacuum chamber having an exterior surface, an underlying interior surface, and defines an enclosed space. At least a portion of the exterior surface can define or surround a subject cavity. The subject cavity is adapted to receive a subject. The subject cavity can be adapted to receive a human, a mouse or a rat, e.g.
The apparatus can further comprise a charged particle beam generator having a proximal and a spaced distal end. The electron beam generator can generate a flat or curved electron sheet. The electron beam generator can have a scanning speed from about 25 Hz to about 50 Hz. The apparatus can comprise a single electron beam generator or a plurality of electron beam generators.
The apparatus can include a focusing mechanism adapted to selectively focus charged particles generated by the charged particle beam generator and a target adapted to generate X-rays upon receipt of charged particles from the charged particle beam generator. If a plurality of electron beam generators are used, the apparatus may comprise a plurality of focusing mechanisms.
The apparatus further comprises a detector or a series of detectors surrounding the target. There are a large variety of detectors that can be used in the disclosed apparatuses, systems and methods. Two representative types are a) thin-film transistors (TFT, .alpha.-Si:H) and b) mono-crystalline silicone CCD/CMOS detectors. Although their quantum efficiency is high, the readout speed of TFT detectors is generally less than 30 frames per second, rarely reaching 100 frames per second. On the other hand, the readout speed of CCD/CMOS detectors can be extremely high, such as 10,000-30,000 frames per second, and are coupled with fiber-optical tapers, resulting in low quantum efficiency. For example, the 1000 Series camera from Spectral Instruments (Tucson, Ariz.), can be used. This camera is compact, measuring 92 by 92 by 168 mm. Two, three, and four-phase architecture CCDs from Fairchild Imaging (Milpitas, Calif.), E2V (Elmsford, N.Y.), Kodak (Rochester, N.Y.), and Atmel (San Jose, Calif.) can be placed in the selected camera. The readout and digitization can use 16-bit digitizer. The pixel readout rate can be varied from 50 kHz to 1 MHz. The gain of the analog processor can be modified under computer control to compensate for the gain change of the dual slope integrator at different readout speeds. The 1000 Series system offers fully programmable readout of sub arrays and independent serial and parallel register binning. In addition, specialized readout modes, such as time delay and integration (TDI) using an internal or external time base can be used. These capabilities allow the readout of only the area of the CCD of interest at variable resolution in order to optimize image signal to noise ratio.
To facilitate a better understanding of the present invention, the following examples of certain aspects of some embodiments are given. In no way should the following examples be read to limit, or define, the scope of the invention.
EXAMPLES
Example 1
To verify and showcase the fast FBP algorithms of the present invention, numerical tests were performed using the Clock phantom. This phantom consists of ellipses, as parameterized in Table IV. In the simulations, the origin of the reconstruction coordinate system was set to the center of each phantom. The spherical phantom support was of 375 mm for the experiment. Three sources were arranged uniformly along a circle with their corresponding detectors on the opposite side. The source-detector distance was 1000 mm. Projections were generated from 1000 view angles while the sources and the detectors were constantly moved along three helixes in one turn. The helix was of 750 mm in radius and 100 mm in pitch. The detector plane consisted of 1300×200 detection elements of 1.0×1.0 mm2.







TABLE IV







Parameters of the Clock Phantom















No.
xc
yc
zc
a
b
c
θ
ρ


















1
0
0
0
1
1
1
0
1


2
0
0.8
0
0.1
0.1
0.1
0
1


3
0.4
0.69
0.01
0.1
0.1
0.1
0
1


4
0.69
0.4
0.02
0.1
0.1
0.1
0
1


5
0.8
0
0.03
0.1
0.1
0.1
0
1


6
0.69
−0.4
0.04
0.1
0.1
0.1
0
1


7
0.4
−0.69
0.05
0.1
0.1
0.1
0
1


8
0
−0.8
0.06
0.1
0.1
0.1
0
1


9
−0.4
−0.69
0.07
0.1
0.1
0.1
0
1


10
−0.69
−0.4
0.08
0.1
0.1
0.1
0
1


11
−0.8
0
0.09
0.1
0.1
0.1
0
1


12
−0.69
0.4
0.1
0.1
0.1
0.1
0
1


13
−0.4
0.69
0.11
0.1
0.1
0.1
0
1


14
0
0.5
0
0.05
0.05
0.05
0
1


15
0.25
0.43
−0.01
0.05
0.05
0.05
0
1


16
0.43
0.25
−0.02
0.05
0.05
0.05
0
1


17
0.5
0
−0.03
0.05
0.05
0.05
0
1


18
0.43
−0.25
−0.04
0.05
0.05
0.05
0
1


19
0.25
−0.43
−0.05
0.05
0.05
0.05
0
1


20
0
−0.5
−0.06
0.05
0.05
0.05
0
1


21
−0.25
−0.43
−0.07
0.05
0.05
0.05
0
1


22
−0.43
−0.25
−0.08
0.05
0.05
0.05
0
1


23
−0.5
0
−0.09
0.05
0.05
0.05
0
1


24
−0.43
0.25
−0.1
0.05
0.05
0.05
0
1


25
−0.25
0.43
−0.11
0.05
0.05
0.05
0
1









The algorithms were coded in MATLAB and executed on a regular PC (Intel Core2 Duo CPU 3.06 GHz, 4 GB RAM). Reconstructed images are shown in FIG. 22. Our numerical results show that in the case of r=0.495 R both two algorithms produced high quality images.
Example 2
Auxiliary Lemmas were used as described below. A point was fixed at xεΩ and its three associated inter-PI lines were found as shown in FIG. 3. Then, a source position was selected as sε(sj3, sje), jε{1, 2, 3} and how the inter-PI lines project onto the corresponding detector plane was determined. For simplicity, in this disclosure the projection of yi(s), jε{1, 2, 3} on a detector plane is denoted by Ŝj.
Lemma 1.
On a detector plane, the slopes of the projected inter-PI lines ŜjsŜjmod3+1e and Ŝ(j+1)mod3+1s Ŝje are always positive, and that of the inter-PI line Ŝjmod3+1sŜ(j+1)mod3+1e, jε{1, 2, 3} is always negative.
Proof of Lemma 1.
Without loss of generality or wishing to be limited by theory, the source position was selected to be y1(s0), s0ε(s1s, s1e). By construction, s2e−s1s<2π, s1e−s3s<2π and s3e−s2s>0. Hence, the projections of s1s, s3e, and s3s were always to the left of those of s2e, s2e, and s1e respectively (FIG. 18). When sε(s1s, s1e) changed, the point {circumflex over (x)}, i.e., the projection of x onto the detector plane, was moved inside the region G:=G1∪G2∪G3. Clearly, {circumflex over (x)} could reach its highest (respectively, lowest) position in the vertical direction when {circumflex over (x)} was at the intersection of Γ1 and Γ+1 (respectively, of Γ, and Γ−1). Also, the vertical coordinates at these points are




  
    
      v
      max
    
    =
    
      
        1.3794
        ⁢
        
          Dh
          
            2
            ⁢
            π
            ⁢
            
              
            
            ⁢
            R
          
        
        ⁢
        
          
        
        ⁢
        and
        ⁢
        
          
        
        ⁢
        
          v
          min
        
      
      =
      
        
          -
          1.3794
        
        ⁢
        
          Dh
          
            2
            ⁢
            π
            ⁢
            
              
            
            ⁢
            R
          
        
      
    
  
  ,




respectively. Moreover, the lowest point on Γ+2 and the highest point on Γ−2 are





  
    v
    min
    ′
  
  =
  
    
      1.3801
      ⁢
      
        Dh
        
          2
          ⁢
          π
          ⁢
          
            
          
          ⁢
          R
        
      
      ⁢
      
        
      
      ⁢
      and
      ⁢
      
        
      
      ⁢
      
        v
        max
        ′
      
    
    =
    
      
        -
        1.3801
      
      ⁢
      
        Dh
        
          2
          ⁢
          π
          ⁢
          
            
          
          ⁢
          R
        
      
    
  




respectively. Evidently, vmax≦vmin′ and vmax′≦vmin. Since ŝ1s and ŝ3s were to the left of ŝ2e and ŝ1e, the slopes of the inter-PI lines ŝ1sŝ2e and ŝ3sŝ1e were positive for all {circumflex over (x)} in G.

The inter-PI line ŝ2sŝ3e was satisfied by EQUATION 18 below:




  
    
      
        
          
            x
            1
          
          =
          
            
              Rt
              ⁢
              
                
              
              ⁢
              cos
              ⁢
              
                
              
              ⁢
              
                s
                2
                s
              
            
            +
            
              
                R
                ⁡
                
                  (
                  
                    1
                    -
                    t
                  
                  )
                
              
              ⁢
              cos
              ⁢
              
                
              
              ⁢
              
                s
                3
                e
              
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            x
            2
          
          =
          
            
              Rt
              ⁢
              
                
              
              ⁢
              sin
              ⁢
              
                
              
              ⁢
              
                s
                2
                s
              
            
            +
            
              
                R
                ⁡
                
                  (
                  
                    1
                    -
                    t
                  
                  )
                
              
              ⁢
              sin
              ⁢
              
                
              
              ⁢
              
                s
                3
                e
              
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            
              x
              3
            
            =
            
              
                
                  h
                  
                    2
                    ⁢
                    π
                  
                
                ⁢
                
                  t
                  ⁡
                  
                    (
                    
                      
                        s
                        2
                        s
                      
                      -
                      
                        
                          2
                          ⁢
                          π
                        
                        3
                      
                    
                    )
                  
                
              
              +
              
                
                  h
                  
                    2
                    ⁢
                    π
                  
                
                ⁢
                
                  (
                  
                    1
                    -
                    t
                  
                  )
                
                ⁢
                
                  (
                  
                    
                      s
                      3
                      e
                    
                    -
                    
                      
                        4
                        ⁢
                        π
                      
                      3
                    
                  
                  )
                
              
            
          
          ,
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          18
        
        )
      
    
  




where tε[0,1] and s2s, s3eε(s0, s0+2π).
By allowing

x1=r0 cos μ0 

x2=r0 sin μ0  (EQUATION 19)

where r0ε[0,0.495 R] and μ0ε[0,2π].
The following was left:




  
    
      
        
          
            cos
            ⁢
            
              
                
                  s
                  3
                  e
                
                -
                
                  s
                  2
                  s
                
              
              2
            
          
          =
          
            
              
                r
                0
              
              ⁢
              
                sin
                ⁡
                
                  (
                  
                    
                      μ
                      0
                    
                    -
                    
                      s
                      2
                      s
                    
                  
                  )
                
              
            
            
              
                
                  R
                  2
                
                +
                
                  r
                  0
                  2
                
                -
                
                  2
                  ⁢
                  
                    Rr
                    0
                  
                  ⁢
                  
                    cos
                    ⁡
                    
                      (
                      
                        
                          μ
                          0
                        
                        -
                        
                          s
                          2
                          s
                        
                      
                      )
                    
                  
                
              
            
          
        
        ,
        and
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          20
        
        )
      
    
  




then, EQUATION 20 was rewritten as EQUATION 21 below:




  
    
      
        
          cos
          ⁢
          
            
              
                s
                3
                e
              
              -
              
                s
                2
                s
              
            
            2
          
        
        =
        
          
            
              sin
              ⁡
              
                (
                
                  
                    μ
                    0
                  
                  -
                  
                    s
                    2
                    s
                  
                
                )
              
            
            
              
                
                  
                    (
                    
                      
                        R
                        
                          r
                          0
                        
                      
                      -
                      
                        cos
                        ⁡
                        
                          (
                          
                            
                              μ
                              0
                            
                            -
                            
                              s
                              2
                              s
                            
                          
                          )
                        
                      
                    
                    )
                  
                  2
                
                +
                
                  
                    sin
                    2
                  
                  ⁡
                  
                    (
                    
                      
                        μ
                        0
                      
                      -
                      
                        s
                        2
                        s
                      
                    
                    )
                  
                
              
            
          
          .
        
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          21
        
        )
      
    
  




When μ0−s2e was fixed and r0 was reduced,




  
  
    cos
    ⁢
    
      
        
          s
          3
          e
        
        -
        
          s
          2
          s
        
      
      2
    
  
  




decreased. Therefore, the right side of EQUATION 21 reached its maximum or minimum when r0 is maximized, i.e. at r0=0.495 R. Those maximum and minimum values were numerically calculated (FIG. 19A), and shown below in EQUATIONS 22 and 23:





  
    
      
        
          -
          0.4949
        
        ≤
        
          cos
          ⁢
          
            
              
                s
                3
                e
              
              -
              
                s
                2
                s
              
            
            2
          
        
        ≤
        0.4949
      
    
    
      
        (
        
          EQUATION
          ⁢
          
            
          
          ⁢
          22
        
        )
      
    
  




or 2.1062<s3e−s2s<4.1770  (EQUATION 23)

Next, it was shown that 0≦s2ss0≦4.1773 implied

v1(s3e−s0)−v2(s2s−s0)>0,  (EQUATION 24)

where




  
    
      v
      1
    
    ⁡
    
      (
      s
      )
    
  
  =
  
    
      
        Dh
        
          2
          ⁢
          π
          ⁢
          
            
          
          ⁢
          R
        
      
      ⁢
      
        
          s
          -
          
            4
            ⁢
            
              π
              /
              3
            
          
        
        
          1
          -
          
            cos
            ⁢
            
              
            
            ⁢
            s
          
        
      
      ⁢
      
        
      
      ⁢
      and
      ⁢
      
        
      
      ⁢
      
        
          v
          2
        
        ⁡
        
          (
          s
          )
        
      
    
    =
    
      
        Dh
        
          2
          ⁢
          π
          ⁢
          
            
          
          ⁢
          R
        
      
      ⁢
      
        
          
            s
            -
            
              2
              ⁢
              
                π
                /
                3
              
            
          
          
            1
            -
            
              cos
              ⁢
              
                
              
              ⁢
              s
            
          
        
        .
      
    
  




FIG. 19B shows the function Φ(s=1−cos s−(s−4/3π) sin s in the range sε(0, 2π), demonstrating that {dot over (v)}1(s) was always positive. Hence v1(s3e−s0) was monotonically increasing. s3e=s2s+2.1062 was fixed and the function) Ψ(s2e)=v1(s3e−s0)−v2(s2s−s0) was plotted in FIG. 19C. Clearly, this function was always positive in the range 0≦s2s−s0≦4.1773. Note that s=4.1773+s0 was the intersection of Γ1 and Γ+1 and ŝ2s could not be to the left of this point (otherwise, {circumflex over (x)} is outside G). EQUATION 24 indicated that ŝ2s was always lower than ŝ3e in the vertical direction. Since ŝ2s was to the right of ŝ3e, the slope of the inter-PI line ŝ2sŝ3e was always negative. Due to symmetry, the other two cases sε(s2s, s2e) and sε(s3s, s3e) was handled similarly. This finishes the proof.
Lemma 2.
A T-curve cannot be tangent to an A-curve at an interior point of T.
Proof of Lemma 2.
The interior point of T-curve can be any point of the T-curve except an endpoint. It has been proved previously that a T-curve is smooth everywhere, except possibly at a cusp. When sie, iε{1, 2, 3} is chosen to be the point where the cusp occurs. It was assumed that a T-curve was tangent to an A-curve at α(s). If s=sicε(sis, sie), where sis, sie were the endpoints of the T-curve, then the osculating plane Πe(x) intersected the helix yi(s) at only one point and it contained one inter-PI line. By construction, Πc(x) intersected the detector plane at the asymptote of the Tam-Danielson window boundary and {circumflex over (x)} belonged to the asymptote. Connecting {circumflex over (x)} and ŝis, {circumflex over (x)} and ŝie we detected two inter-PI lines. Clearly, Πe(x) could not contain any of them. By Lemma 1, the third inter-PI line had a negative slope, thus it would not overlap the asymptote. Hence, Πe(x) could not contain it. Consequently, s≠sie and T-curve was smooth in a neighborhood of α(s). If θ was chosen to be the polar angle for the great circle (x−yi(se))1, saε{sis, sie}, iε{1, 2, 3}, then the A-curve could consist of all the unit vectors αi(θ)ε(x−yi(se))⊥. Clearly, {dot over (α)}1(θ) could be perpendicular to α1(θ) and (x−y1(sa)). By construction, the T-curve was tangent to the A-curve at α(s). Hence, {dot over (α)}(s) was be parallel to {dot over (α)}1(θ). That is, {dot over (α)}(s) was perpendicular to α1(θ) and (x−yi(sa)). Because {dot over (α)}(s) was also perpendicular to (x−yi(s1)), s1ε(sis, sie), (x−yi(s1)) was parallel to (x−yi(sa)) and s1=sa, which contradicted the assumption that si is an inner point of the T-curve. This finishes the proof.
Lemma 3.
Case 1, 2, 4, 6, 7 do not occur for r<0.495 R and Case 5 does not occur for r<0.265 R.
Proof of Lemma 3.
Cases 1 and 2 were impossible because they mean that there can be a plane containing three inter-PI lines or tangent to three inter-PI arcs. In Case 4, there can be one plane Π containing one inter-PI line and tangent to two inter-PI arcs. If one assumes this inter-PI line is s1ss2e, chooses a point s=s1s on yi(s) and denotes L=Π∩DP (s0), then by construction, on the detector plane {circumflex over (x)} is on Γ+1 and it overlaps ŝ2e. Then, L may tangent to Γ+1 and Γ+4 or Γ−1 and Γ−1, see FIG. 20. Because the points of tangency are on the inter-PI arcs, the endpoint ŝ1e is to the left of the tangency for case A and ŝ3e was to the right of the tangency for case B. Connecting ŝ1e and {circumflex over (x)}(or ŝ3s and {circumflex over (x)}) we find that the slope of the inter-PI line ŝ3sŝ1e could be negative. By Lemma 1, these two cases were impossible.
In Case 5, there was one plane containing one inter-PI line and tangent to one inter-PI arc at the inflection point. Thus, the inter-PI line on the detector overlapped the inflection line, i.e., ŝ2sŝ3e overlapped Lit when s2s=su and s3e=2su, where su is difference between s21 and s3e. When looking at EQUATION 21 and μ0−s2e is fixed, the absolute value of




  cos
  ⁢
  
    
      
        s
        3
        e
      
      -
      
        s
        2
        s
      
    
    2
  




was monotonically decreasing when r0 was reduced. Then, the range of s3e−s2e was narrowed. In other words, the difference between see and s2e became closer to π. Case 5 occurred when s3e−s2e=su. If order to exclude Case 5, the range of s3e−s2e could not cover the value su=2.6053. Hence, the minimum range of s3e−s2e is 2.6053<s3e−s2s2π−2.6053. That is,





  
    -
    0.2649
  
  <
  
    cos
    ⁢
    
      
        
          s
          3
          e
        
        -
        
          s
          2
          s
        
      
      2
    
  
  <
  
    0.2649
    ⁢
    
      
    
    ⁢
    and
    ⁢
    
      
    
    ⁢
    cos
    ⁢
    
      
        
          s
          3
          e
        
        -
        
          s
          2
          s
        
      
      2
    
  




reached its extreme when r0=0.265 R. From EQUATIONS 21 and 22 we have s3e−s2s=su only for r≧0.265 R thereby contradicting our condition. Hence, Case 5 was impossible for r<0.265 R.

In Case 6, a T-curve will intersect one A-curve twice before meeting a cusp. Suppose this took place at inter-PI arc . A point s=s1s on y1(s) was chosen and observations of what happens on the detector plane when s moves were taken. By construction, the plane Π containing {dot over (y)}1(s) and x intersected the detector plane at the line L which was parallel to the helix tangent across {circumflex over (x)}. At s=s1s, {circumflex over (x)} was on Γ+1 and Π contained inter-PI line ŝ1sŝ2e. As s moved along y1(s), {circumflex over (x)} moved downwards. Notice that L was parallel to the asymptote of the Tam-Dannielson window, so it would not intersect Γ−2 provided that {circumflex over (x)} moved across the asymptote, at where the cusp occurred. Hence, Π would not contain the inter-PI line ŝ1sŝ2e and Case 6 was impossible. By Lemma 2, Case 7 was impossible. This finishes the proof.
Lemma 4.
The inflection point ŝu(ŝd) is inside the inter-PI arc when {circumflex over (x)} is in G21 (G22).
Proof of Lemma 4.
By Lemma 3, any point in the area r<0.265 R had three cusps in the diagram. Note that there was one IP in each inter-PI arc within D4. Since all three cusps were in D4, an osculating plane of one inter-PI arc intersected two other inter-PI arcs exactly once at one point. Assuming that this osculating plane Πe contained x and considering Πe of the second inter-PI arc (i.e., of y2(s)). s10 was set to be the point where it intersected the first inter-PI arc (i.e., on y1(s)). s was moved along the first inter-PI arc and the results were observed with {circumflex over (x)} on the detector when s=s1s, {circumflex over (x)} entered the Zhao window through Γ+1, and when s=s120, {circumflex over (x)} belongs to Lu. As follows from the diagram, the point su must be inside the second inter-PI arc, i.e., between ŝ2s and s2e. As the point s moved further, the difference s2s−s became smaller, and the point ŝ2e moved to the right of ŝu along Γ+. The inter-PI line ŝ1sŝ2e had a positive slope. Thus, as long as {circumflex over (x)} was inside G21, the point ŝ2e was always to the left of ŝu. The case where {circumflex over (x)} was in G22 can be similarly treated. This proves Lemma 4.
Lemma 5.
An L-curve never intersects an A-curve, for r<0.265 R.
Proof of Lemma 5.
If an L- and A-curve intersect, an L-line through {circumflex over (x)} can overlap the inter-PI line. One point on the inter-PI arc  was chosen and the slope of the inter-PI line was considered on the corresponding detector plane. By construction, an L-curve always started from a cusp of a T-curve and ended on a cusp of another T-curve. For the osculating plane Πe, its intersection with the detector plane was the line tangent to Γ+1 (Γ−1) at ŝu(ŝd). By Lemma 4, the endpoints of the inter-PI arc on Γ+1 (Γ−1) were on different sides of ŝu (ŝd), and one of them on the right (left) side was also an endpoint of the inter-PI line for helices y2(s) and y3(s), and was denoted as ŝ2s(ŝ3e) in FIG. 22.
If {circumflex over (x)} was in G2 and above ŝu, the L-line was formed by connecting {circumflex over (x)}, ŝu, and the inter-PI line was formed connecting {circumflex over (x)}, ŝ2s. If z was in G2 and below ŝd, we formed the L-line by connecting {circumflex over (x)}, ŝd, and the inter-PI line was formed by connecting {circumflex over (x)}, ŝ3e. Clearly, in any case the slope of L-line was between zero and the slope of the inter-PI line. That is, the L-line could not overlap the inter-PI line. For other {circumflex over (x)}, the L-line was parallel to the u axis. By Lemma 1, it was always between two inter-PI lines and could not overlap with any of them. For the point on other inter-PI arcs, the situation was the same. This finishes the proof.
The present invention has been described with reference to particular embodiments having various features. It will be apparent to those skilled in the art that various modifications and variations can be made in the practice of the present invention without departing from the scope or spirit of the invention. One skilled in the art will recognize that these features may be used singularly or in any combination based on the requirements and specifications of a given application or design. Other embodiments of the invention will be apparent to those skilled in the art from consideration of the specification and practice of the invention. It is intended that the specification and examples be considered as exemplary in nature and that variations that do not depart from the essence of the invention are intended to be within the scope of the invention.
Therefore, the present invention is well adapted to attain the ends and advantages mentioned as well as those that are inherent therein. The particular embodiments disclosed above are illustrative only, as the present invention may be modified and practiced in different but equivalent manners apparent to those skilled in the art having the benefit of the teachings herein. Furthermore, no limitations are intended to the details of construction or design herein shown, other than as described in the claims below. It is therefore evident that the particular illustrative embodiments disclosed above may be altered or modified and all such variations are considered within the scope and spirit of the present invention. While compositions and methods are described in terms of “comprising,” “containing,” or “including” various components or steps, the compositions and methods can also “consist essentially of” or “consist of” the various components and steps. All numbers and ranges disclosed above may vary by some amount. Whenever a numerical range with a lower limit and an upper limit is disclosed, any number and any included range falling within the range is specifically disclosed. In particular, every range of values (of the form, “from about a to about b,” or, equivalently, “from approximately a to b,” or, equivalently, “from approximately a-b”) disclosed herein is to be understood to set forth every number and range encompassed within the broader range of values. Also, the terms in the claims have their plain, ordinary meaning unless otherwise explicitly and clearly defined by the patentee. Moreover, the indefinite articles “a” or “an,” as used in the claims, are defined herein to mean one or more than one of the element that it introduces. If there is any conflict in the usages of a word or term in this specification and one or more patent or other documents that may be incorporated herein by reference, the definitions that are consistent with this specification should be adopted.
Throughout this application, various publications are referenced. The disclosures of these publications in their entireties are hereby incorporated by reference into this application in order to more fully describe the features of the invention and/or the state of the art to which this pertains. The references disclosed are also individually and specifically incorporated by reference herein for the material contained in them that is discussed in the portion of this disclosure in which the reference is relied upon.
G. Wang, C. R. Crawford, and W. A. Kalender, “Guest editorial-Multirow detector and cone-beam, spiral/helical CT,” Medical Imaging, IEEE Transactions on, vol. 19, pp. 817-821, 2000.
T. G. Flohr, C. H. McCollough, H. Bruder, M. Petersilka, K. Gruber, C. Süβ, M. Grasruck, K. Stierstorfer, B. Krauss, and R. Raupach, “First performance evaluation of a dual-source CT (DSCT) system,” European Radiology, vol. 16, pp. 256-268, 2006.
M. Vannier and G. Wang, “Spiral CT refines imaging of temporal bone disorders,” Diagnostic imaging, vol. 15, p. 116-121, 1993.
G. Wang, S. Zhao, H. Yu, C. Miller, P. Abbas, B. Gantz, S. Lee, and J. Rubinstein, “Design, analysis and simulation for development of the first clinical micro-CT scanner1,” Academic Radiology, vol. 12, pp. 511-525, 2005.
G. Wang, T.-H. Lin, P.-C. Cheng, and D. M. Shinozaki, “A general cone-beam reconstruction algorithm,” IEEE Transactions on Medical Imaging, vol. 12, p. 486, 1993.
A. Katsevich, “Theoretically exact filtered backprojection-type inversion algorithm for spiral CT,” SIAM Journal on Applied Mathematics, vol. 62, p. 2012, 2002.
S. S. Orlov, “Theory of three-dimensional reconstruction. 1. Conditions of a complete set of projections,” Sov. Phys. Crystallogr, vol. 20, pp. 312-314, 1975.
I. Gel'fand and M. Graev, “Crofton's function and inversion formulas in real integral geometry,” Functional Analysis and Its Applications, vol. 25, pp. 1-5, 1991.
H. Rullg{dot over (a)}rd, “An explicit inversion formula for the exponential Radon transform using data from 180,” Arkiv för Matematik, vol. 42, pp. 353-362, 2004.
G. Wang, Y. Ye, and H. Yu, “Approximate and exact cone-beam reconstruction with standard and non-standard spiral scanning,” Physics in Medicine and Biology; vol. 52, pp. 1-13, 2007.
J. Zhao, M. Jiang, T. Zhuang, and G. Wang, “An exact reconstruction algorithm for triple-source helical cone-beam CT,” Journal of X-Ray Science and Technology, vol. 14, p. 191, 2006.
J. Zhao, Y. Jin, Y. Lu and G. Wang, “A Filtered Backprojection Algorithm for Triple-Source Helical Cone-Beam CT,” Medical Imaging, IEEE Transactions on, vol. 28, pp. 384-393, 2009.

## Claims

1. A computed tomography (CT) imaging method comprising:
scanning an object using triple-source helical cone-beam computed tomography (CBCT) to acquire projection data relating to the object being imaged, where each x-ray source of the CBCT is disposed opposite a detector and has a scanning radius that is a distance R from a rotation axis, and where each detector covers a field of view less than 0.495 R; and
reconstructing the scanned portion of the object into an image by performing a computationally efficient filtered backprojection (FBP) and theoretically exact/quasi-exact algorithm to generate image data.

2. The method of claim 1, further comprising:
supporting the object in a stationary position; and
moving each source of the triple-source and its associated detector about the object at a constant speed to generate three spiral scans with source trajectories y1(s), y2(s), and y3(s) defined as:





  {
  
    
      
        
          
            
              
                y
                1
              
              ⁡
              
                (
                s
                )
              
            
            =
            
              (
              
                
                  R
                  ⁢
                  
                    
                  
                  ⁢
                  cos
                  ⁢
                  
                    
                  
                  ⁢
                  s
                
                ,
                
                  
                
                ⁢
                
                  R
                  ⁢
                  
                    
                  
                  ⁢
                  sin
                  ⁢
                  
                    
                  
                  ⁢
                  s
                
                ,
                
                  
                    h
                    
                      2
                      ⁢
                      π
                    
                  
                  ⁢
                  s
                
              
              )
            
          
        
      
      
        
          
            
              
                y
                2
              
              ⁡
              
                (
                s
                )
              
            
            =
            
              (
              
                
                  R
                  ⁢
                  
                    
                  
                  ⁢
                  
                    cos
                    ⁡
                    
                      (
                      
                        s
                        +
                        
                          
                            2
                            3
                          
                          ⁢
                          π
                        
                      
                      )
                    
                  
                
                ,
                
                  R
                  ⁢
                  
                    
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      
                        s
                        +
                        
                          
                            2
                            3
                          
                          ⁢
                          π
                        
                      
                      )
                    
                  
                
                ,
                
                  
                    h
                    
                      2
                      ⁢
                      π
                    
                  
                  ⁢
                  s
                
              
              )
            
          
        
      
      
        
          
            
              
                y
                3
              
              ⁡
              
                (
                s
                )
              
            
            =
            
              (
              
                
                  R
                  ⁢
                  
                    
                  
                  ⁢
                  
                    cos
                    ⁡
                    
                      (
                      
                        s
                        +
                        
                          
                            4
                            3
                          
                          ⁢
                          π
                        
                      
                      )
                    
                  
                
                ,
                
                  R
                  ⁢
                  
                    
                  
                  ⁢
                  
                    sin
                    ⁡
                    
                      (
                      
                        s
                        +
                        
                          
                            4
                            3
                          
                          ⁢
                          π
                        
                      
                      )
                    
                  
                
                ,
                
                  
                    h
                    
                      2
                      ⁢
                      π
                    
                  
                  ⁢
                  s
                
              
              )
            
          
        
      
    
     
  



where R is the distance from the x-ray source to the rotation axis, h is helical pitch, s is a scan path corresponding to source position.

3. A method of computing images derived from triple-source spiral computed tomography scan with three detectors, comprising the steps of:
(a) collecting cone beam data from three detectors during a scan of an object;
(b) for each source position yj(s), jε{1, 2, 3}, identifying two families of lines on a detector plane DP(s) corresponding to a source position s and containing the corresponding detector and intersecting the cone beam, and two families of lines include:
i. a first family of lines parallel to L0, where
L0 is the projection of the helical tangent at current source position;

ii. a second family of lines tangent to Γ+1 and Γ−1, or parallel to the horizontal axis of the plane DP(s), where
Γ+1 is the projection of the helical turn yjmod3+1(s) defined by s<q<s+2π onto the plane DP(s);
Γ−1 is the projection of the helical turn y(j+1)mod3+1(s) defined by s<q<s+2π onto the plane DP(s);
q is the parameter along the scan path which describes the point being projected;


(c) computing a derivative of the cone beam data with respect to the source position;
(d) performing Hilbert transform of the derivative of the cone beam data along the two families of lines, where the Hilbert transform is a convolution between the derivative of the cone beam data and a kernel function h(t)=1/(πt);
(e) back projecting said filtered data to form a precursor of said image; and
(f) repeating steps a, b, c, d and e to obtain an image.

4. The method of claim 3, wherein identifying the second family of lines includes:
the lines tangent to Γ+1, when the projection of point x onto DP(s) is located in the area bounded by Γl, Lit, and Γ+1;
the lines tangent to Γ−1, when the projection of point x onto DP(s) is located in the area bounded by Γr, Lib and Γ−1;
the lines parallel to the horizontal axis of the plane DP(s), when the projection of x onto DP(s) is located in the area bounded by Γl, Γr, Lit, Lib, Γ+1 and Γ−1, where
Γl and Γr are the projections of the object support limitation r=0.495 R onto DP(s);
Lit is the inflection line of Γ+1;
Lib is the inflection line of Γ−1;
r is the radius of the object support, and
R is the radius of the scanning trajectory.

5. The method of claim 3, wherein the back projection step(e) includes:
(ei) fixing a reconstruction point x, which represents a point inside the object being scanned where it is required to reconstruct the image;
(eii) determining the three inter-PI arcs for x;
(eiii) finding the projection {circumflex over (x)} of x onto a detector plane DP(s);
(eiv) identifying lines from the two families of lines and points on the said lines that are passing through the said projection {circumflex over (x)};
(ev) computing contribution from filtered cone beam data to the image being reconstructed at the point x by multiplying





  
    -
    
      1
      
        4
        ⁢
        
          π
          2
        
        ⁢
        
          
          
            x
            -
            
              y
              ⁡
              
                (
                s
                )
              
            
          
          
        
      
    
  
  ;



(evi) adding the contribution from filtered cone beam data to the image being reconstructed at the point x according to the three inter-PI arcs;
(evii) going to step (ei) and choose a different reconstruction point x.

6. The method of claim 5, wherein the three inter-PI arcs for x are determined according to the following rules:
the endpoints of the inter-PI arc on a first helical turn y1(s) are s=s1s and s=s1e, s1e>s1s;
the endpoints of the inter-PI arc on a second helical turn y2(s) are s=s2s and s=s2e, s2e>s2s;
the endpoints of the inter-PI arc on a third helical turn y3(s) are s=s3s and s=s3e, s3e>s3s;


|s2e−s1s|<2π;



|s3e−s2s|<2π;



|s1e−s3s|<2π;


the line connecting y1(s1e) and y3(s3s) passes through x;
the line connecting y2(s2e) and y1(s1s) passes through x;
the line connecting y3(s3e) and y2(s2e) passes through x.

7. A method of computing images derived from triple-source spiral computed tomography scan with three detectors, comprising the steps of:
(a) collecting cone beam data from three detectors during a scan of an object;
(b) for each source position yj(s), jε{1, 2, 3}, identifying two families of lines on a detector plane DP(s) corresponding to a source position s and containing the corresponding detector and intersecting the cone beam, and two families of lines include:
i. a first family of lines tangent to Γ+2 and Γ−2, where
Γ+2 is the projection of the current helical turn defined by s<q<s+2π onto the plane DP(s);
Γ−2 is the projection of the current helical turn defined by s−2π<q<s onto the plane DP(s);

ii. a second family of lines tangent to Γ+1 and Γ−1, or parallel to the horizontal axis of the plane DP(s), where
Γ+1 is the projection of the helical turn yjmod3+1(s) defined by s<q<s+2π onto the plane DP(s);
Γ−1 is the projection of the helical turn y(j+1)mod3+1(s) defined by s<q<s+2π onto the plane DP(s);


(c) computing the derivative of the cone beam data with respect to the source position;
(d) performing the Hilbert transform of the derivative of the cone beam data along the two families of lines, where the Hilbert transform is a convolution between the derivative of the cone beam data and a kernel function h(t)=1/(πt);
(e) back projecting said filtered data to form a precursor of said image; and
(f) repeating steps a, b, c, d and e until an image of the object is completed.

8. The method of claim 7, wherein identifying the first family of lines includes:
the lines tangent to Γ+2, when the projection of x onto DP(s) is located above L0;
the lines tangent to Γ−2, when the projection of x onto DP(s) is located below L0;
where L0 is the projection of the helical tangent at current source position.

9. The method of claim 7, wherein identifying the second family of lines includes:
the lines tangent to Γ+1, when the projection of x onto DP(s) is located in the area bounded by Γl, Lit and Γ+1;
the lines tangent to Γ−1, when the projection of x onto DP(s) is located in the area bounded by Γr, Lib and Γ−1;
the lines parallel to the horizontal axis of the plane DP(s), when the projection of x onto DP(s) is located in the area bounded by Γl, Γr, Lit, Lib, Γ+1 and Γ−1, where
Γl and Γr are the projections of the object support limitation r=0.495 R onto DP(s);
Lit is the inflection line of Γ+1;
Lib is the inflection line of Γ−1;
r is the radius of the object support, and
R is the radius of the scanning trajectory.

10. The method of claim 7, wherein the back projection step(e) includes:
(ei) fixing a reconstruction point x, which represents a point inside the object being scanned where it is required to reconstruct the image;
(eii) determining the three inter-PI arcs for x;
(eiii) finding the projection {circumflex over (x)} of x onto a detector plane DP(s);
(eiv) identifying lines from the two families of lines and points on the said lines that are passing through the said projection {circumflex over (x)};
(ev) computing contribution from filtered cone beam data to the image being reconstructed at the point x by multiplying





  
    -
    
      1
      
        4
        ⁢
        
          π
          2
        
        ⁢
        
          
          
            x
            -
            
              y
              ⁡
              
                (
                s
                )
              
            
          
          
        
      
    
  
  ;



(evi) adding the contribution from filtered cone beam data to the image being reconstructed at the point v according to the three inter-PI arcs;
(evii) going to step (ei) and choose a different reconstruction point x.

11. The method of claim 7, wherein the three inter-PI arcs for x are determined according to the following rules:
the endpoints of the inter-PI arc on a first helical turn y1(s) are s=s1s and s=s1e, s1e>s1s;
the endpoints of the inter-PI arc on a second helical turn y2(s) are s=s2s and s=s2e, s2e>s2s;
the endpoints of the inter-PI arc on a third helical turn y3(s) are s=s3s and s=s3e, s3e>s3s;


|s2e−s1s|<2π;



|s3e−s2s|<2π;



|s1e−s3s|<2π;


the line connecting y1(s1e) and y3(s3s) passes through x;
the line connecting y2(s2e) and y1(s1s) passes through x;
the line connecting y3(s3e) and y2(s2e) passes through x.