Patent Publication Number: US-8971599-B2

Title: Tomographic iterative reconstruction

Description:
BACKGROUND OF THE INVENTION 
     The present disclosure is directed generally to the reconstruction of tomographic images and, more specifically, to approaches for reconstructing tomographic images using iterative reconstruction. 
     Non-invasive imaging technologies allow images of the internal structures of a patient or object to be obtained without performing an invasive procedure on the patient or object. In particular, technologies such as computed tomography (CT), magnetic resonance imaging (MRI), positron emission tomography (PET), and single positron emission computed tomography (SPECT) use various physical principles to acquire image data and to construct tomographic images (e.g., three-dimensional representations of the interior of the human body or of other imaged structures). 
     The reconstruction of tomographic images may be a laborious and time-consuming process. For example, certain reconstruction processes attempt to utilize an approach where the image reconstruction processes are iterated until some degree of convergence based on some image quality criteria is achieved. Certain of these iterative reconstruction approaches utilize simultaneous update methods that are suitable for parallel implementation or processing. However, these types of iterative approaches may also be subject to slow convergence speed for various high frequency image features. Further, such iterative approaches may be very processor-intensive, requiring large amounts of processor power to implement effectively. 
     BRIEF DESCRIPTION OF THE INVENTION 
     In one embodiment, an iterative image reconstruction method is provided. The method includes the steps of generating an image, generating a set of calculated measurements based on the image, and generating a set of actual measurements. A determination is made whether a completion criterion is satisfied. If the completion criterion is satisfied, the image is provided as a final image. If the completion criterion is not satisfied: an image update of the image is generated, a deconvolution filter or a square root of the deconvolution filter is applied to one or both of the image or the image update, and the steps of generating the set of calculated measurements and determining whether the completion criterion is satisfied is repeated. 
     In a further embodiment, an iterative reconstruction method is provided. In accordance with this method a first image is iteratively reconstructed. The first image, when filtered using a deconvolution filter, corresponds to a current estimate of a second image that represents a scanned object. The deconvolution filter is applied during each iteration of the iterative reconstruction. 
     In another embodiment, an image processing system is provided. The image processing system includes a processing component configured to implement one or more algorithms. The one or more algorithms, when implemented by the processing component: iteratively reconstruct a sequence of images, wherein a deconvolution filter is applied during each iteration of the iterative reconstruction; and terminate the iterative reconstruction upon satisfaction of a completion criterion. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       These and other features, aspects, and advantages of the present invention will become better understood when the following detailed description is read with reference to the accompanying drawings in which like characters represent like parts throughout the drawings, wherein: 
         FIG. 1  is a diagrammatical view of a CT imaging system for use in producing images in accordance with aspects of the present disclosure; 
         FIG. 2  depicts a flow chart describing steps of an algorithm for implementing aspects of the present disclosure; 
         FIG. 3  depicts a flow chart describing steps of an algorithm for implementing further aspects of the present disclosure; 
         FIG. 4  depicts a graph displaying convergence speeds obtained using deconvolution filtering during iterative reconstruction in accordance with aspects of the present disclosure; 
         FIG. 5  depicts a graph displaying convergence speeds obtained using deconvolution filtering in conjunction with ordered sets during iterative reconstruction in accordance with aspects of the present disclosure; and 
         FIG. 6  depicts another graph displaying convergence speeds obtained using deconvolution filtering in conjunction with ordered sets during iterative reconstruction in accordance with aspects of the present disclosure 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     The present approach is directed to the use of iterative reconstruction for tomographic data sets. In particular, certain embodiments are directed to the use of model-based iterative reconstruction approaches that perform simultaneous updates. In certain implementations, an image domain deconvolution filter, such as a ramp filter or other deconvolution filter, is employed to achieve faster convergence in the iterative process. In other implementations, the deconvolution filter may be applied in the projection or other non-image domain. Thus, instead of reconstructing the final image directly, an image (i.e, a pre-image) is reconstructed which, after application of a ramp filter or other deconvolution filter, provides the final image. In certain embodiments, the deconvolution filter may be applied in combination with statistical weight modeling and/or image regularization (i.e., use of a prior). 
     The approaches disclosed herein are suitable for use with a wide range of tomographic reconstruction systems, such as, but not limited to, computed tomography (CT), magnetic resonance imaging (MRI), positron emission tomography (PET), and single positron emission computed tomography (SPECT). To facilitate explanation, the present disclosure will primarily discuss the present iterative reconstruction approaches in the context of a CT system. However, it should be understood that the following discussion is equally applicable to other tomographic reconstruction systems, such as those listed above as well as others. 
     With this in mind, an example of a computer tomography (CT) imaging system  10  designed to acquire X-ray attenuation data at a variety of views around a patient and suitable for tomographic reconstruction is provided in  FIG. 1 . In the embodiment illustrated in  FIG. 1 , imaging system  10  includes a source of X-ray radiation  12  positioned adjacent to a collimator  14 . The X-ray source  12  may be an X-ray tube, a distributed X-ray source (such as a solid-state or thermionic X-ray source) or any other source of X-ray radiation suitable for the acquisition of medical or other images. 
     The collimator  14  permits X-rays  16  to pass into a region in which a patient  18 , is positioned. A portion of the X-ray radiation  20  passes through or around the patient  18  and impacts a detector array, represented generally at reference numeral  22 . Detector elements of the array produce electrical signals that represent the intensity of the incident X-rays  20 . These signals are acquired and processed to reconstruct images of the features within the patient  18 . 
     Source  12  is controlled by a system controller  24 , which furnishes both power, and control signals for CT examination sequences. In the depicted embodiment, the system controller  24  controls the source  12  via an X-ray controller  26  which may be a component of the system controller  24 . In such an embodiment, the X-ray controller  26  may be configured to provide power and timing signals to the X-ray source  12 . 
     Moreover, the detector  22  is coupled to the system controller  24 , which controls acquisition of the signals generated in the detector  22 . In the depicted embodiment, the system controller  24  acquires the signals generated by the detector using a data acquisition system  28 . The data acquisition system  28  receives data collected by readout electronics of the detector  22 . The data acquisition system  28  may receive sampled analog signals from the detector  22  and convert the data to digital signals for subsequent processing by a processor  30  discussed below. Alternatively, in other embodiments the digital-to-analog conversion may be performed by circuitry provided on the detector  22  itself. The system controller  24  may also execute various signal processing and filtration functions with regard to the acquired image signals, such as for initial adjustment of dynamic ranges, interleaving of digital image data, and so forth. 
     In the embodiment illustrated in  FIG. 1 , system controller  24  is coupled to a rotational subsystem  32  and a linear positioning subsystem  34 . The rotational subsystem  32  enables the X-ray source  12 , collimator  14  and the detector  22  to be rotated one or multiple turns around the patient  18 . It should be noted that the rotational subsystem  32  might include a gantry upon which the respective X-ray emission and detection components are disposed. Thus, in such an embodiment, the system controller  24  may be utilized to operate the gantry. The linear positioning subsystem  34  may enable the patient  18 , or more specifically a table supporting the patient, to be displaced within the bore of the CT system  10 . Thus, the table may be linearly moved within the gantry to generate images of particular areas of the patient  18 . In the depicted embodiment, the system controller  24  controls the movement of the rotational subsystem  32  and/or the linear positioning subsystem  34  via a motor controller  36 . 
     In general, system controller  24  commands operation of the imaging system  10  (such as via the operation of the source  12 , detector  22 , and positioning systems described above) to execute examination protocols and to process acquired data. For example, the system controller  24 , via the systems and controllers noted above, may rotate a gantry supporting the source  12  and detector  22  about a subject of interest so that X-ray attenuation data may be obtained at a variety of views relative to the subject. In the present context, system controller  24  may also includes signal processing circuitry, associated memory circuitry for storing programs and routines executed by the computer (such as routines for executing image processing techniques described herein), as well as configuration parameters, image data, and so forth. 
     In the depicted embodiment, the image signals acquired and processed by the system controller  24  are provided to a processing component  30  for reconstruction of images. The processing component  30  may be one or more conventional microprocessors. The data collected by the data acquisition system  28  may be transmitted to the processing component  30  directly or after storage in a memory  38 . Any type of memory suitable for storing data might be utilized by such an exemplary system  10 . For example, the memory  38  may include one or more optical, magnetic, and/or solid state memory storage structures. Moreover, the memory  38  may be located at the acquisition system site and/or may include remote storage devices for storing data, processing parameters, and/or routines for iterative image reconstruction described below. 
     The processing component  30  may be configured to receive commands and scanning parameters from an operator via an operator workstation  40 , typically equipped with a keyboard and/or other input devices. An operator may control the system  10  via the operator workstation  40 . Thus, the operator may observe the reconstructed images and/or otherwise operate the system  10  using the operator workstation  40 . For example, a display  42  coupled to the operator workstation  40  may be utilized to observe the reconstructed images and to control imaging. Additionally, the images may also be printed by a printer  44  which may be coupled to the operator workstation  40 . 
     Further, the processing component  30  and operator workstation  40  may be coupled to other output devices, which may include standard or special purpose computer monitors and associated processing circuitry. One or more operator workstations  40  may be further linked in the system for outputting system parameters, requesting examinations, viewing images, and so forth. In general, displays, printers, workstations, and similar devices supplied within the system may be local to the data acquisition components, or may be remote from these components, such as elsewhere within an institution or hospital, or in an entirely different location, linked to the image acquisition system via one or more configurable networks, such as the Internet, virtual private networks, and so forth. 
     It should be further noted that the operator workstation  40  may also be coupled to a picture archiving and communications system (PACS)  46 . PACS  46  may in turn be coupled to a remote client  48 , radiology department information system (RIS), hospital information system (HIS) or to an internal or external network, so that others at different locations may gain access to the raw or processed image data. 
     While the preceding discussion has treated the various exemplary components of the imaging system  10  separately, these various components may be provided within a common platform or in interconnected platforms. For example, the processing component  30 , memory  38 , and operator workstation  40  may be provided collectively as a general or special purpose computer or workstation configured to operate in accordance with the aspects of the present disclosure. In such embodiments, the general or special purpose computer may be provided as a separate component with respect to the data acquisition components of the system  10  or may be provided in a common platform with such components. Likewise, the system controller  24  may be provided as part of such a computer or workstation or as part of a separate system dedicated to image acquisition. 
     With the foregoing example of one suitable image acquisition system in mind, the reconstruction of image data generated using such a tomographic imaging system is discussed. For example, a tomographic scanner such as system  10  may measure a dataset y that is a function of the object or image dataset x through a system operator A. In the simplest case this is a linear model and can be written as a matrix multiplication:
 
 y=A*x   (1)
 
This equation gets more complicated when including exponential attenuation, spatial-dependent blurring, polychromaticity, and so forth. However, in all cases there is a notion of a system matrix A with elements a ij , which relate an individual measurement y i  to an individual image voxel x j . In different imaging modalitites, the property corresponding to x may vary. For example, x can be the activity in PET or SPECT, the attenuation coefficient or material density in CT, or any other quantity to be tomographically reconstructed. Correspondingly, y is typically a detector signal, representing a detected energy, a number of photons, or an electric current. Thus, x is typically the property being tomographically reconstructed, y is typically the detector signal or measurement at the detector and the element a ij  of matrix A provides the necessary transformation to get from one detector value y i  to a particular voxel x j .
 
     The most direct approach to reconstruct x from y is to use an analytic inversion approach. In the simplest case, this may be represented as:
 
 x=Â   −1   *y   (2)
 
where Â is an estimate of A, and Â −1  is an analytically computed inverse of Â, such as a filtered backprojection (FBP).
 
     The main alternative to the analytic inversion process described above is to use iterative inversion or iterative reconstruction. In such an iterative approach, for every iteration a current estimate of the image x (or a subset of it) is run through a forward model including Â, resulting in an estimated dataset ŷ (or a subset of it). The estimated data set ŷ is compared with the measured dataset y, and the image x is updated based on the residual discrepancy. In some cases the discrepancy between y and ŷ is weighted (implicitly or explicitly) based on the estimated statistical reliability of a measurement i. Approaches where a weighting is used based on some estimate of the statistical reliability of the measurements may be described as statistical iterative reconstruction approaches. 
     One example of a statistical iterative reconstruction is an ordered subsets expectation maximization (OSEM) approach, which is an implicitly-weighted statistical iterative reconstruction algorithm that may be used in emission tomography. In an OSEM algorithm, the entire image is updated every iteration based on a subset of the measurements. Another example, of a statistical iterative reconstruction is an iterative coordinate descent (ICD) approach, which is an explicitly-weighted statistical iterative reconstruction algorithm usually targetted at transmission tomography. In an ICD reconstruction, one pixel of the image is updated every iteration based on all the measurements. While the examples and embodiments discussed herein are generally presented in the context of statistical iterative reconstruction, it should be appreciated that the approaches disclosed herein are suitable for other, non-statistical, iterative reconstruction approaches. One example of a non-statistical iterative reconstruction algorithm is the algebraic reconstruction technique (ART). In an ART reconstruction, all pixels are updated every iteration based on one single measurement. 
     These various approaches have different strengths and weaknesses. For example, single-coordinate approaches like ICD update one pixel at the time, resulting in fast convergence of the high frequencies (i.e., few iterations), but each iteration is slow. Conversely, simultaneous update approaches like OSEM update all pixels at the same time, resulting in slow convergence (i.e., many iterations), but each iteration is relatively fast. 
     Various approaches may be employed to improve the efficiency of the update step in iterative reconstruction approaches that employ simultaneous updates. For example, algebraically, it may be desirable to use a Newton Raphson step to compute an efficient update step:
 
 x   n+1   =x   n   +H   −1 *∇
 
where ∇ is the easy-to-compute gradient and H is the Hessian matrix (i.e., matrix of second derivatives). However, the size of the Hessian matrix is typically about 1 million×1 million in two-dimensional reconstructions, and its size can be orders of magnitude greater in three-dimensional cases. Therefore, it is typically not feasible to directly invert the Hessian matrix. Further, approaches that approximate or simplify this algebraic approach typically result in poor convergence (i.e., many iterations).
 
     With the foregoing in mind, the present approach attempts to optimize a statistical cost function (such as a weighted least squares, point weighted least squares, maximum likelihood, maximum a posteriori probability, or other cost function) efficiently without inverting or directly manipulating the Hessian matrix. In particular, in one embodiment a change of variables is implemented. For example, in one implementation, instead of reconstructing x directly, an image z is instead reconstructed which, when filtered, provides x. In one such embodiment the filter employed is a deconvolution filter F, such as a ramp filter or other high pass filters, such that:
 
 ŷ=A*F*z   (4)
 
or
 
 ŷ=B*z   (5)
 
In two-dimensional reconstruction the deconvolution filter can be provided as a two-dimensional filter kernel while in three-dimensional reconstruction the deconvolution filter can be provided as a three-dimensional filter kernel. In one implementation, the two-dimensional or three-dimensional filter kernels consist of a large positive central element while the remaining elements (i.e., the surrounding or peripheral elements) consist of small negative or small positive values, like a ramp filter. The filter has the effect of a deconvolution and by incorporating the filter in the reconstruction, neighboring pixel values are decoupled, thereby accelerating reconstruction. In one embodiment, the deconvolution filter is implemented as a ramp filter operation in the image domain using fast Fourier transform (FFT) operations.
 
     In certain embodiments the image z may be initially estimated using a simple backprojection. In other embodiments, the image z may be provided by a filtered backprojection reconstruction (FBP). In such an embodiment, the reconstruction used as the image z may or may not be smoothed. In one embodiment in which an iterative three-dimensional cone-beam CT reconstruction is performed, a three dimensional filtered backprojection (such as a backprojection employing Katsevich&#39;s filter) may be employed. 
     Typically iterative processing continues until a specified completion criterion is met which signifies that no additional iterations are to be performed. In certain embodiments, the completion criterion may be a predetermined number of iterations, such as an empirically determined number of iterations corresponding to a desired degree of image enhancement. In other embodiments, the completion criterion may generally reflect the degree of convergence attained between the calculated and measured images or image characteristics. That is, the completion criterion may be based upon a desired degree of convergence between the calculated and measured images or image characteristics such that the degree of similarity is deemed sufficient and/or such that further iterations will not further converge the calculated and measured images or image characteristics to an appreciable or desired extent. In certain implementations, the completion criterion may be implemented as a cost function, such as a statistically weighted cost function C, that may be optimized to assess convergence (and, thus, to determine when the iterative processing is complete) such that:
 
 C ( ŷ ( x ))≡ C ( ŷ ( z )).  (6)
 
A prior or regularizer can be defined as in other iterative reconstruction approaches, either in the x-domain or in the z-domain. Upon optimization of the cost function (generally signifying little or no image improvement would be attained by additional iterations) convergence is achieved and image z may be converted to image space by application of the appropriate filter:
 
 x=F*z.   (7)
 
The change of variables as discussed herein does not change the final image x, but instead makes the update step almost perfectly separable from the remainder of the reconstruction process. That is, the unknowns (e.g., image pixel values) are decoupled, allowing simultaneous update algorithms to converge more quickly. Since every image can be written as the convolution of an image and filter (e.g., a ramp filter), the present approach does not pose any limitation on the reconstructed images.
 
     With the foregoing general discussion and examples in mind, various specific implementations are discussed below to further illustrate the present approach. In one example, the above-described change of variable approach is implemented in the context of a gradient descent algorithm. In this example, in the original image space x, the image update equation can be written as:
 
 x   (n+1)   =x   (n)   +αA   T   w ( y−Ax   (n) ).  (8)
 
After changing variable from x to z, the image update equation becomes
 
 z   (n+1)   =z   (n)   +αF   T   A   T   w ( y−AFz   (n) ).  (9)
 
where F is a suitable deconvolution filter, such as a ramp filter, or the square root of a ramp filter. When effective convergence is reached, we may convert z to the image space by computing x=F*z, as disclosed with respect to equation (7) above.
 
     For example, turning to  FIG. 2 , flowchart  80  depicts steps associated with one implementation in which the square root of the ramp filter is applied at different points in the image reconstruction process. In the depicted example, the square root of the ramp filter is initially applied (block  82 ) to a pre-image  84  (e.g., an image which when ramp-filtered, gives the actual image) to generate image  86 . In addition, the square root of the ramp filter is applied (block  88 ) to the image update generated by the optimization process  90  used to update the pre-image  84  as the iterative process proceeds. 
     Based upon the noise statistics (e.g., the differences between the measured and calculated sinograms) a cost function may be evaluated to evaluate the degree of difference between the simulation  92  and the measurements  94  to determine whether additional iterations will be applied. In one embodiment a suitable cost function used to assess convergence may be a quadratically penalized weighted least squares (WLS) cost function, such as:
 
 C ( x )=½ ∥y−Ax∥   D   2 +½ βx   T   Rx   (10)
 
where D is the statistical weighting, and β is the regularization strength. In general, the relevant completion criteria may be based on a projection domain comparison, a suitable cost function, image domain statistics, a MAP cost function, or iteration number, or any other suitable convergence indication to determine algorithm termination.
 
     If no additional iterations are needed, the image  86  is determined to be the final image. Otherwise, an image update is generated by the optimizer  90  based on the noise statistics  96 . In certain embodiments, the optimizer  90  may be a formula or algorithm implemented on a processor or other computing component to update the image  86 . The updated image may be ramp filtered to approximate inversion of the Hessian matrix and to generate an updated image for processing in the next iteration. As noted above, other suitable deconvolution filters (such as high pass filters) may be applied in other embodiments. In certain embodiments the ramp filter (or other deconvolution filter) is applied in the image domain. However, in other embodiments the ramp filter (or other deconvolution filter) is applied in the projection domain (or other non-image domain, depending on the imaging modality). 
     In yet other embodiments, the ramp filter (or other deconvolution filter) may be applied at different times or in different forms. For example, both sides of Equation (9) can be multiplied by filter F (e.g., a ramp filter or other deconvolution filter) to obtain:
 
 Fz   (n+1)   =Fz   (n)   +αFF   T   A   T   w ( y−AFz   (n) )  (11)
 
which can be re-written as
 
 x   (n+1)   =x   (n)   +αFF   T   A   T   w ( y−Ax   (n) ).  (12)
 
In this implementations, the iteration is now in the original image space. To implement equation (12), an image domain filter FF T  may be applied to the image increment α A T w(y−Ax (n) ). As discussed above, to approximate the inverse of the Hessian matrix, F may be implemented as the square root of a ramp filter, i.e., FF T  becomes a standard ramp filter. Including a prior term in equation (12) gives
 
 x   (n+1)   =x   (n)   +αFF   T   [A   T   w ( y−Ax   (n) )−β∇ R ( X   (n) )]  (13)
 
where ∇R is the gradient of the prior cost function and β is a regularization parameter.
 
     With the foregoing in mind,  FIG. 3  provides a flow chart depicting one example of an implementation of the present approach in which a ramp filter is employed as part of a model-based iterative reconstruction process  100 . In accordance with this example, an image  86  may be derived based on image statistics  102  (e.g., an image prior or regularizer) determined for the image reconstruction. The image  86  may be processed based upon the system optics  104  (e.g., a forward model determined for the corresponding imaging system) to derive a simulation  92  (e.g., a calculated sinogram). Noise statistics  96  (e.g., differences or residual values) may be derived based upon the simulation (e.g., the calculated sinogram) and a set of measurements  94  (e.g., a measured sinogram). Based upon the noise statistics (e.g., the differences between the measured and calculated sinograms) a cost function may be evaluated to determine whether additional iterations will be applied. If no additional iterations are needed, the image  86  is determined to be the final image. Otherwise, an image update is generated by the optimizer  90  based on the noise statistics  96 . In certain embodiments, the optimizer  90  may be a formula or algorithm implemented on a processor or other computing component to update the image  86 . The updated image may be ramp filtered (block  108 ) to approximate inversion of the Hessian matrix and to generate an updated image for processing in the next iteration. Other suitable deconvolution filters (such as high pass filters) may be applied in other embodiments. In certain embodiments the ramp filter (or other deconvolution filter) is applied in the image domain. However, in other embodiments the ramp filter (or other deconvolution filter) is applied in the projection domain (or other non-image domain, depending on the imaging modality). 
     With the flowchart  80  of  FIG. 3  in mind, in one embodiment a gradient descent approach may be employed for minimization, in accordance with:
 
 {circumflex over (x)}   (n+1)   ={circumflex over (x)}   (n) −α( A   T   D ( A{circumflex over (x)}   (n)   −y )+β∇ R{circumflex over (x)}   (n) )  (14)
 
where β is the regularization strength. In one such implementation, adding a positive-definite deconvolution filter (preconditioner) results in:
 
 {circumflex over (x)}   (n+1)   ={circumflex over (x)}   (n)   −αM ( A   T   D ( A{circumflex over (x)}   (n)   −y )+β∇ R{circumflex over (x)}   (n) )  (15)
 
where M is a positive-definite deconvolution filter. In one implementation:
 
 M =( A   T   DA+βR ) −1 .  (16)
 
Further, as discussed herein, the Hessian matrix may be approximated, such as with a spatially-weighted ramp filter, where:
 
( A   T   DA+βR ) −1 ≈( A   T   DA ) −1 ≈(   A   T   A   ) −1 ≈   −1 [ramp filter]   −1   (17)
 
and where:
 
 =diag{κ j }, and  (18)
 
     
       
         
           
             
               
                 
                   
                     κ 
                     j 
                   
                   = 
                   
                     
                       
                         
                           ∑ 
                           i 
                           
                               
                           
                         
                         ⁢ 
                         
                           
                             A 
                             ij 
                             2 
                           
                           ⁢ 
                           
                             
                               D 
                               ii 
                             
                             / 
                             
                               
                                 ∑ 
                                 i 
                                 
                                     
                                 
                               
                               ⁢ 
                               
                                 A 
                                 ij 
                                 2 
                               
                             
                           
                         
                       
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   19 
                   ) 
                 
               
             
           
         
       
     
     Such a ramp-filtered model-based iterative reconstruction may be used to accelerate convergence by approximating the Hessian matrix using a deconvolution filter. In one example of an implementation, a pixelized Shepp-Logan phantom (64×64) was imaged to generate a 64×64 two-dimensional image. A system matrix A was computed and a noiseless projection computed by a forward projection using the system matrix. Gaussian noise was added with a diagonal covariance matrix D −1 . A q-generalized Gaussian Markov random field (qGGMRF) regularization term was employed as the prior and penalized differences between eight neighboring pixels in two-dimension. The regularization strength β was chosen manually. The image data was iteratively reconstructed using five different approaches: iterative coordinate descent (ICD), steepest descent (SD), conjugate gradient (CG), and SD and CG implementations employing ramp-filter approximation of the Hessian gradient, as discussed herein. A filtered backprojection image was employed as the initial image in each reconstruction. The ICD reconstruction employed random update order with the non-negativity constraint disabled. The CG reconstruction was of the Polak-Ribiere form. Results of these reconstructions are depicted in  FIG. 4 . As depicted in  FIG. 4 , the ramp filter based approaches converge more quickly than the non-ramp filter based counter-parts, and with a similar number of iterations as the ICD reconstruction, which only updates a single pixel each iteration as opposed to the other reconstruction approaches which perform simultaneous pixel updates each iteration. 
     To further improve convergence speeds, other iterative optimization techniques (such as line search, conjugate gradient, ordered subsets, Nesterov optimization, and so forth) can be combined with filter operation to provide further acceleration. For example, Nesterov optimal gradient iteration may be performed in certain embodiments in conjunction with application of a deconvolution filter. For example, in conventional gradient descent for L2-error minimization may be performed in accordance with:
 
 x   n+1   =x   n   +τA   T ( y−Ax   n )  (20)
 
However, in implementations in which Nesterov optimization is employed, the Nesterov optimal gradient for L2-error minimization uses two prior iterations to generate the next solution update in accordance with:
 
                       x     n   +   1       =       g   n     +     τ   ⁢           ⁢       A   T     ⁡     (     y   -     Ag   n       )             ⁢     
     ⁢     where   ⁢     :               (   21   )                 g   n     =       x   n     +         n   -   1       n   +   2       ⁢     (       x   n     -     x     n   -   1         )                 (   22   )               
In certain implementations (such as iterative CT reconstruction) it may be useful to utilize edge-preserving priors to recover sharp boundaries and maintain image resolution. However, because edge-preserving priors are space-variant, a ramp filter may not be a suitable approximation to the inverse of the Hessian matrix. In such implementations, the ramp filter may be applied to the data fit term alone, leaving the prior term unfiltered:
 
 x   (n+1)   =x   (n)   +αFF   T   [A   T   w ( y−Ax   (n) )]−β∇R( x   (n) ).  (23)
 
That is, in such embodiments, the filtered likelihood term and the unfiltered prior term are combined to provide edge preservation.
 
     In a further embodiment, it is recognized that one theoretical constraint is that the image domain ramp filter has infinite support, which would require reconstruction of an image with infinite support. To obviate this possible constraint in implementations where it might arise, the order of ramp filter application and backprojection may be switched such that the ramp-filter is implemented in the projection domain:
 
 x   (n+1)   =x   (n)   +α{[A   T   fw ( y−Ax   (n) )]−β∇R( x   (n) )}  (24)
 
where f is a projection domain ramp filter. In one such embodiment, rebinning of projection data (e.g., cone-to-parallel) can be performed so that projection domain filtering is more accurate.
 
     In a further example, ordered subsets (OS) may be used to further accelerate convergence in accordance with the presently disclosed approaches. In such an implementation, instead of using all projections simultaneously to compute an updated image estimate, a subset of projections can be used in each iteration. For instance, an OS version of equation (24) may be represented as
 
 x   (n+1)   =x   (n)   +αFF   T   [A   k   T   f   k   w   k ( y   k   −A   k   x   (n) )−β∇ R ( x   (n) )],  k=n  mod  K  
 
where K is the total number of subsets, and A k , f k , w k , and y k , represent the subset k of A, f, w, and y, respectively. Turning to  FIGS. 5 and 6 , respective graphs are depicted demonstrating the convergence speeds of iterative coordinate descent (ICD), steepest descent (SD), conjugate gradient (CG) reconstructions as well as SD and CG implementations employing ramp-filter approximation of the Hessian matrix, ordered subsets, or both ramp filtering and ordered subsets. In the depicted examples, it can be observed that ramp filtering and ramp filtering with ordered subsets yield superior convergence speeds relative to other reconstruction approaches in both cases.
 
     Technical effects of the invention include iterative reconstruction of tomographic images by applying a deconvolution filter in the image domain or, in other embodiments, in a non-image domain. One technical effect of the invention it to allow iterative reconstruction of a tomographic image by approximating the inversion of the Hessian matrix using a deconvolution filter. Another technical effect of the invention is the iterative reconstruction not of a desired image directly, but instead iteratively reconstructing an image which when deconvolution filtered (e.g., ramp filtered) yields the desired image. 
     This written description uses examples to disclose the invention, including the best mode, and also to enable any person skilled in the art to practice the invention, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the invention is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims if they have structural elements that do not differ from the literal language of the claims, or if they include equivalent structural elements with insubstantial differences from the literal languages of the claims.