The present invention, relates generally to computed tomography (CT) systems, and more particularly to a CT industrial inspection system and method for generating high-resolution CT images using pre-correction techniques.
In present date industrial inspection processes, different types of measurement systems are available such as CT, coordinate measuring machines (CMM) and laser-based profilometry. Each inspection modality has its own advantages and disadvantages associated therewith. Modalities such as CMM and laser-based profilometry can measure external surfaces with high accuracy, but cannot measure internal features unless the part is cut open. To date, CT is one of the more versatile of the measurement/inspection systems for revealing both the internal and external structures of industrial parts in a non-destructive manner. The potential industrial applications of CT include reverse engineering, rapid prototyping, quality assurance, casting simulation & validation, tire development, first article inspection, ceramic porosity inspection, process validation, parts qualification and defect detection, to name a few. However, improved inspection accuracy of industrial CT is desirable, for widespread applications thereof.
For example in the area of reverse engineering, CT has not been optimized for capturing detailed external surface features, which can be crucial for capturing the design intent. The factors affecting CT accuracy in this regard include (among other aspects) beam-hardening, partial volume effect, scattering and off-focal radiation. Thus, in order to improve CT inspection accuracy, more effective methods are needed for removing the effects of these artifacts. In the area of CT image reconstruction, filtered backprojection (FBP) is a common technique because of its fast computation and ease of implementation. However, because FBP oversimplifies the CT data acquisition into an ideal Radon transform (i.e., Fan Beam transform, cone beam transform or any other transform depending on the particular acquisition geometry), the reconstructed image suffers from artifacts such as beam hardening and partial volume as discussed above.
In order to improve image quality, iterative reconstruction techniques have been employed to correct system imperfections such as focal spot size, detector point spread (PSF) function, detector time lag, non-linear partial volume error, scatter, beam hardening etc. Iterative reconstruction techniques are based on different mathematical principles, such as the statistical approach of maximum likelihood, and the least squares approach, for example. These techniques permit the incorporation of a dedicated forward model of the data acquisition. Typically, in an iterative reconstruction approach, the reconstructed image is incrementally updated using the differences between the initial projection measurements and the forward projection model.
Although iterative reconstruction techniques significantly improve image quality, the computational complexity associated with iterative reconstruction is highly intensive, as these techniques require multiple applications of computationally expensive forward and backprojections. Hence, iterative methods are not yet widely used in CT. Accordingly, it is desirable to be able to provide a technique for capturing both internal and external features of an object to be inspected and a technique that improves the image quality of reconstructed images without a significant increase in computation time.
Embodiments of the present technique address this and other needs. In one embodiment, a method for reconstructing image data from measured sinogram data acquired from a CT system is provided. The CT system is configured for industrial imaging. The method includes pre-processing the measured sinogram data. The pre-processing includes performing a beam hardening correction on the measured sinogram data and performing a detector point spread function (PSF) correction and a detector lag correction on the measured sinogram data. The pre-processed sinogram data is reconstructed to generate the image data.
In another embodiment, a CT system for reconstructing image data from measured sinogram data is provided. The CT system is configured for industrial imaging. The system includes an X-ray source configured to project polychromatic X-ray beams through an object and a detector configured to produce electrical signals in response to received X-ray beams from the source. The system further includes a system controller configured to process the electrical signals to generate the measured sinogram data. The system controller is further configured to pre-process the measured sinogram data. Pre-processing the measured sinogram data includes performing a beam hardening correction; a detector point spread function (PSF) correction and a detector lag correction on the measured sinogram data. The system controller then reconstructs the pre-processed sinogram data to generate the image data.
Disclosed herein, is a CT system configured for industrial imaging. As used herein, “industrial imaging” refers to the imaging of both the internal and external structures of industrial parts in a non-destructive manner. As will be appreciated by those skilled in the art, industrial X-ray CT provides rapid three-dimensional measurement for imaging both external and internal features. Disclosed herein, is a CT reconstruction method that provides for improved image quality, in which pre-correction techniques are applied to projection data measurements prior to image reconstruction. Thereby, a more accurately calculated image projection (or sinogram) is used to reconstruct the image of the object. In addition, the reconstruction method, in accordance with the present technique, achieves high resolution while reducing the reconstruction time.
Turning now to the drawings, referring first to
For the exemplary embodiment of
For the exemplary embodiment of
It should be borne in mind that the controllers, and indeed various circuitry described herein, may be defined by hardware circuitry, firmware or software. The particular protocols for imaging sequences, for example, will generally be defined by code executed by the system controllers. Moreover, initial processing, conditioning, filtering, and other operations required on the projection data acquired by the scanner 12 may be performed in one or more of the components depicted in
For the exemplary embodiment of
System controller 22 or operator interface 24, or any remote systems and workstations, may include software for image processing and reconstruction. Therefore, some or all of the image processing may be performed remotely by additional computing resources based upon raw or partially processed image data. As will be appreciated by those skilled in the art, such processing of CT projection data may be performed by a number of mathematical algorithms and techniques. A remote interface 28 may be included in the system for transmitting data from the imaging system to such remote processing stations or memory devices.
Although in the present discussion reference is made to a CT scanning system in which a source and detector rotate on a gantry arrangement, it should be borne in mind that the present technique is not limited to data collected on any particular type of scanner. For example, the technique may be applied to data collected via a scanner in which an X-ray source and a detector are effectively stationary and an object is rotated, or in which the detector is stationary but an X-ray source rotates. Further, the data could originate in a scanner in which both the X-ray source and detector are stationary, as where the X-ray source is distributed and can generate X-rays at different locations as will be described in greater detail with respect to
A number of alternative configurations for emitters or distributed sources may, of course, be envisaged. Moreover, the individual X-ray sources in the distributed source may emit various types and shapes of X-ray beams. These may include, for example, fan-shaped beams, cone-shaped beams, and beams of various cross-sectional geometries. Similarly, the various components comprising the distributed X-ray source may also vary. The emission devices may be one of many available electron emission devices, for example, thermionic emitters, carbon-based emitters, photo emitters, ferroelectric emitters, laser diodes, monolithic semiconductors, etc. Although a distributed source configuration has been discussed in detail here, any combination of one or more rotating-anode, stationary-anode, or distributed X-ray sources may be utilized in the CT system 10.
As will be appreciated by those skilled in the art, conventional X-ray sources for CT imaging typically comprise rotating anode tubes (or X-ray tubes) that possess polychromatic spectrums. That is, the X-ray photons emitted from such X-ray tubes do not all possess the same X-ray energy levels. Moreover, attenuation processes in matter are typically energy dependent. Non-uniform attenuation of different energies results in the preferential depletion of X-rays in energy ranges with higher attenuation coefficients. In general, X-rays in energy ranges that are more easily attenuated are referred to as soft X-rays while those in ranges that are more penetrating are referred to as hard X-rays.
Beam hardening is the process of selective removal of soft X-rays from an X-ray beam. As these X-rays are removed, the beam becomes progressively harder or alternatively more penetrating. The amount of beam hardening typically depends on the initial X-ray spectrum as well as on the composition of the material traversed. In general, for a pre-determined initial X-ray spectrum and tissue type, the process of beam hardening represents a monotonic increase in beam hardness as a function of material thickness traversed. As a result of beam hardening, the effective attenuation coefficient of a material depends on the thickness of material traversed. This effect causes beam-hardening artifacts in CT images. As will be appreciated, beam-hardening artifacts, if uncorrected, create a cupping effect, or a reduction of the reconstructed attenuation coefficient toward the center of the object being imaged.
In order to appreciate the benefits provided by embodiments of the present technique, reference is made to
d(E)=(1−exp(−μ(E)thkness))E (1)
where E is the energy dependent attenuation co-efficient for the detector material, thkness, is the thickness of the detector material and μ is the attenuation. The X-ray source spectrum is scaled by the detector efficiency curve, d(E) to obtain a corrected X-ray source spectrum.
In step 68, an attenuation coefficient value for each material traversed by X-rays, at one or more X-ray energy levels is determined. In particular, the mass attenuation coefficient value, denoted by μ/ρ (where μ denotes the attenuation value and ρ denotes the density) is determined for each material. As will be appreciated by those skilled in the art, the mass attenuation coefficient μ/ρ and the mass energy absorption coefficient, denoted by μen/ρ are basic quantities used to compute the penetration and the energy deposition by photons (which may comprise, for example, X-rays, or bremsstrahlung) in biological, shielding and other such materials. In accordance with this embodiment, to correct beam hardening effects, the mass attenuation coefficients, μ/ρ, for each material is determined. Further, the mass attenuation coefficient μ/ρ at one or more energy levels is pre-computed based upon a material composition and an X-ray source profile as follows.
As will be appreciated by those skilled in the art, the mass attenuation coefficient μ/ρ of a compound or a mixture of elements may be calculated by taking the weighted sum of the mass attenuation coefficients for each element, where a weight is assigned according to the fraction of the element contributing to the compound. Then, the mass absorption co-efficient may be obtained using the formula:
μ/ρ=Σwi(μi/ρi) (2)
where wi is the fraction by weight of the ith atomic constituent. Furthermore, as will be appreciated by those skilled in the art, the values for the mass attenuation co-efficient, μi/ρi, and the mass energy-absorption coefficient, μen/ρ, as a function of photon energy, for various elements may be pre-determined. In addition, the fractions by weight may also be pre-determined. In accordance with the present embodiment, to obtain the μi/ρi values at all the absorption edges of all constituent elements, interpolation is performed. For example, in the case of a constituent element such as liquid water, the composition of Helium (H) and Oxygen (O) have fraction by weights of 0.111898 and 0.888102 respectively. Accordingly, the corresponding mass attenuation co-efficient was water is computed as follows:
(μ/ρ)water=0.111898(μ/ρ)H+0.888102(μ/ρ)o (3)
In step 70, a beam hardening correction curve is analytically derived based on the mass attenuation coefficient at the one or more X-ray energy levels and the corrected X-ray source spectrum obtained in step 66. More specifically, in accordance with the present embodiment, the beam hardening correction curve is analytically derived based on attenuation resulting from beam hardening for an X-ray path length. In particular, a summation of one or more intensity values over the one or more X-ray energy levels is performed, wherein each X-ray energy level has a corresponding mass attenuation coefficient as follows:
I=I0∫Ω(E)−∫μ
where Ω(E) represents the incident X-ray spectrum probability distribution, μM(E) denotes the absorption coefficients for a material M at an energy level E, and I0 and I represent the total incident and transmitted intensities, respectively. The beam hardening correction curve is then analytically generated based on the attenuation resulting from beam hardening and a linearized projection for the X-ray path length S by taking the logarithm of equation (4) as follows:
It may be noted that equation (5) is a non-linear equation. Since, the basic assumption for the filtered backprojection in CT reconstruction is the linear integral equation from Lambert-Beers law, the main purpose of the beam hardening correction is therefore to linearize equation (5), so that filtered backprojection may be applied to the sinogram data for reconstruction without artifacts. In one embodiment of the present invention, this is accomplished by performing a linearization around the effective absorption coefficient μeff, as follows
where μeff=μ(Eeff) is computed for a user defined energy level Eeff. Since, both pbh and p are both functions of the path length for a fixed material, in accordance with embodiments of the present technique, a look up table between pbh and p is created by eliminating the path length S as an argument. Then, the corrected sinogram after beam hardening correction is given by:
p=pcorrected=μeffS=μeffpbh−1(measured sinogram) (7)
As may be observed from equation (7), the effective absorption co-efficient μeff is a scaling factor for the corrected sinogram, thereby enabling the choice of positive numbers for μeff. Referring again to equations (5) and (6), the values of pbh and p are computed for each path length S. The integration in equation (5) is implemented as a summation so that equation (5) is a function of the path length S. After calculating the pbh and p values associated with all possible values of the path length S, either a look-up table or a plot as shown in
Referring again to
As will be appreciated by those skilled in the art, the detector lag and the detector PSF typically cause azimuthal and radial blur in the image domain, respectively. The detector time lag occurs from the fact that if a solid-state detector is exposed to X-ray photons for a period time and is then quickly switched off, the input X-ray and the detector reading does not reach a zero value immediately. Moreover, if the CT gantry speed or rotating table speed is faster than the decay rate of the solid state material, then the detector reading at each projection view is affected from the detector readings from past projection views. This, in turn, causes the azimuthal blur in the reconstruction image.
The detector PSF is typically used to access the spatial resolution of an imaging system. The detector PSF generally originates from the scintillator behavior for the incident X-ray. This results in an X-Ray spread over several detector elements, which is characterized as the detector PSF. Since the PSF is across each detector element at each incident projection view, the PSF is referred to as a view angle invariant characteristic. Mathematically, the measurement from the detector which incorporates this characteristics is denoted by:
y=Hx+n (8)
where y denotes the detector reading, x is the incident X-ray intensity, n denotes the noise and H denotes the convolution operator due to the detector lag and the detector PSF. In accordance with embodiments of the present technique, the convolution operator H is composed of two separable filters in detector view and detector row directions and is indicated as follows:
H=h*g (9)
where h and g denote the blur due to detector lag and detector PSF respectively.
Referring to equation (8), in accordance with the present embodiment, the detector lag and the detector PSF correction is used to obtain the ideal X-ray intensity, x from the detector reading, y by deconvolving, H. In particular, an iterative deconvolution algorithm is applied to the measured sinogram data, to optimally correct the detector PSF and detector time lag. The iterative deconvolution algorithm is formulated as a nonlinear optimization problem as will be described in greater detail below.
Referring again, to equation (8), again, the least square solution for equation (8) is obtained by the following equation:
{circumflex over (X)}=(HTH)+HTy=HT(HHT)+y (10)
where the superscript, T denotes the transpose, and + is the generalized inverse. As is known by those skilled in the art, a disadvantage of the least square solution, obtained using (10) is that the above least square solution does not take into consideration, the non-negative value for the resultant estimate of the intensity {circumflex over (X)}. To obtain a non-negative estimate, a truncation to zero for the negative values may be computed, but such a direct inversion applied to equation (10) may result in an erroneous reconstruction. Therefore, in accordance with the present technique, a novel iterative deconvolution algorithm is applied, to the measured sinogram data. As mentioned above, the iterative deconvolution algorithm optimally corrects the detector PSF and detector time lag. Since, in general, the raw projection data should have a positive intensity, conventional deconvolution techniques that use linear filters (such as a Weiner filter, for example) are not optimal in nature. In accordance with embodiments of the present technique, the deconvolution problem is formulated as a non-linear optimization problem. The optimization problem results in an accurate deconvolution of the detector PSF and detector time lag. The optimization problem is formulated as follows:
subject to the constraint, x≧0
where ∥y−Hx∥2w=(y−Hx)TW(y−Hx) denotes the weighted norm with a weighting matrix, denoted by W. Since, solving the constraint minimization problem presents certain difficulties, an unconstrained minimization problem is formulated to obtain the solution of equation (11). Further, by replacing x=exp(−s), the formulation of equation (11) becomes an unconstrained minimization problem as follows:
C(s)=∥y−H exp(−s)∥2w (12)
wherein
It may be noted that for any real number s, the expression x=exp(−s) cannot be negative. Therefore, in accordance with the present embodiment, the optimization problem depicted in equation (12) is solved without the use of any additional constraints. Also, since the solution to s may be any real number, the optimization problem may be solved using any standard non-linear optimization method such as for example, gradient search. As will be appreciated by those skilled in the art, an optimization problem that is solved without the use of additional constraints has a low associated computational cost. In contrast, more sophisticated algorithms exist that specifically deal with constraint problems, but come with a high associated computational cost.
The embodiments illustrated and described above, provide a technique for reconstructing image data acquired from an imaging system. The beam hardening correction technique described above may be used for a plurality of materials without any tedious and expensive calibration procedures and provides a very accurate correction when combined with detector PSF and detector time lag correction. In addition, the proposed PSF and time lag correction technique, though iterative in nature, has negligible computational complexity as compared to iterative reconstruction techniques, since the iteration is performed in the measurement domain without resorting to computationally expensive forward and backprojections. For example, if 5˜15 iterations are needed to get a reconstruction image by using a certain iterative reconstruction technique, 5˜15 times the number of forward and backprojections would typically need to be performed. Assuming that the time required for one backprojection is Tbp and forward projection costs the same as backprojection, the total time needed for reconstruction will be 10˜30 Tbp using an iterative reconstruction method. On the contrary, the pre-correction step performed in accordance with the present technique, requires one correction step and one back projection step to be performed. Therefore, the total computation time is equal to Tbp+ pre-correction time. In view of the fact that the pre-correction time is known to be smaller than Tbp, the total computation time for pre-correction is lesser than twice the value of Tbp, which is a considerable improvement over the existing iterative reconstruction techniques. Furthermore, the resultant reconstructed CT image with pre-correction in accordance with the present technique achieves high image quality with minimal increments to computational time as compared to conventional filtered backprojection techniques.
While the invention may be susceptible to various modifications and alternative forms, specific embodiments have been shown by way of example in the drawings and have been described in detail herein. However, it should be understood that the invention is not intended to be limited to the particular forms disclosed. Rather, the invention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention as defined by the following appended claims.
Number | Name | Date | Kind |
---|---|---|---|
5099505 | Seppi et al. | Mar 1992 | A |
5692507 | Seppi et al. | Dec 1997 | A |
5774521 | Close et al. | Jun 1998 | A |
5878108 | Baba et al. | Mar 1999 | A |
6067342 | Gordon | May 2000 | A |
6285799 | Dance et al. | Sep 2001 | B1 |
6345113 | Crawford et al. | Feb 2002 | B1 |
6377654 | Willems et al. | Apr 2002 | B1 |
6421409 | Paulus et al. | Jul 2002 | B1 |
6928182 | Chui | Aug 2005 | B1 |
20030007601 | Jaffray et al. | Jan 2003 | A1 |
20030185339 | Heumann et al. | Oct 2003 | A1 |
20050105693 | Zhao et al. | May 2005 | A1 |
20050286749 | De Man et al. | Dec 2005 | A1 |
20060002504 | De Man et al. | Jan 2006 | A1 |
Number | Date | Country |
---|---|---|
0487243 | May 1992 | EP |
1072861 | Jan 2001 | EP |
WO 0157795 | Aug 2001 | WO |
Number | Date | Country | |
---|---|---|---|
20060067461 A1 | Mar 2006 | US |