Abstract:
The present invention provides methods and apparatuses for constructing observed data from a Compton camera in order to provide a three-dimensional image of a radiopharmaceutical source distribution within a patient. No intermediate two-dimensional images are formed. Observed data may be analyzed by a processor in order to construct a three-dimensional image representing the source distribution. An idealized mathematical model may express the observed Compton camera data in terms of an integral over the source distribution. An exact analytic inversion is then found for this idealized model. The new analytic solutions arise from a generalization of the integral equation used to model the Compton camera. The kernel of the integral equation is modified by the introduction of an index p that describes the effect of source distance from the camera.

Description:
This application claims priority to provisional U.S. Application. No. 60/458,131, filed Mar. 27, 2003. 

   FIELD OF THE INVENTION 
   The present invention relates to nuclear medicine for reconstructing a source distribution from Compton camera data in the application of nuclear medicine. In particular, the source distribution corresponds to a radiopharmaceutical distribution within a patient. 
   BACKGROUND OF THE INVENTION 
   In nuclear medicine, a radiopharmaceutical is injected into a patient and accumulates in tissues of physiological interest (e.g., cancer tumors or heart muscle). The gamma-ray emissions from this radiopharmaceutical are then observed and the local density of radiopharmaceutical within the patient reconstructed. Two clinical imaging techniques, single photon emission computed tomography (SPECT) and positron emission tomography (PET), are based on this simple idea and are currently available commercially. Both these imaging techniques use filtered backprojection (FBP) algorithms for the reconstruction of the three-dimensional (3D) distribution of radiopharmaceuticals within the patient. The advantages and limitations of SPECT and PET have been reviewed extensively in the literature and will not be further discussed here. It is generally acknowledged, however, that one of the major limitations of SPECT is the loss of counts due to the collimator in the gamma-ray camera. Typically only one gamma ray is observed for every ten thousand emissions, the vast majority of these gamma rays are absorbed in the collimator. In 1974 Todd, Nightingale, and Everett proposed a new type of gamma camera that would replace the collimator with a method of “electronic collimation” based on Compton scattering. Compton kinematics implies a relationship between the angle of scatter and the energy deposited in a detector. Knowledge of this energy and, therefore, the angle of scatter can be used to estimate the direction of an incident gamma-ray. The success of this strategy is notably demonstrated by the gamma-ray telescopes currently surveying the sky for sources of gamma emissions. In medical applications, however, the patient is located at a finite distance, so that the interpretation of the angular information provided by Compton events is complicated by huge parallax effects. In landmark research during the 1980s and 1990s, Manbir Singh and his group constructed and tested a prototype Compton camera. That research revealed both the promise and difficulties associated with Compton cameras. Many of the technical difficulties associated with Compton cameras (e.g., limited energy resolution, restrictive detector geometries, and events involving multiple scatters) are currently being remedied by advances in detector technology. Among the most significant of these problems is the 3D reconstruction of the radiopharmaceutical distribution from the data provided by a Compton camera. 
   The reconstruction of 3D images from Compton camera data poses a formidable mathematical problem. Unlike collimated cameras that reveal the direction of each gamma-ray, the incident direction of each event in a Compton camera is restricted to a cone. This geometric difference completely changes the nature of the reconstruction problem. In his original studies, Singh in collaboration with Hebert and Leahy implemented a maximum likelihood (ML) reconstruction (and tested other iterative methods). In recent years numerous alternative algorithms have been proposed and studied. Roughly speaking, these algorithms can be divided into iterative and analytic inversion methods. Both types of algorithm require the construction of a mathematical model that describes the Compton camera data in terms of the source distribution of radiopharmaceutical substance. Generally, the iterative methods can accommodate more intricate and accurate models, whereas, the analytic inversions require simpler models that are amenable to direct mathematical analysis. Most experimental groups working on Compton cameras have used iterative algorithms similar to the ML techniques of Singh et al. However, the reason for preferring iterative algorithms was not the accuracy of the mathematical model. After all, FBP algorithms (that are prototypes of analytic inversion) are generally used in clinical SPECT and PET despite the existence of ML algorithms for more accurate models. Basically, the iterative algorithms dominate Compton reconstruction because, until the late 1990s, no analytic algorithm existed for Compton cameras. 
   An algorithm for Compton cameras developed by Cree and Bones restricts the acceptable data as much as would a standard collimator—and is, therefore, unsuitable for practical applications. Shortly thereafter, Basko, Gullberg and Zeng produced and patented an algorithm (for the same mathematical model) that removed all restrictions on the Compton data. This algorithm explicitly rebins the Compton data into a form that leads to a 2D Radon transform (and, therefore, FBP). This rebinning procedure, which is the major drawback of the technique, overcomes the basic problem of Compton camera reconstruction; namely, the cone of incident gamma-ray directions. The rebinning implements a complex summation over spherical harmonics that arises from the cone of acceptance. Other researchers developed analytic-inversion algorithms that explicitly replace integrals over the acceptance cone with a summation over spherical harmonics. For example, Parra and Tomitani and Hirasawa retain the spherical harmonic expansions—truncating the series to a finite number of terms for numerical implementation. Although Parra and Tomitani and Hirasawa use a slightly different mathematical model than Basko et al, they face the same problem; namely, the additional degree of freedom introduced by the acceptance cone of the Compton camera. In reading the papers of Parra and Tomitani and Hirasawa, one is reminded of the initial papers of Cormack in which the Radon transform was rediscovered in terms of circular harmonics. One recalls that, when the problem was properly reformulated, the circular harmonics disappeared and the corresponding FBP algorithms emerged. 
   While these analytic algorithms were being developed, other groups examined alternative strategies. In the process of designing and building the C-SPRINT system, the group at the University of Michigan engaged in a broad-based campaign to establish Compton cameras as a viable imaging technology. Not only were simulations and sophisticated statistical analysis used for camera design, but the reconstruction problem was attacked using a variety of techniques. As a leading research institution in ML and statistical reconstruction, the University of Michigan implemented the most sophisticated iterative algorithms currently available. Furthermore, non-iterative reconstruction algorithms were also examined. Sauve et al studied matrix inversion techniques for a model of their system. The geometrical symmetries of the system were exploited in an attempt to handle the otherwise intractably large set of equations. In another approach, that is relevant to the research reported in this paper, Wilderman et al developed the software for list-mode back-projection of Compton scatter events into the appropriate cone geometries. In this research, they were joined by Rohe et al who also developed backprojection algorithms for the cone geometry of Compton cameras. Unfortunately, simple backprojection onto the appropriate cones (without filtering) does not produce accurate reconstruction of the source distribution. 
   With the current art, reconstruction of a radiopharmaceutical source distribution typically utilizes iterative algorithms. The current art also provides analytic-inversion techniques that may be inherently faster to execute but are typically deficient with respect to relative immunity to statistical noise and the desired accuracy of reconstructing the source distribution. Hence, there is a real need to develop methods and apparatus that utilize analytic-inversion techniques while providing accuracy in reconstructing a radiopharmaceutical source distribution within a patient in order to properly diagnose the patient. 
   BRIEF SUMMARY OF THE INVENTION 
   The present invention provides methods and apparatuses for acquisition and processing of the observed data from a Compton camera in order to construct a three-dimensional image of a radiopharmaceutical source distribution within a patient. No intermediate two-dimensional images are formed. 
   Another aspect of the invention provides apparatuses that include a Compton camera for obtaining observed data from a radiopharmaceutical source distribution. The observed data is analyzed by a processor in order to construct a three-dimensional image representing the source distribution. The three-dimensional image may be displayed for a clinician. Moreover, the clinician may provide input data for configuring the processor for processing the observed data. 
   With another aspect of the invention, an idealized mathematical model expresses the observed Compton camera data in terms of an integral over the source distribution. An exact analytic inversion is then found for this idealized model. The new analytic solutions arise from a generalization of the integral equation used to model the Compton camera. The kernel of the integral equation is modified by the introduction of an index p that describes the effect of source distance from the camera. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
       FIG. 1  shows a basic Compton cones geometry in accordance with an embodiment of the invention. 
       FIG. 2A  shows a rectangular detector in accordance with an embodiment of the invention. 
       FIG. 2B  shows a square prism detector in accordance with an embodiment of the invention. 
       FIG. 2C  shows a cylindrical geometry in accordance with an embodiment of the invention. 
       FIG. 2D  shows a circular geometry in accordance with an embodiment of the invention. 
       FIG. 2E  shows a spherical detector in accordance with an embodiment of the invention. 
       FIG. 3  shows alternative detector D 2  geometries in accordance with an embodiment of the invention. 
       FIG. 4A  shows alternative scanning geometries for a linear scan in accordance with an embodiment of the invention. 
       FIG. 4B  shows alternative scanning geometries for an orbital scan in accordance with an embodiment of the invention. 
       FIG. 5  shows a conceptual flowchart for filtered backprojection (FBP) algorithms in accordance with an embodiment of the invention. 
       FIG. 6  shows a detailed flowchart of filtered backprojection (FBP) algorithms in accordance with an embodiment of the invention. 
       FIG. 7A  illustrates an example for the reconstruction of a three-dimensional source distribution, in which the y=0 slice of the reconstructed helix (for p=½ with the appropriate filter) is shown in accordance with an embodiment of the invention. 
       FIG. 7B  illustrates an example of the reconstruction of a three-dimensional source distribution, in which the y=0 slice of the reconstructed helix (for p=½ with the no filter) is shown in accordance with an embodiment of the invention. 
       FIG. 8  illustrates projections of the reconstructed helix, as shown in  FIGS. 7A and 7B , as the imaging volume is rotated by 90 degrees around the x axis in accordance with an embodiment of the invention. 
       FIG. 9  shows an apparatus for reconstructing a three-dimensional image representing a source distribution in accordance with an embodiment of the invention. 
       FIG. 10  shows a detector (D 1 ) flux in accordance with an embodiment of the invention. 
       FIG. 11  shows a detector (D 1 ) interaction probability in accordance with an embodiment of the invention. 
   

   DETAILED DESCRIPTION OF THE INVENTION 
   The following terminology and notation is referenced in the following discussion.
         BP backprojection operator   FBP filtered backprojection   FT Fourier transformation   LHS left hand side of an equation   RHS right hand side of an equation   voxel The smallest distinguishable box-shaped part of a three-dimensional space. A particular voxel will be identified by the x, y and z coordinates of one of its eight corners, or perhaps its centre. The term is generally used in three dimensional modeling and sometimes generalized to higher dimensions.   δ( ) Dirac delta function. The function is defined by the relation f(0)=∫dx f(x)δ(x) for all integrable functions f(x). An alternative definition is 
         δ   ⁡     (   x   )       =       lim     λ   →   0       ⁢       1         2   ⁢           ⁢   π       ⁢   λ       ⁢           ⁢       exp   ⁡     (     -       x   2       λ   2         )       .             
   {right arrow over (e)} x , {right arrow over (e)} y , {right arrow over (e)} z  unit vectors corresponding to the x, y, and z coordinates, respectively.       

   A mathematical model is proposed for Compton cameras. This model describes the expected data from a Compton camera in terms of the source distribution of radiopharmaceutical in the patient. In this model, a parameter p is introduced that describes the effect of the distance between the source and the primary (scatter) detector. A detailed justification of this model is given in Appendix A. The analysis in Appendix A indicates that the p=0 case corresponds most closely to actual Compton cameras. Backprojection operators (with a distance-dependent index q that is analogous to the index p of the model) are then introduced. The analysis of these backprojection operators requires the evaluation of an integral that is discussed in Appendix B. Appendix C demonstrates how one can select the index q to produce delta functions for p=½ and 3/2, thereby simplifying the task of the inversion and eliminating the complicated summations over spherical harmonics. The backprojection operators are combined with filters (that are dependent on camera geometry) and analytic inversions of the Compton camera models are found. Because the case p=0 is not solved, one must consider these analytic inversions as approximations. The numerical implementation and noise immunity of the algorithms are discussed and exemplary results are presented. The implications of these algorithms for the design of Compton camera geometry are discussed. 
     FIG. 1  shows a basic Compton cones geometry in accordance with an embodiment of the invention. The basic idea of electronic collimation is illustrated in  FIG. 1 . A gamma ray of energy E is emitted at source position (s)  101  (within the patient) and propagates to the primary detector (D 1 )  105  where it undergoes a Compton scattering that deposits energy E 1  at position (P 1 )  111 . The scattered gamma-ray propagates to point (P 2 )  107  in the secondary detector (D 2 )  109  where it is absorbed by a photoelectric interaction and deposits the remaining energy E 2 . (Generally, a detector with low atomic number Z is preferable for primary detector  105  so that Compton scatter is more likely; whereas, a detector with high Z is preferable for secondary detector  109  so that photoelectric absorption is more probable. In addition, primary detector  105  must be sufficiently thick for reasonable probability of Compton interaction, but not so thick that the scattered gamma ray is reabsorbed in primary detector  105 .) In such a Compton camera  100 , an event is registered whenever coincidence counts are observed in the primary and secondary detectors  105 , 109  that satisfy
   E=E   1   +E   2 .  (1) 
   [In reality, the energies can be measured with only limited accuracy, so that events are accepted if the sum (E 1 +E 2 ) is within a window surrounding E.] Thus, for each event, Compton camera  100  provides four parameters: the positions of interaction (P 1 )  111  and (P 2 )  107 , and the energies E 1  and E 2 . Simple relativistic kinematics yields the Compton relation 
                   m   e     ⁢           ⁢     E   1           (       E   1     +     E   2       )     ⁢     E   2         =     (     1   -     cos   ⁢           ⁢   θ       )             (   2   )             
 
where m e  is the electron mass (511 keV) and θ is the angle of the Compton scattering. Two important conclusions can be drawn from equation (2). First, because −1≦cos θ≦1, one finds that 
                     E   1     ≤       E   2       (       2   ⁢   E     +     m   e       )             and               m   e     ⁢           ⁢   E       (       2   ⁢   E     +     m   e       )       ≤     E   2                   (   3   )             
 
or, equivalently, 
                     E   1   2     +     2   ⁢     m   e     ⁢     E   1           -     E   1       ≤     2   ⁢       E   2     .               (   4   )             
 
   Inequality (4) provides an additional constraint on the acceptance of events, even in cases for which the incident energy E is unknown. Second, one observes that equation (2) provides an explicit relationship between the angle of scatter and the observed energies 
               σ   ≡     cos   ⁢           ⁢   θ       =       (       E   2   2     +       E   1     ⁢     E   2       -       m   e     ⁢           ⁢     E   1         )         (       E   1     +     E   2       )     ⁢           ⁢     E   2                 (   5   )             
 
   Compton camera  100  is based on this relationship; namely, the energies E 1  and E 2  determine the angle of scattering (θ)  103 . Thus, knowledge of the energies E 1  and E 2  implies that the incident gamma ray must have originated from one of the directions on the cone defined by σ=cos θ; so that, the source (s)  101  must lie on a cone  113  (as shown in  FIG. 1 ). The basic mathematical problem posed by Compton camera  100  is that individual events are associated with cone  113  and not a unique direction (as in SPECT or PET). Reconstruction of the source distribution from such information has proven significantly more difficult than the standard FBP algorithms for SPECT and PET. 
   In addition to relativistic kinematics, further details of the Compton scattering are required for the reconstruction of the source distribution from the data of Compton camera  100 . The probability of scattering at angle (θ)  103  is proportional to the differential scattering cross section. For a free electron, the Compton scattering cross-section is given by the Klein-Nishina formula: 
                 ⅆ       σ   KN     ⁡     (     θ   ,   E     )           ⅆ   Ω       =         r   e   2     2     ⁢       {         (     1   +       cos   2     ⁢           ⁢   θ       )     ⁡     [     1   +       E     m   e       ⁢     (     1   -     cos   ⁢           ⁢   θ       )         ]       +         E   2       m   e   2       ⁢       (     1   -     cos   ⁢           ⁢   θ       )     2         }         [     1   +       E     m   e       ⁢     (     1   -     cos   ⁢           ⁢   θ       )         ]     3                 (   6   )             
 
where r e  is the classical electron radius. In the following analysis, the weighting function W defined by 
                     W   ⁡     (     θ   ,   E     )       =       ⁢       2     r   e   2       ⁢       ⅆ       σ   KN     ⁡     (     θ   ,   E     )           ⅆ   Ω                     =       ⁢       {         (     1   +       cos   2     ⁢   θ       )     ⁡     [     1   +       E     m   e       ⁢     (     1   -     cos   ⁢           ⁢   θ       )         ]       +         E   2       m   e   2       ⁢       (     1   -     cos   ⁢           ⁢   θ       )     2         }         [     1   +       E     m   e       ⁢     (     1   -     cos   ⁢           ⁢   θ       )         ]     3                     (   7   )             
 
will be used. The function W gives the relative probability of a photon of energy E scattering at angle (θ)  103  and has the units of steradian −1 . Unfortunately, the electrons in detector (D 1 )  105  (from which the gamma-ray scatters) are bound, not free. The bound state of the atomic electrons causes two effects that influence the performance of Compton camera  100 . First, the electrons oscillate within the atoms of the detector leading to Doppler shifts that affect the relativistic kinematics of equation (5). This Doppler shift is characterized by momentum distributions [J(p z )] associated with the various atoms in the detector. This effect has been studied and the resultant errors in the scatter angle characterized. For low-energy gamma-rays, these errors can be significant. For this reason, materials that minimize the Doppler shift should be preferred for detector (D 1 )  105 . Second, scattering at small angles (θ)  103  is suppressed if the energy transferred to the recoiling electron is insufficient for ejection from the atom. The scatter event is forbidden if the excited electron is driven into a previously occupied state, so that the cross section becomes dependent on the occupation of outer electron shells of the atom (or the Fermi surface of a metal). Extensive tables of the incoherent scattering function (which characterizes this suppression for the small angle scattering) are available and can be incorporated into the weighting function W. In general, both these effects are minor if the initial energy E is large (E&gt;m e ). Despite potential problems for low energies (e.g., 140 keV emissions of Tc-99m), both effects are ignored in an embodiment of the invention.
 
   The Compton camera geometry shown in  FIG. 1  with a small primary detector (D 1 )  105  located in front of a large secondary detector (D 2 )  109  is a realistic model of the prototype geometry built by most researchers. The major limitation is the size of the D 1  detector  105  which is generally a small, pixelated solid-state detector. A standard gamma-ray camera is often used for the D 2  detector  109 . Since the first development of Compton cameras, researchers have realized that this geometry is insufficient for 3D reconstruction. The problem can be solved by scanning the camera around the patient, thereby effectively increasing the area of the D 1  detector  105 . Experimentalists tend to rotate the camera around the patient in a circular orbit in emulation of SPECT. Some theoreticians extend the D 1  and D 2  detectors  105 , 109  to infinity and discuss planar detectors. Others surround the patient with spherical or cylindrical banks of detectors. Still others rotate infinite planar detectors around the patient. Because each geometry potentially requires a different reconstruction algorithm, one might expect a plethora of algorithms. However, the filtered backprojection algorithms, in an embodiment of the invention, display a significant invariance. The backprojection operation is unchanged by the detector geometry; whereas, the filter used on the backprojections depends on the geometry of primary detector (D 1 )  105  and its scanning trajectory. Consequently, one can make a list of relevant geometries (or scanning motions) for primary detector (D 1 )  105  and associate a filter with each geometry. This strategy provides a unified structure for development and evaluation of algorithms based on the geometry of primary detector (D 1 )  105 . Five of these geometries are illustrated in  FIGS. 2A–2E . 
     FIG. 2A  shows a rectangular detector in accordance with an embodiment of the invention.  FIG. 2B  shows a square prism detector in accordance with an embodiment of the invention.  FIG. 2C  shows a cylindrical geometry in accordance with an embodiment of the invention.  FIG. 2D  shows a circular geometry in accordance with an embodiment of the invention.  FIG. 2E  shows a spherical detector in accordance with an embodiment of the invention. Although all of the geometries in  FIGS. 2A–2E  have been discussed in the literature, few are currently considered practical.  FIG. 2A  shows one of the more practical designs, namely, a rectangular D 1  detector  201  with dimensions L×Δ.  FIG. 2B  shows a (4-sided) rectilinear box detector  203  with two open ends and each side L×(2R c ).  FIG. 2C  is one of the favorites in the literature: D 1  is a cylinder with radius R c  and axial length L.  FIG. 2D  shows circular D 1  detector  207  with radius L. The limit L→∞ of  FIG. 2D  corresponds to the infinite planar detector that is popular with theorists. Finally,  FIG. 2E  shows a spherical detector consisting of concentric spheres for detectors D 1  and D 2   209 , 211 . 
   Possible D 1  geometries are restricted by the size and shape of the patient. For this work, the patient (and, therefore, the source of radiation) will be restricted to a truncated cylinder of radius B and axial length A. One observes that B&lt;A for human subjects. For some applications (i.e. cardiac or brain imaging), the axial length A may be truncated (to the chest or head region). The D 1  detector is characterized by a radius R c . The D 1  detector is restricted to regions outside a cylinder of radius R c &gt;B for all geometries except spherical. For spherical detectors, the D 1  radius, R c , must be larger than √{square root over (A 2 +B 2 )} where A is the full height of the patient (so that the full patient can fit into the spherical detector). 
   The secondary D 2  detectors are also subject to an important requirement. Not all scatter events in the D 1  detector will enter the D 2  detector. In particular, gamma-rays that are scattered back toward the source will miss the D 2  detector. Obviously, one wants a large D 2  detector (as shown in  FIG. 1 ) so that as many events as possible will be observed. As will be discussed, the equations describing Compton camera performance imply an important relation between the D 1  and D 2  detectors. At each point {right arrow over (z)}εD 1 , the D 2  detector is visible (unobstructed by any intervening material) in a solid angle Ω({square root over (z)}, D 2 ). If this solid angle Ω ({right arrow over (z)}, D 2 ) is a hemisphere for every point {right arrow over (z)}ε=D 1 , then the information is sufficient for the filtered backprojection algorithms in an embodiment of the invention. The specific geometry of the D 2  detector is not important. The only restriction on the D 2  geometry is that, from every point in the D 1  detector, the D 2  detector must be capable of recording events from a complete hemisphere of scatter directions. In  FIG. 1 , this requirement is satisfied only if the planar D 2  detector extends to infinity. The algorithms in an embodiment of the invention can still be implemented if this condition on D 2  geometry is violated, but the resulting reconstructions may be suboptimal. 
     FIG. 3  illustrates three equivalent D 2  geometries  301 , 303 , 305  that satisfy the requirement for a small rectangular D 1  detector  351 . All the detector geometries in  FIG. 2A–2E  satisfy this requirement. 
     FIGS. 4A and 4B  illustrate that the D 1  detector geometries in  FIGS. 2A and 2C  can be accomplished without the construction of special detector arrays. The same geometry can be achieved by scanning with a smaller detector. Furthermore, the geometry of rotating a small Δ×Δ detector around the patient is equivalent to the cylindrical geometry of  FIG. 2C  with L=Δ. 
   An embodiment of the invention provides a set of FBP algorithms that predict the distribution of radiopharmaceutical source distribution within a patient from the data acquired by Compton camera  100 . For mathematical simplicity, the emissions of the radiopharmaceutical source distribution will be assumed monoenergetic with energy E. [This assumption can be dropped without affecting the final imaging equation (11). However, the notation becomes more complicated.] The distribution of radiopharmaceuticals is described by a function f ({right arrow over (x)}) that denotes the density of radiopharmaceuticals in the patient. The function f has dimensions cm − 3 (emissions/cm 3 ) and is a bounded, positive function with support within a cylinder of radius B and axial length A (i.e., the radiopharmaceuticals are contained within the body of the patient). 
   The data generated by Compton camera  100  are a collection of N coincidence events from the primary and secondary detectors  105 , 109 . For SPECT studies, one typically collects N≈10 7  counts. Whereas, for Compton camera  100 , one expects (at least) a one hundred-fold increase in sensitivity, so that N≈10 9 . Each event provides four parameters that describe the coincident detection: the position z of the Compton scattering in primary (D 1 ) detector  105 , the position w of the photoelectric absorption in secondary (D 2 ) detector  109 , and the energies E 1  and E 2  deposited in these detectors. In summary, a coincident event provides two positions and energies: ({right arrow over (z)}, E 1 ; {right arrow over (w)}, E 2 ) where {right arrow over (z)}εD 1  and {right arrow over (w)}εD 2 . For purposes of data analysis, these data will be transformed into a new set of parameters:
 
({right arrow over (z)}, E 1 ; {right arrow over (w)}, E 2 )           ({right arrow over (z)}, {right arrow over (n)}, σ, E)
 
by the algebraic relations 
                 n   -&gt;     =       (       w   -&gt;     -     z   -&gt;       )              w   -&gt;     -     z   -&gt;                ,     E   =       E   1     +     E   2         ,       and   ⁢           ⁢   σ     =       (       E   2   2     +       E   1     ⁢     E   2       -       m   e     ⁢     E   1         )         (       E   1     +     E   2       )     ⁢     E   2                   (   8   )             
 
{right arrow over (z)} represents the position of Compton scattering in the primary detector, {right arrow over (n)} denotes the direction of the γ-ray after scattering, σ corresponds to the cosine of the angle between incident and scattered directions of the γ-ray, and E corresponds to the energy of the incident γ-ray. Each event is characterized by four parameters ({right arrow over (z)}, {right arrow over (n)}, σ, E), such that
 
 {right arrow over (z)}εD   1 , | {right arrow over (n)}|= 1, and −1&lt;σ&lt;1.  (9)

   Thus, the two forms of the data have the same number of dimensions (6) and the parameters E 1  and E 2  are interchangeable with σ and E. However, the vectors {right arrow over (w)} and {right arrow over (n)} are subtly related by the geometry of the detector system; {right arrow over (w)}εD 2  whereas |{right arrow over (n)}|=1. The vector {right arrow over (n)} denotes the direction of the γ-ray propagation after scatter. At each point {right arrow over (z)}εD 1  some directions {right arrow over (n)} will not intersect the detector D 2   109  and, therefore, are not sampled by the detector system. For this reason, Compton camera  100  is often called a “limited angle” imaging system and the data are called “incomplete.” From this argument, many researchers have concluded that backprojection algorithms are not suitable for Compton cameras and that alternative (e.g., maximum likelihood) algorithms are the only appropriate method for reconstruction. One of the important ideas, in an embodiment of the invention, is that the data are complete if the vector {right arrow over (n)} is sampled on a hemisphere at each point {right arrow over (z)}εD 1  [see the discussion following equation (12)]. Many current Compton camera designs meet this criterion. 
   [If the energy E of the incident γ-ray is known, the energy detected in secondary detector D 2   109  becomes redundant. As a result, one could accept scatter events in detector D 2   109 , rather than requiring total absorption in a photoelectric interaction. This strategy is often proposed for high-energy radiation because the probability of photoelectric absorption in the secondary detector drops for higher energies. However, higher rates of accidental coincidence counts may result. In an embodiment of the invention, the reconstruction algorithms require the cosine of the scattering angle and the energy of the incident γ-ray E. The method by which these parameters are determined does not affect the algorithms.] 
   The list of events from the camera is next combined into a histogram of events g ({right arrow over (z)}, {right arrow over (n)}, σ, E) that describes the imaging process. Because {right arrow over (n)} and σ are expressed in dimensionless parameters (steradians and cosines, respectively), the resulting histogram ĝ function has dimensions cm −2  corresponding to the density of observed events at the position {right arrow over (z)} on the 2D surface D 1 . The observed function g ({right arrow over (z)}, {right arrow over (n)}, σ, E) is related to the source distribution f ({right arrow over (x)}) by the integral equation: 
                     g   ⁡     (       z   -&gt;     ,     n   -&gt;     ,   σ   ,   E     )       =       ⁢       W   ⁡     (     σ   ,   E     )       ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     x   -&gt;       ⁢           ⁢     f   ⁡     (     x   -&gt;     )       ⁢           ⁢     1     2   ⁢           ⁢   π   ⁢           ⁢     R   c   p     ⁢              z   -&gt;     -     x   -&gt;              2   -   p                                     ⁢     δ   ⁡     [           n   -&gt;     ·     (       z   -&gt;     -     x   -&gt;       )                z   -&gt;     -     x   -&gt;              -   σ     ]                   =       ⁢       W   ⁡     (     σ   ,   E     )       ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     x   -&gt;       ⁢           ⁢     f   ⁡     (     x   -&gt;     )       ⁢           ⁢     1     2   ⁢           ⁢   π   ⁢           ⁢     R   c   p     ⁢              z   -&gt;     -     x   -&gt;              1   -   p                                     ⁢       δ   ⁡     [         n   -&gt;     ·     (       z   -&gt;     -     x   -&gt;       )       -     σ   ⁢            z   -&gt;     -     x   -&gt;                ]       .                   (   10   )             
 
   The physical justification and limitations of equation (10) are discussed in Appendix A. For all realistic imaging processes, the parameter p should vanish; the case p=0 implies that the radiation dies off as r −2  as D 1  detector  105  moves away from the source. For technical reasons, however, analytic inversion for the p=0 case (and more generally any integer value of p) is difficult. On the other hand, the p=±½ and 3/2 cases yield analytic solutions. The imaging equation (10) is a generalization that includes these analytic inversions. The function W is the dimensionless weighting function [defined in equation (7)] that accounts for the different count rates at various scattering angles. One can immediately absorb the function W into the LHS of equation (10) and write that 
                       g   ⁡     (       z   -&gt;     ,     n   -&gt;     ,   σ   ,   E     )         W   ⁡     (     σ   ,   E     )         =       ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     x   -&gt;       ⁢           ⁢     f   ⁡     (     x   -&gt;     )       ⁢           ⁢     1     2   ⁢           ⁢   π   ⁢           ⁢     R   c   p     ⁢              z   -&gt;     -     x   -&gt;              2   -   p                                   ⁢     δ   ⁡     [           n   -&gt;     ·     (       z   -&gt;     -     x   -&gt;       )                z   -&gt;     -     x   -&gt;              -   σ     ]                   =       ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     x   -&gt;       ⁢           ⁢     f   ⁡     (     x   -&gt;     )       ⁢           ⁢     1     2   ⁢           ⁢   π   ⁢           ⁢     R   c   p     ⁢              z   -&gt;     -     x   -&gt;              1   -   p                                   ⁢     δ   ⁡     [         n   -&gt;     ·     (       z   -&gt;     -     x   -&gt;       )       -     σ   ⁢            z   -&gt;     -     x   -&gt;                ]                   ≡       ⁢           C   ^     p     ⁡     [   f   ]       .                   (   11   )             
 
   All energy dependence in equation (10) arises from the function W; so that, the resulting ratio g/W is independent of the total energy E. Equation (11) will be called the Compton imaging equation with index p and Ĉ p  is the Compton camera operator. This notation is slightly misleading because, once the function W is factored out, equation (11) actually describes an isotropic scattering process. However, the important idea, namely that the angle of scatter is determined by the energies, is retained. Straightforward dimensional analysis confirms that f has units of cm −3  and g has units cm −2 . One should also note the symmetry 
                 g   ⁡     (       z   -&gt;     ,     n   -&gt;     ,   σ   ,   E     )         W   ⁡     (     σ   ,   E     )         =       g   ⁡     (       z   -&gt;     ,     -     n   -&gt;       ,     -   σ     ,   E     )         W   ⁡     (       -   σ     ,   E     )                 (   12   )             
 
   This symmetry is very important because, at any point {right arrow over (z)}, in the primary detector, half of the directions {right arrow over (n)} are generally excluded from measurement. Equation (12) guarantees that the remaining directions {right arrow over (n)} are sufficient for the determination of g over the full domain of {right arrow over (n)}. Thus, equation (12) is the basis of the geometric condition that detector (D 2 )  109  must cover a hemisphere of directions {right arrow over (n)} for each point in detector (D 1 )  105 . Another important property of equation (11) is that the integral on the RHS is a convolution integral in ({right arrow over (z)}−{right arrow over (x)}). Normally, convolution integrals are easily inverted with a Fourier transformation (FT). Unfortunately, the function g is known only on detector (D 1 )  105 , not over the entire 3D range of {right arrow over (z)} and, therefore, cannot be Fourier transformed. Thus, the most obvious inversion technique is not feasible. 
     FIG. 5  shows a conceptual flowchart for filtered backprojection (FBP) algorithms in accordance with an embodiment of the invention. The inverse problem for a Compton camera [of index p] can now be simply stated: one must find the function f in equation (11) that generates the observed function g. The inversion strategy employed, in an embodiment of the invention, is the standard backprojection strategy as illustrated in  FIG. 5 . The source distribution f ({right arrow over (x)})  501  produces the Compton camera data g ({right arrow over (z)}, {right arrow over (n)}, σ, E)  503  by experimental observation [as modeled by equation (10) with index p]. The inversion requires the selection of an appropriate combination of weighted backprojection operator (BP)  505  and filter function (Φ)  507 . An important step in this process is the selection of the correct backprojection operator. For the success of this inversion strategy, the FT of the backprojected data must be proportional to the FT of the original source distribution. The relation between these FTs is an integral equation [which is given in equation (24)] and the strategy can only succeed if the kernel in this integral equation is proportional to a delta function. The introduction of the apparently ad hoc weighting indices p [in equation (10)] and q [in equation (13) below] is required because the correct combination of p and q yields a delta function in equation (24) that renders filtered backprojection feasible. 
   In an embodiment of the invention, the backprojection operator (with index q) of g is defined by the equation: 
                         BP     q     ⁡     [     g   W     ]       =       ⁢       b   q     ⁡     (     s   -&gt;     )                   ≡       ⁢       ∫   D1             ⁢         ⅆ   2     ⁢     z   -&gt;       ⁢           ⁢       ∫            n   -&gt;          =   1               ⁢         ⅆ   2     ⁢     n   -&gt;       ⁢           ⁢       ∫     -   1     1     ⁢       ⅆ   σ     ⁢       g   ⁡     (       z   -&gt;     ,     n   -&gt;     ,   σ   ,   E     )         W   ⁡     (     σ   ,   E     )                                     ⁢       1     2   ⁢           ⁢   π   ⁢           ⁢     R   c     1   +   q       ⁢              z   -&gt;     -     s   -&gt;              2   -   q           ⁢           ⁢     δ   ⁡     [           n   -&gt;     ·     (       z   -&gt;     -     s   -&gt;       )                z   -&gt;     -     s   -&gt;              -   σ     ]                       (   13   )             
 
where D 1  corresponds to the surface of primary detector  105 . The integrals in equation (13) are analogous to the standard backprojection done in SPECT and PET; except that, rather than backprojecting the data onto a line, one backprojects onto a cone (described by the delta function). [Mathematically, the operator {circumflex over (B)}{circumflex over (P)} p  is the adjoint of Ĉ p .] The index q is analogous to the index p in the imaging equation (11) and permits weighting of the backprojection with various powers of |{right arrow over (s)}−{right arrow over (z)}| as one projects the data inward from the point of scattering {right arrow over (z)} in detector D 1   105 . The motivation for the introduction of the index q will become apparent later in the calculation [see equations (29) and (30) and Appendix C]. The domain ({right arrow over (s)}) of the functions b q  is the 3D region corresponding to the possible support of f and the functions b q  have units of cm −3  for all values of q. Thus, all the functions b q  are candidates for the inversion of Compton camera imaging operator Ĉ p . The integral over {right arrow over (n)} in equation (13) requires integration over the entire 2-sphere. As previously noted, the function g is seldom observed for {right arrow over (n)} over the full 2-sphere. However, if the function g is observed over a hemisphere, one can use the symmetry of equation (12) for evaluation of g in the unobservable directions. In practical terms, this means that the integral over {right arrow over (n)} in equation (13) can be evaluated using an integral over the “observed” directions {right arrow over (n)} on a hemisphere; the contribution to equation (13) from the “unobserved” directions merely doubles the value of the resultant integral.
 
   The backprojection operators will be examined using the functions h q,p  defined by 
                       h     q   ,   p       ⁡     (     s   -&gt;     )       ≡       ⁢         BP     q     ⁢           ⁢         C   ^     p     ⁡     [   f   ]                     =       ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     x   -&gt;       ⁢           ⁢     f   ⁡     (     x   -&gt;     )       ⁢           ⁢       ∫   D1             ⁢       ⅆ   2     ⁢     z   -&gt;                                 ⁢     1         (     2   ⁢           ⁢   π     )     2     ⁢     R   c     1   +   q   +   p       ⁢           ⁢              z   -&gt;     -     s   -&gt;              2   -   q       ⁢           ⁢              z   -&gt;     -     x   -&gt;              2   -   p                           ⁢       ∫            n   -&gt;          =   1               ⁢         ⅆ   2     ⁢     n   -&gt;       ⁢           ⁢       ∫     -   1     1     ⁢       ⅆ   σ     ⁢           ⁢     δ   ⁡     [           n   -&gt;     ·     (       z   -&gt;     -     s   -&gt;       )                z   -&gt;     -     s   -&gt;              -   σ     ]               ⁢                           ⁢       δ   ⁡     [           n   -&gt;     ·     (       z   -&gt;     -     x   -&gt;       )                z   -&gt;     -     x   -&gt;              -   σ     ]       .                   (   14   )             
 
   One might naively hope that h q,p =f for some p and q because, in that case, {circumflex over (B)}{circumflex over (P)} q  would be the inverse operator for Ĉ p . Despite the fact that such hopes are not met, the evaluation of the functions h q,p  will lead to useful FBP algorithms in the same manner that filtered backprojection arises in SPECT and PET. Performing the integral over a in equation (14), one finds 
                       h     q   ,   p       ⁡     (     s   -&gt;     )       =       ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     x   -&gt;       ⁢           ⁢     f   ⁡     (     x   -&gt;     )       ⁢           ⁢       ∫   D1             ⁢       ⅆ   2     ⁢     z   -&gt;                                 ⁢     1         (     2   ⁢           ⁢   π     )     2     ⁢     R   c     1   +   q   +   p       ⁢           ⁢              s   -&gt;     -     z   -&gt;              2   -   q       ⁢           ⁢              x   -&gt;     -     z   -&gt;              2   -   p                           ⁢       ∫            n   -&gt;          =   1               ⁢         ⅆ   2     ⁢     n   -&gt;       ⁢           ⁢       δ   ⁡     [           n   -&gt;     ·     (       z   -&gt;     -     s   -&gt;       )                s   -&gt;     -     z   -&gt;              -         n   -&gt;     ·     (       z   -&gt;     -     x   -&gt;       )                x   -&gt;     -     z   -&gt;                ]       ⁢           .                       (   15   )             
 
   Next, the integral over the 2-sphere {right arrow over (n)} can be evaluated by a simple calculation 
                       ∫            n   -&gt;          =   1               ⁢         ⅆ   2     ⁢     n   -&gt;       ⁢           ⁢     δ   ⁡     (       n   -&gt;     ·     v   -&gt;       )           =       ⁢       ∫   0     2   ⁢           ⁢   π       ⁢       ⅆ   φ     ⁢           ⁢       ∫   0   π     ⁢       ⅆ   θ     ⁢           ⁢   sin   ⁢           ⁢   θ   ⁢           ⁢     δ   ⁡     (            v   -&gt;          ⁢           ⁢   cos   ⁢           ⁢   θ     )                           =       ⁢       ∫   0     2   ⁢           ⁢   π       ⁢       ⅆ   φ     ⁢           ⁢       ∫     -   1     1     ⁢       ⅆ   μ     ⁢           ⁢     δ   ⁡     (            v   -&gt;          ⁢           ⁢   μ     )                           =       ⁢     2   ⁢           ⁢   π   ⁢           ⁢       ∫     -   1     1     ⁢       ⅆ   μ     ⁢           ⁢     δ   ⁡     (            v   -&gt;          ⁢           ⁢   μ     )                         =       ⁢           2   ⁢           ⁢   π            v   -&gt;            ⁢       ∫     -   1     1     ⁢           ⁢     ⅆ     μδ   ⁡     (   μ   )             =       2   ⁢           ⁢   π            v   -&gt;                            (   16   )             
 
so that in equation (15) one finds 
                       h     q   ,   p       ⁡     (     s   -&gt;     )       =       ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     x   -&gt;       ⁢           ⁢     f   ⁡     (   x   )       ⁢           ⁢       ∫   D1             ⁢       ⅆ   2     ⁢     z   -&gt;                                 ⁢       1     2   ⁢           ⁢   π   ⁢           ⁢     R   c     1   +   q   +   p       ⁢              s   -&gt;     -     z   -&gt;              2   -   q       ⁢           ⁢              x   -&gt;     -     z   -&gt;              2   -   p       ⁢              (       x   -&gt;     -     z   -&gt;       )              x   -&gt;     -     z   -&gt;              -       (       s   -&gt;     -     z   -&gt;       )              s   -&gt;     -     z   -&gt;                         .                   (   17   )             
 
   Next, the 3D Fourier transformations (F and H q,p ) of the functions f and h q,p  are defined by the equations
 
 F ({right arrow over (ν)})=∫∫∫ d   3   {right arrow over (x)}f ( {right arrow over (x)} )exp [−2 πi ({right arrow over (ν)}· {right arrow over (x)} )]  (18)
 
and  H   q,p ({right arrow over (ω)})=∫∫∫ d   3   {right arrow over (s)}h   q,p ( {right arrow over (s)} )exp [−2 πi ({right arrow over (ω)}· {right arrow over (s)} )]  (19)
 
   Thus, equation (17) can be written in the spatial-frequency domain as 
                       H     q   ,   p       ⁡     (     ω   -&gt;     )       =       ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     v   -&gt;       ⁢           ⁢     F   ⁡     (     v   -&gt;     )       ⁢           ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     x   -&gt;       ⁢     ∫     ∫     ∫       ⅆ   3     ⁢     s   -&gt;                                             ⁢       exp   ⁡     [       2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢       v   -&gt;     ·     x   -&gt;         -     2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢       ω   -&gt;     ·     s   -&gt;           ]       ⁢           ⁢       ∫   D1             ⁢       ⅆ   2     ⁢     z   -&gt;                           ⁢       1     2   ⁢           ⁢   π   ⁢           ⁢     R   c     1   +   q   +   p       ⁢              s   -&gt;     -     z   -&gt;              2   -   q       ⁢           ⁢              x   -&gt;     -     z   -&gt;              2   -   p       ⁢              (       x   -&gt;     -     z   -&gt;       )              x   -&gt;     -     z   -&gt;              -       (       s   -&gt;     -     z   -&gt;       )              s   -&gt;     -     z   -&gt;                         .                   (   20   )             
 
   One can translate the dummy variables of integration {right arrow over (s)} and {right arrow over (x)},
 
{right arrow over (x)}→{right arrow over (x)}+{right arrow over (z)} and {right arrow over (s)}→{right arrow over (s)}+{right arrow over (z)},
 
so that 
                       H     q   ,   p       ⁡     (     ω   -&gt;     )       =       ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     v   -&gt;       ⁢           ⁢     F   ⁡     (     v   -&gt;     )       ⁢           ⁢       ∫   D1             ⁢         ⅆ   2     ⁢     z   -&gt;       ⁢           ⁢     exp   ⁡     [     2   ⁢           ⁢   π   ⁢           ⁢       ⅈ   ⁡     (       v   -&gt;     -     ω   -&gt;       )       ·     z   -&gt;         ]                                   ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     x   -&gt;       ⁢           ⁢     ∫     ∫     ∫       ⅆ   3     ⁢     s   -&gt;                                     ⁢         exp   ⁡     [       2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢       v   -&gt;     ·     x   -&gt;         -     2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢       ω   -&gt;     ·     s   -&gt;           ]         2   ⁢           ⁢   π   ⁢           ⁢     R   c     1   +   q   +   p       ⁢           ⁢            s   -&gt;            2   -   q       ⁢           ⁢            x   -&gt;            2   -   p       ⁢           ⁢              x   -&gt;            x   -&gt;            -       s   -&gt;            s   -&gt;                       .                   (   21   )             
 
   Thus, one finds that the integrand on the RHS of equation (21) contains two independent terms: 
                 M     q   ,   p       ⁡     (       v   -&gt;     ,     ω   -&gt;       )       =     ∫     ∫     ∫         ⅆ   3     ⁢     x   -&gt;       ⁢           ⁢     ∫     ∫     ∫       ⅆ   3     ⁢     s   -&gt;                             (   22   )                       ⁢       exp   ⁡     [       2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢       v   -&gt;     ·     x   -&gt;         -     2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢       ω   -&gt;     ·     s   -&gt;           ]         2   ⁢           ⁢   π   ⁢           ⁢     R   c     1   +   q   +   p       ⁢           ⁢            s   -&gt;            2   -   q       ⁢           ⁢            x   -&gt;            2   -   p       ⁢           ⁢              x   -&gt;            x   -&gt;            -       s   -&gt;            s   -&gt;                                             and                             D   ⁡     (     v   -&gt;     )       =       ∫   D1             ⁢         ⅆ   2     ⁢     z   -&gt;       ⁢           ⁢     exp   ⁡     (       -   2     ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢       v   -&gt;     ·     z   -&gt;         )             ,           (   23   )               such   ⁢           ⁢   that                               H     q   ,   p       ⁡     (     ω   -&gt;     )       =     ∫     ∫     ∫         ⅆ   3     ⁢     v   -&gt;       ⁢           ⁢     F   ⁡     (     v   -&gt;     )       ⁢           ⁢       M     q   ,   p       ⁡     (       v   -&gt;     ,     ω   -&gt;       )       ⁢           ⁢     D   ⁡     (       ω   -&gt;     -     v   -&gt;       )                       (   24   )             
 
   The function M q,p  contains all the important information about the backprojection operation, whereas the function D is a simple Fourier transform over the surface of primary detector D 1   105 . In the subsequent analysis, one finds that the properties of M q,p  determine the weighting index q of the backprojection operator and the function D determines the appropriate filter function. The function M q,p  will be examined first. 
   The function M q,p  is evaluated by explicit integration in Appendix B. Unfortunately, the integral in equation (22) is only convergent if −1&lt;p&lt;1 and −1&lt;q&lt;1. The lower limits on convergence (−1&lt;p and −1&lt;q) arise from the small distance behavior of the integrand (|{right arrow over (x)}|→0 and |{right arrow over (s)}|→0) and are intrinsic to the problem. However, the upper limits (p&lt;1 and q&lt;1) arise from the large distance behavior of the integrand (|{right arrow over (x)}|→∞ and |{right arrow over (s)}|→∞). For all practical applications, both the source distribution f and the backprojections h q,p  have compact support, so that the behavior at extremely large distances should be unimportant. This suggests that the integral can be regularized and the divergences at large distances suppressed. In particular, one can choose 
                       M     q   ,   p       ⁡     (       v   -&gt;     ,     ω   -&gt;       )       =       ⁢       lim     ɛ   →   0       ⁢     ∫     ∫     ∫         ⅆ   3     ⁢     x   -&gt;       ⁢           ⁢     ∫     ∫     ∫       ⅆ   3     ⁢     s   -&gt;                                       ⁢       exp   [       2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢       v   -&gt;     ·     x   -&gt;         -     2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢       ω   -&gt;     ·     s   -&gt;         -     ɛ   ⁢       (              x   -&gt;            +            s   -&gt;              )         R   c             ]       2   ⁢           ⁢   π   ⁢           ⁢     R   c     1   +   q   +   p       ⁢            s   -&gt;            2   -   q       ⁢           ⁢            x   -&gt;            2   -   p       ⁢           ⁢              x   -&gt;            x   -&gt;            -       s   -&gt;            s   -&gt;                                     (   25   )             
 
which is convergent for −1&lt;p and −1&lt;q. For all positive values of             the additional term in the exponential eliminates the divergences in the integral as |{right arrow over (x)}|→∞ and |{right arrow over (s)}|→∞. Furthermore, as shown in Appendix B, the limit as ε→0 is well-defined. The regularized integral in equation (25) will henceforth be used as the definition of M q,p . This regularization is significant because the cases q= 3/2 and 5/2 will be important in the subsequent inversions as will be discussed.

   For mathematical convenience, a unit vector on the 2-sphere will be denoted by 
                       k   -&gt;            k   -&gt;            ≡       ⁢       Ω   -&gt;     k                 =       ⁢       Ω   -&gt;     ⁡     (       θ   k     ,     φ   k       )                   =       ⁢       cos   ⁢           ⁢     θ   k     ⁢           ⁢       e   -&gt;     3       +     sin   ⁢           ⁢     θ   k     ⁢           ⁢   sin   ⁢           ⁢     φ   k     ⁢           ⁢       e   -&gt;     2       +     sin   ⁢           ⁢     θ   k     ⁢           ⁢   cos   ⁢           ⁢     φ   k     ⁢           ⁢         e   -&gt;     1     .                       (   26   )             
 
   The evaluation of the function M q,p  in Appendix B yields that 
                       M     q   ,   p       ⁡     (       v   -&gt;     ,     ω   -&gt;       )       =       ⁢         π     1   -   p   -   q           R   c     1   +   q   +   p       ⁢            ω   -&gt;            q   +   1       ⁢           ⁢            v   -&gt;            p   +   1           ⁢       ∑     n   =   0     ∞     ⁢       2     (       2   ⁢   n     +   1     )       ⁢       Γ   ⁡     (       n   2     +     q   2     +     1   2       )         Γ   ⁡     (       n   2     -     q   2     +   1     )                               ⁢         Γ   ⁡     (       n   2     +     p   2     +     1   2       )         Γ   ⁡     (       n   2     -     p   2     +   1     )         ⁢           ⁢       ∑     m   =     -   n       n     ⁢           Y   n   m     (       Ω   -&gt;     v     )     ⁡     [       Y   n   m     ⁡     (       Ω   -&gt;     ω     )       ]       *                       (   27   )             
 
for −1&lt;p and −1&lt;q, where the functions Y n   m ({right arrow over (Ω)} ν )=Y n   m (θ ν , φ ν ) are the spherical harmonics. Recalling that 
                   ∑     n   =   0     ∞     ⁢       ∑     m   =     -   n       n     ⁢           Y   n   m     ⁡     (       Ω   -&gt;     ν     )       ⁢           [       Y   n   m     ⁡     (       Ω   -&gt;     ω     )       ]     *         =       δ   2     ⁡     (         Ω   -&gt;     ν     -       Ω   -&gt;     ω       )         ,           (   28   )             
 
one observes (Appendix C) that 
                 M       3   /   2     ,     1   /   2         ⁡     (       ν   -&gt;     ,     ω   -&gt;       )       =         δ   2     ⁡     (         Ω   -&gt;     ν     -       Ω   -&gt;     ω       )         2   ⁢           ⁢   π   ⁢           ⁢     R   c   3     ⁢           ⁢            ω   -&gt;            5   /   2       ⁢           ⁢            ν   -&gt;            3   /   2                   (   29   )             and                             M       1   /   2     ,     3   /   2         ⁡     (       ν   -&gt;     ,     ω   -&gt;       )       =           δ   2     ⁡     (         Ω   -&gt;     ν     -       Ω   -&gt;     ω       )         2   ⁢           ⁢   π   ⁢           ⁢     R   c   3     ⁢           ⁢            ω   -&gt;            3   /   2       ⁢           ⁢            ν   -&gt;            5   /   2           .             (   30   )             
 
   The delta functions in equations (29) and (30) are extremely suggestive. The inversion algorithms for the Compton camera equation with p=½ and 3/2 are a direct result of these equations. Indeed, the discovery of functions M q,p  (or combinations of such functions) that produce delta functions of the form δ 2 ({right arrow over (Ω)} ν −{right arrow over (Ω)} ω ) is an important idea in an embodiment of the invention. 
   For the case p=−½, the functions M 5/2,−1/2  and M 1/2,−1/2  prove useful. In particular, one finds that 
                 M       5   /   2     ,       -   1     /   2         ⁡     (       ν   -&gt;     ,     ω   -&gt;       )       =       1     2   ⁢           ⁢   π   ⁢           ⁢     R   c   3     ⁢           ⁢            ω   -&gt;            7   /   2       ⁢           ⁢            ν   -&gt;            1   /   2           ⁢           ⁢       ∑     n   =   0     ∞     ⁢       (       n   2     +   n   -     3   4       )       (       n   2     +   n   +     1   4       )                   (   31   )                       ⁢       ∑     m   =     -   n       n     ⁢           Y   n   m     ⁡     (       Ω   -&gt;     v     )       ⁢           [       Y   n   m     ⁡     (       Ω   -&gt;     ω     )       ]     *                             and                             M       1   /   2     ,       -   1     /   2         ⁡     (       ν   -&gt;     ,     ω   -&gt;       )       =         4   ⁢           ⁢     π   2     ⁢           ⁢     R   c   2         2   ⁢           ⁢   π   ⁢           ⁢     R   c   3     ⁢           ⁢            ω   -&gt;            3   /   2       ⁢           ⁢            ν   -&gt;            1   /   2           ⁢           ⁢       ∑     n   =   0     ∞     ⁢     1     (       n   2     +   n   +     1   4       )                   (   32   )                       ⁢       ∑     m   =     -   n       n     ⁢           Y   n   m     ⁡     (       Ω   -&gt;     v     )       ⁢           [       Y   n   m     ⁡     (       Ω   -&gt;     ω     )       ]     *                               so   ⁢           ⁢   that                                 M       5   /   2     ,       -   1     /   2         ⁡     (       ν   -&gt;     ,     ω   -&gt;       )       +       1       (     2   ⁢           ⁢   π   ⁢           ⁢     R   c     ⁢           ⁢          ω   -&gt;            )     2       ⁢       M       1   /   2     ,       -   1     /   2         ⁡     (       ν   -&gt;     ,     ω   -&gt;       )           =           (   33   )                     δ   2     ⁡     (         Ω   -&gt;     v     -       Ω   -&gt;     ω       )         2   ⁢           ⁢   π   ⁢           ⁢     R   c   3     ⁢           ⁢            ω   -&gt;            7   /   2       ⁢           ⁢            ν   -&gt;            1   /   2           .                         
 
   This last identity provides the basis of inversion for the p=−½ case. 
   The three equations (29), (30) and (33) (with various parameters p and q) all yield the same basic expression that will be used in the subsequent inversions. For convenience, that mathematical expression will be defined as 
                   M   ^     k     ⁡     (       ν   -&gt;     ,     ω   -&gt;       )       ≡         δ   2     ⁡     (         Ω   -&gt;     v     -       Ω   -&gt;     ω       )         2   ⁢           ⁢   π   ⁢           ⁢     R   c   3     ⁢            ω   -&gt;            2   +   k       ⁢           ⁢            ν   -&gt;            2   -   k                   (   34   )             
 
   For those combinations of p and q that produce {circumflex over (M)} k , one can rewrite equation (24) as 
                   H   ^     k     ⁡     (     ω   -&gt;     )       =     ∫     ∫     ∫         ⅆ   3     ⁢     ν   -&gt;       ⁢           ⁢     F   ⁡     (     ν   -&gt;     )       ⁢           ⁢         M   ^     k     ⁡     (       ν   -&gt;     ,     ω   -&gt;       )       ⁢           ⁢     D   ⁡     (       ω   -&gt;     -     ν   -&gt;       )                       (   35   )               so   ⁢           ⁢   that                             4   ⁢           ⁢   π   ⁢           ⁢     R   c   3     ⁢          ω          2   +   k       ⁢           ⁢         H   ^     k     ⁡     (     ω   ⁢           ⁢       Ω   -&gt;     ω       )         =       ∫     -   ∞     ∞     ⁢       ⅆ   ν     ⁢           ⁢          ν        k     ⁢           ⁢     F   ⁡     (     ν   ⁢           ⁢       Ω   -&gt;     ω       )       ⁢           ⁢     D   ⁡     [       (     ω   -   ν     )     ⁢           ⁢       Ω   -&gt;     ω       ]                   (   36   )             
 
   [In equation (36) the identity 
                     ∫     ∫     ∫       ⅆ   3     ⁢     ν   -&gt;             =       ∫   0   ∞     ⁢       ⅆ   ν     ⁢           ⁢          ν        2     ⁢           ⁢     ∫     ∫       ⅆ   2     ⁢       Ω   -&gt;     ν                           =       1   2     ⁢       ∫     -   ∞     ∞     ⁢       ⅆ   ν     ⁢           ⁢          ν        2     ⁢           ⁢     ∫     ∫       ⅆ   2     ⁢       Ω   -&gt;     ν                               (   37   )             
 
is used so that the range of the integral is (−∞, ∞) rather than [0, ∞).] An important observation is that the RHS of equation (36) is a 1D convolution integral. Thus, standard Fourier techniques permit the analytic inversion of equation (36) and the recovery of the function F from the backprojected data, Ĥ k . The function D is obviously important in such an inversion.
 
   The function D is determined by the geometry of the primary detector D 1   105 . The definition of D [equation (25)] suggests that, for many geometries, D will contain delta functions. Such delta functions may simplify the integral in equation (36). However, not all geometries will produce delta functions and care is required in many geometries that do produce delta functions. 
   In an embodiment of the invention, five basic detector geometries are examined ( FIGS. 2A–2E ). In Table 1, each of these detector geometries is described and the function D displayed. 
   
     
       
             
           
             
             
             
             
           
         
             
               TABLE 1 
             
           
           
             
                 
             
             
               The function D(v) for different camera geometries 
             
           
        
         
             
               Geometry 
               D1 Surface 
               D({right arrow over (v)}) 
               
                 
                   
                     
                       
                         
                           
                             
                               
                                 
                                   lim 
                                   
                                     L 
                                     → 
                                     ∞ 
                                   
                                 
                                 ⁢ 
                                 
                                   D 
                                   ( 
                                   
                                     v 
                                     → 
                                   
                                   ) 
                                 
                               
                               ⁢ 
                               
                                 
 
                               
                               ⁢ 
                               
                                 ( 
                                 
                                   L 
                                   ⪢ 
                                   
                                     4 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       max 
                                       ⁡ 
                                       
                                         ( 
                                         
                                           A 
                                           , 
                                           
                                             R 
                                             c 
                                           
                                         
                                         ) 
                                       
                                     
                                   
                                 
                                 ) 
                               
                             
                           
                         
                       
                         
                     
                   
                 
               
             
             
                 
             
             
               1. Rectangular 
               −L/2 &lt; z &lt; L/2, −Δ/2 ≦ y ≦ Δ/2, 
               D 1 ({right arrow over (v)}) = ΔLsinc(πΔ{right arrow over (v)} · {right arrow over (e)} y  )sinc(πL{right arrow over (v)} · {right arrow over (e)} z  ) × 
               D 1 ({right arrow over (v)}) = Δδ({right arrow over (v)} · {right arrow over (e)} y  )sinc(πL{right arrow over (v)} · {right arrow over (e)} z  ) × 
             
             
               Detector 
               x = R c   
               exp(2πiR c {right arrow over (v)} · {right arrow over (e)} x  ) 
               exp(2πiR c {right arrow over (v)} · {right arrow over (e)} x  ) 
             
             
                 
             
             
               2. SquarePrism Detector 
               −L/2 &lt; z &lt; L/2, −R c  ≦ y ≦ R c ,x = ± R c −L/2 &lt; z &lt; L/2, −R c  ≦ x ≦ R c ,y = ± R c   
               
                 
                   
                     
                       
                         
                           
                             
                               
                                 
                                   D 
                                   2 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     v 
                                     → 
                                   
                                   ) 
                                 
                               
                               = 
                                 
                               ⁢ 
                               
                                 4 
                                 ⁢ 
                                 
                                   R 
                                   c 
                                 
                                 ⁢ 
                                 L 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 sin 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   c 
                                   ( 
                                   
                                     π 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     L 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       
                                         v 
                                         → 
                                       
                                       · 
                                       
                                         
                                           e 
                                           → 
                                         
                                         z 
                                       
                                     
                                   
                                   ) 
                                 
                                 × 
                               
                             
                           
                         
                         
                           
                             
                                 
                               ⁢ 
                               
                                 [ 
                                 
                                   
                                     
                                       
                                         
                                           sin 
                                           ⁢ 
                                           
                                               
                                           
                                           ⁢ 
                                           
                                             c 
                                             ( 
                                             
                                               2 
                                               ⁢ 
                                               
                                                 πR 
                                                 c 
                                               
                                               ⁢ 
                                               
                                                 
                                                   v 
                                                   → 
                                                 
                                                 · 
                                                 
                                                   
                                                     e 
                                                     → 
                                                   
                                                   y 
                                                 
                                               
                                             
                                             ) 
                                           
                                           ⁢ 
                                           
                                             cos 
                                             ⁡ 
                                             
                                               ( 
                                               
                                                 2 
                                                 ⁢ 
                                                 
                                                   πR 
                                                   c 
                                                 
                                                 ⁢ 
                                                 
                                                   
                                                     v 
                                                     → 
                                                   
                                                   · 
                                                   
                                                     
                                                       e 
                                                       → 
                                                     
                                                     x 
                                                   
                                                 
                                               
                                               ) 
                                             
                                           
                                         
                                         + 
                                       
                                     
                                   
                                   
                                     
                                       
                                         sin 
                                         ⁢ 
                                         
                                             
                                         
                                         ⁢ 
                                         
                                           c 
                                           ( 
                                           
                                             2 
                                             ⁢ 
                                             
                                               πR 
                                               c 
                                             
                                             ⁢ 
                                             
                                               
                                                 v 
                                                 → 
                                               
                                               · 
                                               
                                                 
                                                   e 
                                                   → 
                                                 
                                                 x 
                                               
                                             
                                           
                                           ) 
                                         
                                         ⁢ 
                                         
                                           cos 
                                           ⁡ 
                                           
                                             ( 
                                             
                                               2 
                                               ⁢ 
                                               
                                                 πR 
                                                 c 
                                               
                                               ⁢ 
                                               
                                                 
                                                   v 
                                                   → 
                                                 
                                                 · 
                                                 
                                                   
                                                     e 
                                                     → 
                                                   
                                                   y 
                                                 
                                               
                                             
                                             ) 
                                           
                                         
                                       
                                     
                                   
                                 
                                 ] 
                               
                             
                           
                         
                       
                         
                     
                   
                 
               
               
                 
                   
                     
                       
                         
                           
                             
                               
                                 
                                   D 
                                   2 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     v 
                                     → 
                                   
                                   ) 
                                 
                               
                               = 
                                 
                               ⁢ 
                               
                                 4 
                                 ⁢ 
                                 
                                   R 
                                   c 
                                 
                                 ⁢ 
                                 
                                   δ 
                                   ( 
                                   
                                       
                                   
                                   ⁢ 
                                   
                                     
                                       v 
                                       → 
                                     
                                     · 
                                     
                                       
                                         e 
                                         → 
                                       
                                       z 
                                     
                                   
                                   ) 
                                 
                                 × 
                               
                             
                           
                         
                         
                           
                             
                                 
                               ⁢ 
                               
                                 [ 
                                 
                                   
                                     
                                       
                                         
                                           sin 
                                           ⁢ 
                                           
                                               
                                           
                                           ⁢ 
                                           
                                             c 
                                             ( 
                                             
                                               2 
                                               ⁢ 
                                               
                                                 πR 
                                                 c 
                                               
                                               ⁢ 
                                               
                                                 
                                                   v 
                                                   → 
                                                 
                                                 · 
                                                 
                                                   
                                                     e 
                                                     → 
                                                   
                                                   y 
                                                 
                                               
                                             
                                             ) 
                                           
                                           ⁢ 
                                           
                                             cos 
                                             ⁡ 
                                             
                                               ( 
                                               
                                                 2 
                                                 ⁢ 
                                                 
                                                   πR 
                                                   c 
                                                 
                                                 ⁢ 
                                                 
                                                   
                                                     v 
                                                     → 
                                                   
                                                   · 
                                                   
                                                     
                                                       e 
                                                       → 
                                                     
                                                     x 
                                                   
                                                 
                                               
                                               ) 
                                             
                                           
                                         
                                         + 
                                       
                                     
                                   
                                   
                                     
                                       
                                         sin 
                                         ⁢ 
                                         
                                             
                                         
                                         ⁢ 
                                         
                                           c 
                                           ( 
                                           
                                             2 
                                             ⁢ 
                                             
                                               πR 
                                               c 
                                             
                                             ⁢ 
                                             
                                               
                                                 v 
                                                 → 
                                               
                                               · 
                                               
                                                 
                                                   e 
                                                   → 
                                                 
                                                 x 
                                               
                                             
                                           
                                           ) 
                                         
                                         ⁢ 
                                         
                                           cos 
                                           ⁡ 
                                           
                                             ( 
                                             
                                               2 
                                               ⁢ 
                                               
                                                 πR 
                                                 c 
                                               
                                               ⁢ 
                                               
                                                 
                                                   v 
                                                   → 
                                                 
                                                 · 
                                                 
                                                   
                                                     e 
                                                     → 
                                                   
                                                   y 
                                                 
                                               
                                             
                                             ) 
                                           
                                         
                                       
                                     
                                   
                                 
                                 ] 
                               
                             
                           
                         
                       
                         
                     
                   
                 
               
             
             
                 
             
             
               3. CylindricalDetector 
               −L/2 &lt; z &lt; L/2, x 2  + y 2  = R c   2   
               
                 
                   
                     
                       
                         
                           
                             
                               
                                 
                                   D 
                                   3 
                                 
                                 ( 
                                 
                                   v 
                                   → 
                                 
                                 ) 
                               
                               = 
                                 
                               ⁢ 
                               
                                 2 
                                 ⁢ 
                                 
                                   πR 
                                   c 
                                 
                                 ⁢ 
                                 L 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 sin 
                                 ⁢ 
                                 
                                     
                                 
                                 ⁢ 
                                 
                                   c 
                                   ( 
                                   
                                     πL 
                                     ⁢ 
                                     
                                       
                                         v 
                                         → 
                                       
                                       · 
                                       
                                         
                                           e 
                                           → 
                                         
                                         z 
                                       
                                     
                                   
                                   ) 
                                 
                                 × 
                               
                             
                           
                         
                         
                           
                             
                                 
                               ⁢ 
                               
                                 
                                   J 
                                   0 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     2 
                                     ⁢ 
                                     
                                       πR 
                                       c 
                                     
                                     ⁢ 
                                     
                                       
                                         
                                           
                                             ( 
                                             
                                               
                                                 v 
                                                 → 
                                               
                                               · 
                                               
                                                 
                                                   e 
                                                   → 
                                                 
                                                 x 
                                               
                                             
                                             ) 
                                           
                                           2 
                                         
                                         + 
                                         
                                           
                                             ( 
                                             
                                               
                                                 v 
                                                 → 
                                               
                                               · 
                                               
                                                 
                                                   e 
                                                   → 
                                                 
                                                 y 
                                               
                                             
                                             ) 
                                           
                                           2 
                                         
                                       
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                         
                       
                         
                     
                   
                 
               
               
                 
                   
                     
                       
                         
                           
                             
                               
                                 
                                   D 
                                   3 
                                 
                                 ( 
                                 
                                   v 
                                   → 
                                 
                                 ) 
                               
                               = 
                                 
                               ⁢ 
                               
                                 2 
                                 ⁢ 
                                 
                                   πR 
                                   c 
                                 
                                 ⁢ 
                                 
                                   δ 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       
                                         v 
                                         → 
                                       
                                       · 
                                       
                                         
                                           e 
                                           → 
                                         
                                         z 
                                       
                                     
                                     ) 
                                   
                                 
                                 × 
                               
                             
                           
                         
                         
                           
                             
                                 
                               ⁢ 
                               
                                 
                                   J 
                                   0 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     2 
                                     ⁢ 
                                     
                                       πR 
                                       c 
                                     
                                     ⁢ 
                                     
                                       
                                         
                                           
                                             ( 
                                             
                                               
                                                 v 
                                                 → 
                                               
                                               · 
                                               
                                                 
                                                   e 
                                                   → 
                                                 
                                                 x 
                                               
                                             
                                             ) 
                                           
                                           2 
                                         
                                         + 
                                         
                                           
                                             ( 
                                             
                                               
                                                 v 
                                                 → 
                                               
                                               · 
                                               
                                                 
                                                   e 
                                                   → 
                                                 
                                                 y 
                                               
                                             
                                             ) 
                                           
                                           2 
                                         
                                       
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                         
                       
                         
                     
                   
                 
               
             
             
                 
             
             
               4. CircularDetectors 
               0 ≦ (y 2  + z 2 ) ≦ L 2 , x = R c   
               
                 
                   
                     
                       
                         
                           
                             
                               
                                 
                                   
                                     
                                       D 
                                       4 
                                     
                                     ( 
                                     
                                       v 
                                       → 
                                     
                                     ) 
                                   
                                   = 
                                     
                                   ⁢ 
                                   
                                     4 
                                     ⁢ 
                                     
                                       πL 
                                       2 
                                     
                                     ⁢ 
                                     
                                       
                                         
                                           J 
                                           1 
                                         
                                         ⁡ 
                                         
                                           ( 
                                           
                                             2 
                                             ⁢ 
                                             πL 
                                             ⁢ 
                                             
                                               
                                                 
                                                   
                                                     ( 
                                                     
                                                       
                                                         v 
                                                         → 
                                                       
                                                       · 
                                                       
                                                         
                                                           e 
                                                           → 
                                                         
                                                         y 
                                                       
                                                     
                                                     ) 
                                                   
                                                   2 
                                                 
                                                 + 
                                                 
                                                   
                                                     ( 
                                                     
                                                       
                                                         v 
                                                         → 
                                                       
                                                       · 
                                                       
                                                         
                                                           e 
                                                           → 
                                                         
                                                         z 
                                                       
                                                     
                                                     ) 
                                                   
                                                   2 
                                                 
                                               
                                             
                                           
                                           ) 
                                         
                                       
                                       
                                         2 
                                         ⁢ 
                                         πL 
                                         ⁢ 
                                         
                                           
                                             
                                               
                                                 ( 
                                                 
                                                   
                                                     v 
                                                     → 
                                                   
                                                   · 
                                                   
                                                     
                                                       e 
                                                       → 
                                                     
                                                     y 
                                                   
                                                 
                                                 ) 
                                               
                                               2 
                                             
                                             + 
                                             
                                               
                                                 ( 
                                                 
                                                   
                                                     v 
                                                     → 
                                                   
                                                   · 
                                                   
                                                     
                                                       e 
                                                       → 
                                                     
                                                     z 
                                                   
                                                 
                                                 ) 
                                               
                                               2 
                                             
                                           
                                         
                                       
                                     
                                   
                                 
                                   
                               
                               × 
                             
                           
                         
                         
                           
                             
                                 
                               ⁢ 
                               
                                 exp 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     2 
                                     ⁢ 
                                     
                                       πiR 
                                       c 
                                     
                                     ⁢ 
                                     
                                       
                                         v 
                                         → 
                                       
                                       · 
                                       
                                         
                                           e 
                                           → 
                                         
                                         x 
                                       
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                         
                       
                         
                     
                   
                 
               
               
                 
                   
                     
                       
                         
                           
                             
                               
                                 
                                   
                                     
                                       D 
                                       4 
                                     
                                     ( 
                                     
                                       v 
                                       → 
                                     
                                     ) 
                                   
                                   = 
                                     
                                   ⁢ 
                                   
                                     4 
                                     ⁢ 
                                     L 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       δ 
                                       ⁡ 
                                       
                                         ( 
                                         
                                           
                                             
                                               
                                                 ( 
                                                 
                                                   
                                                     v 
                                                     → 
                                                   
                                                   · 
                                                   
                                                     
                                                       e 
                                                       → 
                                                     
                                                     y 
                                                   
                                                 
                                                 ) 
                                               
                                               2 
                                             
                                             + 
                                             
                                               
                                                 ( 
                                                 
                                                   
                                                     v 
                                                     → 
                                                   
                                                   · 
                                                   
                                                     
                                                       e 
                                                       → 
                                                     
                                                     z 
                                                   
                                                 
                                                 ) 
                                               
                                               2 
                                             
                                           
                                         
                                         ) 
                                       
                                     
                                   
                                 
                                   
                               
                               × 
                             
                           
                         
                         
                           
                             
                                 
                               ⁢ 
                               
                                 exp 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     2 
                                     ⁢ 
                                     
                                       πiR 
                                       c 
                                     
                                     ⁢ 
                                     
                                       
                                         v 
                                         → 
                                       
                                       · 
                                       
                                         
                                           e 
                                           → 
                                         
                                         x 
                                       
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                         
                       
                         
                     
                   
                 
               
             
             
                 
             
             
               5. SphericalDetector 
               x 2  + y 2  + z 2  = R c   2   
               D 5 ({right arrow over (v)}) = 4πR c   2 sinc(2πR c |{right arrow over (v)}|) 
               D 5 ({right arrow over (v)}) = 4πR c   2 sinc(2πR c |{right arrow over (v)}|) 
             
             
                 
             
           
        
       
     
   
   Detector geometry  1  (first row of Table 1) is a rectangular area with dimensions L×□. The actual detector need not have physical dimensions L×Δ, one could scan a much smaller detector over the region ( FIG. 4A ) and the resulting function D would be the same. Detector geometry  2  (second row of Table 1) is a square prism (an elongated box without top or bottom). This geometry is a simple application of the first geometry; the function D can be evaluated by a summation of the terms for the four sides. Detector geometry  3  (third row of Table 1) is a cylinder of radius R c  and length L. Once again, the primary detector D 1   105  need not be a cylinder if the surface D 1  is sampled by scanning. One might build a ring detector and scan it along the axis or one might build a narrow strip detector  403  and rotate it around the axis (as shown in  FIG. 4B ). Either configuration would yield the same function D. Detector geometry  4  (fourth row of Table 1) is a large circular detector of radius L. Finally, detector geometry  5  (bottom row of Table 1) is a spherical detector with radius R c  surrounding the imaging region. 
   Detector geometries  1  through  4  will be examined in the limit as L→∞. Recalling that 
                 lim     L   →   ∞       ⁢     [     L   ⁢           ⁢   sinc   ⁢           ⁢     (     π   ⁢           ⁢   L   ⁢           ⁢   ν     )       ]       =     δ   ⁡     (   ν   )               (   38   )             and                             lim     L   →   ∞       ⁢     [     π   ⁢           ⁢   L   ⁢           ⁢         J   1     ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢   L   ⁢           ⁢   ν     )         2   ⁢           ⁢   π   ⁢           ⁢   L   ⁢           ⁢   ν         ]       =     δ   ⁡     (   ν   )               (   39   )             
 
one observes that these limits yield the delta functions needed for the simple inversion of equation (36). In reality, of course, no detector can be extended to infinity. However, because the region of source reconstruction has finite dimensions (cylindrical radius R c  and axial length A where generally R c &lt;A), one need not sample the spatial frequency in steps smaller than (Δν)≈1/(2A). The width of the sinc function in equation (38) is approximately (Δν) 2/L. Thus, if L&gt;4A, then L should be sufficiently large that the function D can be approximated as a delta function. [A similar argument applies for the Bessel function in equation (39).] In particular, the last column of Table 1 displays the L→∞ limit of the function D for detector geometries  1 – 4 . For the spherical detector geometry ( 5 ), no limit yields a delta function.
 
   Equation (36) provides the prototypical equation for all the inversion algorithms in an embodiment of the invention. Using the delta function approximations for D i  of detector geometries  1  through  3 , one finds that 
                       D   i     ⁡     [       (     ω   -   ν     )     ⁢           ⁢       Ω   -&gt;     ω       ]       ⁢           =       K   i     ⁢           ⁢     δ   ⁡     [       (     ω   -   ν     )     ⁢           ⁢     (         Ω   -&gt;     ω     ·       e   -&gt;     z       )       ]                     =         K   i     ⁢           ⁢   δ   ⁢           ⁢     (     ω   -   ν     )                  Ω   -&gt;     ω     ·       e   -&gt;     z                            (   40   )             
 
where the constants K i  are given by
 
K 1 =Δ, K 2 =8R c , and K 3 =2πR c .  (41)
 
   Substitution of equation (40) into equation (36) yields 
               4   ⁢           ⁢   π   ⁢           ⁢     R   c   3     ⁢          ω          2   +   k       ⁢           ⁢         H   ^     k     ⁡     (     ω   ⁢           ⁢       Ω   -&gt;     ω       )         =           K   i     ⁢          ω        k                  Ω   -&gt;     ω     ·       e   -&gt;     z              ⁢           ⁢     F   ⁡     (     ω   ⁢           ⁢       Ω   -&gt;     ω       )                 (   42   )             or                           F   ⁡     (     ν   -&gt;     )       =         4   ⁢           ⁢   π   ⁢           ⁢     R   c   3         K   i       ⁢           ⁢          ν   -&gt;          ⁢           ⁢            ν   -&gt;     ·       e   -&gt;     z            ⁢           ⁢         H   ^     k     ⁡     (     ν   -&gt;     )                 (   43   )             
 
   Equation (43) provides a filtered backprojection algorithm for the source distribution F in detector geometries  1 ,  2 , and  3 . [The second equality in equation (40) presumes that |{right arrow over (ω)} ω ·{right arrow over (e)} z |≠0. Consequently, the results in equation (43) are not rigorous on the |{right arrow over (ν)}·{right arrow over (e)} z |=0 plane. Nonetheless, equation (43) is satisfied for {right arrow over (ν)} arbitrarily close to the |{right arrow over (ν)}·{right arrow over (e)} z |=0 plane. Thus, the continuity of the function F suggests that equation (43) will also be valid on the |{right arrow over (ν)}·{right arrow over (e)} z |=0 plane.] Surprisingly, the algorithm does not depend on the index k because the same power appears on both sides of equation (42) and is canceled. For circular detector geometry  4 , the delta function approximation yields a similar result, namely, 
               F   ⁡     (     ν   -&gt;     )       =         π   ⁢           ⁢     R   c   3       L     ⁢           ⁢          ν   -&gt;          ⁢           ⁢           (       ν   -&gt;     ·       e   -&gt;     y       )     2     +       (       ν   -&gt;     ·       e   -&gt;     z       )     2         ⁢           ⁢         H   ^     k     ⁡     (     ν   -&gt;     )                 (   44   )             
 
   Equations (43) and (44) demonstrate that, if the detector geometry produces a function D containing a delta function, a simple filtered backprojection algorithm follows immediately. 
   Unfortunately, the spherical detector geometry (5) does not produce a function D 5  containing a delta function. The substitution of D 5  (from Table 1) in equation (36) gives the 1D convolution equation 
                     2   ⁢     R   c   2     ⁢          ω          2   +   k       ⁢         H   ^     k     ⁡     (     ω   ⁢           ⁢       Ω   -&gt;     ω       )         =       ⁢       ∫     -   ∞     ∞     ⁢       ⅆ   ν     ⁢           ⁢          ν        k     ⁢           ⁢     F   ⁡     (     ν   ⁢           ⁢       Ω   -&gt;     ω       )                           ⁢     2   ⁢           ⁢     R   c     ⁢           ⁢   sin   ⁢           ⁢     c   ⁡     [     2   ⁢           ⁢   π   ⁢           ⁢       R   c     ⁡     (     ω   -   ν     )         ]                       (   45   )             
 
   Fourier analysis provides a straightforward method for inversion of such a convolution equation. However, a property of the sinc function provides an even more powerful method. If one defines 
               s   ⁡     (   ω   )       ≡     2   ⁢     R   c     ⁢           ⁢   sin   ⁢           ⁢     c   ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢     R   c     ⁢   ω     )                 (   46   )             then                           s   ⁡     (   ν   )       =       ∫     -   ∞     ∞     ⁢       ⅆ   ω     ⁢           ⁢     s   ⁡     (     ν   -   ω     )       ⁢           ⁢     s   ⁡     (   ω   )                   (   47   )             
 
   [Equation (47) follows immediately from the fact the Fourier transform of s is equal to unity on the interval (−R c , R c ) and vanishes elsewhere. As a result, the square of the Fourier transform of s is equal to itself. By the convolution theorem, therefore, s must satisfy equation (47).] Thus, by taking a convolution integral with s on both sides of equation (45), one finds that 
                       ∫     -   ∞     ∞     ⁢       ⅆ   ω     ⁢           ⁢   2   ⁢     R   c   2     ⁢           ⁢          ω          2   +   k       ⁢           ⁢         H   ^     k     ⁡     (     ω   ⁢           ⁢       Ω   -&gt;     ω       )       ⁢           ⁢     s   ⁡     (     ν   -   ω     )           =                 ∫     -   ∞     ∞     ⁢       ⅆ   ω     ⁢           ⁢          ω        k     ⁢           ⁢     F   ⁡     (     ω   ⁢           ⁢       Ω   -&gt;     ω       )       ⁢           ⁢     s   ⁡     (     ν   -   ω     )                       (   48   )             
 
   The obvious solution of equation (48) is
 
 F ({right arrow over (ω)})=2 R   c   2 |{right arrow over (ω)}| 2   Ĥ   k ({right arrow over (ω)})  (49)
 
   Of course, solutions of the homogeneous integral equation 
                 ∫     -   ∞     ∞     ⁢       ⅆ   ω     ⁢           ⁢          ω        k     ⁢           ⁢       F   homo     ⁡     (     ω   ⁢           ⁢       Ω   -&gt;     ω       )       ⁢           ⁢     s   ⁡     (     ν   -   ω     )           =   0           (   50   )             
 
can be added to equation (49). The selection of a particular homogeneous solution requires an additional criterion. For all practical purposes, the homogeneous term can be assumed to vanish and equation (49) gives the resulting solution. Equation (48) implies a consistency condition on the function Ĥ k ; namely, 
                            ν          2   +   k       ⁢         H   ^     k     ⁡     (     ν   ⁢           ⁢       Ω   -&gt;     ω       )         =       ⁢       ∫     -   ∞     ∞     ⁢       ⅆ   ω     ⁢           ⁢          ω          2   +   k       ⁢           ⁢         H   ^     k     ⁡     (     ω   ⁢           ⁢       Ω   -&gt;     ω       )                           ⁢     2   ⁢           ⁢     R   c     ⁢           ⁢   sin   ⁢           ⁢       c   ⁡     [     2   ⁢           ⁢   π   ⁢           ⁢       R   c     ⁡     (     ν   -   ω     )         ]       .                     (   51   )             
 
   Filtered backprojection algorithms can now be stated for the Compton imaging equation of indices p=±½ and 3/2. These FBP algorithms are based on five assumptions that have been stated earlier, are implicit in the previous analysis, and are repeated here for completeness. First, the data from the Compton camera take the form of a function g ({right arrow over (z)}, {right arrow over (n)}, σ, E) that is known on the domain defined in equation (9). Second, the primary detector surface D 1  is assumed to be one of the five geometries described in Table 1 ( FIGS. 2A–2E ). Third, the source distribution f ({right arrow over (x)}) has support only within a compact region. For detector geometries 1–4, the support of f is a cylindrical region of radius R c  and axial length A. For spherical detector geometry (5), the support of f is a spherical region of radius R c . Fourth, the parameter L in geometries 1 through 4 must satisfy the condition L&gt;&gt;4 max (A, R c ). Finally, the functions g and f are related by the Compton imaging equation (11) with index p=±½ or 3/2; that is 
                 g   ⁡     (       z   -&gt;     ,     n   -&gt;     ,   σ   ,   E     )         W   ⁡     (     σ   ,   E     )         =         C   ^     p     ⁡     [   f   ]               (   11   )             
 
   Of course, the FBP algorithms will be different for the various indices p and detector geometries. All the FBP algorithms can be divided into four steps. The first step involves the backprojection of the data and is dependent only on the index p. The second step requires the 3D Fourier transformation of the backprojected images. The third step involves the filtering of the backprojection data and depends only on the detector geometry. The final step is a simple 3D Fourier transformation that recovers the source distribution f and is independent of either index p or detector geometry. 
     FIG. 6  shows a detailed flowchart of filtered backprojection (FBP) algorithms in accordance with an embodiment of the invention.  FIG. 6  delineates the steps (steps  603 ,  605 ,  607 , and  609 ) that constitute FBP algorithms in accordance with an embodiment of the invention. 
   Step  1  (step  603 ): Backprojections for Index p 
   In step  603 , the Compton camera data g ({right arrow over (z)}, {right arrow over (n)}, σ, E) (representing input data  601  and generated by Compton imaging equation of index p) are mapped by backprojection into various functions b q  ({right arrow over (x)})={circumflex over (B)}{circumflex over (P)} q [g/W] [equation (13)]. For each index p, different indices q are required for the backprojected functions b q  ({right arrow over (x)}). In particular, the algorithms for the three cases p=±½ and 3/2 require the following backprojections:
 
 p=  3/2 requires  b   1/2 ( {right arrow over (x)} ) [= h   1/2, 3/2 ( {right arrow over (x)} )],
 
 p= ½ requires  b   3/2 ( {right arrow over (x)} ) [= h   3/2, 1/2 ( {right arrow over (x)} )],
 
and
 
 p=− ½ requires both  b   5/2 ( {right arrow over (x)} ) [= h   5/2,−1/2 ( {right arrow over (x)} )]
 
and  b   1/2 ( {right arrow over (x)} )[= h   1/2, −1/2 ( {right arrow over (x)} )].
 
Step  2  (Step  605 ): 3D Fourier Transformation of the Backprojected Data
 
   In step  605 , the backprojected data are Fourier transformed and produce an intermediate function H ({right arrow over (ν)}). For notational convenience, the 3D Fourier transform of the backprojected data is defined by
 
 B   q ({right arrow over (ν)})=∫∫∫ d   3   {right arrow over (x)}b   q ( {right arrow over (x)} )exp(−2 πi{right arrow over (ν)}·{right arrow over (x)} )  (52)
 
   Algorithms for the three cases p=±½ and 3/2 are then given by:
 
 p=  3/2 H ({right arrow over (ν)})= B   1/2 ({right arrow over (ν)}) [= H   1/2, 3/2 ({right arrow over (ν)})],  (53)
 
 p= ½ H ({right arrow over (ν)})= B   3/2 ({right arrow over (ν)}) [= H   3/2, 1/2 ({right arrow over (ν)})].  (54)
 
and,  p=− ½
 
                     H   ⁡     (     ν   -&gt;     )       =       ⁢         B     5   /   2       ⁡     (     ν   -&gt;     )       +       1       (     2   ⁢           ⁢   π   ⁢           ⁢     R   c     ⁢          ν   -&gt;            )     2       ⁢       B     1   /   2       ⁡     (     ν   -&gt;     )                           ⁢       [     =         H       5   /   2     ,       -   1     /   2         ⁡     (     ν   -&gt;     )       +       1       (     2   ⁢           ⁢   π   ⁢           ⁢     R   c     ⁢          ν   -&gt;            )     2       ⁢       H       1   /   2     ,       -   1     /   2         ⁡     (     ν   -&gt;     )             ]     .                   (   55   )             
 
   For each value of p, the function H corresponds to the combination of functions H q,p  that produced {circumflex over (M)} k  in equations (29), (30) and (33). The function H, therefore, satisfies the integral equation (36) as previously discussed. 
   Step  3  (step  607 ): Geometric Filtering 
   In step  607 , the intermediate function H is converted to the Fourier transform of the source distribution, F, by the appropriate filtering (as previously discussed). The filtering is dependent on the D 1  detector geometry. Thus, the function F is given by
 
 F ({right arrow over (ν)})=Φ({right arrow over (ν)}) H ({right arrow over (ν)})  (56)
 
where the function Φ({right arrow over (ν)}) is a filter function determined by the D 1  geometry.
 
For detector geometries (i=1–3) (corresponding to  FIGS. 2A ,  2 B, and  2 C):
 
[L&gt;&gt;4 max (A, R c )], the function Φ is given by 
               Φ   ⁡     (     ν   -&gt;     )       =         4   ⁢           ⁢   π   ⁢           ⁢     R   c   3         K   i       ⁢           ⁢          ν   -&gt;          ⁢           ⁢            ν   -&gt;     ·       e   -&gt;     z                      (   57   )             
 
where K i  are given by equations (41).
 
For detector geometry ( 4 ) (corresponding to  FIG. 2D ):
 
[L&gt;&gt;4 max (A, R c )], the function Φ is given by 
               Φ   ⁡     (     v   -&gt;     )       =         π   ⁢           ⁢     R   c   3       L     ⁢           ⁢          v   -&gt;          ⁢           ⁢           (       v   -&gt;     ·       e   -&gt;     y       )     2     +       (       v   -&gt;     ·       e   -&gt;     z       )     2                   (   58   )             
 
For the spherical detector geometry ( 5 ) (corresponding to  FIG. 2E ),
 
the function Φ is given by
 
Φ({right arrow over (ν)})=2 R   c   2 |{right arrow over (ν)}| 2 .  (59)
 
Step  4  (step  609 ): Inverse 3D Fourier Transformation
 
   In step  609 , the source distribution f is recovered from F by inverse Fourier transformation [of equation (18)]:
 
 f ( {right arrow over (x)} )=∫∫∫ d   3   {right arrow over (ν)}F ({right arrow over (ν)})exp [2 πi ({right arrow over (ν)}·{right arrow over ( x )})]  (60)
 
   Steps  601 – 609  constitute the FBP algorithms in accordance with an embodiment of the invention. In condensed operator notation, these algorithms can be written as 
                       C   ^       3   /   2       -   1       =       ℱ   3     -   1       ⁢       [     Φ   ⁢     (     v   -&gt;     )       ]       ︷   Filter       ⁢           ⁢     ℱ   3     ⁢           ⁢         BP       1   /   2         ︷     Backprojection   Weighted                           C   ^       1   /   2       -   1       =         ℱ   3     -   1       ⁡     [     Φ   ⁡     (     v   -&gt;     )       ]       ⁢           ⁢     ℱ   3     ⁢           ⁢       BP       3   /   2                         C   ^         -   1     /   2       -   1       =           ℱ   3     -   1       ⁡     [     Φ   ⁡     (     v   -&gt;     )       ]       ⁢           ⁢     ℱ   3     ⁢           ⁢       BP       5   /   2         +         ℱ   3     -   1       ⁡     [       Φ   ⁡     (     v   -&gt;     )           (     2   ⁢           ⁢   π   ⁢           ⁢     R   c     ⁢           ⁢          v   -&gt;            )     2       ]       ⁢           ⁢     ℱ   3     ⁢           ⁢       BP       1   /   2                         (   61   )             
 
where              3  denotes the 3D Fourier transformation and Φ({right arrow over (ν)}) is the filter function associated with the D 1  detector geometry.

   Throughout the derivation of these algorithms, Compton camera data  601  are tabulated as a histogram function g ({right arrow over (z)}, {right arrow over (n)}, σ, E). In reality, of course, the data are typically accumulated as a list of N individual (coincident) events of the form
 
( {right arrow over (z)}   i   , {right arrow over (n)}   i , σ i   , E   i ) for  i= 1 , . . . , N 
 
   Organizing these events into a histogram g ({right arrow over (z)}, {right arrow over (n)}, σ, E) is one method of handling the data. An alternative strategy is typically more practical. If one estimates the function g ({right arrow over (n)}, {right arrow over (n)}, σ, E) by 
                 g   ⁡     (       z   -&gt;     ,     n   -&gt;     ,   σ   ,   E     )       =       ∑     i   =   1     N     ⁢       δ   ⁡     (       z   -&gt;     -       z   -&gt;     i       )       ⁢           ⁢     δ   ⁡     (       n   -&gt;     -       n   -&gt;     i       )       ⁢           ⁢     δ   ⁡     (     σ   -     σ   i       )       ⁢           ⁢     δ   ⁡     (     E   -     E   i       )             ,           (   62   )             
 
the (linear) backprojection operation in Step  1  can be performed on each individual event and accumulated into the functions b q ({right arrow over (x)}) as 
                 b   q     ⁡     (     x   -&gt;     )       =       ∑     i   =   1     N     ⁢         BP     q     [         δ   ⁡     (       z   -&gt;     -       z   -&gt;     i       )       ⁢           ⁢     δ   ⁡     (       n   -&gt;     -       n   -&gt;     i       )       ⁢           ⁢     δ   ⁡     (     σ   -     σ   i       )       ⁢           ⁢     δ   ⁡     (     E   -     E   i       )           W   ⁡     (       σ   i     ,     E   i       )         ]               (   63   )             
 
   Thus, one jumps directly from the list mode data to the backprojection function b q ({right arrow over (x)}). The backprojection operations {circumflex over (B)}{circumflex over (P)} q  require much more computation than binning the data into voxels of a histogram. Consequently, if many events fall into the same voxel of g, the accumulation of events into a histogram would reduce the computational burden. Assuming that events are limited to a narrow energy window, the domain of the function g({right arrow over (z)}, {right arrow over (n)}, σ, E 0 ) is five dimensional, whereas, the domain of b q ({right arrow over (x)}) is three dimensional. This difference is important. If one samples the functions g ({right arrow over (z)}, {right arrow over (n)}, σ, E 0 ) and b q  ({right arrow over (x)}) with 128 (2 7 ) subdivisions for each dimension, then g requires 3.4×10 10  (2 35 ) voxels, whereas, b q  requires only 2.1×10 6  (2 21 ) voxels. For a typical acquisition with Compton camera  100 , the current estimate gives the number of accumulated counts as N≈10 7 −10 9 . If 10 9  counts are distributed among 3.4×10 10  voxels, it is unlikely many events will fall into the same voxel. Thus, any computational advantage of rebinning the data into a histogram g is lost. Thus, one concludes that the histogram function g ({right arrow over (z)}, {right arrow over (n)}, σ, E 0 ) should not be explicitly calculated. Instead, the list mode data should be accumulated directly into the appropriate backprojection functions b q  ({right arrow over (x)}) as prescribed by equation (63). [For p=½ or 3/2 only one function b q  is required; whereas, for p=−½ two functions b q  are required.] This strategy requires significantly less computer memory and should not slow the reconstruction. 
   A major consideration in image reconstruction is noise immunity. FBP algorithms for CT, SPECT, and PET are surprisingly immune to noisy data, but the FBP algorithms for Compton camera  100  have worrisome features that suggest noise immunity could be a problem. In the (2D) FBP algorithms for CT, SPECT, and PET, a ramp filter |{right arrow over (ν)}| appears in the reconstruction. This ramp filter enhances high spatial-frequencies and noise. This ramp filter is generally multiplied by a (low-pass) noise filter (e.g., a Butterworth filter) that suppresses the high-frequency noise. In general, the (3D) FBP algorithms for Compton cameras contain a quadratic filter |{right arrow over (ν)}| 2  rather than a linear (ramp) filter and, therefore, enhance the high-frequency noise even more. One might expect, therefore, that a stronger noise filter would be required for Compton camera reconstruction. Inherently, a noise filter causes a loss of image resolution, so that one expects either noisier images or a loss of resolution in the reconstruction. Fortunately, there is an alternative method of handling the quadratic filter that apparently avoids this problem. Noting that
 
             3   −1 [−4π 2 |{right arrow over (ν)}| 2 ]           3   =∇   2 ,  (64)
 
one can write the filtered backprojection reconstruction algorithms [equation (61)] of the previous section as 
                       C   ^       3   /   2       -   1       =         ℱ   3     -   1       ⁡     [       Φ   ⁢     (     v   -&gt;     )         4   ⁢           ⁢     π   2     ⁢           ⁢            v   -&gt;          2         ]       ⁢           ⁢     ℱ   3     ⁢           ⁢     (     -     ∇   2       )     ⁢           ⁢       BP       1   /   2                         C   ^       1   /   2       -   1       =         ℱ   3     -   1       ⁡     [       Φ   ⁡     (     v   -&gt;     )         4   ⁢           ⁢     π   2     ⁢           ⁢            v   -&gt;          2         ]       ⁢           ⁢     ℱ   3     ⁢           ⁢     (     -     ∇   2       )     ⁢           ⁢       BP       3   /   2                         C   ^         -   1     /   2       -   1       =         ℱ   3     -   1       ⁡     [       Φ   ⁡     (     v   -&gt;     )         4   ⁢           ⁢     π   2     ⁢            v   -&gt;          2         ]       ⁢           ⁢     ℱ   3     ⁢           ⁢       {         (     -     ∇   2       )     ⁢           ⁢       BP       5   /   2         +       1     R   c   2       ⁢           ⁢       BP       1   /   2           }     .                     (   65   )             

   The filter function Φ({right arrow over (ν)})/|{right arrow over (ν)}| 2  is well-behaved at large spatial-frequencies and does not amplify noise. A numerical implementation of the Laplacians in equation (65) has been used successfully for data containing fewer than a million counts (so that significant statistical noise would be expected) without any noise filter. 
   Another concern about noise immunity arises from the weighting applied in the backprojection process. The weighted backprojection operators {circumflex over (B)}{circumflex over (P)} q  add different amounts to the images depending on the distance from point of scattering in D 1  detector  105 . If q&lt;2, the backprojection operator amplifies the signal near the point of detection in detector (D 1 )  105  and suppresses the signal far away. On the other hand, if q&gt;2, the backprojection operator suppresses the signal near the point of detection and amplifies it far away. The latter case is bothersome because experience with attenuation correction algorithms implies that amplification of signals far from the point of detection amplifies noise. Only a small region is close to the point of scatter (in D 1 ), whereas most of the universe is far from this point. Thus, amplification of the signal far from the point of detection implies amplification almost everywhere. The algorithm must subtly cancel such large signals in the far field for successful reconstruction. Noise disrupts this subtle cancellation and produces significant artifacts in the reconstruction. Thus, the backprojection operator {circumflex over (B)}{circumflex over (P)}  5/2  that appears in the inversion for p=−½(Ĉ −1/2   −1 ) causes concern. This potential problem is best evaluated with computer studies that simulate statistical noise. 
   Because of these concerns, preliminary computer simulations were run to test whether the algorithms amplified noise. These initial simulations were performed for the spherical detector geometry ( 5 ) with the D 1  detector radius of 6.3 cm. The reconstruction volume was (128×128×128) voxels with 1 mm on each voxel side. Inside this volume, we placed a helix source phantom 6 cm long (z axis) and having diameter 1 cm. This phanton can be visualized as a simple helical wire defined by the affine parameter t:
 
 z= 6 t− 3
 
 x= 0.5 cos(4 πt )
 
 y= 0.5 sin(4 πt )  (66)
         where 0&lt;t&lt;1.       

   The study involved the simulation of 1 million detected scatter events. Each event required the tracking of a photon emitted in a random direction from a point on the helix, then propagated to the D 1  detector  105 , scattered at an angle sampled from the function W, and then propagated to D 2  detector  109  where the detector was assumed thick enough to assure absorption. Thus, the simulation produced data appropriate for the physical (p=0) model of the Compton camera. The scattering was assumed isotropic (W=1, rather than Klein-Nishina) and no Doppler broadening effects were included. (Later simulations using the Klein-Nishina cross-section provided virtually identical results.) The spatial and energy resolution of the detectors were assumed perfect, so that the only sources of error were the statistical noise and the index p of the reconstruction algorithm. Reconstructions were performed for the filtered backprojection algorithms described by p= 1/2, 3/2, and −½. No noise filters were used in the reconstructions. The algorithms all produced comparable results—with the conclusion that p=½ gave the marginally best results, 3/2 next best, and −½ the least successful results. The observed resolutions for all these FBP reconstructions were &lt;1.5 mm—essentially limited by the voxel size (1 mm) of the simulation.  FIGS. 7 and 8  illustrate the results of the FBP reconstruction for p=½. 
     FIGS. 7A and 7B  show the reconstruction results on a slice through the y=0 plane.  FIG. 7A  shows the results of FBP with the appropriate (Laplacian) filter; the four points where the helix intersects the y=0 plane are clearly visible.  FIG. 7B  shows the result of backprojection without the filter; the blurred image is comparable to the results of Rohe and Valentine. Visualization of the 3D distribution reconstructed from the helix source is difficult. 
     FIG. 8  shows a sequence of projections of the reconstructed distribution from various angles as the imaging volume is rotated 90 degrees around the x axis. The helix structure is clearly visible. The “bead” structure observed on the helix is a partial volume effect caused by superimposing the finite voxel size (1 mm) on the infinitely thin helical wire. These images suggest that neither the statistical noise nor the artificial index p pose significant problems for the FBP algorithms proposed in an embodiment of the invention. However, the unphysical assumptions (perfect spatial and energy resolution) leave major questions about the application of the algorithms to realistic data. More extensive simulations are required and currently underway. 
   In an embodiment of the invention, an integral equation (11) was introduced that describes the imaging process in Compton cameras. That integral equation was generalized with an index p that characterizes the response of the camera to sources at different distances from the primary detector. The case p=0 (corresponding to the radiation flux decreasing with distance as r −2 ) is an accurate model of Compton camera performance. In an embodiment of the invention, analytic inversions were found for the cases p=±½ and 3/2 for certain detector geometries. These analytic inversions take the form of FBP algorithms. If the inverse Ĉ 0   −1  were known, the analytic problem would be solved. Because Ĉ 0   −1  is not known, one must hope that Ĉ 1/2   −1 , Ĉ −1/2   −1 , or Ĉ 3/2   −1  will serve as an acceptable approximation of Ĉ 0   −1 . (Attempts to find a FBP algorithm for Ĉ 0   −1  are still under study. However, the most promising result yields an infinite series of different weighted backprojections and filters. Aside from serious questions of convergence, such expressions for Ĉ 0   − 1 are completely unsuitable for practical applications.) If one uses Ĉ p   −1  rather than Ĉ 0   −1  for image reconstruction, the resulting algorithm will produce systematic errors. These errors can be studied either analytically or experimentally with computer simulations. Analytically, one examines the difference
 
 Ê   p   =Ĉ   p   −1   Ĉ   0   −{circumflex over (l)}   (67)
 
where Ê p f gives the error in the reconstruction of a source distribution f. The analytic evaluation of Ê p  is not very difficult—the calculation is basically the same as previously discussed. The calculations of H p,q , M q,p  and D recur and can be evaluated in the same manner. A (chi-squared) norm of the error operator (i.e., ∥Ê p ∥ 2   2 ) can be found, but such norms are notoriously bad indicators of performance in image reconstruction. Experimental computer simulations and observer performance studies are required for evaluation of these algorithms. Such studies are now underway.
 
   In addition to the actual FBP algorithms, the mathematical analysis provides significant insights into the geometrical design of Compton cameras. Typically, the primary detectors of corresponding Compton cameras are small solid-state detectors. These small detectors may be sufficient for proving the concept because the small detector can be moved around the patient (or experimental phantom) in simulation of an extended primary (D 1 ) detector. However, most researchers then discuss moving the camera around the patient in a circular orbit. The D 1  geometry resulting from such an orbital scan is a thin ring around the patient. Unfortunately, a thin ring geometry for D 1  does not produce sufficient data for the FBP algorithms in an embodiment of the invention. In an embodiment of the invention, a simple linear scan down the axis of the patient body does provide sufficient information. This counterintuitive result requires careful explanation. 
   The important idea of FBP algorithms is that a 3D delta function in the spatial frequency domain permits the identification of the Fourier transform of a backprojected image with the Fourier transform of the source distribution. In the preceding derivation two functions M and D were evaluated. The function M produced a 2D delta function by careful selection of the backprojection weighting; whereas, the function D produced a 1D delta function from the detector geometry. Without both delta functions, the FBP strategy generally fails (spherical geometry is an exception to this rule). For the thin ring geometry, no delta function appears in the function D (see Table 1, Geometry  3 ), so the FBP strategy fails. Nonetheless, one might attempt analytic inversion for the thin ring geometry. After all, the 2D delta function from the function M yields the 1D convolution integral of equation (36) that should be invertible by Fourier techniques. This approach requires two additional 3D Fourier transforms and an additional filter function. Unfortunately, the additional filter generally entails dividing by a (Bessel) function that has zeros. The divergences caused by these zeros can only be removed by orbiting the camera at two (or three) different radii and taking appropriate combinations of filters. One concludes that such a reconstruction may be possible with two or three rings and with twice the computer execution time as the FBP algorithms. In contrast, a linear scan along the z axis yields a function D with a delta function, so that a FBP algorithm follows immediately. 
   A Compton camera with simple linear D 1  geometry provides similar data. Of course, two linear scans along the axis would be better than one, and a cylindrical detector would be better yet. However, the preceding analytic analysis demonstrates conclusively that a simple ring detector (or equivalently a small primary detector orbiting the patient) may not be a good configuration. The lesson is simple: the geometry of the primary D 1  detector (or equivalently, the surface scanned by the D 1  detector) is important in the reconstruction. A naive choice of this geometry can doom the best detector electronics and engineering. (On the other hand, the geometry of the secondary detector D 2  is relatively unimportant. If the D 2  detector samples {right arrow over (n)} on a half sphere surrounding each point on D 1 , then sufficient information will be available for FBP.) 
     FIG. 9  shows an apparatus  900  for reconstructing a three-dimensional image representing a source distribution in accordance with an embodiment of the invention. Apparatus  900  comprises Compton camera  100 , processor  901 , input module  903 , output module  905 , and positioning module  907 . Compton camera  100  interfaces with processor  901  and collects observed data. Processor  901  processes the observed data to reconstruct a three-dimensional image of a radiopharmaceutical source distribution in accordance with process  600  as shown in  FIG. 6 . 
   The reconstructed image may be displayed to a clinician on output module  905 . In addition, the clinician may provide processing parameters (e.g., index p) to processor  901  when reconstructing the three-dimensional image. Also, the clinician may instruct Compton camera  100  to be repositioned through positioning module  907 . Alternatively, repositioning module  907  may automatically reposition Compton camera  100  such as by linearly scanning along a patient&#39;s axis in order to provide a three-dimensional image. 
   Appendix A: Physical Justification of the Compton Camera Equation (10) 
   The integral equation (10) represents a model for Compton camera imaging. In this appendix, the physical assumptions underlying that equation are elucidated. The original data from a coincident observation in the Compton camera consists of four parameters ({right arrow over (z)}, E 1 ; {right arrow over (w)}, E 2 ) where {right arrow over (z)}εD 1  and {right arrow over (w)}εD 2 . The original data is converted into the canonical coordinates in an embodiment of the invention
 
({right arrow over (z)}, E 1 ; {right arrow over (w)}, E 2 )           ({right arrow over (z)}, {right arrow over (n)}, σ, E)
 
by the algebraic relations in equation (8): 
                 n   -&gt;     =       (       w   -&gt;     -     z   -&gt;       )              w   -&gt;     -     z   -&gt;                ,     E   =       E   1     +     E   2         ,       and   ⁢           ⁢   σ     =         (       E   2   2     +       E   1     ⁢     E   2       -       m     e   ⁢               ⁢     E   1         )         (       E   1     +     E   2       )     ⁢     E   2         .               (   8   )             

   Equation (10) can be best understood by examination of the special case of a point source of radiation. If one assumes that
 
 f ( {right arrow over (x)} )=δ 3 ( {right arrow over (x)}−{right arrow over (x)}   0 )  (A.1)
 
and that the original energy of emission is E, then, for p=0, equation (10) gives 
                 g   ^     ⁡     (       z   -&gt;     ,     n   -&gt;     ,   σ   ,   E     )       =         W   ⁡     (     σ   ,   E     )         2   ⁢           ⁢   π   ⁢           ⁢              z   -&gt;     -       x   -&gt;     0            2         ⁢           ⁢     δ   ⁡     [           n   -&gt;     ·     (       z   -&gt;     -       x   -&gt;     0       )                z   -&gt;     -       x   -&gt;     0              -   σ     ]                 (     A   ⁢   .2     )             
 
   This expression for ĝ requires physical justification. Because ĝ describes the density of observed events in the six dimensional space ({right arrow over (z)}, {right arrow over (n)}, σ, E), one expects inverse dimensions for each of these parameters. Because the vector {right arrow over (z)} is located on a 2D surface, one expects cm −2  for that parameter. The vector {right arrow over (n)} occurs on a 2-sphere and, therefore, its density is given in steradians −1  (which is technically dimensionless). Likewise, the cosine σ is dimensionless, whereas, the energy is assumed fixed (integrated over a narrow energy window). Thus, the function ĝ has dimensions cm −2  and gives the density of scatter events detected in the primary detector D 1  at point {right arrow over (z)} that are subsequently observed in the secondary detector with parameters {right arrow over (n)}, σ, and E. [In equation (A.2) the term |{right arrow over (z)}−{right arrow over (x)} 0 | −2  gives units cm −2 , W (σ, E) gives the units steradians −1 , and finally, the delta function gives the units cosine −1 .] Heuristically, one expects that ĝ will be proportional to a product of five terms: 
               g   ^     ∝     Φ   ⁢           ⁢     P   1     ⁢           ⁢     P   2     ⁢           ⁢     P   3     ⁢           ⁢     P   4               (     A   ⁢   .3     )             where                           Φ   =     {           γ   ⁢     -     ⁢   ray   ⁢           ⁢   flux   ⁢           ⁢   at   ⁢           ⁢     z   -&gt;     ⁢           ⁢   in   ⁢           ⁢   D1               originating   ⁢           ⁢   from   ⁢           ⁢   a               point   ⁢           ⁢   source   ⁢           ⁢   at   ⁢           ⁢       x   -&gt;     0             }       ,           (     A   ⁢   .4     )                   P   1     =     (           Probability   ⁢           ⁢   of   ⁢           ⁢   an   ⁢           ⁢   incident   ⁢           ⁢   γ   ⁢           ⁢   ray               undergoing   ⁢           ⁢   Compton   ⁢           ⁢   scattering               in   ⁢           ⁢   the   ⁢           ⁢   detector   ⁢           ⁢   D1           )       ,           (     A   ⁢   .5     )                   P   2     =     (           Probability   ⁢           ⁢   of   ⁢           ⁢   the   ⁢           ⁢   γ   ⁢           ⁢   ray               being   ⁢           ⁢   scattered   ⁢           ⁢   in   ⁢           ⁢   the   ⁢           ⁢   direction   ⁢           ⁢     n   -&gt;             )       ,           (     A   ⁢   .6     )                   P   3     =     (           Probability   ⁢           ⁢   of   ⁢           ⁢   the   ⁢           ⁢   scattered   ⁢           ⁢   γ   ⁢           ⁢   ray               propagating   ⁢           ⁢   from   ⁢           ⁢   D1   ⁢           ⁢   to   ⁢           ⁢   D2           )       ,           (     A   ⁢   .7     )             and                           P   4     =       (           Probability   ⁢           ⁢   of   ⁢           ⁢   the   ⁢           ⁢   scattered   ⁢           ⁢   γ   ⁢           ⁢   ray               being   ⁢           ⁢   detected   ⁢           ⁢   in   ⁢           ⁢   D2           )     .             (     A   ⁢   .8     )             
 
     FIG. 10  shows a detector D 1  flux in accordance with an embodiment of the invention. As illustrated in  FIG. 10 , the flux Φ is proportional to the expression 
             Φ   ∝       cos   ⁢           ⁢   φ       4   ⁢           ⁢   π   ⁢           ⁢              z   -&gt;     -       x   -&gt;     0            2                 (     A   ⁢   .9     )               
where cos φ is the cosine of the angle of the incident radiation with the normal to the D 1  detector surface. In particular, from  FIG. 10 , one finds that 
               cos   ⁢           ⁢   φ     =       [       m   -&gt;     ·     (       z   -&gt;     -       x   -&gt;     0       )       ]              z   -&gt;     -       x   -&gt;     0                      (     A   ⁢   .10     )               
where {right arrow over (m)} is a unit normal to the D 1  surface. As illustrated in  FIG. 11 , the probability (P 1 ) of an incident gamma-ray interacting by Compton scattering in the D 1  detector is given by 
               P   1     =       {       σ   Compton       σ   Total       }     ⁢           [     1   -     exp   ⁡     (     -       μ   ⁢           ⁢   T       cos   ⁢           ⁢   φ         )         ]             (     A   ⁢   .11     )               
where μ is the attenuation coefficient associated with gamma-rays of energy E in the D 1  material, T is the thickness of the detector, and σ Compton  and σ Total  are the cross sections for Compton and all interactions, respectively. In equation (A.11), the first expression in curly brackets { } is the probability of a Compton interaction compared with all other interactions; whereas, the second expression in square brackets [ ] gives the probability of the gamma-ray interacting (by any type of interaction) within the detector. The product of these two expressions gives the probability of a Compton interaction within detector D 1 . The probability (P 2 ) of Compton scattering in the direction {right arrow over (n)} is proportional to the differential scattering cross section [as given by W in equation (7)], 
                 P   2     ∝       W   ⁡     (     σ   ,   E     )       ⁢           ⁢     δ   ⁡     [           n   -&gt;     ·     (       z   -&gt;     -       x   -&gt;     0       )                z   -&gt;     -       x   -&gt;     0              -   σ     ]           ,           (     A   ⁢   .12     )               
where the delta function assures that the vector {right arrow over (n)} lies on the cone with the appropriate cosine σ.
 
   The remaining probabilities, P 3  and P 4 , are both assumed constant. In reality, the probability P 3  of the scattered photon propagating from a point in detector D 1  to a point in detector D 2  is extremely complicated. Basically, the only phenomenon inhibiting the propagation of the scattered photon from detector D 1  to detector D 2  is a second interaction within the primary (D 1 ) detector. Whether the scattered photon escapes detector D 1  depends on the thickness of material it must transverse to escape the detector and the energy of the scattered photon. Consequently, the probability of escape depends on the depth of interaction within the detector crystal and the direction of scatter. This probability can best be studied with Monte Carlo simulations. However, for the purposes of the embodiment of the invention (namely, development of reconstruction algorithms), the probability P 3  will be assumed constant. On the other hand, the probability P 4  of detection of a scattered photon in detector D 2  can be accurately estimated as P 4 ≈1. The secondary detector D 2  is generally designed with thick crystals so that all incident photons are absorbed. Thus, the probability of absorption P 4  can be assumed constant without loss of generality. 
   Based on the heuristic model of previous paragraphs, the function ĝ satisfies the proportionality: 
                     g   ^     ∝       ⁢           P   3     ⁢           ⁢     P   4     ⁢           ⁢     σ   Compton         2   ⁢           ⁢     σ   Total         ⁢           ⁢   cos   ⁢           ⁢     φ   ⁢           [     1   -     exp   ⁡     (     -       μ   ⁢           ⁢   T       cos   ⁢           ⁢   φ         )         ]                       ⁢       {         W   ⁡     (     σ   ,   E     )         2   ⁢           ⁢   π   ⁢           ⁢              z   -&gt;     -       x   -&gt;     0            2         ⁢           ⁢     δ   ⁡     [           n   -&gt;     ·     (       z   -&gt;     -       x   -&gt;     0       )                z   -&gt;     -       x   -&gt;     0              -   σ     ]         }     .                   (     A   ⁢   .13     )             
 
   In equation (A.13), the expression in curly brackets { } is the desired result [equation (A.2)]. Of the remaining terms, all are constant except those containing cos φ. Offhand, this cos φ dependence suggests that equation (A.2) does not accurately model the behavior of Compton cameras. However, Compton cameras are generally designed with very thin D 1  detectors (so that the scattered photons are not immediately reabsorbed in D 1 ). As a result, the expression μT in the exponential of equation (A.13) is generally small (μT&lt;&lt;1), so that 
               [     1   -     exp   ⁡     (     -       μ   ⁢           ⁢   T       cos   ⁢           ⁢   φ         )         ]     ≈         μ   ⁢           ⁢   T       cos   ⁢           ⁢   φ       -         (     μ   ⁢           ⁢   T     )     2       2   ⁢           ⁢     cos   2     ⁢           ⁢   φ                 (     A   ⁢   .14     )             or                           cos   ⁢           ⁢     φ   ⁢           [     1   -     exp   ⁡     (     -       μ   ⁢           ⁢   T       cos   ⁢           ⁢   φ         )         ]       ≈     μ   ⁢           ⁢     T   ⁡     [     1   -       μ   ⁢           ⁢   T       2   ⁢           ⁢   cos   ⁢           ⁢   φ         ]                 (     A   ⁢   .15     )             
 
   If cos φ&gt;&gt;μT, one can approximate 
               cos   ⁢           ⁢     φ   ⁢           [     1   -     exp   ⁡     (     -       μ   ⁢           ⁢   T       cos   ⁢           ⁢   φ         )         ]       ≈     μ   ⁢           ⁢     T   .               (     A   ⁢   .16     )             
 
   Thus, the cos φ dependence is only important for extremely obtuse angles (cos φ≈0). For all other angles, one can write 
               g   ^     ∝           P   3     ⁢           ⁢     P   4     ⁢           ⁢     σ     Compton   ⁢               ⁢           ⁢   μ   ⁢           ⁢   T       2   ⁢           ⁢     σ   Total         ⁢       {         W   ⁡     (     σ   ,   E     )         2   ⁢           ⁢   π   ⁢           ⁢              z   -&gt;     -       x   -&gt;     0            2         ⁢           ⁢     δ   ⁡     [           n   -&gt;     ·     (       z   -&gt;     -       x   -&gt;     0       )                z   -&gt;     -       x   -&gt;     0              -   σ     ]         }     .               (     A   ⁢   .17     )             
 
   The terms containing cos φ in equation (A. 13) are approximately constant because the decrease in flux (cos φ) is offset by an increase in the interactions within the D 1  detector [1−exp(−μT/cos φ)] due to the increased path length (as shown in  FIG. 11 ). If one drops the constant coefficients preceding the curly brackets { } in equation (A.17), the resulting expression is just that of equation (A.2). 
   If the above approximation fails, the terms containing cos p will reduce the observed flux ĝ—especially for obtuse angles (cos φ≈0). For a fixed source position, the angle will be smallest for the D 1  pixels closest to the source, whereas, the angle will be most obtuse for pixels farthest from the source. This observation implies that, if the approximation breaks down, the reduction in flux will affect distant pixels more than near pixels. Consequently, one would expect that p&lt;0 would provide a better approximation than p&gt;0. 
   Other researchers have used alternative expressions for Compton camera response. Legitimate arguments can be made for different approximations. The relative merits of these various formulations are unclear at this time. For purposes of comparison, some of these alternative formulations will be expressed in terms of the Compton imaging operators (Ĉ p ) in an embodiment of the invention. One may denote the response of the Compton camera by the function “q” (as referenced in a paper authored by R. Basko, et al., “Application of spherical harmonics to image reconstruction for the Compton camera,”  Physics in Medicine and Biology,  43, 887–894, 1998), which (using the notation of  FIG. 1 ) can be expressed as 
                     q   ⁡     (     P1   ,   P2   ,   θ     )       ∝       ∫     [       Possible     Cone   ⁢           ⁢   of         Emission   Points       ]               ⁢           ⁢       ⅆ   A     ⁢           ⁢     f   ⁡     (     x   -&gt;     )                               ⁢     ∝         1   -     σ   2         ⁢           ⁢           C   ^     1     ⁡     [   f   ]       .                       (     A   ⁢   .18     )             
 
   One can verify the second proportionality in equation (A.18) by assuming that the function f is given by 
               f   ⁡     (     x   -&gt;     )       =     {           f   0             if   ⁢           ⁢   the   ⁢           ⁢   distance   ⁢           ⁢   between   ⁢           ⁢     x   -&gt;     ⁢           ⁢   and   ⁢           ⁢   P1     ≤   R             0           if   ⁢           ⁢   the   ⁢           ⁢   distance   ⁢           ⁢   between   ⁢           ⁢     x   -&gt;     ⁢           ⁢   and   ⁢           ⁢   P1     &gt;   R                     (     A   ⁢   .19     )             
 
for arbitrary R. In that case, one finds that
 
 q ( P   1 ,  P   2 , θ)∝π R   2   f   0  sin θ=π R   2   f   0 √{square root over (1−σ 2 )},   (A.20)
 
whereas, 
                   C   ^     p     ⁡     [   f   ]       ∝         f   0     ⁢           ⁢     R     p   +   1             (     p   +   1     )     ⁢           ⁢     R   c   p                 (     A   ⁢   .21     )             
 
so that p=1 yields the appropriate weighting for the function q. One may use the same mathematical expression for the function “λ” (as referenced in a paper authored by Cree, et al., “Towards Direct Reconstruction from a Gamma Camera Based on Compton Scattering,”  IEEE Transactions on Medical Imaging , MI-13, 398–407, 1994) and, therefore, are also dealing with the p=1 problem. In both cases the sine term (√{square root over (1−σ 2 )}) has no significant affect on the analysis because it can be absorbed in the definition of g in the same manner as W in equation (11). Basically, this term arises because they parameterize the scatter angle in terms of θ rather than σ=cos θ[dσ=−sin θdθ].
 
   One may introduce an imaging kernel “k” (as referenced in a paper authored by Tomitani, et al., “Image reconstruction from limited angle Compton camera data,”  Phys. Med. Biol.,  47, 2129–2145, 2002) that can be written in the notation of this paper as 
                 k   ⁡     (       x   -&gt;     ,       z   -&gt;     ;   θ       )       =       1     2   ⁢           ⁢   π       ⁢           ⁢     δ   ⁡     [           n   -&gt;     ·     (       z   -&gt;     -     x   -&gt;       )                z   -&gt;     -     x   -&gt;              -   σ     ]           ⁢     
     ⁢       which   ⁢           ⁢   implies   ⁢           ⁢   that   ⁢           ⁢   p     =   2.             (     A   ⁢   .22     )             
 
Appendix B: Evaluation of the Function M for General Indices (q,p)
 
   According to the notation of equation (26), the regularized form of M given in equation (25) can be written in spherical coordinates as 
                       M     q   ,   p       ⁡     (       v   -&gt;     ,     ω   -&gt;       )       =       ⁢       lim     ɛ   -&gt;   0       ⁢       ∫   0   ∞     ⁢       ⅆ   x     ⁢           ⁢       ∫   0   ∞     ⁢       ⅆ   s     ⁢           ⁢         s   q     ⁢           ⁢     x   p         2   ⁢           ⁢   π   ⁢           ⁢     R   c     1   +   q   +   p                                     ⁢       exp   ⁡     [       -   ɛ     ⁢       (       x     +     s       )         R   c           ]       ⁢           ⁢     ∫     ∫       ⅆ   2     ⁢       Ω   -&gt;     x                             ⁢     ∫     ∫         ⅆ   2     ⁢       Ω   -&gt;     s       ⁢       exp   [       2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢   x   ⁢           ⁢       v   -&gt;     ·       Ω   -&gt;     x         -     2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢   s   ⁢           ⁢       ω   -&gt;     ·       Ω   -&gt;     s           ]                Ω   -&gt;     x     -       Ω   -&gt;     s                                  (     B   ⁢   .1     )             
 
   The |{right arrow over (Ω)} x −{right arrow over (Ω)} s | term can found from the generating function of the Legendre polynomials [60] as 
                          Ω   -&gt;     x     -       Ω   -&gt;     s              -   1       =         [     2   -     2   ⁢           ⁢         Ω   -&gt;     x     ·       Ω   -&gt;     s           ]       -     1   2         =       ∑     n   =   0     ∞     ⁢         P   n     ⁡     (         Ω   -&gt;     x     ·       Ω   -&gt;     s       )       .                 (     B   ⁢   .2     )             
 
   Furthermore, from the addition theorem of spherical harmonics, one finds that 
                 P   n     ⁡     (         Ω   -&gt;     x     ·       Ω   -&gt;     s       )       =       P   n     ⁡     [       cos   ⁢           ⁢     θ   x     ⁢           ⁢   cos   ⁢           ⁢     θ   s       +     sin   ⁢           ⁢     θ   x     ⁢           ⁢   sin   ⁢           ⁢     θ   s     ⁢           ⁢     cos   ⁡     (       φ   x     -     φ   s       )           ]               (     B   ⁢   .3     )                       ⁢     =         4   ⁢           ⁢   π       (       2   ⁢   n     +   1     )       ⁢       ∑     m   =     -   n       n     ⁢           Y   n   m     ⁡     (       θ   x     ,     φ   x       )       ⁢           [       Y   n   m     ⁡     (       θ   s     ,     φ   s       )       ]     *                                   so   ⁢           ⁢   that                                        Ω   -&gt;     x     -       Ω   -&gt;     s              -   1       =       ∑     n   =   0     ∞     ⁢         4   ⁢           ⁢   π       (       2   ⁢   n     +   1     )       ⁢           ⁢       ∑     m   =     -   n       n     ⁢           Y   n   m     ⁡     (       Ω   -&gt;     x     )       ⁢           [       Y   n   m     ⁡     (       Ω   -&gt;     s     )       ]     *                   (     B   ⁢   .4     )             
 
   The oscillatory Fourier terms in equation (B.1) can be written as 
               exp   ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢   x   ⁢           ⁢                   ⁢   Ω     -&gt;     x     ·     v   -&gt;         )       =       ∑     k   =   0     ∞     ⁢         2   ⁢           ⁢   π   ⁢           ⁢     ⅈ   k                  v   →          ⁢           ⁢   x         ⁢       J     k   +     1   /   2         ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢          v   →          ⁢           ⁢   x     )                   (     B   ⁢   .5     )                       ⁢       ∑     h   =     -   k       k     ⁢           Y   k   h     ⁡     (       Ω   -&gt;     v     )       ⁢           [       Y   k   h     ⁡     (       Ω   -&gt;     x     )       ]     *                             and                           exp   ⁡     (       -   2     ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢   s   ⁢           ⁢         Ω   -&gt;     s     ·     ω   -&gt;         )       =       ∑     u   =   0     ∞     ⁢         2   ⁢           ⁢   π   ⁢           ⁢     ⅈ     -   u                    ω   →          ⁢           ⁢   s         ⁢           ⁢       J     u   +     1   /   2         ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢          ω   →          ⁢   s     )                   (     B   ⁢   .6     )                       ⁢       ∑     v   =     -   u       u     ⁢             Y   u   v     ⁡     (       Ω   -&gt;     s     )       ⁢           [       Y   u   v     ⁡     (       Ω   -&gt;     ω     )       ]     *     .                             
 
   The integrals over {right arrow over (Ω)} x  and {right arrow over (Ω)} s  in equation (B.1) follow immediately from the orthogonality of the spherical harmonics 
               ∫     ∫         ⅆ   2     ⁢     Ω   -&gt;       ⁢           ⁢           Y   n   m     ⁡     (     θ   ,   φ     )       ⁢           [       Y   u   v     ⁡     (     θ   ,   φ     )       ]     *           =       δ   nu     ⁢           ⁢     δ   mv               (     B   ⁢   .7     )               so   ⁢           ⁢   that                             ∫     ∫         ⅆ   2     ⁢       Ω   -&gt;     x       ⁢           ⁢     ∫     ∫         ⅆ   2     ⁢       Ω   -&gt;     s       ⁢       exp   ⁡     [       2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢   x   ⁢           ⁢       v   -&gt;     ·       Ω   -&gt;     x         -     2   ⁢           ⁢   π   ⁢           ⁢   ⅈ   ⁢           ⁢   s   ⁢           ⁢       ω   -&gt;     ·       Ω   -&gt;     s           ]                  Ω   -&gt;     x     -       Ω   -&gt;     s                          =           (     B   ⁢   .8     )                   16   ⁢           ⁢     π   3           x   ⁢           ⁢   s         ⁢     1              ω   →          ⁢           ⁢          v   →                ⁢           ⁢       ∑     n   =   0     ∞     ⁢       1     (       2   ⁢   n     +   1     )       ⁢           ⁢     J     n   +     1   /   2         ⁢           ⁢     (     2   ⁢           ⁢   π   ⁢           ⁢          ω   →          ⁢           ⁢   s     )     ⁢           ⁢       J     n   +     1   /   2         ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢          v   →          ⁢           ⁢   x     )                                     ∑     m   =     -   n       n     ⁢           Y   n   m     ⁡     (       Ω   -&gt;     v     )       ⁢           [       Y   n   m     ⁡     (       Ω   -&gt;     ω     )       ]     *                           
 
   By substitution of equation (B.8) into equation (B.1), one finds that 
                       M     q   ,   p       ⁡     (       v   →     ,     ω   →       )       =       ⁢         8   ⁢           ⁢     π   2                    v   →          ⁢           ⁢          ω   →              ⁢           ⁢     R   c     1   +   q   +   p           ⁢           ⁢       ∑     n   =   0     ∞     ⁢     1     (       2   ⁢   n     +   1     )                           ⁢     {       ∑     m   =     -   n       n     ⁢           Y   n   m     ⁡     (       Ω   -&gt;     v     )       ⁢           [       Y   n   m     ⁡     (       Ω   -&gt;     ω     )       ]     *       }                     ⁢     {       lim     ɛ   →   0       ⁢       ∫   0   ∞     ⁢       ⅆ   x     ⁢           ⁢     x     p   -     1   /   2         ⁢           ⁢       J     n   +     1   /   2         ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢          v   →          ⁢           ⁢   x     )       ⁢           ⁢     exp   ⁡     [       -   ɛ     ⁢       x         R   c           ]             }                     ⁢       {       lim     ɛ   →   0       ⁢       ∫   0   ∞     ⁢       ⅆ   s     ⁢           ⁢     s     q   -     1   /   2         ⁢           ⁢       J     n   +     1   /   2         ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢          ω   →          ⁢           ⁢   s     )       ⁢           ⁢     exp   ⁡     [       -   ɛ     ⁢       s         R   c           ]             }     .                   (     B   ⁢   .9     )             
 
   The remaining integrals over x and S differ only in the index of the power (p or q) and, therefore, only the integral over x will be evaluated. If (−n−1&lt;p&lt;1), no regularization is required; one can set ε=0 and find [61] that 
                       ∫   0   ∞     ⁢       ⅆ   x     ⁢           ⁢     x     p   -     1   /   2         ⁢           ⁢       J     n   +     1   /   2         ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢          v   →          ⁢           ⁢   x     )           =     1     2   ⁢           ⁢     π     p   +     1   /   2         ⁢           ⁢            v   →            p   +     1   /   2                                 ⁢         Γ   ⁡     (       n   2     +     p   2     +     1   2       )         Γ   ⁡     (       n   2     -     p   2     +   1     )         .                   (     B   ⁢   .10     )             
 
   Because the sum in equation (B.9) requires all terms n&gt;0 one concludes that equation (B.10) is only useful for (−1&lt;p&lt;1). However, the subsequent analysis requires p&gt;1, hence the need for the regularization parameter. The regularization used in this analysis was chosen because the integral can be evaluated in terms of well-behaved hypergeometric functions. The expression [which is valid for ε&gt;0, |{right arrow over (ν)}|&gt;0, and −1&lt;p+n] contains four terms involving successively higher powers of ε. Writing out only the lowest order term, one finds that 
                       ∫   0   ∞     ⁢       ⅆ   x     ⁢           ⁢     x     p   -     1   /   2         ⁢           ⁢       J     n   +     1   /   2         ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢          v   →          ⁢           ⁢   x     )       ⁢           ⁢     exp   ⁡     [       -   ɛ     ⁢       x         R   c           ]           =                 1     2   ⁢           ⁢     π     p   +     1   /   2         ⁢           ⁢            v   →            p   +     1   /   2             ⁢       Γ   ⁡     (       n   2     +     p   2     +     1   2       )         Γ   ⁡     (       n   2     -     p   2     +   1     )                         F   3           2         ⁡     (         p   -   n     2     ,         p   +   n   +   1     2     ;     1   4       ,     1   2     ,       3   4     ;     -       ɛ   4         (     16   ⁢           ⁢   π   ⁢           ⁢          v   →          ⁢           ⁢     R   c       )     2             )       +       o   ⁡     (   ɛ   )       .                   (     B   ⁢   .11     )             
 
   Because the limit as ε→0 of the hypergeometric function is one, the limit of equation (B.11) yields 
                       lim     ɛ   →   0       ⁢       ∫   0   ∞     ⁢       ⅆ   x     ⁢           ⁢     x     p   -     1   /   2         ⁢           ⁢       J     n   +     1   /   2         ⁡     (     2   ⁢           ⁢   π   ⁢           ⁢          v   →          ⁢           ⁢   x     )       ⁢           ⁢     exp   ⁡     [       -   ɛ     ⁢       x         R   c           ]             =                 1     2   ⁢           ⁢     π     p   +     1   /   2         ⁢           ⁢            v   →            p   +     1   /   2             ⁢         Γ   ⁡     (       n   2     +     p   2     +     1   2       )         Γ   ⁡     (       n   2     -     p   2     +   1     )         .                   (     B   ⁢   .12     )             
 
   Thus, the regularization does not effect the preliminary results of equation (B.10); however, the range of the index p is extended to p&gt;−1. The application of equation (B. 12) in equation (B.9) gives that 
                       M     q   ,   p       ⁡     (       v   →     ,     ω   →       )       =       2   ⁢           ⁢     π     1   -   p   -   q                    v   →            p   +   1       ⁢           ⁢            ω   →            q   +   1       ⁢           ⁢     R   c     1   +   q   +   p                               ⁢       ∑     n   =   0     ∞     ⁢       1     (       2   ⁢   n     +   1     )       ⁢       Γ   ⁡     (       n   2     +     p   2     +     1   2       )         Γ   ⁡     (       n   2     -     p   2     +   1     )                                 ⁢         Γ   ⁡     (       n   2     +     q   2     +     1   2       )         Γ   ⁡     (       n   2     -     q   2     +   1     )         ⁢           ⁢       ∑     m   =     -   n       n     ⁢           Y   n   m     ⁡     (       Ω   -&gt;     v     )       ⁡     [       Y   n   m     ⁡     (       Ω   -&gt;     ω     )       ]       *                       for   ⁢           ⁢   all   ⁢           ⁢   p     &gt;       -   1     ⁢           ⁢   and   ⁢           ⁢   q     &gt;     -   1.                   (     B   ⁢   .13     )             
 
Appendix C: The Selection of Special Cases of p and q
 
   The astute reader has probably wondered why the indices p and q were introduced and their role in the reconstruction algorithms. The technical answer is that the proper choice of p and q reduces the function M q,p  in equation (27) to a delta function that permits easy inversion of the integral equation (24). The difficulty in equation (27) is the summation over spherical harmonics 
                 M     q   ,   p       (       v   →     ,     ω   →       )     ∝       ∑     n   =   0     ∞     ⁢           ⁢       1     (       2   ⁢   n     +   1     )       ⁢       Γ   ⁡     (       n   2     +     p   2     +     1   2       )         Γ   ⁡     (       n   2     -     p   2     +   1     )         ⁢       Γ   ⁡     (       n   2     +     q   2     +     1   2       )         Γ   ⁡     (       n   2     -     q   2     +   1     )         ⁢       ∑     m   =     -   n       n     ⁢           ⁢             Y   n   m     (       Ω   →     v     )     [       Y   n   m     (       Ω   →     ω     )     ]     *     .                   (   C1   )             
 
   Virtually every paper concerning backprojection of Compton camera data [30–35] contains a similar summation. However, because other researchers chose values of p and q at the onset of their calculation, there is no way to adjust the coefficients in the summation. The ability to make such adjustments is important because the identity 
                 ∑     n   =   0     ∞     ⁢       ∑     m   =     -   n       n     ⁢           ⁢           Y   n   m     (       Ω   →     v     )     [       Y   n   m     (       Ω   →     ω     )     ]     *         =       δ   2     (         Ω   →     v     -       Ω   →     ω       )             (   C2   )             
 
simplifies the inversion of integral equation (24). The choice of p and q is, therefore, determined by a simple condition: the coefficient 
                 K     p   ,   q       ⁡     (   n   )       =         1     (       2   ⁢   n     +   1     )       ⁢       Γ   ⁡     (       n   2     +     p   2     +     1   2       )         Γ   ⁡     (       n   2     -     p   2     +   1     )         ⁢       Γ   ⁡     (       n   2     +     q   2     +     1   2       )         Γ   ⁡     (       n   2     -     q   2     +   1     )           =       K     q   ,   p       ⁡     (   n   )                 (   C3   )             
 
in equation (C1) must be a non-vanishing constant (that is, independent of n). The symmetry of equation (C3) permits the interchange of p and q (p           q).

   For large n, Stirling&#39;s formula for the gamma functions yields 
                 lim     n   →   ∞       ⁢           ⁢       K     p   ,   q       ⁡     (   n   )         ∝       n     p   +   q   -   2       .             (   C4   )             
 
From equation (C4), one concludes that p+q≦2 (C5)
 
where the equality (p+q=2) must hold if K p,q  is to be a non-vanishing constant.
 
   The special case p=½, q= 3/2 yields 
                 K       1   /   2     ,     3   /   2         ⁡     (   n   )       =         1     (       2   ⁢   n     +   1     )       ⁢       Γ   ⁡     (       n   2     +     4   3       )         Γ   ⁡     (       n   2     +     3   4       )         ⁢       Γ   ⁡     (       n   2     +     5   4       )         Γ   ⁡     (       n   2     +     1   4       )           =         1     (       2   ⁢   n     +   1     )       ⁢       Γ   ⁡     (       n   2     +     5   4       )         Γ   ⁡     (       n   2     +     1   4       )           =     1   4                 (   C6   )             
 
which provides the desired behavior. Because of the symmetry under exchange (p           q), the case p= 3/2, q=½ provides the same behavior. Thus, both the cases p=½, q= 3/2 and p= 3/2, q=½ yield the delta function of equation (C2).

   APPLICABILITY OF EMBODIMENTS 
   As can be appreciated by one skilled in the art, a computer system with an associated computer-readable medium containing instructions for controlling the computer system can be utilized to implement the exemplary embodiments that are disclosed herein. The computer system may include at least one computer such as a microprocessor, digital signal processor, and associated peripheral electronic circuitry. 
   While the invention has been described with respect to specific examples including presently preferred modes of carrying out the invention, those skilled in the art will appreciate that there are numerous variations and permutations of the above described systems and techniques that fall within the spirit and scope of the invention as set forth in the appended claims.