Regularization Technique for Densitometric Correction

1. Field of the Invention

The present invention relates to scatter-glare correction for x-ray desitometry using deconvolution techniques, and in particular to deconvolution techniques using a calibration phantom and a regularization technique.

2. Description of the Prior Art

X-ray imaging is widely used in medical diagnostics and other fields of application. However, an x-ray image obtained from a subject such as a patient is degraded by many factors resulting in artifacts in the image. For example, an x-ray image of a patient is degraded by x-ray scattering within the patient. As a result of this scattering, the image obtained using an x-ray detector is a superposition of a primary component, or the image of the subject, and a scattered component associated with the primary component. In x-ray fluoroscopy where an image intensifier is used in front of the x-ray detector, images are further degraded by veiling glare in the image intensifier. Veiling glare arises from the scatter of electrons and photons within the x-ray image intensifier and causes the electrons and photons produced from the brighter areas of the image to be spread into the darker areas, resulting in a blurring of the image.

These artifacts must be corrected before accurate quantitative information can be derived from the images. Various post-processing techniques have been developed to remove the scatter and glare effects from recorded x-ray images, such as deconvolution. Deconvolution techniques employ a convolution model to describe the scatter-glare component. Mathematically, the scatter-glare component of an image can be modeled by a spatial convolution of the primary component with a two-dimensional response function. When the response function is known, an estimate of the primary component of the image can be calculated from the measured image intensities by a deconvolution technique (also known as Wiener filtering or constrained restoration). Scatter-glare correction using deconvolution techniques is described in, for example, J. A. Seibert et al., Removal of Image Intensifier Veiling Glare by Mathematical Deconvolution Techniques, MED. PHYS., 12(3), 281-288 (1985); S. Naimuddin et al., Scatter-Glare Correction Using a Convolution Algorithm with Variable Weighting, MED. PHYS., 14(3), 330-334 (1987); L. A. Love et al., Scatter Estimation for a Digital Radiographic System Using Convolution Filtering, MED. PHYS., vol. 14(2), 178-185 (1987); L. A. Seibert et al., X-ray Scatter Removal by Deconvolution, MED. PHYS., 15(4), 567-575 (1988); and Scatter-Glare Corrections in Quantitative Dual-Energy Fluoroscopy, MED. PHYS., 15(3), 289-297 (1988). Various functions may be used for the response function, such as exponential or Gaussian functions. The parameters of a response function are usually determined empirically. Ideally, because the response function for each image is dependent upon various acquisition parameters such as beam energy, imaging geometry and patient thickness, which are variable from patient to patient, the precise response function parameters should be determined for each individual patient to achieve accurate scatter-glare removal.

Various methods for determining the response function parameters have been developed. In the methods described in the above references, the response function parameters are estimated by using an object such as a Lucite slab to mimic a patient, placing a series of lead objects before the Lucite object to completely block the direct x-rays in parts of the image, and measuring the x-ray intensities in the blocked areas of the image which presumably consist solely of the scatter-glare component. These methods only give the estimated parameters for an average or typical patient.

To account for dependence of the response function parameters on acquisition parameters, one deconvolution method allows variation of the response function parameters based on empirical relationships between the acquisition parameters and the response function parameters. However, this method still does not allow response function parameters to be optimized for each individual patient. In another method, an array of beam stops is placed between the x-ray source and the patient when the patient is being imaged to block the direct x-ray, and the response function parameters for the particular patient are determined by comparing the image intensities behind the beam stops (which presumably would be scatter-glare only) with scatter-glare intensities calculated from the convolution model. A disadvantage of this method is a loss of information due to complete blocking of parts of the image. Other techniques involve taking multiple images of the patient, such as with and without beam stops, to estimate the parameters. These methods suffer the disadvantage of increased patient radiation dose and data acquisition time.

Thus, each of these methods suffers one or more of the following disadvantages: (1) the response function parameters are not optimized for each individual patient being imaged, thus reducing the accuracy of the scatter-glare removal; (2) more than one image must be taken to obtain one corrected image, thus increasing patient radiation dose and the length of the x-ray process, and reducing throughput in the x-ray process; (3) the image of the patient is completely blocked in some areas by beam stops or other lead objects placed in the beam. These and other disadvantages often make these techniques unsuitable for clinical applications.

Therefore, it is an object of the present invention to provide a method and apparatus for scatter-glare correction for x-ray images suitable for clinical applications. It is another object of the present invention to provide a method and apparatus for scatter-glare correction without requirements of specific prior knowledge about the subject. It is a further objective of the present invention to provide a post-processing method for scatter-glare correction for a single x-ray image without the need for multiple images. It is yet a further objective of the present invention to provide a method and apparatus for scatter-glare correction without using beam stops or other objects that block parts of the subject image.

These and other objectives are achieved by the preferred embodiments of the present invention. According to a preferred embodiment, a reference object is placed between the x-ray source and the subject and imaged in superposition with the subject. The density distribution of the reference object is designed to have very small spatial correlation with the density distribution of the subject. The scatter-glare component is modeled by a convolution of the primary component with a parameterized response function. During post-processing of the image, the unknown parameters of the response function are adjusted to minimize a cost function based on the expectation that the correlation between the reference object and the subject should be small. As a result of the post-processing, the scatter-glare component of the subject and the reference object, as well as the reference object image, are removed from the raw image and an image of the subject free of scatter-glare effects is obtained.

More particularly, a densitometric calibration phantom having a known spatial density distribution is placed in the x-ray beam while the x-ray image of a subject is being taken. The calibration phantom produces on the recorded x-ray image a primary image of the phantom and an associated scatter-glare component, which are superimposed on the primary image and the scatter-glare component produced by the subject. The densities of the calibration phantom is sufficiently low so that it does not completely block any part of the subject image. The density distribution of the phantom is designed such that the spatial correlation of the gird density and the subject density is very small. The x-ray image is digitized, and the digital image is processed numerically using a regularization technique as follows.

The scatter-glare component in the observed x-ray image is modeled by a spatial convolution of a parameterized response function with a primary component, which represents the density of the x-ray absorbing entities in the x-ray path, such as the subject and the calibration phantom. The parameterized response function used in the convolution has a known form, but contains parameters that are not known prior to processing of the x-ray image. If the statistical noise in the x-ray image is neglected, an estimate of the subject density may be calculated from the measured image intensities, the known calibration phantom density, and a set of estimated response function parameters. By adjusting the values of the estimated parameters to minimize a cost function, the optimum response function parameters can be determined. The scatter-glare component of the subject and the calibration phantom can then be calculated and removed along with the primary image of the calibration phantom to produce a subject image free of scatter-glare components.

Further, other artifacts that affect the x-ray image quality, such as beam hardening, and image blurring due to finite x-ray aperture size, can also be corrected using the regularization technique. This may be done by modeling these effects and including the model parameters in the cost function to be optimized.

The regularization technique may be numerically implemented in computer programs in the disclosed manner. A computer may perform minimization of the cost function, obtain the optimum parameter values, and subtract the scatter-glare component of the subject as well as the primary and the scatter-glare components of the phantom from the digitized image. The digital image so processed may be reproduced graphically on a CRT display or as a photographic image for visualization, or displayed in other formats for further processing.

FIG. 1 is a schematic diagram of the apparatus of a preferred embodiment of the present invention.

FIG. 2 is a block diagram of the regularization technique according to a preferred embodiment of the present invention.

FIG. 1 is a schematic illustration of a method for recording and correcting an x-ray image of a subject using a reference object according to a preferred embodiment of the present invention. An x-ray source 11 forms an x-ray image of a subject 12 on an x-ray detector 13. The x-ray source, such as an x-ray tube, and the x-ray detector, such as a TV camera or a CCD camera coupled with an x-ray image intensifier, are well known in the art. A calibration phantom 15 is placed between the x-ray source and the detector as a reference object. The calibration phantom absorbs x-rays and produces an image on the x-ray detector, but does not completely block any area of the subject image. The phantom may be constructed from a plate of appropriate x-ray absorbing material, preferably tin, by drilling a depression of controlled depth at each pixel. The x-ray image obtained on the x-ray detector contains an image of the subject 12 and the calibration phantom 15, as well as blurred components generated by x-ray scattering and veiling glare in the image intensifier. The x-ray image is digitized using conventional technology to produce a digital representation of the image 16, which is subsequently processed using a novel regularization technique to remove the effects of scatter and glare and to produce a restored image 17 of the subject.

The processing of the digitized x-ray image requires an understanding of the effects of the system on the digitized image. Empirically, the observed x-ray image 16 includes a primary component, a scatter-glare component, an offset, and random statistical noises. The primary component P is produced by x-ray attenuation by the subject and the calibration phantom as follows:

where I.sub.0 is the incident x-ray intensity from the x-ray source, and .rho..sub.s and .rho..sub.c are the spatial x-ray density distributions of the subject and the calibration phantom, respectively. X-ray density is defined here as the negative logarithm of x-ray attenuation by the respective objects. The density distribution of an object is a two-dimensional distribution defined in a desired plane substantially perpendicular to the x-ray beams, where the density is integrated in the third dimension, or the direction of the x-ray beams. The density distributions .rho..sub.s and .rho..sub.c used in equation (1) are defined in the image plane, so that they include any relative magnification or blurring of the actual objects due to the finite x-ray aperture size. The effects of magnification and blurring will be discussed in more detail below. While the x-ray density of the object .rho..sub.c is known, .rho..sub.s is generally unknown and must be derived.

The scatter-glare component of the image is modeled as a convolution of the primary component with a response function H parameterized by one or more parameters .sigma..sub.0. The response function H is preferably Gaussian or exponential. The total observed image intensity I(x,y) is thus given by :

where .alpha..sub.0 is a proportional constant, .beta..sub.0 is the offset component, and N is the random statistical noise. The offset may arise from internal noise in the x-ray detector, such as dark currents. The statistical noise arises from the fact that the observed x-ray intensities are produced by a finite number of photons striking on the detector. The nature of these noises are well understood by those skilled in the art. By taking the spatial Fourier transform of equation (2), solving for P by taking the inverse Fourier transform, substituting equation (1) for P, and taking the natural logarithm of the resulting equation, we have: ##EQU1## where the primes (') refer to the spatial Fourier transform of the respective functions of the image intensity, noise and the response function, k is the spatial frequency, .delta. represents a delta function which is the Fourier transform of a constant function, and F.sup.-1 represents the inverse Fourier transform. Equation (3) is exact if the values of parameters .alpha..sub.0, .beta..sub.0 and .sigma..sub.0 are known.

Generally, the values of these parameters are unknown and are approximated by a set of values .alpha., .beta. and .sigma., and the estimate of the subject density is:

where ##EQU2## is a quantity calculated from the measured x-ray intensities l, noise N and estimated parameter values .alpha.,.beta. and .sigma.. Since the incident x-ray intensity I.sub.0 is assumed to be spatially uniform but unknown in value, only the spatially non-uniform component of the subject density .rho..sub.s can be estimated:

where the tilde symbol () indicates the difference from the mean value. The constant term In(I.sub.0) in equation (4) drops out from equation (6).

When the incident x-ray intensity is sufficiently strong, the statistical noise N will be much smaller than the image intensities I, and may thus be neglected in the regularization technique of the preferred embodiment. When the noise term N' is neglected, the above equation (6) gives the estimate of the non-uniform component of the subject density .rho..sub.s calculated from the measured image intensities I, the known phantom density .rho..sub.c, and a set of estimated parameters .alpha.,.beta. and .sigma..

The convolution theory for the scatter-glare effects has been described above. To correct for such effects, the spatial correlation (or the lack of thereof) between the subject density and the calibration phantom is used.

The correlation coefficient between the density distributions of the subject .rho..sub.s and the calibration phantom .rho..sub.c is defined as ##EQU3## where the summation is over all the pixels in the interested region of the image.

In the preferred embodiment, the calibration phantom is designed such that the phantom density .rho..sub.c is spatially uncorrelated to the subject density .rho..sub.s ; i.e., C(.rho..sub.s,.rho..sub.c).apprxeq.O. Preferably, the phantom density is such that it is uncorrelated to any subject density typically seen in actual medical applications. Various designs of such phantom densities are possible. In the preferred embodiment, the phantom density is a constant term (based upon the uniformly thick layer of tin) plus uncorrelated Gaussian distributed random noise (based on the depressions drilled in the tin) having a flat power spectrum of spatial frequencies (up to the highest frequency which is the inverse pixel size) and nearly uniform local variance. Thus, the noise at each pixel of the phantom is preferably a random number having a Gaussian distribution and a variance that is uniform for all pixels.

The ensemble mean correlation coefficient for such calibration phantoms with any signal density is zero as a result of the distribution being Gaussian, and the ensemble variance is ##EQU4## where n is the number of pixels in the image region of interest. If a large number of pixels are considered, then this variance is extremely small.

The correlation condition in equation (8) is incorporated in the regularization technique of the preferred embodiment as a part of the cost function used for minimizing scatter effects. The theoretical mean expected value of .vertline.C.vertline..sup.2 (the square of the correlation coefficient between the calibration phantom and the estimated subject density) is 1/n. Empirically, this value may be determined for a particular calibration phantom by averaging measured correlations between the phantom and typical subject images in which scatter-glare component of the subject has been removed using other techniques, such as by using beam stops. Since the mean expected value of the square of the correlation coefficient .vertline.C.vertline..sup.2 for a particular calibration phantom is a property of the phantom, this value only needs to be determined once for each phantom and is then used in the regularization process for subsequent images.

In the convolution model for image formation described above, equation (6) relates .rho..sub.s, the estimated non-uniform component of the subject density distribution, to the known object density .rho..sub.c, the measured image intensity I and the estimated parameters .alpha.,.beta., and .sigma., the random statistical noise N having been neglected. Since this convolution constraint applies at each of the n pixels, equation (6) gives n constraints. Another constraint that relates the estimated subject density to the known and measurable quantities is the correlation constraint, or equation (8). These n+1 constraints are to be satisfied by varying a total of n+3 unknowns, namely, .alpha.,.beta.,.sigma. and .rho..sub.s. The solution is therefore not unique, and some form of regularization is required.

In the preferred embodiment, regularization is numerically performed by minimizing a cost function with respect to parameters .alpha.,.beta., and .sigma., such as the following cost function: ##EQU5## where .alpha..sub.1,.beta..sub.1 and .sigma..sub.1 are initial parameter values, i represent the pixels, and R represents subregions of the image as described in more details below. The partial derivatives of .rho..sub.s to the parameters are calculated from equations (6) and (5), and are evaluated at the initial parameter values and fixed during minimization. The partial derivatives may be approximated by finite differences. .lambda. is a coefficient representing the relative weight of the cost terms. The parameter .lambda. is set at a fixed value during minimization, and is adjusted until the minimization yields a solution .rho..sub.s such that the sum .SIGMA..vertline.C.sub.R .vertline..sup.2 falls within an acceptable range as described in more details below.

The first three terms in the cost function associated with .alpha.,.beta., and .sigma. represent a requirement that the parameters determined by the minimization not deviate significantly from their initial values. The initial values .alpha..sub.1,.beta..sub.1 , and .sigma..sub.1 used in the cost function preferably reflect the expected physical properties of the subject relevant to the scatter-glare effects, and the properties of the imaging system, such as imaging geometry. Such values used may be typical values for these parameters as determined by other correction techniques which use the same convolution model for the scatter-glare effects, such as by using beam stops. They may be obtained from a typical patient or by averaging over several patients. These terms as written in partial derivatives to the parameters prevent linear degeneracies from allowing the solution .alpha.,.beta., and .sigma. to wander far from physically sensible values. Linear degeneracies reflect the fact that the effects of parameters .alpha.,.beta., and .sigma. may not be entirely independent, and that each parameter may wander far away from its initial value (which is the physically sensible value) yet the combined effect of all parameters still produce a low correlation between the subject and the phantom.

The last term of the cost function of equation (9) represents the requirement that the correlation between the calibration phantom and the subject density be preferably small. In the regularization process of the preferred embodiment, the image is divided into a plurality of subregions of the image as represented by R in equation (9), and the correlation coefficient for each individual subregion are calculated, squared, and summed up. The number of subregions may be four (two by two), or other numbers. Subdividing the image and calculating the correlation in each region separately reduces the possibility that subregions of the image have chance correlations with the phantom which cancel each other when the correlation coefficient of the entire image is calculated.

The mean expected value for the sum .SIGMA..vertline.C.sub.R .vertline..sup.2, denoted <.SIGMA..vertline.C.sub.R .vertline..sup.2 >, is determined in advance for a particular phantom and a particular division of the image into subregions. For example, if the image is divided into R subregions of equal sizes, the theoretical mean expected value <.SIGMA..vertline.C.sub.R .vertline..sup.2 >is R.sup.2 /n for a calibration phantom having uncorrelated Gaussian distributed random noise as described hereinbefore. Empirically, the mean expected value of the sum for the particular phantom and the particular image division may be determined by calculating the average correlation between the phantom and subject images that have been corrected for scatter-glare effects using other techniques, as described above.

FIG. 2 is a flow chart of the regularization process according to the preferred embodiment. The parameters .alpha.,.beta., and .sigma. are initialized in step 21, followed by the initialization of the weighting constant .lambda. in step 22. The initial values .alpha..sub.1,.beta..sub.1, and .sigma..sub.1 are set at the values for a typical or average patient, as described above. The initial value of .lambda. is not critical, and may be set at, for example, 1. After initialization, the cost function is minimized with respect to parameters .alpha.,.beta., and .sigma. in step 23. Minimization is performed numerically using a standard technique well known in the art, such as the Davidon-Fletcher-Powell algorithm as described in, for example, Press, Flannery, Teukolsky and Vetterling, NUMERICAL RECIPES (Cambridge University Press, 1986). The cost function is evaluated using equations (5) and (6), when the partial derivatives to the parameters are approximated by finite differences and are evaluated at the initial parameter values .alpha..sub.1,.beta..sub.1 and .sigma..sub.1. In evaluating the cost function, standard two-dimensional fast Fourier transform (FFT) and inverse FFT techniques are used.

After a minimum is reached, the sum of square correlations .SIGMA..vertline.C.sub.R .vertline..sup.2 is evaluated in step 24 to determined whether it falls within an acceptable range around the mean expected value <.SIGMA..vertline.C.sub.R .vertline..sup.2 >as established above. If the sum .SIGMA..vertline.C.sub.R .vertline..sup.2 falls outside of the acceptable range, such as 5% on either side of the mean expected value, then the value of the weighting constant .lambda. is adjusted (in step 26) to a new value .lambda..sub.new, for example,

where .lambda..sub.old is the old .lambda. value used in the previous minimization process. Other .lambda. adjustment methods may be used; however, the adjustment must be in the appropriate direction determined by the value of the sum .SIGMA..vertline.C.sub.R .vertline..sup.2. Generally, if the sum falls outside of the acceptable range and is greater than the mean expected value, .lambda. is usually increased to give the .SIGMA..vertline.C.sub.R .sup.2 term more weight in the cost function. If the sum falls outside of the acceptable range and is less than the mean expected value, .lambda. is usually decreased to give the .SIGMA..vertline.C.sub.R .vertline..sup.2 term less weight in the cost function.

The minimization procedure is then repeated in step 23 using the new .lambda. value in the cost function. The process of adjusting .lambda. and minimizing the cost function is repeated until a .lambda. value is found such that the sum .SIGMA..vertline.C.sub.R .vertline..sup.2 falls within the acceptable range when the cost function is minimized. When this occurs, the minimization process is terminated, and the subject density distribution .rho..sub.s calculated in the last cycle of minimization is saved as the optimum subject density distribution, together with the optimum parameters, in step 25. This density distribution gives the restored image (the non-uniform component) from which scatter-glare component of the subject and the calibration phantom, as well as the primary image of the phantom, have been removed.

Because the above regularization technique adjusts the weighting coefficient .lambda. and minimizes the cost function while holding .lambda. fixed, it is possible that not a single .lambda., but .lambda. values within a range will produce the result that the sum .SIGMA..vertline.C.sub.R .vertline..sup.2 falls within the acceptable range when the cost function is minimized for that .lambda. value. This corresponds to a plurality of sets of parameter values .alpha.,.beta. and .sigma., all of which satisfying the above regularization requirements. These sets of .alpha.,.GAMMA. and .sigma. values are all considered acceptable and are not discriminated.

The regularization technique described above may also be used to optimize other acquisition parameters, such as beam hardening and blurring due to finite x-ray aperture size. Beam hardening is an artifact arising from the wavelength dependence of x-ray absorption coefficient of matter. X-rays of longer wavelengths (softer x-rays) are absorbed more by the same material than harder x-rays are, causing the relationship between the thickness h of the absorbing object and its x-ray absorption density .rho. to deviate from a linear relationship. This effect may be empirically modeled by including a linear variation with the primary density in the proportionality constant .mu. between the thickness h and the x-ray density .rho.:

and

where .mu..sub.0 and .eta. are constants. The parameters .mu..sub.0 and .eta. may be determined empirically, or may be optimized in the regularization process. If these parameters are included in the regularization, then the regularization is based on the phantom thickness h.sub.c and effective subject thickness h.sub.s =.rho..sub.s /.mu., rather than x-ray densities .rho..sub.c and .rho..sub.s. Further, the calibration phantom should be constructed from a material with absorption characteristics similar to that of the x-ray contrast material used in the patient. For example, if iodine is used as the x-ray contrast material, tin is preferably used to construct the calibration phantom.

An alternative method for reducing the beam hardening artifact is to construct the calibration phantom from opaque specks narrower than the x-ray aperture but thicker in the x-ray direction than the half-value thickness (the thickness of a material which attenuates the x-ray beam to a half of its incident intensity) for diagnostic energy photons, which is less than 0.17 mm Pb for 70 keV photons. The specks are preferably elongated rods or wires and aligned with the fan beam of x-rays. Such a phantom would attenuate x-rays without changing the x-ray energy spectrum.

Effects of image magnification and blurring may also be corrected using the regularization process. Magnification of an object refers to the relationship between the actual size of the object and the size of the image on the x-ray detector. Blurring refers to half-shadows of an object formed on the image due to the finite size of the aperture of the x-ray source. The amounts of magnification and blurring are dependent upon the distances from the x-ray aperture to the object (subject or phantom) and from the object to the x-ray detector. The calibration phantom and the subject (patient) are typically placed at different distances from the x-ray source and the detector, causing their respective images to have different magnification and blurring. The effects of magnification and blurring of the calibration phantom may also be modeled, and the model parameters, such as the x-ray aperture size and the source-to-image-intensifier distance, may be optimized in the regularization process. When the magnification and blurring effects are included in the regularization, the regularization is based on the density distribution of the calibration phantom defined at the position of the phantom, rather than the density distribution defined at the x-ray image plane. The subject density distribution may be defined either at the image plane or the subject location. Thus, the regularization technique described here may be carried out using either the density distributions defined at the image plane, or the density distributions defined at the respective locations of the subject and the phantom.

Although a preferred embodiment is described in detail, it should be apparent to those skilled in the art that variations in apparatus setup and modeling techniques are possible. For example, the calibration phantom may be physically constructed in a variety of ways to give a desired x-ray density distribution. Further, other forms of response function, cost function, or minimization techniques than described herein may also be employed in the regularization process.

In addition, other techniques may be utilized to improve the performance of scatter-glare removal using the regularization technique of the preferred embodiment. For example, during diagnostics, multiple temporal images are often obtained from a patient. Since the response function parameters and other acquisition parameters should remain the same, they may be optimized globally using all of the temporal images while giving multiple temporal images free of scatter-glare effects. This improves the statistics for the parameter optimization. Further, the calibration phantom may be shifted and/or rotated between successive temporal images to minimize chance correlation between the calibration phantom and the subject.

Moreover, the calibration phantom may be placed at various positions along the x-ray beam path depending on what effects are to be corrected. To correct for both scatter and glare, the calibration phantom is preferably placed between the x-ray source and the patient, so that the image of the phantom is subject to the scatter effect due to the patient. To correct for veiling glare only, the phantom is preferably placed in front of the image intensifier but behind the subject so that only veiling glare effect is present for the phantom. Furthermore, two phantoms may be used to correct separately for scatter and veiling glare. One phantom placed in front of the image intensifier may be used to perform veiling glare correction independent of scatter, and a second phantom may be placed between the x-ray source and the subject to correct subsequently for scatter. If scatter and glare are separately corrected, the response functions used are preferably exponential for veiling glare and Gaussian for scatter.

Although a calibration phantom density distribution having simulated random noise is described, other calibration phantoms having highly structured x-ray density distribution may also be used. One example would be a phantom comprising pie-shaped segments, each segment having a density distribution proportional to the radius with different proportional constants for different segments. Another example would be a phantom comprising a grid, for example a 7 by 7 grid, where each rectangular segment of the grid has a uniform density.

These and other variations should be apparent to those skilled in the art and within the scope of the present invention, which should be determined by the claims.