Computer Axial Tomography (CAT), sometimes known generally as Computerized Tomography (CT), is used in many applications, especially medical radiology, to obtain two or three dimensional views of the interior of three dimensional bodies (CT or CAT Scans). The technique involves subjecting a three dimensional body to radiation that enters the body from many different angles. The amount of radiation that is scattered or reflected by the body is then detected as a function of the angle of scattering. The scattered data is then analyzed to construct an image of the interior of the body. A two dimensional “slice” of the interior can be “reconstructed”, for example on a screen, and viewed. The slice can be reconstructed for any desired angle of intersection with the body.
While computerized tomography is well known, various specific mathematical reconstruction algorithms have been proposed to construct the image from the scattered radiation. However, it is becoming more challenging for known conventional reconstruction methods to meet the stringent constraints of current imaging applications. For example, the rates at which the impinging radiation beam scans the body has increased dramatically over the years, and the impinging radiation dosage has dropped significantly, especially in medical applications, because of patient safety concerns.
Two of the most common reconstruction algorithms that were developed to meet the stringent constraints of current imaging applications are the Filtered Back Projection method (FBP, sometimes written Filtered Backprojection method) and the Iterative Reconstruction Method (IR). In FBP, as a radiation detector rotates around the body being exposed to radiation, it detects and creates a series of planar images. Each of these images is a projection of the entire body mass for a particular angle, much as a common medical X-ray image. At each angle that the camera makes with the body, only radiation moving perpendicular to the camera is detected. As much of this radiation originates from various depths in the body, the result is an overlapping of the images of all elements of the body along the specific path, much in the same manner that a medical X-ray radiograph is a superposition of all anatomical structures from three dimensions into two dimensions.
An FBP study consists of many of such planar images acquired at various angles. After all the projections are acquired, they are subdivided by taking all the projections for a single, thin slice of the body at a time. All the projections for each slice are then ordered into an image called a sinogram. It represents the projection of the radiation distribution in the body for that single slice for every angle of the acquisition. The aim of the process is to retrieve the spatial distribution of the radiation from the projection data, thereby obtaining a clearer image of the body.
The quality of the image is further improved by filtering the data, which is most easily performed in the frequency domain rather than the spatial domain using well know Fast Fourier Transform techniques (FFT). The filtered data is than transformed back to the spatial domain.
The main reconstruction step involves a process known as back projection. As the original data was collected by only allowing radiation emitted perpendicular to the camera face to enter the camera, back projection smears the data from the filtered sinogram back along the same lines from where the radiation was emitted. Regions where back projection lines from different angles intersect represent areas which are of particular interest.
FBP and FFT-based reconstruction methods are popular in commercial CT scanners for their efficiency and accurate results when enough projection data is available. However, to lower the X-ray dosage (i.e., in low-dose CT imaging), the signal-to-noise ratio and the number of projections are reduced, which often results in undesirable “streaking artifacts”. There are typically two main sources of streaking artifacts First, the filter may amplify the imaging noise in the sinogram, which often produces streaking lines during back projection. Second, strong sharp edges, caused by, e.g., metal or bones, are back projected along the projection direction, causing streaking artifacts, especially when the number of projections is limited. Streaking noise is signal dependent and spatial variant, often making it extremely difficult to model this artifact accurately and distinguish the noise from true edges in the image. In the field of medical radiology, this typically results in difficult and/or inconsistent clinical interpretation(s).
The second reconstruction technique mentioned above, the Iterative Reconstruction (IR) method, can often produce a better reconstruction when FBP does not yield an acceptable solution. There are a large variety of IR algorithms, but each starts with an assumed image, computes projections from the image, compares the original projection data and updates the image based upon the difference between the calculated and the actual projections. Although conceptually this approach is typically simpler than FBP, for medical applications, it has traditionally lacked the speed of implementation and accuracy. This is due to the slow convergence of the algorithm and high computational demands. For these reasons, it was superseded by the FBP method in the early development of CT.
Therefore, there remains a need to more efficiently and accurately reconstruct an image using FBP, especially in low-dose CT.
In accordance with an embodiment of the present invention, an FBP method, referred to below as Compounded FBP (CFBP), has been developed to improve the reconstruction quality for low-dose CT scans where signal-to-noise (SNR) ratio is low and the number of projections may be restricted. The inventive method exploits complementary information from projection subsets to discriminate between angle dependent streaking noise and angle invariant true edges. This improves the image quality of the filtered back projection. The CFBP method produces image reconstruction of a quality that is typically similar to the quality of images produced by IR methods, but without the computational disadvantages of IR.
In one aspect of the invention, a method and system for constructing, from a computerized tomography (CT) scans an image relating to a physical structure involves obtaining projection data associated with the physical structure. The projection data is then divided into a plurality of subsets. FBP is then applied to each subset in the plurality of subsets. In one embodiment, an estimate of streaking noise for each subset is obtained. The image is then constructed after application of the FBP to each subset in the plurality of subsets. The estimate of streaking noise for each subset can be compared with at least some of the image after construction. In one embodiment, dependency between estimates of streaking noise for different subsets are analyzed. An independence constraint between subsets in the plurality of subsets can also be regularized.
These and other advantages of the invention will be apparent to those of ordinary skill in the art by reference to the following detailed description and the accompanying drawings
a) is a prior art projection image;
b) shows a prior art back projection image;
a) is an image of FBP reconstruction from all projections;
b) shows an image illustrating a regularization term in the first iteration in accordance with an embodiment of the present invention;
c) shows an image of the regularized CFBP reconstruction result in accordance with an embodiment of the present invention;
a) is a reconstructed image from all projections in accordance with an embodiment of the invention;
b) is an enlarged image of the FBP result with streaking noise;
c) is a regularized CFBP reconstruction result in accordance with an embodiment of the present invention;
d) is an enlarged image of the CFBP reconstruction result in accordance with an embodiment of the present invention;
a) is an original image for reconstruction;
b) is an FBP reconstruction image from all projections in accordance with an embodiment of the present invention; and
c) is a CFBP reconstruction result in accordance with an embodiment of the present invention.
In one aspect of the present invention, a new FBP method results in improved reconstruction quality. The invention may be used beneficially for low dose CT scans where signal-to-noise (SNR) ratio is low and the number of projections may be restricted. As noted, this improved method is referred to herein as Compounded FBP (CFBP). The inventive method exploits complementary information from projection subsets to discriminate between angle dependent streaking noise and angle invariant true edges. To aid in understanding the invention, we will discuss below: 1. Filtered Back Projection (FBP); 2. Compounded Filtered Back Projection (CFBP); 3. CFBP Background; 4. FBP Derived from MAP Framework; and 5. CFBP Based on Cross Group Regularization.
1. Filtered Back Projection
As discussed above, FBP involves the construction of an image by analyzing projection images obtained from many different viewing angles. The technique is shown schematically in
For parallel beam tomography, the projections can be expressed as a Radon transform of the object that is to be reconstructed. The Radon transform is defined as:
g(p, θ)=R(f)=∫−∞∞∫−∞∞f(x, y)δ(x cos θ+y sin θ−p)dxdy (2)
i.e., the line integral along a line (a tomography beam) at an angle θ from the y-axis and at a distance |s| from the origin. By rotating the coordinate system,
p=x cos θ+y sin θ
u=−x sin θ+y cos θ
g(p,θ) can be expressed as
g(p,θ)=∫−∞∞f(p cos θ−u sin θ, p sin θ+u cos θ)du (3)
In the context of tomography, the Radon transform is often called a sinogram because the Radon transform of a delta function is the characteristic function of the graph of a sine wave. Consequently, the Radon transform of a number of small objects appears graphically as a number of blurred sine waves with different amplitudes and phases.
The sinogram is a collection of raw projection data that appear as a set of sinusoidal lines and bands. A sinogram typically contains information that is not available from its reconstructed image. For instance, a full CT scanning field of view is recorded in a sinogram, whereas a small subsection or region of interest of the full field of view is usually reconstructed in an image.
In accordance with an embodiment of the present invention, compounded filtered back projection (CFBP) divides projection data into several subsets. Each subset reconstructs the same structure from a slightly different angle relative to FBP. Conventional FBP typically uses the average of these subsets.
As described in more detail below, if we define the reconstruction process as finding the most probable estimation of the true object structure being scanned, FBP can be derived from Maximum a posteriors (MAP) (i.e., an energy minimization) framework based on statistical noise/signal models. Restoration from conventional MAP does not, however, always guarantee optimal results. Better results may often be achieved by fully exploiting the complementary information provided by multiple images. In FBP, the results from different subsets of projections typically produce streaking artifacts at different orientations while the true object edges remain static. By compounding the subsets, the results of the FBP can be improved by reducing the streaking artifacts while preserving the structures (i.e., the edges). Thus, quality similar to iterative methods can effectively be obtained.
2. Compounded Filtered Back Projection
In CT imaging, cross-sectional images I (x,y) are reconstructed from their X-ray projection measurements. Without loss of generality, the projection measurements (i.e., data) (referred to above as a sinogram (i.e., g(p,θ))), are obtained by integrals along parallel rays, where θ is the angular orientation of the projection rays and p is the distance of the ray to the origin. The image reconstruction is equivalent to solving the inverse Radon transform based on the given sinogram.
Filtered back projection (FBP) is derived from signal processing theory. In one embodiment, the image can be reconstructed by first filtering each projection with a ramp filter kernel and then back projecting the P filtered projections (i.e., gf(p,θ)) into the 2D image space, as shown in the following equation.
As described above, for low-dose CT imaging, where the signal-to-noise ratio is low and the number of projection is restricted, FBP reconstruction typically suffers from streaking artifacts (i.e., noise).
3. Compounded Filtered Back Projection (CFBP) Background
According to an embodiment of the present invention, a total of P projection measurements are down-sampled by a factor of M, forming M projection subsets. The conventional FBP can be reformulated using the linearity of the Radon transform.
where Im is the FBP reconstruction from the mth subset of the projections.
Since each subset of projections contains slightly different viewing angles, the streaking noise in the FBP results (i.e., Im) for each subset is slightly rotated as well. Hence, the different Im are reconstructions reflecting the same underlying objects, but corrupted by differently oriented streaking noise. Therefore, the reconstruction problem can be viewed as a multi-image based restoration problem.
Unlike CFBP, conventional FBP sums all Im, which is equivalent to averaging (e.g., up to a normalization factor). While the averaging process may preserve the underlying objects and reduce independent streaking noise, the reconstruction quality under low-dose imaging may still be insufficient for accurate diagnostic interpretation.
4. FBP Derived from MAP Framework
Assume that the projections are divided into two subsets (e.g., odd and even projections). A Bayesian approach (i.e., an approach that uses a statistical estimator resulting with a penalty functional that is minimized by a numerical optimization algorithm that yields the restored signal) may be used to recover the true signal based on the two sub-reconstructions (i.e., I1, I2) which reflect the same underlying signal but are corrupted by independent streaking noise N1 and N2. The problem can be formulated as follows:
P(Î|I1,I2))=c·P(I1,I2|Î)P(Î)=c·PN
where c is a normalization constant and PN
The cost function (i.e., the number of projections require to generate the final reconstruction) guarantees that the final reconstruction looks like the sub-reconstructions. The MAP estimation can be obtained by setting the derivative of the objection function to zero.
The MAP result is the average of the sub-reconstructions, which is the same as the traditional FBP reconstruction (as shown in Eq. 2 above).
The MAP framework does not always guarantee the optimal results. For example, the MAP framework can sometimes generate results that violate some underlying constraints, especially the independency between the streaking noise terms (i.e., the streaking noise on different sub-reconstructions does not correlate with each other due to different back projections angles). With respect to the above description, the final reconstruction is estimated as Î=(I1+I2)/2. If Î is the correct reconstruction, the streaking noise on the sub-reconstruction is
{circumflex over (N)}1=I1−Î=(I1−I2)/2=−{circumflex over (N)}2 (6)
As described above, the streaking noise on different sub-reconstructions does not correlate with each other due to different back projection angles. But if the average of the sub-reconstructions is used in the final reconstruction (i.e., following the traditional FBP method), the streaking noise becomes fully correlated (negative to each other). The violation of the independency constraint indicates reconstruction errors in the FBP method.
5. CFBP Based on Cross Group Regularization
By dividing the projection into sub-groups and applying FBP separately, an estimate of the streaking noise {circumflex over (N)}i in each sub-reconstruction Ii can be obtained. This estimate can further be compared with the final reconstruction result Î. Combined with the fact that each sub-reconstruction is the true reconstruction (i.e., I) corrupted with noise (i.e., Ii=I+Ni), the following equation results:
{circumflex over (N)}i=Ii−Î=Ni+(I−Î) (7)
If the reconstruction Î is perfect (i.e., Î=I), the estimated streaking noise {circumflex over (N)}i is also perfect and is independent across the sub-reconstructions. However, if a reconstruction error exists, the reconstruction error is typically equally reflected in all estimated streaking noises and can cause correlation, such as in Eq. (6). Therefore, analyzing the dependency between the estimated streaking noises can help identify and correct reconstruction errors.
In accordance with an embodiment of the present invention, the compounded filtered back projection (CFBP) method can explicitly enforce independence between the estimated streaking noise on each sub-reconstruction.
In one embodiment, evaluating the strict independency (i.e., p({circumflex over (N)}1, {circumflex over (N)}2)=p({circumflex over (N)}1)p({circumflex over (N)}2)) is very expensive in terms of processing time (i.e., when evaluating p(N1, N2, . . . Nk), joint histogram needs to be estimated, which increases exponentially when the number of images increases). Therefore, one of the properties of the independent random variables is instead relied upon:
E(h1({circumflex over (N)}1)h2({circumflex over (N)}2))=E(h1({circumflex over (N)}1))E(h2({circumflex over (N)}2)) (8)
where E is the expectation, and h1( ) and h2( ) are functions of {circumflex over (N)}1 and {circumflex over (N)}2, respectively. This equation holds for arbitrary h1( ) and h2( ), as long as the random variables {circumflex over (N)}1 and {circumflex over (N)}2 are independent.
For efficient computation, the functions are chosen to be a weighted sum of different order moments. For example, if h1({circumflex over (N)}1)={circumflex over (N)}1 and h2({circumflex over (N)}2)={circumflex over (N)}2 are chosen, the independent constraint to the uncorrelated constraint is reduced.
When the random variables are Gaussian distributed, independence and uncorrelated are equivalent. However, {circumflex over (N)}1 and {circumflex over (N)}2 typically do not follow a simple distribution because of the existence of restoration errors in Eq. (7). Combining higher order moments can better approximate the independence. One advantage of this approximation, based on order statistics, is that the expectations using simple summation operations can be used, which reduces the computational cost.
Therefore, a new energy term can be included for the cross group regularization to integrate with the previous MAP objective function defined in Eq. (4). The new energy term that regularizes the independence constraint between the residuals ({circumflex over (N)}i=Ii−Î and {circumflex over (N)}j=Ij−Î) of sub-reconstruction i and j can be defined as:
ei,j({circumflex over (N)}i, {circumflex over (N)}j)=∥E(h1({circumflex over (N)}1)h2({circumflex over (N)}j))−E(h1({circumflex over (N)}i))E(h2({circumflex over (N)}j))∥2 (9)
This regularization term is also extendable based on the available computing power and the expected accuracy of the independence evaluation. For more accurate evaluation of the independence, several terms can be added based on different choices of h1( ) and h2( ).
If the projections are divided into more sub-groups, the sum of the pairwise independence can be used to approximate the joint independence across the sub-groups. With the new cross groups regularization terms, the previous MAP framework for the reconstruction can be adapted into the following new form:
The optimal solution for this objective function is not FBP any more. It can be found by setting the derivative of the objective function defined in Eq. (10) to zero. Since the new regularization terms ei,j is dependent on the final reconstruction result I, the optimal solution is found through several iterations.
The iteration is initialized with a loss-pass filtered FBP result Î0 (i.e., the average of all sub-group reconstructions followed by a low-pass filtering). In one embodiment, two types of regions are inspected. First, for the smooth region where there are no sharp structures, the low-pass filtered FBP result provides streaking, noise-free reconstruction. The cross group regularization terms are almost at their minimum (i.e., estimated streaking noise is independent) and there is no need to refine these regions.
Second, for regions that have sharp boundaries and streaking noise, low-pass filtering removes the streaking noise and also blurs the sharp boundaries. The blurred boundaries are restoration errors and typically increase the cross group regularization cost. In one embodiment, to reduce the cost, the common part in all the sub-reconstructions are preserved while the uncorrelated part in the final reconstruction are suppressed.
The description herewith describes the present invention in terms of the processing steps required to implement an embodiment of the invention. These steps may be performed by an appropriately programmed computer, the configuration of which is well known in the art. An appropriate computer may be implemented, for example, using well known computer processors, memory units, storage devices, computer software, and other components. A high level block diagram of such a computer is shown in
One skilled in the art will recognize that an implementation of an actual computer will contain other components as well, and that
a)-4(c) illustrate FBP reconstruction in accordance with an embodiment of the present invention. In particular,
In one embodiment, the approximated independence evaluation based on order statistics and summation operation reduces the computational cost of the regularization term. This allows efficient minimization of the objective function iteratively.
a) and 5(b) illustrate an entire reconstructed image as well as an enlarged part of the FBP result with streaking noise. In one embodiment, the compounded filtered back projection method is applied to several low dose CT scans. For example, in a first data set, a sinogram with 100 projections and 367 detectors per projection is used. The dimensions of the reconstructed image are 258 by 258. Filtered back projection of the original sinogram yields a reconstruction with apparent streaking noise as shown in
The CFBP method is then applied, as shown in Eq. (10), with cross group regularization constraints. The projections are divided into two subsets and two images are reconstructed using FBP on each of the subsets. To initialize the iterative optimization scheme, the initial reconstruction result is assigned to be the averaging of the two sub-group FBP results followed by uniform low-pass filtering. The streaking noise is reduced by the low-pass filtering, which is equivalent to the FBP modulated by a strong low-pass windowing function. The sharp boundaries, however, may also be blurred.
In accordance with an embodiment of the present invention, the regularization term of the CFBP method captures the reconstruction error around the sharp boundaries in the initial guess. This is shown in
The regularization term is estimated based on the pixels within a neighborhood of each pixel location and hence cover a (e.g., slightly) larger region than the true boundaries. This does not, however, typically affect the sharpness of the reconstruction result. The final reconstruction, shown in image 420 of
In another embodiment, the CFBP method is applied to another low-dose CT data set consisting of 580 projections and 1344 detectors per projection. The dimensions of the reconstructed image are 949 by 949.
a) illustrates the reconstructed image 504 from all projections.
A quantitative analysis can also be performed.
The foregoing Detailed Description is to be understood as being in every respect illustrative and exemplary, but not restrictive, and the scope of the invention disclosed herein is not to be determined from the Detailed Description, but rather from the claims as interpreted according to the full breadth permitted by the patent laws. It is to be understood that the embodiments shown and described herein are only illustrative of the principles of the present invention and that various modifications may be implemented by those skilled in the art without departing from the scope and spirit of the invention. Those skilled in the art could implement various other feature combinations without departing from the scope and spirit of the invention.
This application claims the benefit of U.S. Provisional Application No. 60/725,945, filed Oct. 12, 2005, which is incorporated herein by reference.
Number | Name | Date | Kind |
---|---|---|---|
4709333 | Crawford | Nov 1987 | A |
4729100 | Tsujii | Mar 1988 | A |
5416815 | Hsieh | May 1995 | A |
5438602 | Crawford et al. | Aug 1995 | A |
5933471 | Kalvin | Aug 1999 | A |
6028908 | Taguchi | Feb 2000 | A |
6035012 | Hsieh | Mar 2000 | A |
6038278 | Hsieh et al. | Mar 2000 | A |
6101236 | Wang et al. | Aug 2000 | A |
6282257 | Basu et al. | Aug 2001 | B1 |
6332035 | Basu et al. | Dec 2001 | B1 |
6351548 | Basu et al. | Feb 2002 | B1 |
6426988 | Yamada et al. | Jul 2002 | B2 |
6507633 | Elbakri et al. | Jan 2003 | B1 |
6721387 | Naidu et al. | Apr 2004 | B1 |
6754298 | Fessler | Jun 2004 | B2 |
6768782 | Hsieh et al. | Jul 2004 | B1 |
6775352 | Toth et al. | Aug 2004 | B2 |
6845142 | Ohishi | Jan 2005 | B2 |
6850586 | Cahill | Feb 2005 | B2 |
6950492 | Besson | Sep 2005 | B2 |
7023951 | Man | Apr 2006 | B2 |
7113569 | Okumura et al. | Sep 2006 | B2 |
20030161443 | Xiao et al. | Aug 2003 | A1 |
20040264625 | Basu et al. | Dec 2004 | A1 |
20050058240 | Claus | Mar 2005 | A1 |
20050105693 | Zhao et al. | May 2005 | A1 |
20060115040 | Chen | Jun 2006 | A1 |
Number | Date | Country | |
---|---|---|---|
20070140407 A1 | Jun 2007 | US |
Number | Date | Country | |
---|---|---|---|
60725945 | Oct 2005 | US |