Abstract:
A technique for simplifying the processing of tomographic imagery (e.g., CT and PET) via statistical smoothing is disclosed. An image smoothness criterion is imposed in the projection domain while identifying data that matches observed data within a desired tolerance. The present invention thus amounts to a computation entirely in the projection domain. As a result, the standard (and somewhat expensive) practice of numerically iterating between the image and projection domains is avoided. Also, the fundamental system of equations having on the order of m p unknowns is decoupled into p independent systems of m unknowns each, where there are p projections of m measurements each. The present invention thus provides a processing technique defined in a single domain, which may be carried out via parallel processing.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS  
       [0001]    The present invention claims priority from U.S. Provisional Application No. 60/357,151 entitled “Simplifying Tomography With Decoupled Regularization” filed Feb. 14, 2002, the entirety of which is hereby incorporated by reference. 
     
    
     
       BACKGROUND OF INVENTION  
         [0002]    Tomographic images are produced by converting observed projections (data) into an image. For example, in x-ray CT imaging, x-ray beams are directed at an object and the beams are attenuated by various amounts due to varying structures within the object. On the other side of the object, the attenuated beams are measured with detectors. Such projections are produced at many different angles around the object. Not only are these measurements noisy, the relative noise level depends upon the amount of attenuation. Projections through dense materials such as bone and especially metal have lower signal-to-noise ratios than projections through flesh, water, or other less dense materials. Coping with the large and spatially varying fluctuations in the number of detected photons often requires a statistical smoothing technique, also known as regularization, to improve the image.  
           [0003]    In the absence of noise, a basic problem of CT scans amounts to determining an unknown image f=f(x, y) from its projections, or Radon transform Rf where (Rf)(t, θ)=∫∫f(xy)δ(t−x cos θ−y sin θ) dxdy, where (x, y) are planar coordinates, t is the location along each projection, and θ∈[0, π] is the orientation of the projection. Unfortunately, because measurements are never perfect, what is actually observed is the noisy projection data g=g(t, θ). To emphasize the essentials of the tomography problem, the unknown f=f (x,y) and the observed g=g(t, θ) are viewed as functions, although the implementation is discrete. Next, assume the standard independence and locality assumptions that the likelihood P(g |f) or conditional distribution of g given f equals Πp(g(t, θ) |(Rf)(t θ,)), where the product is over all (t, θ). This product may be nonzero if a finite number of factors is taken. Here, the observed number of x-ray photons b exp(−g(t, θ)) is Poisson distributed with mean w(t, θ):=b exp(−(Rf)(t, θ)), where b is the mean number of incident photons and f(x, y) is the attenuation coefficient at (x, y). The task is to infer the image f given noisy projections g.  
           [0004]    Because this inverse problem is ill-posed, one typically imposes extra constraints on f. In regularization, a type of penalized maximum likelihood, inferring f amounts to finding that f which minimizes—1n P(g |f)+p(f)g, where p(f) characterizes the extra constraint on f Here we impose smoothness via the gradient  
         ∇   f     =     (         ∂   f       ∂   x       ,       ∂   f       ∂   y         )                           
 
           [0005]    with the quadratic penalty p(f)=β∥∇∥ 2 , where β&gt;0 and ∥f∥ 2 =(f,f), using the inner product (f 1 ,f 2 )=∫f 1 f 2 . Following Saver et al., “A local update strategy for iterative reconstruction from projections”, IEEE Transactions of signal Processing, vol. 41, no.2, pp. 534-548 (1993), and to simplify the presentation, I approximate—1n P(g |f) with the quadratic form ∥g−Rf∥ w   2  where ∥g∥ w   2  (g, Wg) is a weighted norm with (diagonal) weight operator W satisfying (Wg)(t, θ)=w(t, θ)g(t, θ). The weight w(t, θ) is small for those rays passing through dense materials such as bone or metal, and larger otherwise. This formulation of tomography requires solving the following “hard” optimization problem.  
           [0006]    Problem I (Coupled Regularization). Given the projection data g, find the image f that minimizes  
           ∥g−Rf∥ w   2 +⊕∥∇f∥ 2 .  
           [0007]    (Because w depends on Rf, which is unknown, we use g instead to obtain a quadratic functional above. In general, however, one could solve the full nonlinear problem and only make this substitution as a first guess of w. One can compute f iteratively, using the current best estimate of Rf to set w.) To proceed, let Δ denote  
             ∂   2       ∂     x   2         +       ∂   2       ∂     y   2           ,                         
 
           [0008]    the Laplacian in the plane. Using integration by parts and zero boundary conditions, recall that ∥∇f∥ 2 =(f, −Δf). Then the (linear) Euler-Lagrange equation for Problem 1 is  
                     R   *        WRf     -       β        (         ∂   2       ∂     x   2         +       ∂   2       ∂     y   2           )          f       =       R   *        Wg       ,           (   1   )                               
 
           [0009]    where f is unknown and A* denotes the adjoint (or transpose) of linear operator A (R* is also known as the backprojection operator). By examining equation (1), we see that Problem 1 is hard in two related ways.  
           [0010]    First, the problem constraints occur in two different domains. Fidelity to the data (∥g−Rf∥ w   2 ) is enforced in the Radon domain {(t, θ)}, while smoothness (∥□f∥ 2 ) is imposed in the image domain {(x, y)}. Thus in equation (1) the operators R and R* are for shuffling back and forth between these domains; iterative solution techniques typically compute these forward and backprojections explicitly and often at great expense.  
           [0011]    Second, observe that equation (1) is a coupled equation in the two-variable functions, i.e., in the large set of variables {f (x, y), for all x, y} under some reasonable discretization of x and y. The coupling arises first because both x and y derivatives are present; in addition, R and R* are integral operators, and so are not even local. The computational difficulty in solving equation (1) and related tomographic problems (e.g., emission) has spawned a great deal of work in optimization.  
           [0012]    Thus, the need exists for a processing technique which can be formed in a single domain.  
         SUMMARY OF THE INVENTION  
         [0013]    A method of reconstructing an image is comprised of scanning an object to obtain projection data. The projection data is decomposed by projection angle. For each projection angle, a solution is obtained to an inference problem which is defined in a single domain. The solution may include, for example, finding projections h=Rf minimizing:  
                  g   -   h          w   2     +       β   ′          (     h   ,       ℐ     -   3          h       )         ,     where                   β   ′          :          =     β     4                 π       .                             
 
           [0014]    The corresponding Euler-Lagrange equation is:  
             Wh+β′I   −3   h=Wg   (2)  
           [0015]    An image is reconstructed from the obtained solution by using the standard method of filtered backprojection, which computes f=R −1 h.  
           [0016]    The present invention is directed to a technique for simplifying the processing of tomographic imagery (e.g., CT and PET) via statistical smoothing that decouples the fundamental system of equations having on the order of m p unknowns into p independent systems of m unknowns each, where there are p projections of m measurements each. That allows for more rapid processing of medical imagery for clinical use because (1) solving a system of equations having on the order of m p unknowns typically requires much more computation and memory than solving p systems of equations in m unknowns each and (2) the p systems of equations can be solved on separate processors in parallel because they are independent systems of equations.  
           [0017]    The basic idea of regularization is to infer a smooth image whose simulated projections approximate the observed but noisy data. The difficulty is that the standard regularization specifies the two basic constraints in different domains. More specifically, one constraint on a property related to the smoothness of the image must be specified in the image domain while a second constraint on the closeness to the observed data is specified in the projection domain. Enforcing these two constraints computationally requires a constant numerical conversion between the two domains during processing which is difficult and computationally expensive. The present invention (1) imposes an image smoothness criterion in the projection domain before any computations begin, and (2) identifies data that matches observed data within a desired tolerance. The present invention thus amounts to a computation entirely in the projection domain. As a result, the standard (and somewhat expensive) practice of numerically iterating between the image and projection domains is avoided. An added benefit of the analytical conversion of the smoothness constraint to the projection domain is the decoupling of a large system of regularization equations into many small systems of simpler regularization equations. Those advantages and benefits, and others, will be apparent from the description of the invention below. 
       
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0018]    For the present invention to be readily understood and easily practiced, the present invention will now be described, for purposes of illustration and not limitation, in connection with preferred embodiments in which:  
         [0019]    [0019]FIG. 1 and FIG. 2 illustrate typical hardware for generating data on which the method of the present invention may be practiced;  
         [0020]    [0020]FIG. 3 is a flow chart illustrating the basic steps of the method of the present invention;  
         [0021]    [0021]FIG. 4 illustrates streaking artifacts due to the presence of metal in an x-ray CT slice of a hip;  
         [0022]    [0022]FIG. 5 illustrates a simulated projection (solid curve, g (−, θ)) oriented at θ≈45° which is noisy particularly in the 150 position due to the presence of metal; the projection is adaptively smoothed by the decoupled regularization of the present invention (dotted curve, h (−, θ)), because there the weights are the lowest;  
         [0023]    [0023]FIG. 6A illustrates simulated projections g (t, θ), or Radon transform, for the cropped region of  
         [0024]    [0024]FIG. 7A which shows a bright band due to the metal;  
         [0025]    [0025]FIG. 6B illustrates the band of FIG. 6A after being smoothed by decoupled regularization, resulting in h (t, θ);  
         [0026]    [0026]FIG. 6C illustrates the smoothing being localized on the metal band as a result of performing I −1  (a kind of differentiation and the first step in filtered backprojection) on the difference g−h; and  
         [0027]    [0027]FIGS. 7A and 7B illustrate that the streaking artifacts in the cropped image region around the metal (FIG. 7A) are significantly reduced by the application of the decoupled regularization method of the present invention; note that the smoothing effect is primarily near the site of the metal while the rest of the image remains sharp. 
     
    
     DESCRIPTION OF THE INVENTION  
       [0028]    Referring to FIGS. 1 and 2, a computed tomographic (CT) imaging system  10  is shown as including a gantry  12 . Gantry  12  has an x-ray source  14  that projects a beam of x-rays  16  toward a detector array  18  on the opposite side of gantry  12 . Detector array  18  is formed by detector elements  20  which together sense the projected x-rays that pass through an object  22 , for example a medical patient. Each detector element  20  produces an electrical signal that represents the intensity of an impinging x-ray beam and hence the attenuation of the beam as it passes through patient  22 . During a scan to acquire x-ray projection data, gantry  12  and the components mounted thereon rotate about a center of rotation  24 .  
         [0029]    Rotation of gantry  12  and the operation of x-ray source  14  are governed by a control mechanism  26  of CT system  10 . Control mechanism  26  includes an x-ray controller  28  that provides power and timing signals to x-ray source  14  and a gantry motor controller  30  that controls the rotational speed and position of gantry  12 . A data acquisition system  32  in control mechanism  26  samples analog data from detector elements  20  and converts that data to digital signals for subsequent processing by a computer  36  which then stores the image in a mass storage device  38 .  
         [0030]    Computer  36  also receives commands and scanning parameters from an operator via console  40  that has a keyboard. An associated display  42  allows the operator to observe the reconstructed image and other data from computer  36 . The operator supplied commands and parameters are used by computer  36  to provide control signals and information to data acquisition system  32 , x-ray controller  28  and gantry motor controller  30 . In addition, computer  36  operates a table motor controller  44  which controls a motorized table  46  to position patient  22  in gantry opening  48 .  
         [0031]    In a multi-slice imaging system  10 , detector array  18  comprises a plurality of parallel detector rows, wherein each row comprises a plurality of individual detector elements  20 . An imaging system  10  having a multi-slice detector array  18  is capable of providing a plurality of images representative of a volume of object  22 . Each image of the plurality of images corresponds to a separate “slice” of the volume. The “thickness” or aperture of the slice is dependent upon the thickness of the detector rows.  
         [0032]    The hardware described in conjunction with FIGS. 1 and 2 is standard, commercially available hardware. Additionally, because the likelihood function may depend upon the particular imaging modality, I have described the present invention in connection with an x-ray CT scanning for concreteness, although my technique applies more broadly, including other tomographic imaging modalities such as PET. Furthermore, in the following discussion, I have focused on the two-dimensional problem with standard parallel geometry, but the ideas readily extend to three-dimensions and other scanning geometries, e.g. fan geometries.  
         [0033]    This invention formulates the entire inference problem in a single domain. As demonstrated below, working solely in the Radon domain decouples the large joint estimation problem into many smaller ones.  
         [0034]    Recall that in regularized tomography the standard technique for inverting the Radon transform is filtered backprojection (FBP). Let the Fourier transform of g=g(t, θ) with respect to the first variable t be denoted (F 1 g)(r, θ(2π) −1/2 ∫g(t, θ)e −irt dt, where r is the spatial frequency along a projection (for each fixed θ).  
         [0035]    Given function h=h(t, θ), the Riesz potential is the linear operator I α  satisfying (F 1 I α h)(t, θ)=|τ| −α (F 1 h)(t, θ)  
         [0036]    Given noise-free observations h=Rf, one can solve for the unknown f by directly implementing the following classical formula for the inverse of the Radon transform:  
         [0037]    Fact 1 (Filtered Backprojection).  
              -   1       =       1     4                 π            R   *            ℐ     -   1       .                             
 
         [0038]    By defining □ as the Laplacian ∂ 2 /∂t 2  each projection, we can also use the Riesz potential to define the “square root” of −□(Clearly I α1 I α2 =I α1+α2 .)  
         [0039]    Fact 2. I −2=−□   
         [0040]    The proof is that the Fourier transform of −□ is τ 2 =|τ| 2  Now recall the Fourier slice theorem, which says that the two-dimensional Fourier transform off, evaluated at polar coordinates (r, θ), is (F 1 Rf)(r, θ). Using this theorem, one can relate the two Laplacians Δ and □ because the two-dimensional Fourier transform of −Δ is u 2 +v 2  (where u and v are spatial frequencies for x and y, respectively), or |r| 2  in polar coordinates, which is the one-dimensional Fourier transform of −□. For details and the extension to higher dimensions, see Natterer, “The Mathematics of Computerized Tomography”, Wiley, Chirchester (1986) Chapters 1 and 7. This is called “intertwining”.  
         [0041]    Fact 3. The Radon transform R intertwines Δ and □, i.e.,  
         RΔ=□R.  
         [0042]    By applying these facts to Problem 1, we need only solve an equivalent set of smaller equations, thus partially “decoupling” the regularization problem. The idea is to reformulate Problem 1 in terms of h=Rf Then smoothness constraint ∥∇f∥ 2 =&lt;f, −Δf&gt; becomes &lt;R −1 h, −ΔR −1 &gt;=&lt;h, −ΔR −1* ΔR −1 . This constraint can be simplified using the following Lemma, which exploits the “intertwining” Fact 3 to analytically “shuffle” smoothness from the image domain to the Radon domain. Lemma 1 (Decoupling).  
           -          -     1   *              Δ                        -   1         =       1     4      π              ℐ     -   3       .                             
 
         [0043]    The proof begins by observing that R −1*  equals (4 π) −1 I −1 R, using Fact 1 and the symmetry of I α  But then −R −1*  ΔR −1  equals −(4π) −1RΔR   −1 =−(4π) −1 I −1 □RR −1 , using Fact 3. Lemma 1 follows using Fact 2 and because RR −1 s the identity operator.  
         [0044]    Thus Problem 1 can be posed in an equivalent, “easy” form, as follows:  
         [0045]    Problem 2 (Decoupled Regularization). Given observed projections g, find projections h=Rf minimizing  
                  g   -   h          w   2     +       β   ′          (     h   ,       ℐ     -   3          h       )         ,     where                   β   ′          :          =     β     4                 π       .                             
 
         [0046]    The corresponding Euler-Lagrange equation,  
           Wh+β′I   −3   h=WG   (2)  
         [0047]    where h is unknown, is easy exactly where equation (I) is hard. First, the forward and backprojections are eliminated from the optimization; backprojection need only be done once to determine f from solution h. Second, equation (2) is really a decoupled set of systems of equations, where each system corresponds to the unknowns {h(t, θ), for all t}, at each fixed 0. This follows because operator W is pointwise multiplied by a scalar and I −3  acts only along t for each fixed θ (using the definition of the Riesz potential). Thus for each fixed θ, we have an integral equation in the unknown single variable function h=Rf. Because h Rf, we can find f using filtered backprojection. Thus my approach to processing CT data is to solve equation (2).  
         [0048]    The decoupled regularization equation (2) is related to direct algebraic reconstruction tomography (DART) by emphasizing optimization entirely in the Radon domain. However, DART was not developed to deal with the Poisson noise model considered in many statistical approaches to tomography, and in effect assumed an invariance constraint regarding the weights w(t, θ), which is often not the case in real images. Thus decoupled regularization gains the benefits of space-varying filtering while maintaining the simplicity of DART. Although the proposed approach does not impose an image support constraint, the results obtained show strong improvement over the prior art.  
         [0049]    To solve equation (2), W+β′I −3  is a positive definite operator (if W&gt;θ), and thus the conjugate gradient method can be applied. To discretize the equations, I sampled in t and θ uniformly. The operator W was implemented by restriction to the sample locations. The Riesz potential was implemented by taking 1-dimensional FFTS. For greater fidelity to the Poisson noise model, one can use the said conjugate gradient method as a linear solver in an inexact Newton&#39;s method.  
         [0050]    Turning now to FIG. 3, after the object  22  (e.g. patient) is positioned within gantry opening  48 , the method illustrated in FIG. 3 may be carried out under control of the computer  36 . In FIG. 3, at step  52 , the object  22  is scanned to obtain projection data. At step  54 , the projection data is decomposed by slice and by projection angle. Thereafter, at step  56 , for each slice, equation (2) above is solved. Any known technique for solving linear equations, including iterative techniques, may be used to solve equation (2). Thereafter, at step  58 , an image of the slice is reconstructed from the processed projection angle data. The reconstructed image may be displayed at step  60  or may be stored for future use, analysis, or transmission to another location. The method of FIG. 3 (e.g. steps  54 ,  56 ,  58  and  60 ) is preferably implemented in software which may be carried as a set of instructions in a computer readable medium ( 50  in FIG. 2) which may take any known form.  
         [0051]    An application of the present invention is in medical tomographic imaging. PET scans, for example, are very noisy when processed using the standard filtered backprojection (FBP) technique, and therefore, computationally expensive regularization is applied. X-ray CT images can be badly degraded by the presence of metal in the scanned patient, giving rise to offensive streaking artifacts. Metal is often present in dental CT scans (in tooth fillings), in hip scans after hip replacement surgery, in gun shot victims, and in luggage for security inspection at an airport. Bone and other dense materials also create streaking artifacts in CT images as well. The disclosed method allows for much more rapid processing of these degraded images to significantly improve their quality as illustrated in FIGS.  4 - 7 .  
         [0052]    [0052]FIG. 4 illustrates streaking artifacts due to the presence of metal in an x-ray CT slice of a hip. For this sample computation, I cropped out an image region around the metal component to produce the cropped image shown in FIG. 7A. I then simulated the projection data (which were unavailable) by taking the (numerical) Radon transform in Matlab (see FIG. 6A). The bright band results from the metal, and is noisier than elsewhere; this is more evident in the single projection in the solid curve of FIG. 5. For each fixed projection orientation θ, the decoupled regularization equation (2) was solved to produce the nonhomogenously smoothed projection shown in FIG. 5. After this smoothing of each projection independently, a set of processed projections is obtained as shown in FIG. 6B. FIG. 6C illustrates the smoothing being localized on the metal band as a result of performing T* (a kind of differentiation and the first step in filtered backprojection) on the difference g−h. Using filtered backprojection, the final result is obtained as shown in FIG. 7B, and shows reduced streaking artifacts.  
         [0053]    By studying regularized tomography in the continuous domain, I was able to decouple a linear equation in a two-variable function into a one-parameter family of linear equations in single-variable functions. Although I emphasized a quadratic approximation of the x-ray CT likelihood, my result is based upon Lemma 1, which is independent of the imaging modality. For example, one can set up decoupled nonlinear equations to fully capture the Poisson likelihood. This technique can also be applied in three dimensions by noting that Facts 1-3 can be extended to higher dimensions. The gradient smoothness term can be replaced with related higher-order derivative penalties while maintaining the decoupling of the regularization; or one could combine such smoothness or other image constraints provided essentially that one can obtain the analogue of the above intertwining of operators for the constraint in question. The standard regularization approach, which the present invention replaces, implicitly ensures that the inferred image will have nonzero values only in a limited region which includes the patient. This is known as a consistency constraint. The present invention does not specifically enforce this constraint, as doing so would recouple the systems of equations. Not only are these equations smaller and thus easier to solve than standard tomographic equations, but their solutions for different angles can be computed in parallel on separate processors with no communication beyond distributing the workloads.  
         [0054]    This method readily applies in more general situations as well. This method applies to other tomographic imaging modalities such as PET, by altering the noise model. Essentially one must substitute the first term in Problem 1 with the appropriate negative log-likelihood for the modality in question. The realism of the model even with CT can be improved as well, while maintaining the benefits of decoupling. Specifically, the weight matrix need not be completely decoupled; instead, it could allow for intra-projection couplings. In addition, to allow for the reality that real scanners do not exactly follow the Radon transform, one can replace R with BR in Problem 1, where B is a linear operator which blurs within each projection. This captures the fact that x-rays are not always in perfect focus.  
         [0055]    One could also abandon the decoupling aspect and only take advantage of the transformation of the image smoothness or regularity constraint to the projection domain. This would allow more general blur matrices B and image smoothness constraints, at the cost of having to solve a fully coupled equation; the advantage of avoiding repeated forward and backprojections would be maintained.  
         [0056]    While the present invention has been described in connection with exemplary embodiments thereof, those of ordinary skill in the art will recognize that many modifications and variations are possible. The present invention is intended to cover such modifications and variations and is to be limited only by the scope of the following claims.