Methods and apparatus for correcting scatter

[0003] This invention relates generally to X-ray image processing and more particularly to methods and apparatus for performing a scatter correction algorithm.

[0004] At least some known imaging systems use a moving grid of parallel metal plates to facilitate a scatter-signal suppression. The grid rejects a majority of photons that are not essentially parallel to a beam emanating from an X-ray source in a direction of a detector, based on the scatter grid geometry. Since many scattered photons are not parallel to the X-ray beam, a significant proportion of scatter events are absorbed by the grid. In addition to absorbing scatter events, the grid also absorbs a significant portion of the non-scattered, or primary photons. When a scatter grid is used, approximately 20% of the remaining signal is scatter, and approximately 30% of the primary events are rejected. Scatter grid assemblies are also expensive and mechanically complex.

[0005] In one aspect, a method for removing scatter in an image is provided. The method includes acquiring data of an object of interest, and using an iterative equation including a thickness-dependent kernel modulation factor to reconstruct an image of the object.

[0006] In another aspect, a method for removing scatter in an image is provided. The method includes acquiring data of an object of interest, and using an iterative equation to reconstruct an image of the object when a scatter fraction greater than approximately 0.5.

[0007] In a further aspect, a computer readable medium encoded with a program executable by a computer for removing scatter in an image is provided. The program is configured to instruct the computer to acquire data of an object of interest, and use an iterative equation including a thickness-dependent kernel modulation factor to reconstruct an image of the object.

[0008] In a still another aspect, a medical imaging system for removing scatter in an image is provided. The medical imaging system includes a detector array, at least one radiation source, and a computer coupled to said detector array and radiation source and configured to acquire data of an object of interest, and use an iterative equation to reconstruct an image of the object when a scatter fraction is greater than approximately 0.5.

[0009] FIG. 1 is a pictorial view of an imaging system.

[0010] FIG. 2 is a flow diagram of a method for removing scatter.

[0011] FIG. 3 is an illustration of scatter distribution using a typical imaging system shown.

[0012] FIG. 4 is an image acquired with a slit collimator using a CIRS breast phantom.

[0013] FIG. 5 is a detector profile underneath the slit collimator without correcting for kernel asymmetry.

[0014] FIG. 6 is a detector profile underneath the slit collimator correcting for kernel asymmetry.

[0015] FIG. 7 illustrates the performance of the scatter-correction algorithm described herein.

[0016] FIG. 8 illustrates the performance of the scatter-correction algorithm described herein.

[0017] FIG. 9 illustrates the performance of the scatter-correction algorithm described herein.

[0018] FIG. 1 is a pictorial view of an imaging system 10, such as a X-ray mammography imaging system and a tomosynthesis imaging system. System 10 includes at least one radiation source 14, such as an X-ray source, and at least one detector array 16 for collecting views from a single or a plurality of projection angles 18. Specifically, system 10 includes a radiation source 14 which projects a cone-shaped beam of X-rays which pass through object 12 and impinge on detector array 16. Detector array 16 is fabricated in a panel configuration having a plurality of pixels (not shown) arranged in rows and columns such that an image is generated for an entire object 12 of interest. In one embodiment, for example for imaging of the lung, detector array 16 is a chest detector array 16 and object 12 is a patient's chest 12. In an exemplary embodiment, imaging system 10 is used to generate a three-dimensional dataset representative of an imaged object 12. The views obtained at each angle 18 may be used to reconstruct a plurality of slices, i.e., images representative of structures located in planes 20 which are parallel to detector 16. In another embodiment, only a single projection radiograph of the object of interest is acquired.

[0019] Each pixel includes a photosensor, such as a photodiode (not shown), that is coupled via a switching transistor (not shown) to two separate address lines (not shown). In one embodiment the two lines are a scan line and a data line. The radiation incident on a scintillator material and the pixel photosensors measure, by way of change in the charge across the diode, an amount of light generated by X-ray interaction with the scintillator. More specifically, each pixel produces an electronic signal that represents an intensity and/or scatter, after attenuation by object 12, of an X-ray beam impinging on detector array 16. In one embodiment, detector array 16 is approximately 20 centimeters (cm) by 20 cm and is configured to produce views for an entire object 12 of interest. Alternatively, detector array 16 is variably sized depending on the intended use. Additionally, a size of the individual pixels on detector array 16 is selected based on the intended use of detector array 16.

[0020] In the exemplary embodiment, alternative detector technology is used, such that views in digital form are generated by detector 16. In another embodiment, the reconstructed three-dimensional dataset is not arranged in slices corresponding to planes that are parallel to detector 16, but in a more general fashion. In yet another embodiment, detector 16 is a shape other than planar.

[0021] In the exemplary embodiment, radiation source 14 and detector array 16 are moveable relative to the object 12 and each other. More specifically, radiation source 14 and detector array 16 are movable such that the projection angle 18 of the imaged volume is altered. Radiation source 14 and detector array 16 are movable such that projection angle 18 may be any acute or oblique projection angle.

[0022] The operation of radiation source 14 is governed by a control mechanism 28 of imaging system 10. Control mechanism 28 includes a radiation controller 30 that provides power and timing signals to radiation source 14 and a motor controller 32 that controls a respective movement speed and position of radiation source 14 and detector array 16. A data acquisition system (DAS) 34 in control mechanism 28 samples digital data from detector 16 for subsequent processing. A portion of the data acquired includes a compressed breast thickness, which is included along with a detector cover thickness and a compression paddle thickness, which are then used as parameters in the scatter correction process. The measured image and breast thickness information are processed to form a “breast thickness map” which is included in the scatter correction process. An image processor 36 receives a sampled and digitized projection dataset from DAS 34 and performs image processing and scatter correction, as described herein. In one embodiment, image processor 36 also performs 3D reconstruction of imaged object 12 using scatter-corrected projection images. The scatter-corrected images, and/or the reconstructed three-dimensional dataset, representative of imaged object 12, is applied as an input to a computer 38 which stores the datasets in a mass storage device 40. Image processor 36 is programmed to perform functions described herein, and, as used herein, the term image processor refers to computers, processors, microcontrollers, microcomputers, programmable logic controllers, application specific integrated circuits, and other programmable circuits.

[0023] Computer 38 also receives commands and scanning parameters from an operator via console 42 that has an input device. A display 44, such as a cathode ray tube, and a liquid crystal display (LCD), allows the operator to observe the reconstructed three-dimensional dataset and other data from computer 38. The operator supplied commands and parameters are used by computer 38 to provide control signals and information to DAS 34, motor controller 32, and radiation controller 30. In one embodiment, the compressed breast thickness is not acquired by system 10, but entered by an operator via console 42.

[0024] In one embodiment, a patient is positioned such that the object of interest is within the field of view of system 10, i.e., object 12 is positioned within the imaged volume extending between radiation source 14 and detector array 16. Views of object 12, are then acquired from at least two projection angles 18 to generate a projection dataset of the volume of interest. The plurality-of-views represent the tomosynthesis projection dataset. The collected projection dataset is then utilized to generate a three-dimensional dataset, i.e., a plurality of slices for scanned object 12, such as a breast 12, wherein the three-dimensional dataset is representative of the three-dimensional radiographic representation of imaged breast 12. After enabling radiation source 14 so that the radiation beam is emitted at first projection angle 50, a view is collected using detector array 16. Projection angle 18 of system 10 is then altered by translating the position of source 14 so that central axis 48 of the radiation beam is altered to a second projection angle 52 and position of detector array 16 is altered such that breast 12 remains within the field of view of system 10. Radiation source 14 is again enabled and a view is collected for second projection angle 52. The same procedure is then repeated for any number of subsequent projection angles 18.

[0025] Some known imaging systems 10 rely on a differential absorption of X-rays by a plurality of tissue types. For example, X-ray mammography relies on differential absorption of tissue types found in a human breast 12. X-ray photons that are directly transmitted (i.e. primary photons) through object 12 are detected by a film, or by digital detector 16.

[0026] In addition to interacting with matter through absorption, X-rays also undergo scattering with the electronic constituents of object 12 being imaged and other components of X-ray imaging system 10. For example, in X-ray mammography, a majority of photon scatter occurs in breast 12, a compression paddle (not shown) and a plurality of detector materials (not shown). Scattered photons produce a low-frequency “haze” superimposed on the image produced by primary X-ray photons due partially to their nearly isotropic scattering cross section. A signal from scattered photons contributes to the overall noise of the image and contributes to a decrease in an image contrast. In some known imaging systems 10, approximately 40% of the signal in a typical mammogram, i.e. using a breast 12 approximately 4 centimeters (cm) in compressed breast thickness, is composed of scattered events, if a scatter grid (not shown) is not used. For breasts 12 thicker than approximately 6 cm, the scatter fraction (scatter/total signal) may be over 0.5.

[0027] In one embodiment, an acquired x-ray image is processed such that the contribution to the image due to scattered x-ray photons is essentially removed. In yet another embodiment, the scatter-corrected images may be used to quantitatively link the image signal to radiographic characteristics of the constituent tissue components of the imaged object.

[0028] The scatter-correction algorithm described herein facilitates eliminating the need for using the scatter grid in X-ray mammographic imaging. Further, scatter-corrected images according to the present invention show a mean absolute deviation of approximately 3% with respect to scatter-free images, even in the boundary region of breast 12. Therefore, resulting in superior contrast with respect to a Smit-Roentgen scatter grid mammography result, wherein approximately 20% of the remaining signal is due to scatter. In one embodiment, the algorithm described herein may be applied to scatter correction of tomosynthesis projection datasets where a standard scatter grid cannot be used. Further, the scatter correction algorithm facilitates scatter-correction in other X-ray applications, such as, but not limited to, a tissue density decomposition for a plurality of two-dimensional (2D) and three-dimensional (3D) applications.

[0029] In one embodiment, an image y is formed by combining primary events and scatter events. Primary events are photons which travel directly from source 14 to detector 16, and scatter events are photons that scatter from object 12 being imaged or any other object in the path of the X-ray beam into detector 16. An image y may be defined according to:

[0030] y=b+s,  Equation 1

[0031] where b is an image formed by unattenuated photons, i.e. primary events, and s is an image formed by photons scattered into detector 16, i.e. scatter events.

[0032] To recover a direct image b, term s may be modeled as a convolution between a scatter kernel p and the direct image b, according to:

[0033] s=p*b,  Equation 2

[0034] where * is a convolution operator.

[0035] Using Equations 1 and 2, image y may be defined according to:

[0036] y=b+p*b  Equation 3

[0037] For any given pixel j in image y, Equation 3 may be expressed as: $\begin{array}{cc}{y}_j=b_j+\sum _{k=1}^n{p}_{\mathrm{jk}}{b}_k& \mathrm{Equation}4\end{array}$

[0038] where an index k spans all pixels of detector 16, i.e. pixel 1 through pixel n, and a factor, or kernel, p_j,krelates the number of direct events detected in pixel j to a contribution of photons originally going to pixel k.

[0039] In one embodiment, the factor, or kernel p_j,kis spatially invariant such that p_j,k=p_j-k, where j-k includes distance and directional information. For a spatially invariant kernel, a convolution translates into a simple product in Fourier space, which may be numerically calculated through Fourier Transforms, such as, but not limited to, a Fast Fourier Transforms (FFT). For example, the second term in Equation 4 may be calculated through a FFT, if the kernel p_j,kis spatially invariant. In another embodiment, p_j,kis not spatially invariant, and the second term in Equation 4 is not calculated using a FFT. The assumption of a spatially-invariant scatter kernel is a good approximation in the case of scatter within a single slab or a plurality of slabs, wherein the slabs include a homogeneous material and a constant thickness, and extend approximately horizontally beyond the X-ray beam emanating from radiation source 14 and impinging on detector array 16, and wherein the slabs are resting on detector 16 surface with X-ray photons incident approximately perpendicular to the slab's surface. For the slab geometry, the scatter kernel is mainly dependent on the thickness of the material, and only weakly dependent on the material's composition for fatty through glandular-equivalent breast material, and on an X-ray energy spectrum for typical X-ray techniques used in mammography.

[0040] The assumption of space invariance in the kernel may introduce inaccuracies in the calculation of the scatter signal when there are large variations of thickness in object 12 being imaged, such as in mammography near the edge of breast 12. Also, there is a significant contribution of scatter events from the detector materials which are irradiated directly by the X-ray beam, i.e. a region outside the edges of breast 12 in a mammography image.

[0041] FIG. 2 is a flow diagram of a method 60 for removing scatter in an image. Method 60 includes acquiring 62 data of object 12 of interest, and using 64 an iterative equation including a thickness-dependent kernel modulation factor to reconstruct an image of the object.

[0042] In the exemplary embodiment, a thickness-dependent kernel modulation factor α_k, is defined as: $\begin{array}{cc}{a}_k=\frac{{\stackrel{\_}{p}}_k}{{\stackrel{\_}{p}}_{k_0}},& \mathrm{Equation}5\end{array}$

[0043] where {overscore (p)}_kis a norm of scatter kernel p, for each pixel k, of detector 16, which is a function of a thickness of material d_k, i.e. breast thickness map, traversed by the primary ray incident on the pixel k. Additionally, {overscore (p)}_jis the norm of the kernel for the compressed thickness of breast 12, i.e., the distance between a compression paddle and detector 16 surface, and any additional thickness of material due to detector 16 and the compression plate assembly d_k, which can be read out using imaging system 10 or estimated by an alternative means.

[0044] The image y, as defined in Equation 3, is re-defined according to:

[0045] y=b+p*(ab)  Equation 6

[0046] where the kernel p is calculated assuming a thickness d_k.

[0047] A scatter signal may then be defined according to: $\begin{array}{cc}{s}_j=\sum _k{p}_{\mathrm{jk}}{a}_k{b}_k,& \mathrm{Equation}7\end{array}$

[0048] where each contribution of the scatter events from pixel k to pixel j is re-normalized to account for a change in a detection probability of the scatter event based on a change of thickness of breast 12 material traversed by a ray k. In another embodiment, modulation factor α_kis also dependent on pixel j.

[0049] As described previously herein, the norm of the scatter kernel is sensitive to the thickness of material in a slab arrangement, but the shape, as characterized, e.g., by a full-width half-maximum or other parameters, of the scatter kernel is relatively insensitive to the thickness of the uniform slab within a thickness range of approximately 0.3 centimeters (cm) to approximately 10 cm. Therefore, the thickness dependent kernel modulation factor α_kfacilitates compensating for a majority of the change in detection probability of a scatter event due to a variable thickness of object 12, and enables the use of space-invariant scatter kernels for accurate scatter correction.

[0050] FIG. 3 is an illustration of scatter distribution using a typical imaging system 10. In one embodiment, the scatter kernel p_j,k, i.e., a space-invariant scatter kernel, is modeled by a radially symmetric function. Approximating scatter using a radially symmetric function is used when the scatter signal is calculated away from the boundaries of breast 12. For example, and referring to FIG. 3, a kernel asymmetry is introduced near breast 12 boundary by photons scattered outside breast 12 that travel through less breast material before being detected. Therefore, in case 1 p_j,k=p_j′k, but in case 2 p_j,k<p_j′kwhere in case 1 the distance traversed by a scattered X-ray photon through tissue is the same for pixels j and j′, but in case 2 the distance traveled through tissue to pixel j is longer than for pixel j′ (pixel j is inside the breast and pixel j′ is outside the breast). Therefore in case 2 the X-ray photon scattered into pixel j′ is attenuated less than in case 1. Taking this effect into account directly in a kernel calculation facilitates eliminating the spatial invariance of the kernel, but prohibits the use of FFTs in the convolution operation to calculate the scatter signal.

[0051] In the exemplary embodiment, a re-nonmalization map is based on the measured image, which is then processed to find the boundary of breast 12. The pixel distance away from the boundary of breast 12 is used to approximate the kernel asymmetry effect without introducing space-variance in the kernel. The re-normalization map N(m) is estimated for every pixel m outside of breast 12 boundary by calculating a distance d_b(m) between pixel m and a closest pixel belonging to breast 12 boundary.

[0052] The re-normalization map N(m) may be approximated according to:

[0053] N(m)˜e^μd(m),  Equation 8

[0054] where μ is a mean attenuation coefficient for the X-ray photon spectrum in 50% glandular/50% fat breast equivalent material. In one embodiment, N(mn)=1 for pixels inside breast 12, since most of the X-ray photons scattered within the breast region are, on average, subject to the same attenuation if they travel the same distance from a scattering point to where they are detected.

[0055] In use, the re-normalization map N(m), facilitates accounting for the increased probability of detection of an event scattered from breast 12 into detector 16 region outside of breast 12, due to a decrease of path-length through the breast attenuation medium, i.e. the breast tissue.

[0056] The re-normalization map is then incorporated in an iterative equation b⁽ⁿ⁾according to: $\begin{array}{cc}{b}^{\left(n\right)}=\frac{y}{2^l}+\left(\frac{2^l-1}{2^l}\right)b^{\left(n-1\right)}-\frac{Np*{\mathrm{ab}}_1^{\left(n-1\right)}}{2^l}-\frac{p*{\mathrm{ab}}_2^{\left(n-1\right)}}{2^l},& \mathrm{Equation}9\end{array}$

[0057] where the subscript 1 refers to direct events outside breast 12 boundary, subscript 2 refers to direct events inside breast 12 boundary, and l is an integer that satisfies the condition {overscore (p)}<2^l. In one embodiment, l is a quantity of sub-iterations.

[0058] FIG. 4 is an image acquired with a slit collimator using a breast phantom wherein the breast phantom is approximately 5 cm thick and includes approximately 70% glandular 30% fat composition. Data was acquired in open geometry, i.e. without using a slit collimator, and with a slit collimator directly on top of the phantom. The signal directly beneath the slit is almost exclusively due to direct events, i.e., the signal includes no contribution due to scatter. The open-geometry image was processed using the scatter-correction algorithm as described herein, with and without scatter-asymmetry correction.

[0059] FIG. 5 is a detector profile underneath the slit collimator without correcting for kernel asymmetry. FIG. 6 is a detector profile underneath the slit collimator correcting for kernel asymmetry. FIGS. 5 and 6 show the profiles for pixels underneath the slit collimator for open geometry, scatter-corrected open geometry, and scatter-free images, i.e. slit collimator in use. In use, a mean absolute deviation between scatter-corrected and scatter free images is approximately 3% using the kernel asymmetry correction, and approximately 4% without using the kernel asymmetry correction, for the profile up to breast 12 boundary. In breast 12 boundary region itself, where the intensity varies rapidly between channels, i.e., approximately 800 and approximately 1064 (region demarcated by the ellipse on FIG. 5), the mean absolute deviations were approximately 3% and approximately 5% respectively.

[0060] In another exemplary embodiment, a scatter-correction algorithm facilitates removing scatter from an X-ray image, regardless of a relative magnitude of the scatter signal with respect to the direct-events signal. A scatter fraction is defined according to: $\begin{array}{cc}{f}_s=\frac{s}{b+s}=\frac{p*b}{b+p*b},& \mathrm{Equation}10\end{array}$

[0061] wherein the scatter fraction f_srepresents a proportion of scatter events to total events throughout the image. For a given image y, a scatter fraction greater than 0.5 indicates that more scatter events are detected than direct events in a given pixel in the image y, and consequently the norm {overscore (p)} of the scatter kernel is greater than 1 for that pixel.

[0062] Assuming that p is known for a given image acquisition system 10 configuration and object 12 being imaged, and in one embodiment, the scatter kernel p is represented by a sum of Gaussian functions. An initial estimation of the scatter signal can defined as: $\begin{array}{cc}{s}^{\left(0\right)}=\frac{p*y}{1+\stackrel{\_}{p}}=\frac{p*\left(b+p*b\right)}{1+\stackrel{\_}{p}}=\frac{p*b*\left(\delta +p\right)}{1+\stackrel{\_}{p}},& \mathrm{Equation}11\end{array}$

[0063] where δ denotes the Kronecker delta function.

[0064] Equation 12 illustrates redefining a direct image according to:

[0065] b^(o)=y−s^colEquation 12

[0066] An iterative solution is then defined according to:

[0067] b⁽ⁿ⁾=y−p*b⁽ⁿ⁻¹⁾.  Equation 13

[0068] Equation 13 reduces to

[0069] b⁽ⁿ⁾=b−(−l)ⁿpⁿ*(b⁽⁰⁾−b).  Equation 14

[0070] For the case of a scatter fraction less than 0.5, i.e. less than 50% of the measured signal is scatter, the norm {overscore (p)} of the scatter kernel is less than 1.0 and Equation 14 converges toward the direct image b with increasing values of n. However, for {overscore (p)}≧1 Equation 14 diverges for increasing values of n and the iterative approach fails to converge toward the direct image b.

[0071] In the exemplary embodiment, an initial direct image estimation b is defined as:

[0072] b⁽⁰⁾=y−s⁽⁰⁾=b+p*b−s⁽⁰⁾.  Equation 15

[0073] A first interim solution i is calculated according to:

[0074] i⁽¹⁾⁽⁰⁾=y−p*b⁽⁰⁾,  Equation 16

[0075] First interim solution i is averaged with b⁽⁰⁾according to: $\begin{array}{cc}{i}^{\left(1\right)\left(1\right)}=\frac{i^{\left(1\right)\left(0\right)+b^{\left(0\right)}}}{2}& \mathrm{Equation}17\end{array}$

[0076] Equation 17 is repeated l times to obtain a first iterative solution defined as: $\begin{array}{cc}{b}^{\left(1\right)}=i^{\left(1\right)\left(l\right)}=\frac{i^{\left(1\right)\left(l-1\right)}+b^{\left(0\right)}}{2}.& \mathrm{Equation}18\end{array}$

[0077] A second iteration is performed according to:

[0078] i⁽²⁾⁽⁰⁾=y−p*b⁽¹⁾=b+p*b−p*b⁽¹⁾,  Equation 19

[0079] followed by: $\begin{array}{cc}{i}^{\left(2\right)\left(1\right)}=\frac{i^{\left(2\right)\left(0\right)+b^{\left(1\right)}}}{2}.& \mathrm{Equation}20\end{array}$

[0080] Averaging l times, a second iterative solution is defined as: $\begin{array}{cc}{b}^{\left(2\right)}=i^{\left(2\right)\left(l\right)}=\frac{i^{\left(2\right)\left(l-1\right)}+b^{\left(1\right)}}{2}.& \mathrm{Equation}21\end{array}$

[0081] Using equations 15 through 21, an initial direct estimation is defined according to: $\begin{array}{cc}{b}^{\left(n\right)}=\frac{b}{2^l}+\frac{p*b}{2^l}+\left(\frac{2^l-1}{2^l}\delta -\frac{p}{2^l}\right)*b^{\left(n-1\right)},& \mathrm{Equation}22\end{array}$

[0082] where l is a quantity of averaging operations used in each iteration n, and the first two terms can be expressed as $\frac{y}{2^l}.$

[0083] Equation 22 is expanded, and the terms collected, according to: $\begin{array}{cc}{b}^{\left(n\right)}=\left(\delta +p\right)*\left(\delta +\sum _{m=1}^{n-1}{\alpha }^m\right)*\frac{b}{2^l}+{\alpha }^n*b^{\left(0\right)},& \mathrm{Equation}23\end{array}$

[0084] where: $\begin{array}{cc}\alpha =\left(\left(\frac{2^l-1}{2^l}\right)\delta -\frac{p}{2^l}\right).& \mathrm{Equation}24\end{array}$

[0085] In use, a solution converging toward the direct image b may be acquired by applying Equations 11, 12 and 22 , and selecting l such that the scatter kernel norm satisfies {overscore (p)}<2^l.

[0086] In one embodiment, the iterative algorithm for scatter correction is of the form b⁽⁰⁾y−s⁽⁰⁾, as initial estimate, with b⁽ⁿ⁾=α·b⁽ⁿ⁻¹⁾+(1−α)·(y−p*b⁽ⁿ⁻¹⁾), as iterative update, where α is given by $\alpha =1-\frac{2}{{\hat{p}}_{\min }+{\hat{p}}_{\max }+2},$

[0087] and s⁽⁰⁾=p*y, or some other suitable initial estimate of the scatter contribution in the image. The specific choice of α, where {circumflex over (p)}_max={overscore (p)} and {circumflex over (P)}_mindenote the maximum and minimum values, respectively, of {circumflex over (p)}, the Fourier transform of p, gives optimal convergence. However, any value of α that satisfies the condition $\frac{\stackrel{\_}{p}-1}{\stackrel{\_}{p}+1}<\alpha <1$

[0088] can be used. In use, these equations may be applied to obtain a solution converging toward the direct image. This approach may also be combined with the use of the modulation factor a_k, and the use of a re-normalization map as described herein.

[0089] In one embodiment, an initial estimate of the direct image, b⁽⁰⁾=y−s⁽⁰⁾, is used with an iterative improvement, b⁽ⁿ⁾=y−p*b⁽ⁿ⁻¹⁾. Combining these equations yields b⁽ⁿ⁾=b+(−1)ⁿp⁺¹*b+(−1)ⁿ⁺¹pⁿ*s⁽⁰⁾=b+(−1)ⁿ(pⁿ⁺¹*b−pⁿ*s⁽⁰⁾). As discussed previously, this approach does not converge if {overscore (p)}≧1. However, when adding two intermediate results of this iteration, with suitable weights that add up to one, α·b⁽ⁿ⁾+(1−α)·bⁿ⁺¹=b+(−1)ⁿ·(αδ−(1−α)p)*((pⁿ⁺¹*b−pⁿs⁽⁰⁾) can be obtained. Further, by choosing a appropriately, the behavior of the term (αδ−(1−α)p) one can controlled. In particular its zeros in Fourier space can pre-determined. The undesired term in this last expression is the convolution (or product, in Fourier space) of two terms, the second term being identical to the undesired term in the expression for b⁽ⁿ⁾Therefore, terms of the form b⁽ⁿ⁾can be successively combined, as well as suitable linear combinations thereof, to obtain a sequence of images that converges to the direct image. To obtain a sequence of images that converges to the direct image, the weights are chosen such that the zeros in Fourier space are at appropriate locations, and the sum of weights is one.

[0090] In one embodiment, a scatter signal s⁽ⁿ⁾is calculated. The scatter signal estimate s⁽ⁿ⁾is then used by at least one other scatter-correction algorithm to estimate a direct image b⁽ⁿ⁾.

[0091] FIGS. 7 and 8 illustrate the performance of the scatter-correction algorithm described herein on an X-ray phantom image. The phantom includes an acrylic cylinder approximately 5 cm in height with a paraffin core. A discrepancy between the scatter-corrected image and the scatter-free “direct events” data is less than approximately 3% throughout most of the phantom image.

[0092] FIG. 9 illustrates the performance of the scatter-correction algorithm for a breast-tissue equivalent phantom approximately 9 cm thick, wherein the scatter fraction is greater than 0.5, wherein the scatter signal remaining after correction is approximately 5% of the scatter-free images.

[0093] While the invention has been described in terms of various specific embodiments, those skilled in the art will recognize that the invention may be practiced with modification within the spirit and scope of the claims.