# Patent

## Title
System and method reconstructing a nuclear medicine image using deformed attenuation image

## Abstract

According to some embodiments, an emission tomography scanner may acquire emission scan data. One or more anatomical images may be generated using an anatomical imaging system, and the anatomical images may be processed to obtain an initial attenuation image. An emission image and a corrected attenuation image may be jointly reconstructed from the acquired emission scan data, the corrected attenuation image representing a deformation of the initial attenuation image. A final reconstructed emission image may then be calculated based on the reconstructed emission image and/or the corrected attenuation image. The final reconstructed emission image may then be stored in a data storage system and/or displayed on a display system.

## Background
BACKGROUND
The invention relates generally to tomographic imaging and, more particularly, to methods and systems for correcting for motion-induced attenuation mismatch for use in tomographic imaging.
Tomographic imaging has become an integral part of healthcare services. Examples of tomographic imaging include Positron Emission Tomography (“PET”) imaging, Single Photon Emission Computed Tomography (“SPECT”) imaging, X-ray Computed Tomography (“CT”) imaging, and magnetic resonance imaging (“MRI”). In CT imaging, X-rays are propagated through the body and are detected on the other side of the body. The X-rays are attenuated to different degrees depending on encountered bodily structures, resulting in an image showing the structural features of the body. MRI also creates images of the internal structures of the human body by exploiting the property of nuclear magnetic resonance. MRI provides good contrast resolution by using various optimized radio frequency pulse sequences. Note, however, that CT imaging or MRI is not particularly sensitive to biological processes and functions.
On the other hand, PET or SPECT imaging produces images of various biological processes and functions. In PET imaging, a solution including a tracer is injected into a subject or patient to be scanned. The tracer is a pharmaceutical compound including a radioisotope with a relatively short half-life, such as 18F-Fluoro-2-Deoxyglucose (“FDG”), which is a type of sugar that includes radioactive fluorine. The tracer can be adapted such that it is attracted to sites such as lesions within the subject where specific biological or biochemical processes occur. The tracer moves to and is typically taken up in one or more organs of the subject in which these biological and biochemical processes occur. For example, cancer cells may metabolize the tracer, allowing the PET scanner to create an image illuminating the cancerous region. When the radioisotope decays, it emits a positron, which travels a short distance before annihilating with an electron. The short distance, also referred to as the positron range, is typically of the order of 1 mm for FDG in common subjects. The annihilation produces two high-energy photons propagating in substantially opposite directions.
PET imaging uses a photon detector array arranged around a scanning area, usually in a ring-shaped pattern, in which the subject or at least the part of interest of the subject is arranged. When the detector array detects two photons within a short timing window, a so-called “coincidence” is recorded. The line connecting the two detectors that received the photons is called the Line Of Response (“LOR”). The reconstruction of the image is based on the premise that the decayed radioisotope is located somewhere on the LOR. The relatively short positron range may be neglected or may be compensated for in the reconstruction. Each coincidence may be recorded in a list by three entries: two entries representing the two detectors and one entry representing the time of detection. The coincidences in the list may be grouped in one or more sinograms. A sinogram is typically processed using image reconstruction algorithms to obtain volumetric medical images of the subject. However, PET imaging does not generally provide structural details as well as other types of scanners such as CT and MRI scanners.
A PET-CT scanner includes both a CT scanner and a PET scanner installed around a single patient bore. A PET-CT scanner creates a fused image including a PET image spatially registered to a CT image. PET-CT scanners provide the advantage that the functional and biological features shown by the PET scan may be precisely located with respect to the structure illuminated by the CT scan. In a typical PET-CT scan, the patient first undergoes a CT scan, and then the patient undergoes a PET scan before exiting the scanner. After the CT and PET data have been acquired, the PET-CT scanner processes the data and generates a fused PET-CT image.
Patient motion, such as motion due to respiration, may be a significant factor in degrading the quantitative integrity of PET images. Respiratory motion may result in artifacts and/or contrast dilution of lesions from motion blurring. Respiratory-gated acquisition of PET data can reduce motion blur. In a respiratory-gated acquisition, the data is partitioned during each respiratory cycle to produce independent PET images for each partition or gate. Each of these PET images may have reduced motion blur compared to the un-gated image. However, the reduction in blur may come at the expense of increased image noise because each gate has fewer counts than the un-gated image. A PET image, nearly free of respiratory motion, with an improved signal to noise ratio may be generated by an RRA (Reconstruct, Register and Average) procedure wherein the PET images reconstructed independently for each gate are registered and then averaged.
For accurate PET quantitation, it is important to correct for both photon attenuation and respiratory motion. Although the RRA approach using respiratory gated PET data may address respiratory motion blur, a respiratory phase mismatch may exist between an attenuation image that is used in a PET image reconstruction process and each respiratory gated PET data set since the attenuation image may be obtained from a CT image, which is usually taken while the subject is holding his or her breath. This mismatch may generate undesirable artifacts in reconstructed PET images for each gate or the output PET image of the RRA approach, preventing accurate PET quantitation.
It would therefore be desirable to enhance quantitative accuracy in tomographic imaging by providing both attenuation and motion correction in an accurate manner.
BRIEF DESCRIPTION
According to some embodiments, an emission tomography scanner may acquire emission scan data. One or more anatomical images may be generated using an anatomical imaging system, and the anatomical images may be processed to obtain an initial attenuation image. An emission image and a corrected attenuation image may be reconstructed from the acquired emission scan data, the corrected attenuation image representing a deformation of the initial attenuation image. A final reconstructed emission image may then be calculated based on the reconstructed emission image and/or the corrected attenuation image. The final reconstructed emission image may then be stored in a data storage system and/or displayed on a display system.
Other embodiments are associated with systems and/or computer-readable medium storing instructions to perform any of the methods described herein.



DRAWINGS
FIG. 1 is a schematic illustration of a PET-CT imaging system.
FIG. 2 illustrates an embodiment of CT system architecture.
FIG. 3 illustrates an embodiment of PET system architecture.
FIG. 4 is a flow chart of a method for correcting for a motion-induced attenuation mismatch in accordance with aspects disclosed herein.
FIG. 5 represents a data model in accordance with some embodiments.
FIG. 6 illustrates a joint estimation approach that might be implemented according to some embodiments.
FIG. 7 is a flow chart of a method utilizing a cost function and a parameterized deformation model in accordance with aspects disclosed herein.



DETAILED DESCRIPTION
Embodiments disclosed herein include an imaging method for providing motion correction in tomographic imaging. Some embodiments are associated with reconstruction of respiratory-gated PET images using a CT image for PET attenuation correction. Note that patient motion may lead to a mismatch between PET and CT images and result in artifacts in reconstructed PET images. This mismatch may, according to some embodiments, be adjusted between gated PET and CT images. The mismatch adjustment may be performed, for example, in the CT-based attenuation image to correct for motion in a region of interest. An attenuation mismatch may be modeled as a deformation of the initially given CT-based attenuation image and the deformation as well as the PET image for each gate may be jointly estimated from respiratory gated PET data. The resulting estimated deformation may yield a mismatch corrected attenuation image, which may result in a PET image free of attenuation mismatch induced artifacts. As used herein, singular forms such as “a,” “an,” and “the” include plural referents unless the context clearly dictates otherwise.
FIG. 1 illustrates an exemplary embodiment of a PET-CT scanner 100. The PET-CT scanner 100 may include a CT system 200 and a PET system 300 mounted around a bore in a housing 120. The PET-CT scanner 100 may also include a patient table 113, a table bed 114, a processing unit 150, and a control station 115. A patient table controller (not shown) may move the table bed 114 into the bore in response to commands received from the control station 115. The control station 115 may include a display and one or more input devices such as a keyboard, a mouse, or other similar input/controller device. Through the keyboard and associated input devices, the operator may control the operation of the PET-CT scanner 100 and the display of the resulting image on the display.
The processing unit 150 may include one or more processors, one or more memories, and other associated electronics for image processing. The processing unit 150 may process the data acquired by the CT system 200 and the PET system 300 under control of an operator operating the control station 115.
FIG. 2 depicts major components of the CT system 200 of the PET-CT system 100 according to an exemplary embodiment of the present invention. The components of the CT system 200 may be housed both in the housing 120 and in the processing unit 150 shown in FIG. 1. The CT system 200 conducts X-ray tomographic imaging by illuminating a subject 20 with an x-ray beam 204 substantially transverse to an axis through the subject 20. The axis may generally be centered on an object 22 of interest, such as an organ or other tissue structure. The subject 20 may be located on the table bed 114 (as shown in FIG. 1) that translates along the direction of the axis, thereby enabling illumination of a volumetric portion of the subject 20 by the x-ray beam 204.
The CT system 200 may include a source-detector assembly, which in an exemplary embodiment may comprise a gantry 212 rotatable about the axis. An x-ray source 214, such as an x-ray tube, may be mounted on the gantry 212 and may rotate with rotation of the gantry 212. The x-ray source 214, which may comprise a collimating element (not shown), may project the beam 204 of x-rays toward a detector array 216 disposed opposite the source 214 relative to the gantry 212.
The detector array 216 may include numerous individual detector elements 218. Detector elements 218 may together provide information regarding the internal structures of the subject 20, such as the object 22. In one embodiment, each detector element 218 may generate an electrical signal indicating the intensity of a portion of the x-ray beam 204 impinging thereupon.
The signals from detector elements 218 may indicate a degree of attenuation of the beam 204 as the x-rays traverse the material or substance of the subject 20. In one embodiment, the source 214 may be rotated around the subject 20 to execute a scan operation whereby the CT system 200 acquires x-ray data. In another embodiment, the gantry 212, with source 214 attached to a side portion thereof, may rotate about the axis of the subject 20 to acquire x-ray data from numerous different illumination angles or “view angles.”
The rotation operation for the source 214 may be controlled by a control/interface system 220. The control/interface system 220 may include a server computer residing in the processing unit 150 and the operator may interact with the control/interface system 220 by means of the control station 115 and/or other input devices. The control/interface system 220 may provide control for positioning of the gantry 212 relative to the subject 20, such as controlling speed of rotation about the axis and control of relative positions of the table 113 and the gantry 212. A controls section 222 may also provide control over x-ray generation (power and timing) of the source 214. The control/interface system 220 may also include a Data Acquisition System (“DAS”) 224 that samples the detector signals generated from the detector elements 218 and converts the sampled signals into digital data for further processing.
A reconstruction engine 230 may also be housed in the processing unit 150 and may receive the sampled and digitized data (sometimes referred to as “projection data”) from the DAS 224 and may perform image reconstruction to generate CT images. In one embodiment, the reconstruction engine 230 may include a separate processor 232 and/or memory 234. Various algorithms may be utilized for reconstructing a CT image from projection data comprising a plurality of projection views. Generally, the CT image may be generated in a format compatible with the Digital Imaging and Communications in Medicine “DICOM”) standard. The DICOM standard specifies the network protocol by which two DICOM-compatible systems communicate.
In one embodiment, the reconstruction engine 230 may send the reconstructed CT image to, for example, a system management computer 240, which may also reside in the processing unit 150, for storage or further processing. The computer 240 may include a CPU (a processor) 242 and/or at least one memory 244.
FIG. 3 illustrates an embodiment of the PET system 300 of the PET-CT imaging system. The PET system 300 may include detector ring assembly 311 disposed about the patient bore. The detector ring assembly 311 may include multiple detector rings that are spaced along the central axis to form a cylindrical detector ring assembly. Each detector ring of the detector ring assembly 311 may be formed of detector modules 320. Each detector module 320 may include an array (e.g., a 6×6 array) of individual detector crystals that may be formed of bismuth germanate (BGO), for example. Other various detector crystals or materials may also be provided. The detector crystals may detect gamma rays emitted from the patient and in response produce photons. In one embodiment, the array of detector crystals may be positioned in front of four PhotoMultiplier Tubes (“PMTs”). The PMTs may produce analog signals when a scintillation event occurs at one of the detector crystals, e.g., when a gamma ray emitted from the patient is received by one of the detector crystals. Other photon detectors such as solid-state photomultipliers or avalanche photodiodes may be used in another embodiment. A set of acquisition circuits 325 may be mounted within the housing 120 (shown in FIG. 1) to receive these signals and produce digital signals indicating the event coordinates (e.g., the location of the detected gamma ray) and the total energy of the gamma ray. These may be sent through a cable 326 to an event locator circuit 327. In another embodiment, each acquisition circuit 325 may also produce an Event Detection Pulse (“EDP”), which indicates the time the scintillation event took place.
The event locator circuits 327 may form part of a data acquisition processor 330, which periodically samples the signals produced by the acquisition circuits 325. The processor 330 may have an acquisition CPU 329, which controls communications on the local area network 318 and a backplane bus 331. The event locator circuits 327 may assemble the information regarding each valid event into a set of digital numbers that indicate precisely when the event took place and the position of the detector crystal, which detected the event. This event data packet may be conveyed to a coincidence detector 332, which is also part of the data acquisition processor 330.
The coincidence detector 332 may accept the event data packets from the event locator circuits 327 and may determine whether any two of them are in coincidence. In this example, coincidence may be determined by a number of factors. First, the time markers in each event data packet may be required to be within a certain time period of each other, e.g., 12.5 nanoseconds. Second, the locations indicated by the two event data packets may be required to lie on a straight line which passes through the Field Of View (“FOV”) in the patient bore. For a detailed description of the coincidence detector 332, reference is made to U.S. Pat. No. 5,241,181 entitled “Coincidence Detector For A PET Scanner,” which is hereby incorporated by reference in its entirety. Coincidence event pairs may be located and recorded as a coincidence data packet that is conveyed through a link 333 to a storage subsystem 350. In the storage subsystem 350, a sorter 334 may use a lookup table to sort the coincidence events in a 3D projection plane format. For a detailed description of the sorter 334, reference is made to U.S. Pat. No. 5,272,343 entitled “Sorter For Coincidence timing Calibration In A PET Scanner,” which is hereby incorporated by reference in its entirety.
The detected events may be stored in a dynamic histogram memory (histogrammer 335) where the events are ordered by radial positions and projection angles and other parameters. For example, in Time-of-Flight (“TOF”) PET scanners, the difference in arrival time of the two photons may also be recorded. In addition, the information on the energy of the photons may also be used. The PET data for a particular frame may be written to a raw data disk 336. TOF PET imaging enables time-difference measurement, e.g., determines the amount of time between the recording of one event by one of the detectors and the recording of the other event by the other detector. Therefore, if an event occurs at the midpoint between these two detectors, the difference in time would be zero. If the event occurs closer to one detector, there will be a delay before the other detector sees it. Thus, TOF makes it possible for a point of origination of annihilation to be more accurately predicted, which leads to more accurate imaging. Ultimately, improved event localization reduces noise in image data, resulting in higher image quality, shorter imaging times, and lower dose to the patient.
The PET system 300, as shown in FIG. 3, may include one or more additional processors 345 such as, for example, a Prospective Reconstruction Manager (“PRM”), a Compute Job Manager (“CJM”), and a PET Image Processor (“PET IP”). The processors 345 may interact with an array processor 337 in the storage subsystem 350 to process the projection plane format PET data into attenuation corrected PET images.
The PET system 300 may also include a Computed Tomography Attenuation Correction (“CTAC”) server 342. The CTAC server 342 may execute an independent process that runs in the processing unit 150. The CTAC process may receive CT image data from the CT system 200 and convert that CT image data into CTAC data. For example, the CTAC process may receive a request from the CT system and perform a bi-linear or other algorithm to convert the data from CT image units (HU) to a PET 511 keV attenuation coefficient (cm−1), which produces the CTAC correction for PET data from the CT images. Once the CT images are converted to CTAC data, the CTAC server 342 may write the CTAC data to the raw data disk 336 in the storage subsystem 350. At the same time, a record may be transmitted to the PET database 348 to create a data link (CTAC record) to the CTAC data.
The PET system 300 may also include a PET-CT image processor 410 for receiving CT images and PET images. The CT images and the PET images may be spatially registered to each other because the patient undergoes both scans while remaining in the same position on the table bed 114. Registration may be achieved by detecting and estimating patient movement. The PET-CT image processor 410 generates a fused PET-CT image using the input CT and PET images.
It should be appreciated that the arrangement depicted in FIGS. 1 through 3 is exemplary. For instance, the PET-CT scanner 100 may include different configurations or number of processors, memories, and/or other hardware, to perform various additional functions, and these components may be located at other locations such as the control station 115, or at another server or processing unit. It should also be appreciated that the PET-CT system 100 may be further configured or customized to suit various scanning needs.
FIG. 4 illustrates a flow chart of an embodiment of an imaging method 400 for correcting for a motion-induced mismatch in accordance with some embodiments of the present invention. The method 400 may provide attenuation mismatch correction in regions of interest and may be employed with a PET-CT system, such as described in reference to FIGS. 1 through 3. The flow charts described herein do not imply a fixed order to the steps, and embodiments of the present invention may be practiced in any order that is practicable. Note that any of the methods described herein may be performed by hardware, software, or any combination of these approaches. For example, a computer-readable storage medium may store thereon instructions that when executed by a machine result in performance according to any of the embodiments described herein.
At S410, emission scan data may be acquired using an emission tomography scanner. Note that acquiring emission scan data might comprise acquiring PET emission scan data using a PET scanner or acquiring TOF PET emission scan data using a TOF PET scanner. According to other embodiments, emission scan data may be acquired by acquiring SPECT emission scan data using a SPECT scanner. According to some embodiments, the acquisition of emission scan data may involve respiratory gating and acquiring one or more gates of the respiratory gated emission scan data. According to other embodiments, the acquisition of emission scan data may be associated with cardiac gating and acquiring one or more gates of the cardiac gated emission scan data. Thus, the emission scan data may be respiratory and/or cardiac gated, and the method can be applied to each gate independently. Therefore, explicit gate dependence may be ignored here without loss of generality. Note that emission scan data may be acquired in sinogram format. Also note, however, that the method can also be applied to emission scan data acquired in list-mode format.
The emission scan data acquired at S410 may be associated with a statistical model. In particular, PET emission scan data yi may be modeled as independent Poisson random variables whose means are given by

 yi=e−li(μ)si(λ)+ri  (1)

with

li(μ)=[Gμ]i  (2)

si(λ)=[Aλ]i  (3)

for i=1, . . . , Nb where i denotes the sinogram bin index, μ=[μ1, . . . , μNv]′ and λ=[λ1, . . . , λNv]′ denote column vectors representing an attenuation image and an emission image, respectively, ′ denotes the matrix or vector transpose, Gε represents the geometric forward projector, Aε represents the system matrix including geometric forward projection, normalization and detector blurring, ri is the mean contribution of randoms and scatters, Nb and Nv are the number of sinogram bins and image voxels, respectively, and [•]i denotes the ith element of a vector. Here the numbers of image voxels for attenuation and emission images are assumed to be the same without loss of generality. The si terms represent the unattenuated projection of the emission image as if there were no attenuation, and the e−li terms represent the survival probabilities for attenuation. Here, ri includes scatter estimates, which may be estimated using an initially given attenuation image or furthermore may be iteratively updating using corrected attenuation images obtained by the methods described herein. Alternatively, the PET scan data may be modeled as independent Gaussian random variables.

At S420, one or more anatomical images may be generated using an anatomical imaging system. The step may be associated with, for example, generating one or more X-ray CT images from a helical CT scan, a cone beam CT scan, a step-and-shoot CT scan, or a cine CT scan with breath-hold or free breathing using an X-ray CT system and/or generating one or more Magnetic Resonance Imaging (“MRI”) images, which may include a proton density weighted MRI image, a T1-weighted MRI image, a T2-weighted MRI image, a T2*-weighted MRI image, MRI images obtained using the Dixon's method, or an MRI image acquired using ultra-short or zero echo time radio frequency sequences, using an MRI system. According to some embodiments, the generation of the one or more anatomical images may involve respiratory gating and providing one or more gates of the respiratory gated anatomical images. For example, an X-ray CT image may be acquired using quiescent period respiratory gating. According to other embodiments, the generation of the one or more anatomical images may be associated with cardiac gating and providing one or more gates of the cardiac gated anatomical images. According to alternative embodiments, the generation of the one or more anatomical images may involve dual cardiac-respiratory gating and providing one or more gates of the dual cardiac-respiratory gated anatomical images. The one or more anatomical image generated in S420 and the emission scan in S410 may be mismatched due to patient motion such as respiratory and/or cardiac motion, or due to respiratory and/or cardiac phase mismatches.
At S430, the anatomical images may be processed to obtain an initial attenuation image. The processing of the anatomical images to obtain the initial attenuation image may involve segmenting the anatomical images and assigning a pre-determined attenuation coefficient to each segmented region, organ or tissue. The processing of the anatomical might also be associated with an application of a bilinear or trilinear transformation to an X-ray CT image to convert CT image units (HU) to a 511 keV PET attenuation coefficient (cm−1), and/or an application of atlas-based, machine-learning, and/or pattern-recognition methods to one or more MRI images.
At S440, an emission image and a corrected attenuation image may be automatically reconstructed from the acquired emission scan data. The corrected attenuation image may represent, according to some embodiments, a deformation of the initial attenuation image. At S450, a final reconstructed emission image may be calculated based on the reconstructed emission image and/or the corrected attenuation image.
Suppose μ0=[μ10, . . . , μNv0]′ denotes an initial attenuation image, which is obtained from an anatomical image such as x-ray CT or MRI images. In the initial attenuation image there may exist a respiratory and/or cardiac motion induced mismatch between the emission scan and the anatomical imaging scan. Such a mismatch may produce artifacts in reconstructed emission images.
A corrected attenuation image μ where the mismatch is corrected for may be modeled as a deformation of the initial attenuation image. To facilitate interpolation, continuous functions {tilde over (μ)}:  and {tilde over (μ)}0:  may be introduced, whose discrete-space samples are μ and μ0, respectively. That is,

μj={tilde over (μ)}(xj)  (4)

μj0={tilde over (μ)}0(xj)  (5)

where xj=(x1j, x2j, x3j) denotes the Cartesian coordinate of voxel j in the 3D Euclidean space for j=1, . . . , Nv.

The continuous-space version of the initial attenuation image {tilde over (μ)}0 may be modeled as a linear combination of basis functions vl: , that is,

{tilde over (μ)}0(x)=ΣNl=1Nvulvl(x)  (6)

for xε where ul are linear coefficients. For example, the tensor product of cubic B-splines may be used for vl, that is,





  
    
      
        
          
            
              v
              l
            
            ⁡
            
              (
              
                
                  x
                  1
                
                ,
                
                  x
                  2
                
                ,
                
                  x
                  3
                
              
              )
            
          
          =
          
            
              
                β
                3
              
              ⁡
              
                (
                
                  
                    x
                    1
                  
                  -
                  
                    x
                    1
                    l
                  
                
                )
              
            
            ⁢
            
              
                β
                3
              
              ⁡
              
                (
                
                  
                    x
                    2
                  
                  -
                  
                    x
                    2
                    l
                  
                
                )
              
            
            ⁢
            
              
                β
                3
              
              ⁡
              
                (
                
                  
                    x
                    3
                  
                  -
                  
                    x
                    3
                    l
                  
                
                )
              
            
          
        
        ⁢
        
          

        
        ⁢
        
          
            
              for
              ⁢
              
                
              
              ⁢
              l
            
            =
            1
          
          ,
          ...
          ⁢
          
            
          
          ,
          
            
              N
              v
            
            ⁢
            
              
            
            ⁢
            where
          
        
      
    
    
      
        (
        7
        )
      
    
  
  
    
      
        
          
            β
            3
          
          ⁡
          
            (
            t
            )
          
        
        =
        
          {
          
            
              
                
                  
                    
                      
                        
                          
                            
                            t
                            
                          
                          3
                        
                        2
                      
                      -
                      
                        
                          
                          t
                          
                        
                        2
                      
                      +
                      
                        2
                        3
                      
                    
                    ,
                  
                
                
                  
                    0
                    ≤
                    
                      
                      t
                      
                    
                    ≤
                    1
                  
                
              
              
                
                  
                    
                      -
                      
                        
                          
                            (
                            
                              
                                
                                t
                                
                              
                              -
                              2
                            
                            )
                          
                          3
                        
                        6
                      
                    
                    ,
                  
                
                
                  
                    1
                    ≤
                    
                      
                      t
                      
                    
                    ≤
                    2
                  
                
              
              
                
                  
                    0
                    ,
                  
                
                
                  
                    
                      
                      t
                      
                    
                    ≥
                    2
                  
                
              
            
            .
          
        
      
    
    
      
        (
        8
        )
      
    
  




Given μ0, to calculate u=[u1, . . . , uNv]′, the equation

μj0=ΣNl=1Nvulvl(xj),  (9)

for j=1, . . . , Nv, may be solved by deconvolution in the Fourier space.

The continuous-space version of the corrected attenuation image {tilde over (μ)} may be modeled as a deformation of {tilde over (μ)}0 using the following form:

{tilde over (μ)}(x)={tilde over (μ)}0(x1+{tilde over (d)}1(x),x2+{tilde over (d)}2(x),x3+{tilde over (d)}3(x))  (10)

for x=(x1, x2, x3)ε where ({tilde over (d)}1, {tilde over (d)}2, {tilde over (d)}3):  represents a deformation field or a displacement field. Note that embodiments described herein may apply to cases where the deformation is not only a re-interpolation but also takes density changes into account (e.g., by adding a multiplicative term to equation 10). The deformation field may be modeled as a linear combination of basis functions:

{tilde over (d)}q(x)=Σk=1Npθkq{tilde over (b)}k(x)  (11)

for q=1, 2, 3 where {tilde over (b)}k are basis functions, θkq are linear coefficients, and Np is the number of the basis functions. Here, for {tilde over (b)}k, the tensor product of cubic B-splines may be used as





  
    
      
        
          
            
              b
              ~
            
            k
          
          ⁡
          
            (
            
              
                x
                1
              
              ,
              
                x
                2
              
              ,
              
                x
                3
              
            
            )
          
        
        =
        
          
            
              β
              3
            
            (
            
              
                
                  x
                  1
                
                -
                
                  t
                  1
                  k
                
              
              Δ
            
            )
          
          ⁢
          
            
              β
              3
            
            (
            
              
                
                  x
                  2
                
                -
                
                  t
                  2
                  k
                
              
              Δ
            
            )
          
          ⁢
          
            
              β
              3
            
            (
            
              
                
                  x
                  3
                
                -
                
                  t
                  3
                  k
                
              
              Δ
            
            )
          
        
      
    
    
      
        (
        12
        )
      
    
  




where (t1k, t2k, t3k) represents the Cartesian coordinates of the control nodes and A is the node spacing. Alternative basis functions such as DCT (discrete cosine transform) basis functions or polynomial functions may be used for {tilde over (b)}k. Let djq be discrete-space samples of {tilde over (d)}q such that

djq={tilde over (d)}q(xj)=Σk=1Npθkq{tilde over (b)}k(xj)=[Bθq]j  (13)

for q=1, 2, 3 and j=1, . . . , Nv, where [B]jk={tilde over (b)}k(xj), θq=[θ1q, . . . , θNpq]′, and [•]jk denotes the (j, k)th element of a matrix.

Combining (4), (6), (7), (10) and (11) yields

μj=Σl=1NvulΠq=13β3(xqj+xql+Σk=1Np{tilde over (b)}k(xj)θkq).  (14)

This can be rewritten as

μ(θ)=Wβ3β3β3(θ)u  (15)

where θ=[θ1, θ2, θ3]′ and the Nv by Nv matrix W is defined as

[Wfgh]jl=ƒ(x1j−x1l+dj1)g(x2j−x2l+dj2)h(x3j−x3l+dj3)  (16)

for functions ƒ: , g:  and h:  with djq being given in (13). From (5)-(7), it can be seen that

μ0=Wβ3β3β3(0)u  (17)

where 0=[0, . . . , 0]′ε3Np. Therefore, (15) can be rewritten as

μ(θ)=Wβ3β3β3(θ)(Wβ3β3β3(0))−1μ0.  (18)

The above equation represents the corrected attenuation image μ as a deformation of the initial attenuation image μ0 where the deformation is parameterized by the deformation parameters θ.

With respect to the reconstruction of the corrected attenuation image and the emission image, the goal is to estimate, from the PET emission scan data yi, the emission image λ and the deformation parameters θ, which parameterize the corrected attenuation image μ as a deformation of the initial attenuation image μ0, as given in (18). The emission image λ and the deformation parameters θ may be calculated by minimizing a cost function Φ:

({circumflex over (λ)},{circumflex over (θ)})=arg minλ≧0,θΦ(λ,θ)  (19)

Φ(λ,θ)=L(λ,θ)+R(λ,θ)  (20)

where L is the negative Poisson log-likelihood,

L(λ,θ)=Σi=1Nb−yi log(e−li(μ(θ))si(λ)+ri)+(e−li(μ(θ))si(λ)+ri)  (21)

and R is the regularization function or the regularizing penalty function. Instead of the negative Poisson log-likelihood, least squares or weighted least squares may be alternatively used. If the log-likelihood is solely used without any regularization function, the calculation in (19) amounts to maximum likelihood estimation. If the regularization function is used, the calculation in (19) is called penalized likelihood estimation or Maximum A Posteriori (“MAP”) estimation. An example of the regularization function is the quadratic penalty, which encourages smoothness in the emission image and the deformation field:





  
    
      
        
          R
          ⁡
          
            (
            
              λ
              ,
              θ
            
            )
          
        
        =
        
          
            
              γ
              λ
            
            ⁢
            
              
                ∑
                
                  j
                  =
                  1
                
                
                  N
                  v
                
              
              ⁢
              
                
              
              ⁢
              
                
                  ∑
                  
                    k
                    ∈
                    
                      N
                      j
                      λ
                    
                  
                  
                    
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    
                      
                        w
                        jk
                        λ
                      
                      ⁡
                      
                        (
                        
                          
                            λ
                            j
                          
                          -
                          
                            λ
                            k
                          
                        
                        )
                      
                    
                    2
                  
                  2
                
              
            
          
          +
          
            
              γ
              θ
            
            ⁢
            
              
                ∑
                
                  j
                  =
                  1
                
                
                  N
                  p
                
              
              ⁢
              
                
              
              ⁢
              
                
                  ∑
                  
                    k
                    ∈
                    
                      N
                      j
                      θ
                    
                  
                  
                    
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    
                      
                        w
                        jk
                        θ
                      
                      ⁡
                      
                        (
                        
                          
                            θ
                            j
                          
                          -
                          
                            θ
                            k
                          
                        
                        )
                      
                    
                    2
                  
                  2
                
              
            
          
        
      
    
    
      
        (
        22
        )
      
    
  




where yλ and yθ are regularization parameters,  and denote a set of neighbors for voxel j and basis function j′, respectively, and wjkλ and wjkθ are pre-determined non-negative weights. Alternatively, nonquadratic penalties or other types of penalties such as Relative Difference Prior (“RDP”) may be used for the regularization function.

Any numerical optimization algorithm may be used to calculate the emission image λ and the deformation parameters θ by minimizing the cost function Φ in (19). For example, an alternating update approach may be used where the emission image λ and the deformation parameters θ are alternatively updated:

{circumflex over (λ)}n+1=arg minλ≧0Φ(λ,θn)  (23)

{circumflex over (θ)}n+1=arg minθΦ({circumflex over (λ)}n+1,θ)  (24)

for n=0, 1, . . . where n denotes the iteration index. For initialization, {circumflex over (θ)}0=0 may be used. In (23) and (24), accurate minimization is not necessary but a sufficient decrease in the cost function is enough. Iterating between (23) and (24) may be terminated if a certain criterion is satisfied, for example, a pre-determined iteration number is reached, the difference between the previous and the current iterates is smaller than a certain threshold, or the gradient of the cost function is smaller in norm than a certain threshold.

For the update of the emission image in (23), any algorithm for emission image reconstruction may be used including gradient based methods, De Pierro's modified EM, Separable Paraboloidal Surrogate (“SPS”), preconditioned conjugate gradient, and Block Sequential Regularized Expectation Maximization (“BSREM”).
For the update of the deformation parameters in (24), gradient based algorithms may be used utilizing the gradient of the cost function with respect to the deformation parameters: ∇θΦ(λ,θ)=∇θL(λ,θ)+∇θR(λ,θ). The gradient of the negative log-likelihood with respect to the deformation parameters can be calculated as

∇θL(λ,θ)=[∇θ1L(λ,θ),∇θ2L(λ,θ),∇θ3L(λ,θ)]′  (25)

with

∇θ1L(λ,θ)=B′diag{W{dot over (β)}3β3β3(θ)u}∇μL(λ,μ)  (26)

∇θ2L(λ,θ)=B′diag{Wβ3{dot over (β)}3β3(θ)u}∇μL(λ,μ)  (27)

∇θ3L(λ,θ)=B′diag{Wβ3β3{dot over (β)}3(θ)u}∇μL(λ,μ)  (28)

where diag{•} denotes a diagonal matrix,





  
    
      
        
          
            β
            .
          
          3
        
        =
        
          {
          
            
              
                
                  
                    
                      
                        (
                        
                          
                            
                              3
                              ⁢
                              
                                
                                  
                                  t
                                  
                                
                                3
                              
                            
                            -
                            
                              4
                              ⁢
                              
                                
                                t
                                
                              
                            
                          
                          2
                        
                        )
                      
                      ⁢
                      
                        sgn
                        ⁡
                        
                          (
                          t
                          )
                        
                      
                    
                    ,
                  
                
                
                  
                    0
                    ≤
                    
                      
                      t
                      
                    
                    ≤
                    1
                  
                
              
              
                
                  
                    
                      
                        -
                        
                          
                            
                              (
                              
                                
                                  
                                  t
                                  
                                
                                -
                                2
                              
                              )
                            
                            2
                          
                          2
                        
                      
                      ⁢
                      
                        sgn
                        ⁡
                        
                          (
                          t
                          )
                        
                      
                    
                    ,
                  
                
                
                  
                    1
                    ≤
                    
                      
                      t
                      
                    
                    ≤
                    2
                  
                
              
              
                
                  
                    0
                    ,
                  
                
                
                  
                    
                      
                      t
                      
                    
                    ≥
                    2
                  
                
              
            
            ,
          
        
      
    
    
      
        (
        29
        )
      
    
  




and ∇μL denotes the gradient of L with respect to μ given by





  
    
      
        
          
            ∇
            μ
          
          ⁢
          
            L
            ⁡
            
              (
              
                λ
                ,
                μ
              
              )
            
          
        
        =
        
          
            
              G
              ′
            
            ⁡
            
              [
              
                
                  
                    ⋮
                  
                
                
                  
                    
                      
                        (
                        
                          
                            
                              y
                              i
                            
                            
                              
                                y
                                _
                              
                              i
                            
                          
                          -
                          1
                        
                        )
                      
                      ⁢
                      
                        e
                        
                          -
                          
                            
                              l
                              i
                            
                            ⁡
                            
                              (
                              μ
                              )
                            
                          
                        
                      
                      ⁢
                      
                        
                          s
                          i
                        
                        ⁡
                        
                          (
                          λ
                          )
                        
                      
                    
                  
                
                
                  
                    ⋮
                  
                
              
              ]
            
          
          .
        
      
    
    
      
        (
        30
        )
      
    
  




On terminating the updates in (23) and (24) at iteration N, the emission image iterate {circumflex over (λ)}N may be used as a final reconstructed emission image, or the corrected attenuation image μ({circumflex over (θ)}N)=Wβ3β3β3({circumflex over (θ)}N)(Wβ3β3β3(0))−1μ0 may be used for reconstructing a final emission image, for example, by Ordered Subsets Expectation Maximization (“OSEM”) algorithm or penalized-likelihood through BSREM, conjugate gradient or SPS.
At S460, the final reconstructed emission image is output. The final reconstructed emission may be “output,” for example, by storing the final reconstructed emission image in a data storage system and/or displaying the final reconstructed emission image on a display system.
Thus, according to some embodiments, the respiratory phase mismatch between a CT-based attenuation image and a respiratory gated PET emission data set may be corrected for by jointly estimating the mismatch and PET emission images. Specifically, an attenuation image for each respiratory gate may be modeled as a deformation of a given CT-based attenuation image and the deformation and a PET emission image are jointly estimated from each gate of the respiratory gated PET emission data. The output of this method is the corrected attenuation image and the emission image for each respiratory gate, which best explain the gated emission data. The reconstructed emission images may be used to generate a respiratory motion free image by the RRA procedure.
Alternatively, the corrected attenuation image may be subsequently used for each gate for reconstructing a PET emission image with different reconstruction parameters including smoothing parameters and iteration numbers. In addition, some embodiments may help solve a problem introduced when different slices of a given CT-based attenuation image correspond to different respiratory phases, by deforming the given attenuation image so that the measured sinogram data are best explained. Some embodiments may also be applied to a more complicated procedure to generate a motion free emission image where the emission images and their deformation parameters are simultaneously estimated from the entire gates of the gated data set. Note that a joint estimation technique with a CT-based attenuation image may be used in order to address a respiratory phase mismatch between the CT-based attenuation image and respiratory gated PET data.
According to some embodiments, when a CT-based attenuation image and respiratory gated PET emission data are given, a respiratory phase mismatch between the CT-based attenuation image and the gated PET data may be reduced. The CT-based attenuation image, which may be obtained from a CT image usually taken with breath-holding, contains a respiratory motion induced mismatch. According to some embodiments, a respiratory phase matched attenuation image for each gate is modeled as a deformation of the given CT-based attenuation image. For example, FIG. 5 represents a data model 500 in accordance with some embodiments. In particular, in a patient's body 510 a first state of the lungs 520 may be associated with a given CT-based attenuation image μ0. A deformation element 530 may be associated with an unknown deformation parameter θg for gate g. Moreover, a patient's second state of the lungs 522 may be associated with an unknown corrected attenuation image μg for gate g. A data formation element 540 may be associated with an unknown emission image λg for gate g and generate gated emission data yg for gate g. The deformation can be modeled using basis functions such as B-splines, whose coefficients are deformation parameters to be estimated. For each gate, measured PET emission data may be modeled as an attenuated forward projection (analogous to the attenuated Radon transform) of an unknown gate-dependent PET emission image where attenuation factors involve a deformation of the CT-based attenuation image. In other words, the gated PET data set for each gate g may represent a function of the unknown gate-dependent PET emission image and also the deformation parameters, which represent a deformation of the given CT-based attenuation image.
The gate-dependent PET emission image and deformation parameters can be jointly estimated from the corresponding gated PET data. For example, FIG. 6 illustrates a joint estimation approach 600 that might be implemented according to some embodiments. In particular, gated measured emission data yg for gate g may be provided to a joint estimation element 630 along with given CT-based attenuation image μ0 for a patient's body 610 and the first state of the lungs 620. The joint estimation element 630 may output a reconstructed emission image {circumflex over (λ)}g for gate g along with estimated deformation parameters {circumflex over (θ)}g for gate g. Moreover, the estimated deformation parameters {circumflex over (θ)}g for gate g may be used as input to a deformation element 640 which may generate an estimated corrected attenuation image {circumflex over (μ)}g for gate g. Thus, from the estimated deformation parameters, the corresponding attenuation image for each gate can be obtained by warping the CT-based attenuation image accordingly. For each gate g, the deformation parameters and the emission image can be jointly estimated by a Penalized Likelihood (“PL”) estimation that calculates an emission image {circumflex over (λ)}g and deformation parameters {circumflex over (θ)}g for gate g by maximizing a PL objective function as previously described.
According to some embodiments, the deformation parameters and the respiratory motion free emission image are jointly estimated in connection with an objective function:




  Φ
  =
  
    
      
        ∑
        g
        
          
        
      
      ⁢
      
        
      
      ⁢
      
        
          L
          
            (
            g
            )
          
        
        ⁡
        
          (
          
            
              λ
              
                (
                g
                )
              
            
            ,
            
              μ
              
                (
                g
                )
              
            
          
          )
        
      
    
    +
    
      
        D
        1
      
      ⁡
      
        (
        
          
            λ
            
              (
              g
              )
            
          
          ,
          
            
              W
              ⁡
              
                (
                
                  θ
                  
                    (
                    g
                    )
                  
                
                )
              
            
            ⁢
            λ
          
        
        )
      
    
    +
    
      
        D
        2
      
      ⁡
      
        (
        
          
            μ
            
              (
              g
              )
            
          
          ,
          
            
              W
              ⁡
              
                (
                
                  θ
                  
                    (
                    g
                    )
                  
                
                )
              
            
            ⁢
            
              μ
              0
            
          
        
        )
      
    
    +
    R
  




where λ(g) is an emission image for gate g, which is modeled as a deformation of a respiratory motion free emission image λ, W(θ(g)) is a warping matrix representing the deformation for gate g, which is parameterized by deformation parameter θ(g), μ(g) is an attenuation image for gate g, which is modeled as a deformation of the given CT-based attenuation image μ0, D1 and D2 are dissimilarity or distance functions, R is a regularization function, and L(g) is the negative Poisson log-likelihood given by L(g)(λ(g),μ(g))=Σi−yi(g) log(e−[Gμ(g)]i[Aλ(g)]i+ri)+(e−[Gμ(g)]i[Aλ(g)]i+ri).

Here, yi(g) represents gated measured PET emission data for sinogram bin i and gate g, G denotes a forward projection matrix representing line-integration, A denotes a forward projection operator including geometric projection and normalization, and ri represents the mean background contribution of scatters and randoms. The above objective function can be minimized with respect to λ(g), μ(g), θ(g) and λ, for example, by alternative minimization, in order to calculate the gate-dependent respiratory-phase matched attenuation image μ(g) with the respiratory motion free emission image λ.

FIG. 7 is a flow chart of a method 700 utilizing a cost function and a parameterized deformation model in accordance with aspects disclosed herein. At S710, a cost function may be determined based on a statistical model for emission scan data as a function of an emission image and a corrected attenuation image, the corrected attenuation image representing a deformation of the initial attenuation image. Note that the cost function might include, for example, a log-likelihood term, a weighted least squares term, and/or a regularizing penalty term.
At S720, a deformation model that is parameterized by deformation parameters may be provided to represent the deformation of the initial attenuation image. The parameterized deformation model may include, according to some embodiments, a linear combination of deformation basis functions to represent a deformation field or a displacement field with linear coefficients representing deformation parameters. Moreover, the deformation basis functions might include B-splines, Discrete Cosine Transform (“DCT”) basis functions, and/or polynomial functions.
At S730, the emission image and the deformation parameters may be calculated to minimize the cost function by using a numerical algorithm. Note that calculating the emission image and the deformation parameters might include generating intermediate emission images and deformation parameters by iteratively applying a numerical optimization algorithm until a pre-determined criterion is satisfied. Moreover, the numerical optimization algorithm may include alternative minimization algorithms and/or gradient based algorithms. According to some embodiments, the pre-determined criterion includes a pre-determined number of iterations, a criterion on the gradient of the cost function, a criterion on the decrease of the cost function value, or a criterion on the difference of previous and current intermediate emission images and/or intermediate deformation parameters. At S740, the corrected attenuation image may be generated by using these calculated deformation parameters.
Thus, respiratory gated data may be used to generate an improved PET image for accurate PET quantitation, for example, by an RRA procedure with the motion induced mismatch corrected attenuation images. While specific references have been made to an X-ray CT imaging system and PET system, it should be appreciated that embodiments described herein are not intended to be limited to these or to any specific type of imaging system or modality. In general, the embodiments described herein may be used for providing local motion correction for other tomography imaging combinations such as, for example, CT-SPECT and CT-MRI.
It is to be understood that not necessarily all such objects or advantages described above may be achieved in accordance with any particular embodiment. Thus, for example, those skilled in the art will recognize that the systems and techniques described herein may be embodied or carried out in a manner that achieves or optimizes one advantage or group of advantages as taught herein without necessarily achieving other objects or advantages as may be taught or suggested herein.
While only certain features of the invention have been illustrated and described herein, many modifications and changes will occur to those skilled in the art. It is, therefore, to be understood that the appended claims are intended to cover all such modifications and changes as fall within the true spirit of the invention.

## Claims

1. A method for reconstructing a nuclear medicine image, comprising:
acquiring emission scan data using an emission tomography scanner;
generating one or more anatomical images using an anatomical imaging system;
processing the anatomical images to obtain an initial attenuation image;
automatically reconstructing, by a processor, an emission image and a corrected attenuation image from the acquired emission scan data, the corrected attenuation image representing a deformation of the initial attenuation image, wherein automatically reconstructing the emission image and the corrected attenuation image from the acquired emission scan data comprises:
determining a cost function based on a statistical model for emission scan data as a function of an emission image and a corrected attenuation image, the corrected attenuation image representing a deformation of the initial attenuation image,
providing a deformation model that is parameterized by deformation parameters to represent the deformation of the initial attenuation image,
calculating the emission image and the deformation parameters to minimize the cost function by using a numerical algorithm, and
generating the corrected attenuation image by using the calculated deformation parameters;

calculating a final reconstructed emission image based on at least one of: (i) the automatically reconstructed emission image, and (ii) the corrected attenuation image; and
outputting the final reconstructed emission image by performing at least one of: (i) storing the final reconstructed emission image in a data storage system, or (ii) displaying the final reconstructed emission image on a display system.

2. The method of claim 1, wherein acquiring emission scan data comprises at least one of: (i) acquiring Positron Emission Tomography (“PET”) emission scan data using a PET scanner, (ii) acquiring Time-Of-Flight (“TOF”) PET emission scan data using a TOF PET scanner, and (iii) acquiring Single Photon Emission Computed Tomography (“SPECT”) emission scan data using a SPECT scanner.

3. The method of claim 1, wherein acquiring emission scan data comprises at least one of: (i) respiratory gating and acquiring one or more gates of the respiratory gated emission scan data, (ii) cardiac gating and acquiring one or more gates of the cardiac gated emission scan data, and (iii) dual cardiac-respiratory gating and acquiring one or more gates of the dual cardiac-respiratory gated emission scan data.

4. The method of claim 1, wherein generating one or more anatomical images comprises at least one of: (i) generating one or more x-ray Computed Tomography (“CT”) images using an x-ray CT system, and (ii) generating one or more Magnetic Resonance Imaging (“MRI”) images using an MRI system.

5. The method of claim 1, wherein generating one or more anatomical images comprises at least one of: (i) respiratory gating and providing one or more gates of the respiratory gated anatomical images, (ii) cardiac gating and providing one or more gates of the cardiac gated anatomical images, and (iii) dual cardiac-respiratory gating and providing one or more gates of the dual cardiac-respiratory gated anatomical images.

6. The method of claim 1, wherein processing the anatomical images to obtain the initial attenuation image comprises at least one of: (i) segmenting the anatomical images and assigning a pre-determined attenuation coefficient to each segmented region, organ or tissue, (ii) applying a bilinear or trilinear transformation to an x-ray CT image, and (iii) applying atlas-based, machine-learning and/or pattern-recognition methods to one or more MRI images.

7. The method of claim 1, wherein the cost function includes at least one of: (i) a log-likelihood term, (ii) a weighted least squares term, and (iii) a regularizing penalty term.

8. The method of claim 1, wherein the deformation model includes a linear combination of basis functions to represent a deformation field with linear coefficients representing deformation parameters.

9. The method of claim 8, wherein the basis functions include B-splines, Discrete Cosine Transform (“DCT”) basis functions, or polynomial functions.

10. The method of claim 1, wherein calculating the emission image and the deformation parameters comprises generating intermediate emission images and deformation parameters by iteratively applying a numerical optimization algorithm until a pre-determined criterion is satisfied.

11. The method of claim 10, wherein the numerical optimization algorithm includes alternative minimization algorithms, gradient based algorithms, conjugate gradient algorithms, coordinate descent algorithms, minorize-minimize algorithms, Levenberg-Marquardt algorithms, Newton's methods, Gauss-Newton algorithms, or ordered subsets algorithms.

12. The method of claim 10, wherein the pre-determined criterion includes the pre-determined number of iterations, a criterion on the gradient of the cost function, a criterion on the decrease of the cost function value, or a criterion on the difference of previous and current intermediate emission images and/or intermediate deformation parameters.

13. The method of claim 1, wherein calculating the final reconstructed emission image comprises taking the automatically reconstructed emission image or reconstructing an emission image from the acquired emission scan data utilizing the corrected attenuation image.

14. A non-transitory, computer-readable medium storing instructions that, when executed by a computer processor, cause the computer processor to perform a method associated with reconstructing a nuclear medicine image, the method comprising:
acquiring emission scan data using an emission tomography scanner;
generating one or more anatomical images using an anatomical imaging system;
processing the anatomical images to obtain an initial attenuation image;
automatically reconstructing, by a processor, an emission image and a corrected attenuation image from the acquired emission scan data, the corrected attenuation image representing a deformation of the initial attenuation image, wherein automatically reconstructing the emission image and the corrected attenuation image from the acquired emission scan data comprises:
determining a cost function based on a statistical model for emission scan data as a function of an emission image and a corrected attenuation image, the corrected attenuation image representing a deformation of the initial attenuation image,
providing a deformation model that is parameterized by deformation parameters to represent the deformation of the initial attenuation image,
calculating the emission image and the deformation parameters to minimize the cost function by using a numerical algorithm, and
generating the corrected attenuation image by using the calculated deformation parameters;

calculating a final reconstructed emission image based on at least one of: (i) the automatically reconstructed emission image, and (ii) the corrected attenuation image; and
outputting the final reconstructed emission image by performing at least one of: (i) storing the final reconstructed emission image in a data storage system, or (ii) displaying the final reconstructed emission image on a display system.

15. The medium of claim 14, wherein acquiring emission scan data comprises at least one of: (i) acquiring Positron Emission Tomography (“PET”) emission scan data using a PET scanner, (ii) acquiring Time-Of-Flight (“TOF”) PET emission scan data using a TOF PET scanner, and (iii) acquiring Single Photon Emission Computed Tomography (“SPECT”) emission scan data using a SPECT scanner.

16. The medium of claim 14, wherein generating one or more anatomical images comprises at least one of: (i) generating one or more x-ray Computed Tomography (“CT”) images using an x-ray CT system, and (ii) generating one or more Magnetic Resonance Imaging (“MRI”) images using an MRI system.

17. A system associated with reconstructing a nuclear medicine image, comprising:
an emission tomography scanner to acquire emission scan data;
an anatomical imaging system to generate one or more anatomical images; and
a computer system coupled to at least one of the emission tomography scanner and the anatomical imaging system to perform a method, the method comprising:
processing the anatomical images to obtain an initial attenuation image,
automatically reconstructing, by a processor, an emission image and a corrected attenuation image from the acquired emission scan data, the corrected attenuation image representing a deformation of the initial attenuation image, wherein automatically reconstructing the emission image and the corrected attenuation image from the acquired emission scan data comprises:
determining a cost function based on a statistical model for emission scan data as a function of an emission image and a corrected attenuation image, the corrected attenuation image representing a deformation of the initial attenuation image,
providing a deformation model that is parameterized by deformation parameters to represent the deformation of the initial attenuation image,
calculating the emission image and the deformation parameters to minimize the cost function by using a numerical algorithm, and
generating the corrected attenuation image by using the calculated deformation parameters,

calculating a final reconstructed emission image based on at least one of: (i) the automatically reconstructed emission image, and (ii) the corrected attenuation image, and
outputting the final reconstructed emission image by performing at least one of: (i) storing the final reconstructed emission image in a data storage system, or (ii) displaying the final reconstructed emission image on a display system.

18. The system of claim 17, wherein acquiring emission scan data comprises at least one of: (i) respiratory gating and acquiring one or more gates of the respiratory gated emission scan data, (ii) cardiac gating and acquiring one or more gates of the cardiac gated emission scan data, and (iii) dual cardiac-respiratory gating and acquiring one or more gates of the dual cardiac-respiratory gated emission scan data.

19. A method for reconstructing a nuclear medicine image, comprising:
acquiring emission scan data using an emission tomography scanner;
generating one or more anatomical images using an anatomical imaging system;
processing the anatomical images to obtain an initial attenuation image;
automatically reconstructing, by a processor, one or more emission images and one or more corrected attenuation images from the acquired emission scan data, wherein automatically reconstructing the one or more emission images and the one or more corrected attenuation images comprises:
determining a cost function based on a statistical model for emission scan data as a function of one or more emission images and one or more corrected attenuation images, at least one of the corrected attenuation images representing a deformation of the initial attenuation image,
providing a deformation model that is parameterized by deformation parameters to represent the one or more deformations of the initial attenuation image,
calculating the one or more emission images, the one or more corrected attenuation images and the deformation parameters to minimize the cost function by using a numerical algorithm, and
generating one or more final corrected attenuation images based on at least one of: (i) the calculated one or more corrected attenuation images, or (ii) the calculated deformation parameters;

calculating one or more final reconstructed emission images based on at least one of: (i) the automatically reconstructed one or more emission images, or (ii) the one or more final corrected attenuation images; and
outputting the one or more final reconstructed emission images by performing at least one of: (i) storing the one or more final reconstructed emission images in a data storage system, or (ii) displaying the one or more final reconstructed emission images on a display system.

20. The method of claim 19, wherein the one or more corrected attenuation images include a corrected attenuation image for each gate of one of: (i) respiratory gates, (ii) cardiac gates, or (iii) dual cardiac-respiratory gates.

21. The method of claim 19, wherein the one or more emission images include an emission image for each gate of one of: (i) respiratory gates, (ii) cardiac gates, or (iii) dual cardiac-respiratory gates.

22. The method of claim 19, wherein automatically reconstructing the one or more emission images and the one or more corrected attenuation images further comprises:
determining a cost function based on a statistical model for emission scan data as a function of one or more emission images and one or more corrected attenuation images, the one or more emission images including a motion-corrected emission image, at least one of the emission images representing a deformation of the motion-corrected emission image, at least one of the corrected attenuation images representing a deformation of the initial attenuation image,
providing a deformation model that is parameterized by deformation parameters to represent the deformations of the motion-corrected emission image and the initial attenuation image,
calculating the one or more emission images, the one or more corrected attenuation images and the deformation parameters to minimize the cost function by using a numerical algorithm, and
generating one or more final corrected attenuation images based on at least one of: (i) the calculated one or more corrected attenuation images, or (ii) the calculated deformation parameters.

23. The method of claim 19, wherein one or more final reconstructed emission images includes the calculated motion-corrected emission image.

24. The method of claim 19, wherein the cost function includes a dissimilarity function that quantifies the dissimilarity between the motion-corrected emission image and an emission image for a gate of one of: (i) respiratory gates, (ii) cardiac gates, or (iii) dual cardiac-respiratory gates.

25. The method of claim 19, wherein the cost function includes a dissimilarity function that quantifies the dissimilarity between the corrected attenuation image and an attenuation image for a gate of one of: (i) respiratory gates, (ii) cardiac gates, or (iii) dual cardiac-respiratory gates.