Suppression of residual gradients in the flat-field corrected images

Flat field correction, based on using the open beam and the dark field images, is a basic method used for correction of the projection images acquired in computed tomography process. The correction suppresses the stationary noise of the detector, caused by non-uniform gain and offset of each pixel, as well as the non-uniformity of intensity caused by the non-flat intensity profile of the incident X-rays beam. In an ideal case, the corrected image should have a constant intensity in the open-beam areas, i.e., everywhere outside the investigated object. However, this is not always true due to several factors, such as drifts of the spot position, wearing out of the target (pitting), etc. In such a case, the background in the flat-field corrected image is not really “flat” or constant. Particularly in radiography, this is an undesired effect, because the intensity values are also influenced in the object area, leading thus to a possible misinterpretation of the image (even for a constant thickness area, the intensity distribution looks like if the thickness was varying). In this paper, a method is described how to compensate this effect in rather an effective way.


Introduction
The basic idea of the flat-field correction is rather simple: the value of each pixel in the image to be corrected is divided by the value of the corresponding pixel in the open-beam image, which is obtained iluminating the detector without any object in the field of view and using the same X-ray and detector settings.Hence, the pixel value in the corrected image is a number between 0 (completely dark) and 1 (full intensity).For technical reasons, more concretely the methods of the representation of numbers in computers, this value is usually multiplied by a suitable number and the result is saved in an integer format (the multiplication constant is equal to the maximum encodable number in the given format in an ideal case, so for instance 65535 for commonly used 16 bit unsigned integer format).If the used detector exhibits offset (which is usually different in individual pixels), it is necessary to generate the offset image, i.e., with the X-ray off, and subtract this image from the projection as well as for the open beam image before their division.This is very common in case of scincillation detectors, but it might not be needed in case of photon count detectors.In an ideal case, the described procedure supresses the stationary noise of the detector and the non-homogeneity of the intensity accross the detector.However, there are several factors that make the correction of the images imperfect in practice.Thus, several more advanced methods have been reported for improvement of the correction quality [1][2][3][4].A very common practice is to use some area which is outside the investigated object in all the images and use it for the equalization of the image intensities, which can be influences by the variations of the intensity of the beam.Even when using this kind of equalization, we observed gradients in the flat-field corrected images in some of our measurements.This effect is caused by the change in the intensity profile or pattern in the field of view over the measurement time, i.e., the open beam image acquired before the measurement and after the measurement are not identical.This behaviour can be due to the mechanical and electrical temperature drifts, leading to spot movement, and in case of higher power, also to the wearing out effects on the target.The correction open-beam images, however, have to be acquired at once before or after measurement, as their acquisition during the tomography would be too complicated.As a consequence, the gradients in the corrected images are probable for those acquisitions that were performed a significant time period before or after the acquisition of the open-beam images.

Mathematical model
Probably the most correct way how to compensate the described effect would be to analyze the shift in the intensity profile based on the comparison of regions outside the object in particular projections with respect to the open beam image, subsequently to shift the intensity profile and make the flat field correction.However, even neglecting the computational demands of such a procedure, there is a big obstacle caused by the stationary noise of the detector.Therefore, it is much more plausible to perform the compensation of the gradients in the (flat-field) corrected images, where the stationary noise is already suppressed.
The intensity of a (monochromatical) X-ray radiation behind a (homogeneous) object is given by the Beer-Lambert law as Here, I 0 is the intensity of the X-ray before entering the object, d is the thickness of the object, or the attenuation path,  is the attenuation coefficient of the corresponding material.
In a detector pixel, the radiation causes a response V given as where c is the conversion coefficient, V DF is the offset of the pixel value (DF standing for "dark field").
For the case of illumination without the object ("OB" for open beam image), the pixel value is obviously The recovery of the relative I value, or the ratio I/I 0 = e -µd , which is important, as it directly reflects the attenuation and thickness in the corresponding ray path, is then possible knowing V DF and V OB calculating Operation described by ( 4) and performed in each pixel is so-called flat-field correction and is massively used as a basic image correction method in radiography and computed tomography.
It is important to mention that this operation leads to a significant suppression of the stationary noise, or image noise caused by non-uniform response of the individual detector pixels, as can be shown by the following consideration.Let there be two adjacent pixels P 1 , P 2 on the detector with identical incident intensity I and the conversion coefficients c 1 , c 2 and offsets V DF1 , V DF2 , respectively.Using (2), their values will be Now it is very important to realize that even for the same incident intensity of the X-ray, the responses of the pixels can be different, depending on the particular values of the conversion coefficients c 1,2 and the offsets V DF1,2 .This fact applies for the whole detector area and manifests itself as an image noise, often called "stationary noise".This noise is induced by the detector and it does not behave as a random noise over time (the noise patterns are rather characteristic for the particular detector)see Figure 1.
The OB values in the pixels P 1 , P 2 will be, using (3), Applying (4) to the both pixels gives From ( 7) and (8), it can be seen that with the knowledge of appropriate offset and open beam value, the stationary noise can be suppressed and the relative intensity value can be easily recovered.
In practice, two correction images are needed to be exposed, the offset or dark field image and the open beam image.Usually, these correction images are made as an average out of tens to hundreds of expositions to reduce the X-ray related time-random noise.Technically, the correction images are exposed before or after expositions of the object, as for the open-beam images, the object has to be removed.Under the assumption that the offset and open-beam images are stable over time, it does not make any difference if they are exposed before or after the exposition of the set of images with object and the above described procedure can be used to recover reliably the relative intensity, being e -µd behind the object and 1 around it (in the air).However, we have observed that whilst the offset images are very stable over time, the open-beam images can differ significantly.In the particular case of a reflective X-ray tube, higher power and long time measurement (large amount of expositions of the object for CT etc.), the intensity pattern of the beam can drift considerably.The open-beam image provided by a cone-beam tube onto a scintillation detector exhibits a notable center of the maximum pixel value, with a non-linear decay with increasing distance from this "intensity center".We have observed that the intensity center of the beam can vary with the time, see Figure 1.We suppose that several causes for the intensity center drift may apply: the centering currents drift, mechanical drifts in the tube geometry and pitting effects of the target.In any case, the open beam image that has different intensity distribution than the one which would be visible during the expositions of the object does not yield desired results when used as a correction image, because it leads to the gradients in the background of the corrected images, see whereas the pixel value in the shifted-intensity open beam image, which is after the measurement used for the correction, is When the flat-feld correction is performed after (4), with d = 0, it gives in the investigated pixel: Equation (11) bears a very important finding that in an image corrected with a shifted open-beam image, the pixels outside the object do not have the correct value of 1, but at the same time, their value is not dependent on the offset and conversion constant of the individual pixel, i.e., even with a shifted open-beam image, performing the flat-field correction leads to the suppression of the stationary noise.Moreover, it is possible to determine the value of the k coefficient in each pixel knowing its flat-field corrected value, which equals 1/k (if corrected with unshifted open-beam, the flat-field value has to be 1).
In the area of the object, the situation gets more complicated.The intensity behind the object is given by (1).Considering the correct value of the intensity I 0 in one, i.e., the value of unshifted open-beam, which is I 0,c , the pixel value behind the object is Performing the flat-field correction with shifted-intensity open-beam image gives Thus, the pixel value in the object area is not dependent on the pixel offset or conversion constant, ie., the stationary noise is suppressed in the object projection, and the intensity value is again influenced by the coefficient k.However, in this case it is not possible to easily calculate the value of the k coefficient, as the − value is not known.
Nevertheless, the k value can be estimated, if large-enough open-beam areas exist around the investigated object.In such a case, the k values in a particular row and outside the object can be calculated using (11).Two sets of k-values will be obtained in this way, one set for the open-beam area to the left of the object and one set for the open-beam area to the right of the object.We can assume that the function of the k coefficient vs. pixel in the row behaves reasonably, i.e., it is smooth and does not exhibit any steps.A real form of the flat-field corrected values multiplied by 30000 in an image, which was corrected with intensity-shifted open-beam image, can be seen in Figure 2 (right).The k-coefficient is the reciprocal value of the flat-field corrected value, so that Figure 2 shows that our assumption regarding the k-coefficient function in a row is reasonable.Using the two sets of k-values outside the object, we can fit the k-values function over the whole row with a polynomial function using the least-square method (LSM) and a reasonable order of the fitting polynomial.In practice, we observed that 2 nd order polynomial is high-enough for a good performance of the method.Repeating the calculation for every row, we get a corrected image without a notable gradient and with recovered intensity values in the object area.

Experimental verification
Figure 2 shows an image after the flat-field correction, with a notable gradient in the background.The situation gets even more obvious looking at the intensity profile in one row (Figure 2, right).As described above, the proposed method for elimination of this effect requires fitting of the row values outside the measured object with a polynom.In this way, a model of the gradient in the specific row is created, including the part behind the object.In the above text, this corresponds to finding the k coefficients, corresponding to the ratio between the pixel value in the open-beam areas of a flat-field corrected image and the value which would be reached with a correct open-beam image.Hence, ratio of the required open beam constant value (30000 in our case) and the polynom value is calculated in each pixel of the row.The intenstity value in each pixel is then multiplied by this coefficient for every pixel in the row.Consequently, the intensity value outside the object gets close to the required open beam value (with the precision depending on the noise and the order of the fitting polynom) and the intensity values of the pixels in the region of the object are recalculated to the values that they would exhibit if the flat field image would not have the gradient effects.This procedure is performed for each row in the image.Figure 3 shows figure 1 after the reported correction using a 2 nd order fitting polynom.The improvement of the uniformity of the intensity profile outside the object is clearly visible.

Conclusion
Flat-field correction is a common method used for correction of the images in radiography and tomography.The offset, or dark-field, image and an open-beam image are used during the corrections of the set of exposed images.However, the intensity of the beam can vary during the acquisition of the projection images.For this reason, the equalization of the set of images is often made using a defined area, which lies outside the investigated object in all projections.Yet another effect can be observed which is not compensable in the standard way.Not only the intensity itself, but also the intensity pattern of the beam can vary over time.This can lead to undesired gradients in the background of the image, which is supposed to be uniform, and also affects the pixel values in the object.In this paper, we proposed and verified a method, which compensates these effects.It is required that there are two open-beam areas around the investigated object, one on its left side and one on its right side.This requirement is easily met in the case of the tomographical projections.The background is modelled by a fitting polynom in each particular row using the pixel values in the open beam areas.After that, correction coefficients are calculated for each pixel in the row and each pixel is corrected.Consequently, the method not only compensates the influence of the variations in beam intensity over time, but also the effects of the variations of the intensity pattern.The method was successfully experimentally verified and it has been involved in an in-house image processing software.

Figure 2 .
Let us introduce a coefficient k as the ratio of the incident intensity I 0,s in the open-beam image with the shifted intensity pattern and the I 0,c in the correct open-beam image; I 0,s = k.I 0,c .Obviously k is different in cone-beam area.Now, let us consider pixels in the area outside the object, i.e., where the flat-field corrected pixel value should be 1.The pixel value with the correct open-beam image, or the image acquired during the exposition with the object, is

Figure 1 :
Figure 1: Two open-beam images made with the same tube in different time.Notice the stationary noise patterns induced by the detector technology.Left: maximum intensity near to the center of the beam; right: maximum intensity shifted significantly to the left.

Figure 2 :
Figure 2: Left: gradient in the flat-field corrected image; right: the intensity profile in the row 500

Figure 3 :
Figure 3: Left: flat-field corrected image with equalized gradient (2 nd order polynomial); right: the intensity profile in the row 500 of the equalized image