Ultrafast X-ray computed tomography (CT) is an imaging technique with high potential for the investigation of the hydrodynamics in multiphase flows. For correct determination of the phase distribution of such flows, a high accuracy of the reconstructed image data is essential. In X-ray CT, radiation scatter may cause disturbing artefacts. As the scattering is not considered in standard reconstruction algorithms, additional methods are necessary to correct the detector readings or to prevent the detection of scattered photons. In this paper, we present an analysis of the scattering background for the ultrafast X-ray CT imaging system ROFEX at the Helmholtz-Zentrum Dresden-Rossendorf and propose a correction technique based on collimation and deterministic simulation of first-order scattering.
In many industrial fields such as chemical engineering, power engineering or oil and gas industry, multiphase flows play a major role in the efficiency and safety of the respective processes. Therefore, characterization and understanding of the multiphase fluid dynamics is an essential task for process analysis and optimization. The combination of sophisticated simulation techniques with advanced measurement techniques is most promising for a deeper analysis of multiphase flows. There is a wide range of invasive measurement techniques, such as optical and electrical probes and wire-mesh sensors , as well as non-invasive techniques, such as ultrasonic measurements, electrical tomography, gamma-ray tomography  and X-ray tomography . The advantage of non-invasive methods is—beside the fact that the flow is not influenced—the applicability at high temperatures and high pressures as well as for chemically aggressive substances. Further, a high spatial resolution and a good material penetration capability are provided in general.
X-ray tomography is based on measuring the attenuation of X-rays with initial intensity N0 along determined paths L. If one assumes monoenergetic radiation, the intensity N of an X-ray beam is linearly attenuated according to 1.1 where μ(l) is the distribution of the linear attenuation coefficient of the material along the beam path. The sum of the total attenuation of X-rays is caused by four different physical effects: the photoelectric effect, incoherent scattering (Compton scattering), coherent scattering (Rayleigh scattering) and pair production. The latter effect only appears for high energies far beyond the theoretical threshold of 1.022 MeV. The contribution of each effect to the total attenuation depends on the photon energy, on the atomic number of the atom interacting with the photon as well as on the atomic density. The scattering effect is a deflection of photons from their original path. Coherently scattered photons keep their energy, whereas incoherently scattered ones lose an amount of their energy during the scattering event. In tomography systems, scattered photons may be registered by the detector. Standard reconstruction algorithms such as ART and filtered back-projection  only consider attenuation along straight beam paths, as given by equation (1.1). Scattering produces an additional offset at each detector, which depends on geometry. This leads to an underestimation of the real attenuation value. Reconstruction of uncorrected data then leads to significant artefacts in the images such as streaks and cuppings , which hamper quantitative data analysis and may even lead to misinterpretations.
There is a wide range of methods which are able to reduce or even fully correct scattering artefacts. The most common hardware-based method is application of radiation collimation. Placing strongly attenuating materials as boundaries of the beam path between the X-ray sources or detectors prevents to some degree the registering of photons that originate from a scattering event outside the scanning area. Use of beam stop arrays [6,7] and beam hole arrays [8,9] as scatter estimation tools are effective approaches, as well. Due to geometrical circumstances, hardware-based methods cannot prevent the complete scattering, in general. Therefore, several algorithm-based correction methods were developed. Most of these algorithms take a priori knowledge about the material distribution of the scanned object into account. Monte Carlo methods  are very common for scatter correction because of their accuracy. These algorithms simulate the paths and interactions of a sufficiently large number of photons stochastically. Besides, there is a range of other stochastic , deterministic  and empirical  methods.
For our ultrafast X-ray tomography system, we chose a hybrid method consisting of collimation and deterministic simulation. In §2, the tomography system is introduced as well as the implementation of a detector collimation. In §3, a deterministic scattering simulation algorithm is presented, which is based on first-order ray-tracing. In §4, the performance of both methods for a test case of a water-filled column is analysed.
2. Ultrafast X-ray tomography system and experimental procedure
The ultrafast X-ray tomography system ROFEX at Helmholtz-Zentrum Dresden-Rossendorf, introduced in , is a two-dimensional imaging system with two simultaneously sampled horizontal measurement planes. It consists of an electron beam gun, which focuses an electron beam with an acceleration voltage of 150 kV onto an X-ray-producing target, which has circular shape and is positioned around the object to be measured. The detector system currently in use consists of two circular rings. Each ring has 432 CdTe detector pixels with a square active area of 1.3×1.3 mm. The detector pixels integrate the total photon flux without any further energy discrimination (current mode detection). The detector is arranged concentrically within the target with a small axial offset between the focal spot path on the target and the corresponding detector ring. The tomography scanning is realized by deflecting the electron beam in order to generate a rotating focal spot on the target. A sketch of a horizontal cut of the geometry is shown in figure 1a. In two-plane measurement mode, the focal spot alternates between the two target rings after every revolution. Cross-sectional images can be recorded with a frame rate of up to 8000 fps in single plane mode or 4000 fps in dual plane mode. The spatial resolution is around 1 mm. This allows a detailed imaging of a wide range of multiphase-flow scenarios.
Figure 1b shows the axial arrangement of the target and detector. The X-rays are emitted from the target in an isotropic way. In order to shield X-rays which are not directed into the measurement plane, three circular steel plates were installed below and above the detector rings. Due to the design of the scanner, the distance between the plates must be at least 5 mm in order to ensure an unrestricted space for the direct beam paths between target and detector ring. Note that the inner diameter of the steel plates influences the axial shielding. For large objects, the inner diameter of the steel plates must be reduced, which leads to a loss of the axial collimation efficiency.
Typical applications for X-ray tomography are investigations of multiphase flows in vertical columns or pipes . At low gas hold-up scattering and beam-hardening within the liquid phase cause significant cupping artefacts. Therefore, we chose a water-filled column as test case for the scattering analysis methods. The column has a wall thickness of 5 mm and is made of poly(methyl methacrylate). The inner diameter is 100 mm. The inner diameter of the collimator plates is 116 mm. All measurements presented in §4 were performed with a scan rate of 1000 fps. Each frame is reconstructed from a sinogram consisting of 500 projections. For image reconstruction, the filtered back-projection algorithm in connection with the Shepp–Logan filter is applied. The reconstruction program is self-implemented in C++. The resulting slices consist of 280×280 image pixels with a spatial resolution of 0.5 mm per pixel.
3. In-plane scattering simulation
As the collimation presented in §2 cannot prevent scattering within the measurement plane, we chose a simulation-based algorithm for in-plane correction. In , Freud et al. introduced a deterministic first-order scatter simulation algorithm. Based on a discretization of the object of interest, rays are tracked from a source spot to every voxel of the object. Then, the number of scattered photons towards each detector pixel is calculated. In this way, every possible path of photons which are scattered once is taken into account. Basically, simulation of higher order scattering is possible but causes an exponential increase of the computational effort. Further, the amount of higher order scattering is relatively small as shown in . Therefore, we limit the analysis to first-order scattering. Freud et al. validated their method with the Monte Carlo code GEANT4 . The deterministic algorithm provided quantitatively congruent results at a significantly reduced calculation time. Another advantage of this method is the high flexibility regarding the object geometry and discretization.
In the following, we outline the main parts of this method adapted to our tomography system. Based on an initially reconstructed slice without any corrections, a supposed material distribution is set. As the algorithm simulates the three-dimensional scattering, the image pixels of the slice are assumed as voxels, whose height should equal the height of the measurement plane, defined by the focal spot diameter and the height of the detector pixels. According to the numerical studies presented in , the voxel size should be of the same order of magnitude as the detector pixel size. Each of the resulting image voxels is assigned to a chemical element or compound. Then, the number Nsca,pm(Esca) of photons of energy Esca scattered with angle θ in voxel p and registered at detector pixel m is calculated as 3.1 with Φp(E) being the number of incident photons of energy E per unit area at the distance of the voxel from the source. This number is multiplied by the number of atoms Nat,p within voxel p and the differential cross section dσsca(θ,E)/dΩ, which depends on the kind of scattering. Up to here, the product specifies the number of photons scattered towards pixel m. By multiplication with the attenuation of all Nmat materials with attenuation coefficients μn and length ln along the path from the voxel centre to the detector pixel, one obtains the number of scattered photons. Finally, multiplication with the unit solid angle dΩpm between the voxel centre p and detector pixel m gives the number of registered scattered photons in pixel m. For incoherent scattering 3.2 with dσKN(θ,E)/dΩ being the Klein–Nishina differential cross section and S(x, Z) being the incoherent scattering function. For coherent scattering 3.3 with dσT(θ)/dΩ being the Thompson differential cross section and F(x, Z) being the atomic form factor function. The differential cross sections depend on the scattering angle and the additional functions depend on the atomic numbers of the elements contained in the voxel. The formulae and tabulated values are given in . Furthermore, Esca in (3.1) is given as Esca=E in the case of coherent scattering and Esca=ECompton≤E in the case of incoherent scattering. The energy ECompton of an incoherent scattered photon is given as 3.4 where mc2 is the electron rest energy.
The integrated energy of all photons scattered in voxel p and registered at detector pixel m is calculated as 3.5 with Emax being the maximum energy contained in the X-ray spectrum and p(Esca)∈[0,1] being the detection efficiency which is a detector-specific multiplier. For a complete scattering simulation of a single slice, the number of calculations of (3.1) and (3.5) is the product of the number of projections, the number of voxels, the number of detector pixels and the number of discrete energies of the X-ray spectrum.
(a) Analysis of the collimation
In figure 2a,b, the reconstructed grey value images are shown for measurements with and without collimation. The cupping artefact is not immediately visible but becomes obvious when analysing the centre profiles of the images (figure 2c). In order to quantify the cupping artefact without being disturbed by noise, the shown profiles are averaged over 360 equiangular distributed cuts. The grey values v are scaled such that the value of the column wall is vwall=1 and for the surrounding air vair=0. In an ideal image without artefacts, the water value would be videal=0.905. This value is determined as the grey value at the position of maximum curvature of the phase transition from column to water. The local extreme value of the water without collimation is at vnon-collimated=0.797 and with collimation at vcollimated=0.856. This gives relative maximum deviations regarding the air–water contrast of 4.1 and 4.2 This means that the collimation reduces the cupping artefact to 45% compared with the non-collimated case.
In further measurements, the water level was slowly lowered with constant velocity such that the column within the scanning plane was filled at the beginning of the measurements and emptied at the end. In figure 3, the axial cuts of the reconstructed image stack are shown for the non-collimated case (figure 3a) as well as for the collimated case (figure 3b). As the constant liquid draining velocity can be determined by the readings of both detector rings, it is possible to map the axial time scale of the stack onto a spatial scale. In figure 3c, the average grey values of the cross-sectional images are plotted over the height. The reason for the blurred transition from water to air is the decrease of scatter-causing volume with lowering surface. This also explains the non-symmetrical shape of the curves. In the collimated case, the transition from water to air is considerably narrower. The grey value level of water is entirely reached when the surface is 8 mm above the scanning plane. The level of air is reached when the surface is 4 mm below the scanning plane. In the non-collimated case, both levels are reached above and below 10 mm. This means that the collimation limits the scatter-causing volume of the object to a height of 12 mm.
(b) Analysis of the in-plane simulation
For the simulation of the in-plane scattering, the method presented in §3 was applied. The focal spot was assumed to be an isotropic source emitting a constant number of photons per unit solid angle. The measured spectrum of the X-rays (figure 4) was taken into account as distribution proportional to the emitted number of photons. Further, the simulated detection efficiency of the CdTe detector pixels which depends on the geometry and material of the detector pixels was considered as multiplier for the photons to be detected. We simulated the detector as the integrated energy over the spectrum calculated by (4.1). The initial image, on which the selection of the materials was based, was taken from the collimated measurement of the water-filled column (figure 2b). The height of the voxels was set to 1 mm. We simulated 12 voxel planes such that the space 4 mm below and 8 mm above the scanning plane was covered with voxels, according to the results of the previous section. In figure 5a, the integrated energy of scattered photons registered by the detector is plotted over the detector pixels, which are shadowed by the object. The amount caused by coherent scattering is slightly higher compared with incoherent scattering. The incoherent scattering influences a wider range on the detector. The ratio between scattered and non-scattered photons is shown in figure 5b. Behind the object, the ratio has its maximum at about 3.6%. A direct correction of the real detector readings with the simulated energy is not possible because of different sensitivities of the detector pixels, which are corrected by reference measurements. Therefore, we applied the scatter-to-signal ratio to the detector readings. This allowed the correction of detector readings irrespective of the true detector pixel properties.
In figure 6, the grey value profiles of the reconstructed image (cf. figure 2) with and without correction are shown. The extremum of the remaining artefact is at vsimulated=0.872. This gives a ratio of 4.3 compared with 5.4% in the non-simulated collimated case. The relative improvement is 33.3%. In summary, the collimation in connection with the deterministic simulation of the first-order scattering reduces the cupping artefact by 8.3% in relation to the water–air contrast. It is expected that the main contribution of the small remaining artefact comes from inaccurate estimation of the geometrical impact of the collimation. Further, the artefact also consists of contributions of higher order scattering and beam-hardening. The latter leads also to cupping, but was not taken into account by the simulations.
In order to reduce to the computational effort, it is advisable to simulate the scattering based on a reconstructed image with decreased pixel resolution as the calculation time is proportional to the square of the resolution. For this, we repeated the simulations with different pixel resolutions Ni=30, 50, 80, 120, 180, 280 and calculated the relative deviations of the scattering as 4.4 with fsca,i being the detector readings of the different types of scattering for pixel resolution Ni. Figure 7 shows that the deviations are below 2% for simulations with resolutions higher than 120 pixels. A simulation based on a single slice with 120×120 image pixels takes 221 s on an Athlon 2.80 GHz processor. For hydrodynamic studies with scans of several thousand images, it is advisable to simulate reference cases such as the water-filled column for measurements of two-phase flow experiments with low gas hold-ups. As the simulation consists of a huge amount of loop iterations, which are independent from each other, the algorithm is suitable for a parallelized implementation.
In this work, we present an analysis and correction technique of cupping artefacts in ultrafast X-ray tomography caused by scattering. The application of a collimator reduced the artefact significantly. Test measurements with lowering water surface allowed an estimation of the scatter-causing volume of the object limited by the collimator. This information was taken as input for a deterministic simulation of the first-order scatter within the non-collimated object space. As a result, we were able to reduce the cupping artefact from 11.9 to 3.6% of the water–air contrast. For further works, it is conceivable to handle the beam-hardening during simulation and accelerate the calculation with parallelized implementation.
Original reconstructed data sets are available at Zenodo (doi:10.5281/zenodo.16550).
The authors acknowledge the Helmholtz Association for support of the research within the frame of the Helmholtz Energy Alliance ‘Energy Efficient Chemical Multiphase Processes’.
One contribution of 11 to a theme issue ‘X-ray tomographic reconstruction for materials science’.
- Accepted March 13, 2015.
- © 2015 The Author(s) Published by the Royal Society. All rights reserved.