# Polychromatic CT Data Improvement with One-Parameter Power Correction

Authors: Anastasia S. Ingacheva; Marina V. Chukalina
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Standard approaches to tomography reconstruction of the projection data registered with polychromatic emission lead to the appearance of cupping artifacts and irrelevant lines between regions of strong absorption. The main reason for their appearance is the fact that part of the emission with low energy is being absorbed entirely by high absorbing objects. This fact is known as beam hardening (BH). The procedure of processing projection data collected in polychromatic mode is presented; it reduces artifacts relevant to BH and does not require additional calibration experiments. The procedure consists of two steps: the first is to linearize the projection data with one-parameter power correction, and the second is to reconstruct the images from linearized data. Automatic parameter adjustment is the main advantage of the procedure. The optimization problem is formulated. The system flowchart is presented. The reconstruction with different powers of correction is considered to evaluate the quality reconstruction.


---

## Body

## 1. Introduction

X-ray computed tomography (CT) is a nondestructive technique for studying and visualizing interior features of an object without its physical destruction. Nowadays CT has become one of the main noninvasive diagnostic methods in medicine, scientific research, and industrial applications. X-ray computed tomography in medicine is used for 3D visualization of the inner structure of patient’s organs for studying diseases and their course controlling [1, 2]. Dental tomography is used for planning implantation and prosthetics [3, 4]. In industry, CT is used for quality control of products [5, 6]. In progressive scientific research, computed tomography serves as the method for examination of pressure retaining nuclear plant components [7].The task of CT consists of measuring projection data from multiple orientations of an object placed on a tomographic set-up and subsequent processing of the projection data in order to obtain a solution of the inverse tomography problem. Many tomography reconstruction algorithms are based on the assumption that attenuation of the X-rays in a material linearly depends on the length path of probing radiation. The main, becoming classical, types of the algorithms are based on integral or algebraic approaches like Fourier Back Transform (FBP) or Algebraic Reconstruction Technique (ART) [8]. However, nowadays algorithm development work is in progress due to the fact that the objects with nonstandard attenuation properties [9] have to be measured or nonstandard experimental conditions are required (low-dose tomography) [10]. The works are focused also on an improvement of the algorithms properties like convergence [11, 12] or the quality of the reconstruction [13].Projection data which satisfy linear model can be obtained by usage of monochromatic X-ray probe only. However, X-ray laboratory sources produce polychromatic beams; therefore, data collected in lab conditions are not linear. Classical reconstruction algorithms to the polychromatic data lead to different artifacts and distortions such as “cupping artifact” or artificial strips and lines between regions with strong absorption [14]. The main reason for their appearance is that the part of the radiation with low energy is fully absorbed by the regions containing high-Z materials. This effect is known as beam hardening (BH) and was described by Brooks and DiChiro [15]. The authors show that data nonlinearity in tomographic reconstruction can exceed 10% of average value, therefore noticeably limiting the usage of CT in tasks where we need to estimate not only the global structure but also the absorption value in every point. The suppressing of artifacts due to the polychromatic nature of the probing emission was the target for scientists for decades. As a result, nowadays there are a lot of methods for artifacts suppression. These methods can be generally divided into five groups: hardware filtering, usage of dual-energy scaling, preprocessing of projection data, reconstruction with iterative methods using the model of polychromatic radiation, and reconstruction result postprocessing.The first methods group is based on hardware filtration. The idea is to use metal objects to partially absorb the low energy radiation. Material and shape of the filter are highly dependent on the material of anode, current, and peak voltage of the X-ray tube. The main target of researches in this group of methods is to incorporate new filters [16, 17]. This hardware filtration reduces artifacts caused by beam hardening effects and also reduces the number of photons, which leads to bad signal/noise ratio and causes CT image degradation.The second group includes methods for tomography set-ups with two X-ray tubes [18–20]. The usage of two energies allows obtaining two sets of the reconstruction results. If a device with two energies is used, then two reconstruction results are calculated: one for high energy and one for low energy. Thus, it is possible to select the image highly dependent on Z which gives us a clue on the sample material. However, this technique has its own drawbacks such as high complexity of realization and strong sensitivity to noise. Another drawback is that it is required to use X-ray radiation twice for the object which increases the scanning time and radiation dose received by the sample.The methods of the third group use preprocessing of the registered projection data until they transfer to reconstruction. The processing is based on the calibration curves usage. The test samples are measured to draw the calibration curves. First, [21] suggested the use of a polynomial function of power 2 and 3 with one parameter to approximate the curves. Then it was shown that if the sample has high absorbing inclusion, then a polynomial function of power eight or greater has to be used for approximation [22].Next group consists of iterative methods including the correction procedure in the core of reconstruction [23–25]. These methods meet different handicaps and limitations such as requirements for knowledge about the sample or beam spectrum. But the main drawback of these methods is the time of reconstruction.The methods of the last group propose to include the postprocessing of reconstructed data. Methods from this group do not change projection data and use classical reconstruction algorithms. They are focused on the removal of artifacts from the reconstructed images. One of these methods is based on the blind deconvolution method which is widely known in computer vision for correction of blurred images [26]. Another method from this group evaluates missing data using discrete cosine transform and antialiasing filter which preserves borders [27]. The main drawback of methods from this group is that all of them use image processing methods and do not take into account image origin. Therefore these methods can remove important information from the reconstructed image as well as adding nonexistent information to them.In the paper, we present the procedure of processing projection data collected in the polychromatic mode which does not require additional calibration experiments and reduces artifacts relevant to BH. The sections are organized as follows. After the optimization problem statement, we describe the software implementation. In the “Results” the proposed algorithm output is given. The verification of the obtained results is presented in the “Discussion”. In our research, we used the projections, collected with the tomography set-up at FSRC “Crystallography and Photonics” RAS.


## 2. Optimization Problem Statement

If we use X-ray monochromatic radiation with initial intensityI0, then signal P collected by detector cell s (S is the entire set of detector cells, and therefore s∈S) can be described with the law of Beer-Lambert as (1)Ps,α=I0e-∫Lμldl,where α is a rotation angle of the object and L is the straight path of the probe from the source to the detector through the object (dl is a step). Two coordinates (s,α) fully specify the position of the probe with zero thickness and the probing direction. Before the reconstruction of the projection data, there is a linearization procedure:(2)Rs,α=ln⁡I0Ps,α=∫Lμldl.Rs,α is an integral of function μ(l) along the direction L, which is orthogonal to vector n→=(cos⁡α,sin⁡α), so R is function that describes Radon transform. The Radon transform data is often called a sinogram because the Radon transform of an off-center point source is a sinusoid. Most of the reconstruction algorithms work with sinograms. Let us consider total intensity of all detector cells for rotation angle α:(3)Tα=∑s∈SRs,α.Tα from (3) equals summation of the function μ(l) over the entire probing area since the position of the probing beam depends on coordinate s, which takes all values from its range. Therefore the value of Tα is constant for every α∈Α, preserving the property of sinogram invariance relatively to the rotation angle. This is the property of Radon invariant.However, in laboratory tomography scanners typically use X-ray tubes with polychromatic beam spectrum (i.e., wide energy spectrum from the range [0, Emax]). Equation (4) describes the value registered by the detector in this case:(4)Ps,α=∫0EmaxI0EXEe-∫LμE,ldldE,where X(E) is a sensitivity function of the detector. Therefore, (1) has been changed with additional integral through energies. Since the local absorption function of the object μE,l also depends on energy E, the data registered in polychromatic mode cannot be represented as exact Radon transform of any object. Usage of the standard reconstruction algorithms on the data of this type leads to the appearance of the different artifacts. For decreasing the artifacts one has to find the transformation which will linearize data.(5)fPs,α=-∫LμldlWe use one power approximation [28].(6)f=ln⁡∫0EmaxI0EXEdEPs,αγ∫0EmaxI0(E)X(E)dE is easily estimated from the values in the “empty” detector cells containing nonattenuated initial beam.If the power valueγ is evaluated using the property of Radon invariant preserving, then the optimization problem is formulated as(7)1Α∑αϵΑ∑s∈SfPs,αM∑s∈SfPs,α-12→fmin,where M[·] is an operator of a mean value.


## 3. Software Implementation

Figure1 shows the flowchart of our proposed method of polychromatic data processing.Figure 1
Flowchart of our proposed processing method for polychromatic data.The suggested method consists of three steps: with the first step, we find optimal parameter value with adjustment parameter for linearization procedure (APLP) algorithm. With the second step, we correct polychromatic data to make it closer to linear model with linearization procedure (LP) algorithm. After the projection data is corrected, the reconstruction algorithm is executed as the third step. The detailed discussion of each of these three algorithms is provided below.To find optimal parameter value we execute APLP algorithm. The APLP algorithm takes projection dataP as an input. Projection data have three dimensions (N×M×K): N is the number of detector rows, M is the number of detector columns, K is the number of rotation angles. Pseudocode of this procedure is shown in Pseudocode 1.Pseudocode 1:Pseudocode of adjustment parameter for linearization procedure algorithm.
Algorithm 1 Adjustment parameter for linearization procedure (APLP) (1) procedure APLP(P) (2) Input: P is projection data size of N×M×K. (3) Output: γ~ is adjustment beam hardening correction parameter. (4)P′={pijz∈P:i=[0.45N,0.55N];j=[1,M];z=[1,K]}; (5) az←∀z∈[1,K]; (6) mi←0∀i∈[0.45N,0.55N]; (7) di←0∀i∈[0.45N,0.55N]; (8) Dγ∀γ∈[0,500]; (9)for  All i∈[0.45N,0.55N]  do (10)Pi′={pi′jz′∈P′:i′=i,j=[1,M];z=[1,K]} (11)for  All γ∈[0,500]   do (12)Sγ←LC(Pi′,1+γ/100); (13)for  All z∈[1,K]  do (14)az←∑j=1M‍sijz; (15)mi←1/K∑z=1K‍az; (16)di←1/K∑z=1K‍(az/mi-1)2; (17)Dγ+=diγ; (18)Dγ←Dγ/0.1N; (19)γ~←1+argminγDγ/100.On the initial step of the algorithm, we select part of the projection dataP′ which is within 5% distance range from the central projection. If our tomography task is two-dimensional, then the algorithm is executed for all projection data (one line). We search the parameter value γ~ in the range [1.00;6.00] with step equal to 0.01. Every line from the part of the projection data P′ is being corrected with current power value 1+γ/100. For this we used LP algorithm, which will be described in the next section. Then we evaluate Radon invariant for every angle mi, average Radon invariant, and the root-mean-square of the ratio between Radon invariant for every angle di and average Radon invariant mi. Optimal value is found using brute force search over the range Dγ. We use brute force search since there are no obvious reasons to think that dependence between the value and the Radon invariant is differentiable and has only one minimum. As a result of the APLP algorithm, we obtain an optimal parameter for correction projection data.LP algorithm is used for polychromatic projection data correction. Pseudocode of LP algorithm is shown in Pseudocode2. Input parameters for the LP algorithm are polychromatic projection data P and the parameter γ (power value for data correction). Each pixel of the projection data is corrected independently from the others; i.e., correction is performed pixel by pixel.Pseudocode 2:Pseudocode of linearization correction algorithm.
Algorithm 2 Linearization procedure (LP) (1) procedure LP(P,γ) (2) Input: P is projection data size of N×M×K, (3) γ is beam hardening correction factor. (4) Output: S is corrected sinogram size of N×M×K. (5) for  All pijz∈P  do (6)p~ijz←-ln⁡pijz; (7)sijz←(p~ijz)γ.We used the FBP algorithm for data reconstruction since it is well studied and widely used in both medical and industrial CT scanners. It can be noticed that the suggested polychromatic data correction allows usage of any reconstruction algorithm which works in the paradigm of linearized data. The suggested method will complement well the modification reconstruction algorithms focused on suppressing another type of artifacts, for example, artifacts arising from disturbances in the geometry of the CT set-up or the presence of high-level noise.


## 4. Results

Linear dimensions of sample are7×7×12mm. For data acquisition, we used the micro tomograph which was assembled in the “Crystallography and Photonics” Federal Research Center of the Russian Academy of Science [29]. We used an X-ray tube with Mo anode as a radiation source and a vacuum gauge for intensity stability reasons. The projections were collected with the XIMEA-xiRay 11 Mpix detector. The measurement mode parameters for the tube were 40 kV voltage, 20 mA current, and 5 s acquisition time per frame. A source to object distance was 1.2 m and object to detector distance was 0.05 m. Every rotational scan consisted of 400 projections with an angular constant step size of 0.5 deg. The detector pixel resolution was 9 μm. The geometrical configuration of the experiment is such that the X-ray beam can be considered as parallel since its angle of divergence is less than 1 deg. [30]. No filter or other optics was used to produce a monochromatic beam, so the work mode was polychromatic.A CT projection of the sample is shown in the top of Figure2. The blue line indicates the cross-sections chosen for the 2D reconstruction. The sinogram pertaining to the cross-sections is shown at the bottom of Figure 2.Figure 2
Sample projection (top) and sinogram (bottom).We applied the method of polychromatic data processing described above to the registered data. Figure3 shows the behavior of criteria for different correction parameter values. It can be seen in Figure 3 that the optimal value of the parameter is 1.83. The found value agrees with the world publications on this topic, which suggest using the value 2.0 when weak absorption coefficient inclusions are presented in an investigated object and raising the value in case of presence of highly absorbing inclusions. With the found value, we corrected the projection data, thereby preparing its usage of the reconstruction algorithm.Figure 3
NRMSD dynamics depending on the parameter value.We used the FBP algorithm to the projection data. The software implementation for python language is presented in the Astra Toolbox library [31].Figures4(a) and 5(a) show the reconstruction result without correction while Figures 4(b) and 5(b) show reconstruction after the suggested data linearization. Related cross-sections are shown in Figures 4(c) and 5(c).Figure 4
Reconstruction result: (a) without the BH correction, (b) with the BH correction, (c) the respective cross-sections on row 520.
(a) (b) (c)Figure 5
Reconstruction result: (a) without the BH correction, (b) with the BH correction, (c) the respective cross-sections on row 780.
(a) (b) (c)From the cross-section on row 520 (in Figure4(c)), it can be seen that the hollow part inside the tooth has an absorption coefficient greater than the space outside the tooth in the reconstruction without correction (red line). Cross-section of the corrected data reconstruction (blue line) shows the absorption coefficient of the space outside the tooth, and the hollow part inside the tooth is the same; moreover, the absorption coefficient of the tooth itself remains unchanged.It can be seen by cross-section on row 780 (in Figure5(c)) that cupping effect artifact is brightly expressed in the reconstruction result without correction (red line). Cross-section of the corrected data reconstruction (blue line) shows that the cupping artifact completely suppressed, while the level of the absorption coefficient of the inner part remains unchanged.We used the simulation results to test the procedure under high absorption conditions. We have calculated the spectrum of 100 kV voltage X-ray tube with W anode as a radiation source (Figure6) by TASMIP Spectra Calculator [32].Figure 6
Simulated X-ray spectra.Cross-section of the Fe phantom used to calculate the projections in parallel scheme is presented in Figure7. The Astra Toolbox library [32] was used to calculate the projections.Figure 7
Cross-section of the Fe phantom.Phantom size is 250 pixels. Pixel size is 28μm. Number of evenly spaced projection angles is 180. 10 is a rotation angle value. Calculated sinogram is presented in Figure 8.Figure 8
Phantom sinogram.The behavior of criteria for different correction parameter values is presented in Figure9. The minimum of functional (7) is achieved at 1.29.Figure 9
NRMSD dynamics depending on the parameter value.We have used FBP for reconstruction. Reconstruction results are shown in Figure10. Figure 10(a) presents results without linearization procedure. The reconstruction result obtained after linearization procedure with γ=1.29 is presented in Figure 10(b). Two 125-line cross-sections for both cases are presented in Figure 10(c).Figure 10
Reconstruction result: (a) without the BH correction, (b) with the BH correction, (c) the respective cross-sections on row 125.
(a) (b) (c)After reliability inspection under high absorption conditions we used the procedure for data collected during the tomographic experiment. Steel alloy with nickel plated coating (test object) was scanned in cone beam geometry [33]. X-ray laboratory source with W anode was used. A source to object distance was 128.5 mm and source to detector distance was 228 mm. Every rotational scan consisted of 360 projections with an angular constant step size of 1 deg. The detector pixel size was 0.055 mm. Tomographic projection (top) and sinogram of central cross-section (bottom) are presented in Figure 11.Figure 11
Projection (top) and sinogram (bottom) of the central cross-section.We applied the method of polychromatic data processing described above to the registered data. Figure12 shows the behavior of criteria for different correction parameter values. It can be seen in Figure 12 that the optimal value of the parameter is 1.95.Figure 12
NRMSD dynamics depending on the parameter value.We have used FBP for reconstruction. Figure13(a) presents the reconstruction result without correction while Figure 13(b) shows reconstruction after the suggested data linearization. Related cross-sections are shown in Figure 13(c). As the object has nickel plated coating, two areas in Figure 13(c) illustrate the result of BH correction. The first one is near-boundary area (about 400 pixels and about 1200 pixels) and center area.Figure 13
Reconstruction result: (a) without the BH correction, (b) with the BH correction, (c) the respective cross-sections on row 800.
(a) (b) (c)Ring inside the reconstructed area is explained by defects in some pixels of the used detector [34], which is created by the group that provided the data.


## 5. Discussion

To verify the optimality of the found value, we suggest the following analysis. We illustrate the steps of the analysis procedure using data from the baby’s tooth as an example (Figure5). With the first step, we select the region of interest (ROI) on the reconstructed image. The value inside ROI should be homogeneous (see the created ROI mask in Figure 14). With the second step, we calculate the standard deviation and mean value within the region of interest in the reconstruction result.Figure 14
Mask.We made 2D reconstruction of the selected slice with parameters of the correction function in the range [1.0,2.8]. Then we calculated the standard deviation in ROI. The results are shown in Figure 15. Red point in Figure 7 shows the standard deviation for the optimal parameter value.Figure 15
STD of reconstructed value in mask.The analysis shows that the found parameter value is optimal. Experiments with real data have confirmed that our approach is applicable. It should be noted that the suggested method is sensitive to distortions that violate the property of conservation of the Radon invariant. This group includes such distortions as the thermal drift of the center of the X-ray spot or a weak signal-to-noise ratio. On the other hand, the method is resistant to such artifacts as the displacement of the center of rotation.The other way to evaluate the manifestation rate of the cupping artifact on the reconstructed ROI is to estimate the error based seminorm in theW21 space:(8)μfitW21=∫abμfit′2tdt1/2in ROI. To calculate the derivative we fitted the data to exclude salt-and-pepper noise from consideration. The colors used in Figure 16 correspond to the colors in Figure 5. To evaluate the potential gain of the cupping effect reduction, we set several different filter types (Figure 15) to substitute the standard ramp filter in FBP reconstruction [31].Figure 16
ROI for evaluation of the cupping artifact manifestation rate: (left-right, top-down) filters Hann, Hamming, Ram-Lak, Shepp-Logan, and cosine.Effect reduction of cupping artifact is around 3 times for each of them.


## 6. Conclusions

To demonstrate the decrease of the artifacts relevant to BH, we apply the proposed algorithm to the data collected with the tomography set-up at FSRC “Crystallography and Photonics” RAS, data provided by our colleagues from Joint Institute for Nuclear Research, Dubna, and simulated data. First set of data are tomographic projections of baby’s tooth. Object selection is driven by several reasons. The main reason is related to the fact that the presence of artifacts on CT images of teeth [35] is a serious problem when choosing a treatment regimen and when planning implantation procedures. The search is constantly in progress to find new methods to reduce artifacts. The second reason is due to the fact that the baby’s tooth has a cavity formed by nonresorbed roots. It enabled us to evaluate the result of the algorithms operation both in the cross-section passing through the cavity (Figure 4) and in the cross-section passing through the dense part of the object. The second set is the projection data of test steel alloy with nickel plated coating. Its high absorption is a main reason to include the specimen in the consideration. The obtained results allow for the conclusion that the proposed approach, which does not require additional calibration experiments, allows reducing BH artifacts on the reconstructed image. We made reconstruction by FBP for data corrected with different power values to verify the developed system. The solution of the optimization problem gets in line with minimum STD position.



---

*Source: 1405365-2019-04-11.xml*