The present invention relates to an apparatus and a corresponding method for artifact correction of a data set of X-ray projections of an object for generation of a reconstruction image of said object. The invention relates further to an apparatus and a corresponding method for generating a reconstruction image from a data set of X-ray projections of an object. Still further the invention relates to a computer program for implementing said methods on a computer.

Scattered radiation constitutes one of the main problems in cone -beam computed tomography. Especially for system geometries with large cone angle and therefore a large irradiated area, such as C-arm based volume imaging, scattered radiation produces a significant, spatially slowly varying background that is added to the desired detected signal. As a consequence, reconstructed volumes suffer from cupping and streak artifacts or, more generally, from artifacts causing slowly (locally) varying inhomogenities due to scatter, impeding the reporting of absolute Hounsfield units.

Mechanical anti- scatter grids have been designed to prevent detection of scattered radiation, but they have been shown to be ineffective for typical system geometries for volume imaging because they lead to degradation of the signal-to-noise ratio. Therefore, different algorithms for a posteriori software-based scatter compensation have been proposed (e.g. in Maher K. P., Malone J. F., "Computerized scatter correction in diagnostic radiology", Contemporary Physics, vol. 38, no.2, pp. 131-148, 1997) or are currently developed. However, though such methods have the potential to accurately estimate the shape of the spatial distribution of scatter within the projected views, accurate quantitative scatter estimation is difficult to achieve. As a consequence, the absolute local amount of scatter in the projected views is often under- or overestimated, leading to suboptimal reconstruction results.

There are other sources of artifacts in an X-ray projection that also cause spatially slowly varying inhomogenities in a reconstruction image which are, for instance, an incomplete data set used for the reconstruction due to the use of a detector which is smaller

than the object of interest. It will then be desired to complete the data set to avoid the appearance of such artifacts. Standard algorithms (such as described e.g. in R. M. Lewitt, "Processing of incomplete measurement data in computed tomography", Med. Phys., vol. 6, no. 5, pp. 412-417, 1979) require the determination of a projection extension factor. Still further, it is often required to determine a gain factor for normalization of projection data before use in a reconstruction. Sometimes the gain images for normalization are only known with an unknown global factor. Normalizing the measured projections with gain images including a wrong global factor would cause again cupping shaped spatially slowly varying inhomogenities in the reconstructed image.

It is an object of the present invention to provide an apparatus and a corresponding method for artifact correction of a data set of X-ray projections of an object, in particular for correction of artifacts causing cupping or inverse cupping (called capping) shaped spatially slowly varying inhomogenities caused by e.g. scatter, a wrong truncation extension factor or a wrong gain factor. It is a further object to provide an apparatus and a corresponding method for generating a reconstruction image from a data set of X-ray projections of an object including less or no artifacts.

The object is achieved according to the present invention by an apparatus as claimed in claim 1, comprising: an estimation unit for estimating in an X-ray projection the amount of artifact present in said X-ray projection using at least one estimation parameter, a correction unit for correcting said artifact present in the X-ray projection by use of said estimate, - a reconstruction unit for generating an intermediate reconstruction image by use of said data set of X-ray projections including said corrected X-ray projection, and an evaluation unit for evaluating said correction by determining a quantitative measure of inhomogeneity in said intermediate reconstruction image and for optimizing said correction by iteratively repeating said correction using an adjusted estimation parameter determined by use of said quantitative measure until a predetermined stop criterion has been reached.

A reconstruction apparatus according to the invention is defined in claim comprising:

Corresponding methods are defined in claims 15 and 16. The invention relates also to a computer program which may be stored on a record carrier as defined in claim 17. Preferred embodiments of the invention are defined in the dependent claims.

The invention is based on the idea to iteratively optimize the performance of arbitrary (preferably software-based) methods for artifact estimation, in particular scatter estimation. This can be achieved by defining a quantitative measure of inhomogeneity caused by artifacts, in particular (in case of artifacts caused by scatter) a quantitative measure of residual cupping, in the reconstructed volume, and using this measure for repeated evaluation of the efficiency of artifact correction, each time followed by a corresponding adjustment of a parameter in the artifact compensation process.

To construct such a measure that quantifies the degree of performance of the employed compensation method, different options are available. In one embodiment of the invention, for volume images of the human body, the main part of the voxels shows attenuation values similar to that of water, and therefore, inhomogenities, such as cupping which means slowly varying inhomogenities leading to artifacts which are most disturbing in the center of projection images, is most easily visible at medium contrast levels around the attenuation value of water. Hence, starting point is preferably the selection of a contrast window of interest that comprehends the cupping. Then, upper and lower thresholding of the attenuation values in the reconstructed volume is performed. For the quantification of the inhomogenities such as cupping, statistical properties of the voxels belonging to the chosen window, such as standard deviation, may be used. Alternatively, a polynomial or another function may first be fitted to the data within the selected window, or strong low-pass filtering may be applied, followed by another step where the standard deviation or an appropriate curvature measure is deduced from the resulting profile. Having defined the quality measure for artifact correction, a simplex algorithm can be used to optimize the performance of any artifact compensation method. Either one of the internal parameters of the employed compensation method may be used for optimization, or the result of the method (e.g., the estimated spatial profile of artifacts in each projected view) may be globally scaled with an additional factor, which is then adjusted such that the quality measure is maximized. Such a factor would then account for errors of the artifact estimation, which may be either of principal (physical simplifications of the correction method) or algorithmic (implementation, discretization) origin.

The proposed optimization procedure can be fully automated, not requiring any user interaction. Iterative reconstructions can be performed with coarse resolution to

keep computational effort reasonably low. The method is principally not only applicable to C-arm based volume imaging but also to multi-line spiral CT where, however, the amount of scatter is lower and anti-scatter grids are more effective.

In a preferred embodiment the apparatus is particularly adapted for scatter correction, wherein: said estimation unit is adapted for estimating in an X-ray projection the amount of scatter present in said X-ray projection using at least one estimation parameter, said correction unit is adapted for correcting said X-ray projection by subtracting the estimated scatter from the projection data of said X-ray projection, and - said evaluation unit is adapted for evaluating said correction by determining a quantitative measure of residual cupping in said intermediate reconstruction image and for optimizing said correction by iteratively repeating said correction using an adjusted estimation parameter for scatter estimation determined by use of said quantitative measure until a predetermined stop criterion has been reached. It is to be noted that the suggested, although mainly proposed for improvement of scatter correction, is not limited to that application. Alternatively, it can instead be used to optimize performance of computerized beam hardening compensation or truncation correction, or to optimize gain normalization, because beam hardening, truncations and a wrong gain factor cause artifacts such as cupping, too. Therefore, in case of multiple software-based corrections, the correction that is iteratively optimized should be carried out last.

Thus, in a particular embodiment the apparatus is adapted for correction of a truncated projection, wherein: said estimation unit is adapted for estimating in an X-ray projection the degree of truncation present in said X-ray projection using at least one estimation parameter, said correction unit is adapted for correcting said X-ray projection by extending said X-ray projection using an extension factor or another method making use of said truncation estimate, and said evaluation unit is adapted for evaluating said correction by determining a quantitative measure of cupping or capping in said intermediate reconstruction image and for optimizing said correction by iteratively repeating said correction using an adjusted estimation parameter for truncation estimation determined by use of said quantitative measure until a predetermined stop criterion has been reached.

Further, in another particular embodiment the apparatus is adapted for correction of gain normalisation of projection data, wherein: said estimation unit is adapted for estimating a suitable gain factor in an X-ray projection using at least one estimation parameter, - said correction unit is adapted for correcting said X-ray projection by normalizing said X-ray projection using said gain factor, said evaluation unit is adapted for evaluating said correction by determining a quantitative measure of cupping or capping in said intermediate reconstruction image and for optimizing said correction by iteratively repeating said correction using an adjusted gain factor for gain normalization determined by use of said quantitative measure until a predetermined stop criterion has been reached.

The proposed iterative optimization of computerized artifact compensation by means of quantification of inhomogenities such as residual cupping, or by using a similar figure of merit, guarantees optimal image quality for the reconstructed volume as defined by the quality measure. The suggested solution will significantly improve removal of inhomogenities such as cupping caused by artifacts in reconstructed volumes and therefore allow for enhanced simultaneous visibility of low-contrast objects.

The invention will now be explained in more detail by use of exemplary embodiments illustrated in the accompanying drawings in which Fig. 1 illustrates the impact of scatter,

Fig. 2 shows a block diagram of a reconstruction apparatus according to the present invention, Fig. 3 schematically illustrates an artifact correction apparatus according to the present invention,

Figs. 4 to 6 illustrate an embodiment for estimating the effects of artifacts according to the present invention, and

Fig. 7 illustrates the results obtained by use of the method according to the present invention.

Before the invention will be explained in more detail by way of embodiments the impact of scatter and the generation of cupping artifacts caused by scattered radiation

shall be illustrated by way of Fig. 1. While the theory of computed tomography (CT) reconstruction assumes that all photons are either absorbed in an examined object or reach the detector directly, the largest amount of attenuation is, in fact, not caused by absorption but scatter. Therefore, a considerable amount of scattered photons reaches the detector on a non-straight way as can be seen in Fig. Ia.

As shown in Fig. Ib the background signal caused by scattered radiation is generally relatively homogeneous, i.e. especially slowly varying, but its amount is particularly significant. The portion of the total signal intensity caused by scattered radiation can - without anti-scatter grids - amount up to 50% or more. As can be seen from the profiles shown in Fig. Ib the relative error is largest for the total signal in the middle of the attenuation signal. Consequently, the relative error is also largest in the middle of the reconstructed object as shown in Fig. Ic where at the bottom the typical effect of cupping can be seen. For instance, for the head deviations up to -150HU below the correct grey value can be found. Thus, the problems caused by scatter induced artifacts are that scatter impedes the absolute quantification (HU), affects the visibility of low contrast structures and creates problems for further image processing.

Fig. 2 schematically shows the general layout of a reconstruction apparatus according to the present invention. By use of a data acquisition unit 2, for instance a CT or X- ray device, a data set of X-ray projections of an object 1, i.e. a patient's head, is acquired. The acquired data set is generally stored in a memory such as a hard disc of a server in a clinical network or another kind of storage unit of the work station further processing the acquired protection data. Before high-resolution reconstruction images are generated by a reconstruction unit 5 it is foreseen according to the present invention that an artifact correction is carried out by use of an artifact correction apparatus 4 which will be explained in more detail below. The corrected X-ray projections are then used for reconstructing a high resolution reconstruction image for subsequent display on a display unit 6.

Fig. 3 schematically illustrates the layout and the function of an artifact correction apparatus as proposed according to the present invention. In this figure more details of the artifact correction unit 4 shown in Fig. 2 will be illustrated.

The idea of iterative scatter correction is to first apply a generally ambiguous scatter correction algorithm on the available X-ray projections. Thus, the input to the scatter estimation and correction unit 41 can be a data set of original X-ray projections 10 or, in

order to save time required for the subsequent iteration, a reduced data set 11 of less X-ray projections obtained by spatial and/or angular subsampling of the original projections 10.

The estimation and correction carried out by unit 41 is generally incomplete, or the parameters of the correction are incorrectly set in the initial run of the iteration. For instance, in a preferred embodiment, a scatter correction is applied based on the scatter fraction (SF) relative to the minimum measured detection value per projection. Said embodiment can be implemented as follows: search for the minimum measured value per projection (optionally after strong low pass filtering); - a fixed percentage (for instance initially 50%) is assumed to be the minimum value of constant scatter background per projection direction; said constant scatter background per projection direction is subtracted from the projection data of that projection; and the above mentioned fixed percentage, which is called a scatter fraction (SF), is the parameter of this embodiment of the scatter correction which shall be optimised.

The corrected projections are then used for reconstructing an intermediate reconstruction image, for instance using the Feldkamp-David-Kress algorithm for cone-beam filtered back projection in reconstruction unit 42. Said reconstruction can be very coarse having a low resolution in order to be very fast. In an evaluation unit 43 the intermediate reconstruction image will then be evaluated in relation to the scatter-cupping artifact to be removed. Said evaluation will be explained in more detail below with reference to the subsequent figures. Depending on the result of said cupping evaluation the one or more parameters of the scatter correction (in the above explained example the scatter fraction SF) will be adapted and the next iteration starts. After a stop criterion has been reached in the iteration, for instance if a predetermined number of iterations has been run or if no further optimisation has been achieved during the last run(s), the final scatter correction using the adapted parameters (in this example the optimised scatter fraction) is finally applied on the complete data set of high resolution projections in reconstruction unit 5, and a time consuming high resolution final reconstruction is carried out to obtain final full resolution reconstruction image for display on a monitor 6.

An embodiment of a preferred cupping evaluation method will now be explained in more detail. The preferably proposed cupping evaluation method comprises two main steps:

a) selection of appropriate voxels, which can be used for measuring/quantifying cupping and b) calculation of a cupping measure based on said selection.

A large part of the human body consists of water or water-like tissue so that cupping can be best evaluated by use of grey value changes of this water-like tissue caused by scatter. This shall be achieved in the following based on threshold values. The threshold values themselves shall be determined based on histograms. For this purpose, initially a histogram over all voxels of the reconstructed volume is prepared as shown in Fig. 4a. For this histogram it is sufficient to select the range of the values from e.g. 50% up to 200% of the expected attenuation value of water. In Fig. 4b 40 of 64 slices of the reconstructed coarse volume to be evaluated are shown (displayed with a relatively tight HU window).

For the determined histogram values a Gaussian-distribution is fitted as also shown in Fig. 4a. Said Gaussian-distribution comprises an underlying uniform distribution according to the method of least squares using a non-linear least-square-optimisation algorithm (e.g. Levenberg-Marquardt algorithm). Thus, a mean value and a standard deviation of the Gaussian-distribution, a scaling factor for adapting the amplitude and optionally the amplitude of the uniform distribution are obtained.

Based on the result of this fit two pairs of thresholds are now selected: a) for a small range: the mean value +1 standard deviation; b) for a wide range: mean value +2 standard deviations.

Based on these ranges the voxels, which shall be used in the evaluation, are selected according to the following scheme: a) all voxels of the small range are used directly; b) all voxels of the wide range are used, as long as non of their direct neighbours (preferred in a neighbourhood of 6 neighbours) are not in said wider range.

The intermediate result is shown in Fig. 5a. The figure shows again the reconstructed coarse volume, where to all voxels, which have not been selected, the darkest "black" value has been assigned. Due to the partial volume effect there is a still a number of outliers at the edges to water and bones. For the elimination of outliers a polynomial will then be fitted to the selected voxels. Preferably, a second order polynomial is fitted in three coordinates, however without mixed terms (xy, xz, yz). Thus, coefficients for x ² , y ² , z ² , x, y, z and one constant have to be determined. This is done using the method of general linear least squares.

Such a polynomial fit is illustrated as an example in the graph shown in Fig. 6a. In Fig. 6b all selected ("valid water") voxels, which can be about 50.000 voxels, have been replaced by the result of this fit. The image is thus very smooth; black areas are voxels, which have not been selected. It has been found that this fit is reliable and reproducible for the same object.

Subsequently the cupping measure has to be determined. In a preferred embodiment the standard deviation of the result of said polynomial fit, evaluated at all selected voxels (as described above), is used as cupping measure. In further embodiments other cupping measures can be used, such as: - the maximum of the result of said polynomial fit, evaluated at all selected voxels, minus the respective minimum; the respective maximum standard deviation of all slices; or the respective maximum difference between maximum and minimum of all slices. All usable measures are thus reliable and reproducible.

The main idea of the present invention is therefore, to use any kind of scatter correction scheme, which may also be simple and imperfect, and to iteratively adapt said scatter correction scheme to achieve a homogenous reconstructed volume. Besides the above described scatter fraction based scatter correction scheme other scatter correction schemes can be used as well. One class of scatter correction schemes are so-called "self standing methods". Such "self standing methods" are methods, which exclusively use parameters which are adapted by the iterative optimisation method. Such parameters are, for instance: a global constant (i.e. the identical scatter background is subtracted from all projections); - a global scatter fraction (as described above in detail); or generally the one or more parameters of a polynomial which describes scatter in the projections.

Another class of scatter correction schemes are so-called "complementation of single scatter methods". Using models or coarse voxel reconstructions it is possible with reasonable efforts to calculate the single scatter, i.e. quanta which have only been scattered once. However, multiple scatter is still missing, which amounts up to 50% and more of the total scatter. Using the iterative scatter correction proposed according to the present invention, parameters can be optimised which allow to estimate the multiple scatter from the calculated single scatter.

Figs. 7a and 7b show the result of scatter correction achieved by the present invention. As can be seen from the reconstruction image in Fig. 7a the influence of scatter is so strong that the main part of the head can not be recognized any longer since the grey values are no longer in the viewed grey value range but have lower grey values. This is no longer the case after application of the scatter correction method according to the present invention. The reconstruction achieved therefrom is shown in Fig. 7b. In the reconstruction image cupping is hardly visible, the image is homogeneous and can easily be viewed in the desired grey value range. Thus, the mean difference to an ideal image in the tissue region is less than 20HU. While the invention is mainly applied for scatter correction, other applications of the general idea of the invention are possible. For instance, the invention can be applied for adaption of a projection extension factor or for adaption of a gain factor. The method used for adaption of a projection extension factor uses essentially the same steps as described above with reference to Fig. 3 (where the "scatter estimation and correction" is now replaced by a "projection extension, in particular according to "Lewitt et al." (see above) by using an extension factor".

The adaption of a global gain factor on the projection values equals an adaption of a global summand on the line integral values (after taking the logarithm). Thus, in the method illustrated in Fig. 3, the steps of which are essentially used for gain factor adaption the "scatter estimation and correction" is replaced by "applying a global constant" which is to be adapted.