Correction for Polychromatic X-Ray Distortion in CT Images

This invention relates to methods and apparatus for computed tomography. More specifically, this invention relates to image processing apparatus which corrects for polychromatic distortion in images produced by the method of computed tomography.

Machines for producing transverse images of body sections by the methods of computed tomography are known, for example, for U.S. Pat. Nos. 3,778,614 and 3,924,129 which are incorporated herein, by reference, as background material. In such apparatus one or more beams of penetrating radiation, typically X-rays, are projected through the body in a plurality of directions and are measured, typically with electronic radiation detectors, to yield a multiplicity of projections of internal body structures. The projections are then combined typically in a digital computer using, for example, a convolution-backprojection technique, to generate images of transverse sections through the body.

Early methods for computing the transverse image from its projections generally assumed a linear relationship between the lengths of the various constituents and the total attentuation of these constituents so that the integrated tissue density along the path was equal to the logarithm of the ratio of the radiation intensity entering and exiting the body. This assumption, although generally true for a monochromatic radiation source, produces aberrated images if utilized in a scanner having a polychromatic radiation spectrum in conjunction with body constituents having attenuation coefficients which vary with radiation energy. Prior art scanners have included filters for hardening the X-ray beam (to reduce its low energy spectral content) to partially eliminate polychromatic effects. Many prior art scanners have also attempted to compensate for polychromatic effects by effectively assuming a single attenuation function for all body tissues and applying that function, in conjunction with a known spectrum from the X-ray source, as a first order compensation in the image reconstruction calculations (single spectrum or one dimensional corrections).

Virtually all human body tissues are found to have energy dependent X-ray attenuation characteristics which are dominated by the characteristics of water (soft tissues) and bone and can be approximated by a combination of these characteristics. The energy attenuation spectra of water and bone are, however, substantially different. A polychromatic radiation beam propagating through a body which comprises a mixture of bone and soft tissue (on either a macroscopic or microscopic level) will necessarily be influenced by the combined spectra of calcium and water, which interact in a non-linear fashion to distort X-ray intensity values in the measured projections.

The method and apparatus in the present invention operate on an original CT image which has no correction for polychromatic X-ray distortion or only a rough pre-reconstruction correction (i.e. a single spectrum correction) applied to the original projection data. The method determines an error image based on information extracted from the original distorted image. The error image is subtracted from the original image to obtain a corrected image. The following steps are involved in obtaining the corrected image:

Estimates of the projections of the various biological tissues (e.g. bone and soft tissue) are obtained artifically from the original digitized image to the extent that these various tissues can be distinguished by means of their grey levels and by a priori geometric and other structural knowledge of the anatomical section that corresponds to the image;

Error projections are calculated from the projections of the various biological tissue using a precalculated polynomial in as many variables as there are distinguishable biological tissues. The precalculated polynomial is determined by using the X-ray energy spectrum of the X-ray source in the scanning apparatus at the particular kilovoltage at which the original projection measurements were made and the linear attenuation coefficients of the biological tissues as a function of energy in such a way as to enforce a multidimensional linear relationship between the integrated attenuation and the equivalent lengths of the distinguishable tissues through which the X-ray beam passes. If a pre-reconstruction correction was made on the original projection data, then the precalculated multidimensional polynomial is modified to take this pre-reconstruction correction, if any, into account;

The error projections are then filtered to remove ripple that is contributed by the projection of a digitized image;

From the error projections an error image is reconstructed by means of either the same reconstruction process that was used to produce the original image or by some other reconstruction process of sufficient accuracy;

The original image and the error image are then subtracted, pixel by pixel, to obtain a corrected image.

It is, therefore, an object of the invention to correct polychromatic distortion in computed tomographic images.

Apparatus for performing the methods of the invention comprises:

means for analyzing the values in an original image array and for assigning to each pixel element a specific proportion of the attenuation coefficient of two or more constituent tissues;

Means for combining the proportions of said attenuation coefficients to generate a multiplicity of constituent projections for each of said constituents;

Means for combining said constituent projections to generate an error image wherein the value of each pixel represents the difference between a pixel value in the original image and that pixel value in a corrected image; and

Means for subtracting the error image from the original image to generate a corrected image.

The method and apparatus of the present invention thus operate to produce an error image which is subtracted from an original image rather than reconstructing a corrected image from corrected projections. Quantization noise and reconstruction artifacts which might otherwise affect the quality of a corrected image reconstructed from corrected projections are thus reduced.

FIG. 1 illustrates the attenuation of a polychromatic X-ray beam in a homogeneous material;

FIG. 2 is an image correction system of the present invention; and

FIG. 3 illustrates a method for projecting pixel elements.

It is known that for an object consisting of a single homogeneous material, such as calcium or water, the X-ray attenuation function has a shape in the form of curve A in FIG. 1. This function depends on the incident radiation spectrum and on the type of homogeneous material. Its departure from linearity results from the shift and change of the shape of the photon energy spectrum as photons of different energy undergo different attenuation in passing through the same length of material. In FIG. 1, the initial slope B of the curve may be considered as the incremental linear X-ray attenuation coefficient associated with the polychromatic X-ray beam. When such an X-ray beam passes through a composite material made up of two or more distinct homogeneous materials the attenuation is a non-linear function of the lengths of the materials and is not a sum of individual functions of single variables.

Experience indicates that the X-ray attenuation of human body structures may be characterized by the attenuation of a heterogeneous structure of soft tissues (having a water-like energy attenuation spectrum) and bone (having a "compact bone"-like spectrum). Additional body structures may, in fact, comprise air or other gases, but the attenuation of such structures is so low, compared with bone and soft tissue, that its energy spectrum has an insignificant effect on image calculation. The non-linear interaction of bone and soft tissue in the attenuation of a polychromatic X-ray beam may be expressed by a power series of the form

The interaction may be expressed in the form

where s.sub.

is a two-dimensional polynomial of degree N, and L(s.sub.1,s.sub.2) is the linear part of the function U(s.sub.1 s.sub.2). By choosing the criterion for approximation properly, and the degree N of the polynomial and sufficiently high, .epsilon.(s.sub.1,s.sub.s) may be made sufficiently small for calculation purposes, so that T(s.sub.1,s.sub.2) is the correction which must be applied to compensate for the non-linear part of the interaction. Experience indicates that a two-dimensional cubic correction

The interaction of an X-ray beam with a heterogeneous material may be approximated by numerical integration using known spectra for the various constituents of the material and measured energy spectrum data for a particular X-ray source operating at a particular voltage. Attenuation coefficients for water and compact bone are, for example, tabulated in the publication "Photon Cross-sections, Attenuation Coefficients, and Energy Coefficients from 10 KeV to 100 GeV" by J. H. Hubbell, National Bureau of Standards, National Standard Reference Data Series NSRDS-NBS 29, Issued August 1969. Energy spectrum data for a particular X-ray source is normally obtained by direct measurement of each type source at its expected operating voltages. The coefficients of the cubic approximation are calculated using any of the well known approximation algorithms. Typically a set of coefficients will be calculated in advance for each X-ray source and operating voltage and stored for later use with raw images measured at the same spectral parameters. By way of example, Table I tabulates the measured energy spectrum area J(E).DELTA.E for a Tomoscan 200 CT scanner (manufactured by Philips Medical Systems, Incorporated of Shelton, Connecticut, which utilizes a Philips beryllium window X-ray tube) measured at 150 KVP with a three millimeter thick aluminum filter inserted in the beam. Corresponding attenuation coefficients for compact bone (.mu..sub.CB) and for water (.mu..sub.H.sbsb.2.sub.O) are also tabulated. Table II is a listing of a Fortran IV computer program for calculating the polynomial coefficients by a suitable least squares approximation, and Table III are the corresponding coefficients calculated thereby from the data of Table I. The coefficients listed in Table III are utilized in the further examples of image correction methods set forth below.

FIG. 2 is apparatus for correcting images in accordance with the invention. A computerized tomographic scanner 10 which includes an X-ray source 12, a detector bank 14 and an image reconstruction computer 16 functions, in accordance with the methods of the prior art, to project X-rays through a body 18 along a plurality of beam paths to measure and record a series of X-ray projection data taken through the body 18 from a plurality of directions and to subsequently combine those projections, using any of the known image reconstruction algorithms, to produce a matrix of discrete element of a transverse image of the body wherein the numerical values of the elements represent the intensity in corresponding pixels of the transverse image. The matrix of image elements is stored in an image storage device 20, which may, for example, comprise core memory or disc storage. The raw image may be directly displayed, as in prior art scanners, on a display device 22.

In accordance with the present invention the raw transverse image matrix produced by the computer 16 and stored in the image storage 20 is processed in an image correction processor 24 to compensate for polychromatic aberration. A projection generator 26 functions to assign relative proportions of the attenuation coefficient to soft tissue and compact bone in each pixel element represented in the image storage 20. The mixture of soft tissue and bone represented in each pixel of the raw image may correspond to a macroscopic combination of bone and soft tissue structures lying within the pixel area or may, alternately, represent an intimate mixture as in varying bone or cartilage structures. The assignment of a proportion of the attenuation coefficient to bone and soft tissue in each pixel element may be based on a pattern recognition process and known structural details of the raw image, but is most readily accomplished by a multiple thresholding process which assigns a percentage of the attenuation coefficient to compact bone, soft tissue, or contrast media in each pixel element. For example, experience indicates that all pixels having a grey scale level L greater than 100 Hounsfield units may be assumed to contain bone and that the percentage of the attenuation coefficient due to bone and soft tissue in such elements may be approximated by a linear interpolation of the grey scale value, relative to the upper and lower thresholds for soft tissue and compact bone, respectively.

A value of grey level due to soft tissue and bone content is thus assigned to each pixel of the raw image and is used to generate sets of separate projections of soft tissue and bone from the raw image data. The projections thus generated correspond to a decomposition of the projections which were measured by the scanner 10 and were utilized for the original image reconstruction, and the process of generating projections of soft tissue and bone from the raw image in the image storage 20 is the mathematical adjoint of the operation of backprojection used to generate the raw image from the scanner convolved projection data in the computer 16. There are, of course, many algorithms and methods for generating images from projections and it is not necessary that the process for generating bone and soft tissue projection sets from the raw image correspond to the exact adjoint in the algorithm used in the scanner 10 to generate the raw image from the scanner projection data.

FIG. 3 illustrates a preferred method for generating projections of bone or soft tissue from an image matrix and is related to the so-called strip method employed in itterative reconstruction algorithms. It can be implemented in a manner similar to that of backprojection either in a general purpose digital computer or in a dedicated hardware array processor. For each projection angle .phi. a series of equally spaced rays are assumed through the picture matrix and each pixel is assigned to that ray, m, nearest to its center (x,y). The values of the pixel elements assigned to each ray are then summed, the set of sums being the projection at the angle .phi. corresponding to the ray direction. Other projection methods, for example direct projection or Fourier transform projection, are also suitable. By way of example, Tables IV and V are preferred embodiments of machine language computer programs, for operation on PDP 11 series computers, which function to threshold raw image data and generate soft tissue and bone projections, respectively, from that data.

The projection generator 26 thus produces two sets of projection data. A first set corresponds to a plurality of projections, at different angles through the image plane, of the bone or calcium structures in the raw image and is stored in a bone projection storage 28 which may, for example, comprise core memory or disk storage. A second set describes corresponding projections of the soft tissue structures in the raw image and is stored in a second storage area 30.

The bone projections stored in device 28 and the soft tissue projections stored in device 30 are then combined in an error projection generator 32 which utilizes a precalculated polynomial, determined in the manner described above from the X-ray spectrum of the source 12 and the linear attenuation coefficients of the biological tissues, as a function of energy, to calculate projections of polychromatic aberration errors in the raw image data.

If a single-spectrum type pre-reconstruction correction for energy spectrum effects was made during the calculation of the raw image data in the computer 16 the pre-calculated polynomial is modified to take into account this pre-reconstruction correction.

The error projection generator 32 may comprise a dedicated hardware processor or may comprise a general purpose digital computer programmed to calculate the error projections from the soft tissue projection and the bone projection data. By way of example, Table VI is a Fortran language computer program which performs the error projection generator function. The calibration coefficient and system coefficient at lines 20 and 25 of the program are scaling factors related to the particular scanner 10 utilized to generate the raw data and have values of 5,000 and 614 respectively for a Tomoscan 200 scanner. The coefficients at lines 12-18 of the correspond to the polynomial coefficients in Table III calculated for the Tomoscan 200 operating at 150 KVP. Lines 31-38 of the program compensate for the single-spectrum pre-reconstruction correction applied in the Tomoscan 200 scanner. The actual computation of the polynomial value is accomplished at line 66.

The error projections produced by the projection generator 32 are filtered in a digital filter 34 to remove noise which inherently results from the projection of a quantized image. The digital filter 34 is, ideally, tuned to the projection generator 26. A preferred embodiment for use with a projection generator described above comprises a three point averaging filter in cascade with an interpretive filter. The interpretive filter functions, for each data point in the projection, to take the average value of increasingly large sets of points surrounding the data point (i.e. three points, five points, seven points . . . ) until the difference between the data point value and the surrounding average value is less than a predetermined threshold. The filter will not, however, increase or decrease the number of points in the averaging set by more than one point for adjacent data points.

The digital filters described above may be implemented as dedicated hardware units or as program modules in a general purpose digital computer. By way of example, Table VII is a Fortran lanugage program for the three point averaging filter described above while Table VIII is a Fortran language program for the interpretive filter. The filtered error projections from the digital filter 34 are then combined in an image reconstruction computer 36 to produce an error image data set which corresponds, on an element by element basis, to the polychromatic distortion error in the raw image in the image storage 20. The image reconstruction computer 36 may be functionally identical to the image reconstruction function in the computer 16 which computes the raw image from the X-ray projections measured by the scanner 10 and may, thus, comprise any of the hardware or software image computers which are known and described in the prior art.

The error image produced by the image reconstruction computer 36 is then subtracted, on a point by point basis, from the raw image held in image storage 20; the function being preformed in an image subtractor 38. The corrected image thus produced is fully compensated for polychromatic distortion and is held in a corrected image storage device 40 for subsequent display on the display device 22. By way of example, Table IX is a Fortran language computer program which may be utilized to preform the function of the image subtractor 38. As will be recognized by those skilled in the art, the image subtractor 38 may, alternately, comprise a hardware digital subtractor.

Although the preferred embodiments of the invention have been described herein with individual components corresponding to program modules for execution in a general purpose digital computer, it should be recognized that, in a given dedicated system, increases in speed and efficiency may be derived by constructing some or all of the individual components as dedicated digital hardware. It will likewise be recognized that the specific construction of these individual components is necessarily highly dependent on the nature and organization of other computing and data storage components in the system but that the methods for producing such hardware from the software embodiments set forth herein are well known. Further, although the present system utilizes a two-dimensional polynomial to compensate for two tissue constituents, a higher dimensional polynomial may similarly be utilized to compensate for other tissue constituents or contrast media.