Cupping artifacts correction for polychromatic X-ray cone-beam computed tomography based on projection compensation and hardening behavior

Biomedical Signal Processing and Control 57 (2020) 101823

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Cupping artifacts correction for polychromatic X-ray cone-beam computed tomography based on projection compensation and hardening behavior Fuqiang Yang a,b,∗ , Dinghua Zhang a,b , Hua Zhang c , Kuidong Huang a,b,∗ a Key Laboratory of High Performance Manufacturing for Aero Engine (Northwestern Polytechnical University), Ministry of Industry and Information Technology, Xi’an 710072, China b Engineering Research Center of Advanced Manufacturing Technology for Aero Engine (Northwestern Polytechnical University), Ministry of Education, Xi’an 710072, China c School of Mechanical Engineering, Northwestern Polytechnical University, Xi’an 710072, China

a r t i c l e

i n f o

Article history: Received 2 July 2019 Received in revised form 24 October 2019 Accepted 8 December 2019 Keywords: Cupping artifacts Polychromatic Cone-beam computed tomography Hardening behavior

a b s t r a c t Based on the beam hardening behavior of polychromatic X-ray, the aim of this paper is to address and test an improved experimental and theoretical investigations of the distorted polychromatic projection data to generate high qualified imaging by reducing the cupping artifacts and noise in Cone-Beam computed tomography (CBCT). Since the hardening information has a large relationship with the span that the ray passing through the object, a method of obtaining the lengths that the ray across the binary image is proposed. As the knowledge, a new weighted and compensated correction model is constructed to primarily eliminate cupping artifacts. Finally, studies on polychromatic projection are verified by using pieces on CBCT system. Study results show effective suppression of cupping artifacts of polychromatic projection. For experimental, the root mean square error (RMSE) of the proposed method is reduced by 31.6 %, peak signal to noise ratio (PSNR), and universal quality index (UQI) are increased by 16.8 % and 12.8 % for phantom #Bushing, respectively; which also indicate better uniformity of the results for phantom #Wheel, where the RMSE is reduced by 30.6 %, PSNR and UQI were increased by 19.4 % and 15.9 %, respectively. The proposed algorithm is anticipated to find its utility in industrial nondestructive testing wherein the image reconstruction for the complex geometry and high energy X-ray imaging. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Cone-beam computed tomography (CBCT) plays a crucial role in both medical imaging and industrial non-destructive testing (NDT) with high efficiency and precision [1,2], which is equipped with Xray beams that emit photons of different energy and frequency. X-ray computed tomography became an important instrument for quality assurance in industry products as a non-destructive testing tool for inspection, evaluation, analysis and dimensional metrology. Thus, a high-quality image is required. Due to the polychromatic nature of X-ray energy in X-ray computed tomography,

∗ Corresponding authors at: Key Laboratory of High Performance Manufacturing for Aero Engine (Northwestern Polytechnical University), Ministry of Industry and Information Technology; Engineering Research Center of Advanced Manufacturing Technology for Aero Engine (Northwestern Polytechnical University), Ministry of Education, Xi’an 710072, China. E-mail addresses: [email protected] (F. Yang), [email protected] (D. Zhang), [email protected] (H. Zhang), [email protected] (K. Huang). https://doi.org/10.1016/j.bspc.2019.101823 1746-8094/© 2019 Elsevier Ltd. All rights reserved.

this leads to errors in the attenuation coefficient. This leads to a distortion or blurring-like cupping and streak in the reconstruction images, where a significant decrease in imaging quality is observed [3]. However, most reconstruction approaches do not properly examine the non-linear nature of the polychromatic Xrays. When polychromatic X-ray beams used during CT imaging pass through a patient, soft and low energy X-rays [4], which are of little importance in image formation, are preferentially absorbed to a great extent compared to high-energy photons. The consequence of this selective absorption phenomenon is the increase in the patient’s absorbed dose and non-linear increase in the beam’s average energy, which is the beam hardening effect [5,6]. As a result, the total attenuation of X-rays and therefore the associated logprocessed transmission data will no longer be a linear function of tissue thickness. Since the low-energy photons are attenuated more than the high-energy ones, the more transmission, the higher the average energy of photons when interacts with objects [7,8]. The most widely used Filtered back-projection (FBP) algorithms in CT recon-

2

F. Yang, D. Zhang, H. Zhang et al. / Biomedical Signal Processing and Control 57 (2020) 101823

struction assume a linear propagation model for the detected photons and as such fail to consider the non-linear beam hardening effect [9]. Consequently, the reconstructed images exhibit cupping artifacts and reduced CT numbers behind bony structures and streak artifacts around metallic objects [10]. These artifacts may mimic or obscure pathologic lesions, rendering the interpretation of CT images complicated and leading to equivocal findings. This leads to severe noise and other artifacts [11] and inhomogeneity which degrades image quality and makes diagnosis difficult or even impossible [12–14]. The poly-chromaticity of the X-ray beam therefore not only increases the patient dose but also induces beam hardening artifacts. Hence, there is a need for taking appropriate measures toward concomitant patient dose reduction and beam hardening correction. As a focused problem to be solved, several approaches for beam hardening have been proposed, which can be classified as dual-energy methods [15,16], statistical iterative methods [17–20] and sinogram inpainting methods [21–24]. Furthermore, specially designed filters known as beam shapers are used in CT scanners to shape the x-ray beam to deliver dose with the appropriate spatial distribution. As a consequence of pre-filtering, both beam hardening effect and patient’s absorbed dose are intrinsically reduced, however, as a compromise, statistical noise or quantum mottle is increased, which in turn impairs image quality and low-contrast detectability. Despite some good preliminary results that have been obtained, the dual-energy correction requires lengthy postprocessing and the high dose of radiation is a limitation. While, the iterative method requires prior information about the energy spectrum of the incident X-ray and energy-dependent attenuation coefficients of the materials [25,26]. Brabant et al. [27], employed an iterative reconstruction algorithm and took the effect of beam hardening into account during the reconstruction. To solve the issue of sinogram-consistency, Cao et al., [28] presented a beam hardening correction regression model to increase the dimensional measurement accuracy, which is trained with the simulated data afterward. Romano et al., [29] using a linearization procedure of the beam hardening curve applied after the reconstruction process to provide accurate results. Zhao et al., [30] developed a fast and accurate beam hardening correction method by modeling physical interactions between X-ray photons and materials for computed tomography (CT) imaging. The nonlinear attenuation process of the X-ray projection is modeled by reprojecting a template image with the estimated polychromatic spectrum. Sarkar et al., [31] utilizes discrete cosine transform-based missing data estimation and an edge-preserving smoothing filter for the suppression of beam hardening artifacts. Park et al., [32] proposes a sinogram-consistency learning method. to repair inconsistent sinogram by removing the primary metal-induced beam hardening factors along the metal trace in the sinogram in polychromatic computerized tomography (CT). The beam hardening is modeled and incorporated in the forward projector of the simultaneous algebraic reconstruction technique (SART). Various inpainting-based artifact reduction methods had been suggested, such as interpolation [33,34], Poisson in painting [35], wavelet [36], and total variation [37]. But preliminary studies had found limited performance, including the inaccurate interpolation who had introduced additional artifacts to the image and had deteriorated morphological information. Kyriakou et al. [38] proposed a beam hardening correction which replaced the corrupted projection data by interpolating, it had effectively reduced the artifacts, but the edge information lost. The traditional method uses a wedge-shaped phantom to obtain the curve. Huang et al. [8] constructed a new fitting function with satisfactory stability of curve shape through simplifying beam hardening data based on histogram statistics of traversing lengths of the ray. Zou et al. [39] presented algorithm utilizes image segmentation to obtain

the hardening object, and subsequently dilates segmented image to suppress dependence of the exactness of segmentation to correction result based on the modified TV algorithm. The previous methods are designed to reduce the higher-order artifacts induced by beam hardening in require material correction using calibration for global cupping correction. Inspired by the process and the properties of non-linear beam hardening of different materials, we aim to address the improved experimental and theoretical investigations of the distorted polychromatic projection data to generate high qualified imaging by reducing the cupping artifacts and noise in Cone-Beam computed tomography (CBCT) in this work. The method is evaluated using CT images, and the experiments show that the proposed method is attractive as a preferred quality for polychromatic X-rays imaging. The main contributions of the approach are presented as follows: 1) A new effective cupping artifacts correction for a polychromatic X-ray is proposed. The idea is to estimate monochromatic projections from polychromatic projections using a polynomial model and the prior knowledge on the relationship between the ray length and corresponding projection. 2) The artifact corrector is based on projection compensation and an analysis of the hardening characteristics. To reveal the hardening effects of different materials, the introduction and adaptive estimation of the tradeoff factor ␭ is exploited to the model for showing the characteristics and details of the correction from polychromatic to monochromatic X-ray in this work. The remainder of the paper is organized as follows. Section 2 provides a relevant description of notions and the approach for artifacts correction algorithm. Section 3 presents the experimental results and the performance comparison with the state-of-the-art methods. Finally, Section 4 remarks on the conclusion of this paper. 2. Methodology 2.1. Cupping artifacts and correction scheme As has been mentioned that when the polychromatic X-ray beams used during CT imaging pass through the body or the object, soft and low energy X-rays are preferentially absorbed to a great extent compared to high-energy photons. With the X-ray path increased, the proportion of high-energy photons increase and the average energy of ray arise [40]. Then, the attenuation coefficient  is no longer a constant, which shows a decreasing curve as the path gets longer. Clearly, this phenomenon leads to great defects that there is a great gap between the polychromatic projection and monochromatic projection within the same path length (As is showed in Fig. 1(a)). Therefore, the image reconstructed from distorted projection does not reflect the attenuation coefficient of the material correctly (As is showed in Fig. 1(b)), then the attenuation coefficient of the slice at the center is smaller than the edge, that is, a cupping artifact appears (As is showed in Fig. 1(c)∼(d)). The general cupping artifacts correction procedure for polychromatic X-ray based on projection compensation and hardening behavior is illustrated in Fig. 2. We first scan the object and obtain the polychromatic projections using the CBCT system. Projections are stored and filtered, and the projection noise is estimated by subtraction between the original and the filtered. In order to obtain the behavior of beam hardening in advance, the hardening curve that varies with the polychromatic projections is obtained based on the intersection of the ray within the voxel along the path. Then we use the prior knowledge to adjust the detected hardening model previous. Since the linear model has the problem of overcorrected, which would not reflect the correction effect completely,

F. Yang, D. Zhang, H. Zhang et al. / Biomedical Signal Processing and Control 57 (2020) 101823

3

Fig. 1. Relationship of attenuation coefficient, projection with the path length. (a) is the tendency of the attenuation coefficient based on path length, (b) is the tendency of the path length based on projection, (c) shows the reconstruction slice, (d) is the profile along the line.

Fig. 2. The workflow of the beam hardening for cupping artifacts.

the projection is modified by adding weighted tradeoff  to achieve accurate correction for different objects. When the polychromatic projections are corrected by the improved hardening model, the corrected projections are reconstructed to the slice. Details of the process of the cupping artifacts correction are discussed below. 2.2. The proposed algorithm As mentioned in Section above, the key point of the “hardening correction algorithm” is to seek and obtain the relation modeling between the attenuation and the length of the ray passing through the object. The pre-reconstructed slice which is binarized has been useful for exploring the related representations in images, so the sum of the lengths that the X-ray traversed in the voxels along different paths is computed. A notable characteristic of this approach is that the contribution of the pixel to the ray is able to simulate arbitrary relationships between the inner representations of 2D visual features and 3D poses. Since the length of a ray in a voxel appears as a matrix coefficient, it has been computed as the contribution of the pixel to the ray [41]. Fig. 3 shows the intersection between discrete pixels and rays in the 2D case. For the fitting data that does not require information in all ray paths, the X-ray paths in a few different angles are needed based on the location information from the binary image. As is showed in Fig. 4, the preliminary image is binarized first, here we use the projection distance minimization (PDM) method [42] to search for the optimal threshold. Then the ray paths are selected in sparse to obtain coefficient in the matrix, which is equivalent to the length. 2.2.1. Linear hardening correction (LHC) algorithm For polychromatic X-rays, beam hardening is an important factor that affects the quality of CT images [43]. In this part, the hardening behavior is discussed. The relationship between the polychromatic projection and the attenuation changes with the length of the ray that crossing the object, so the polychromatic projection is the line integral of path length and attenuation coefficient. While the attenuation coefficient is the function of the energy of Xray, which also is related to the length of the ray. Therefore, we use the capital P to represent the projections sampled on the detector with the Logarithmic transformation and exploit subscripts to distinguish the correction projections before and after. For exam-

ple, polychromatic projection Ppolychromatic is simplified into Pp , and monochromatic projection Pmonochromatic is simplified into Pm . It is obvious that the polychromatic projection can be expressed as: Ppolychromatic = L · a (E, L)

(1)

As has been mentioned above, the hardening correction (HC) algorithm in reference [19] gives an exponential model that is used to fit the data between the polychromatic projection and the path length, based on the least-squares method. For example, the path length is the function of the hardening projection (the polychromatic projection on the detector) with the constrained parameter: L = ˛ · Pp · exp(ˇ · (Pp ) )

(2)

Generally, the monochromatic projection shows the law of linear growth. It is the product of path length and attenuation coefficient, where the line traverses the coordinate origin [44], so the linear hardening correction (LHC) employs a corrected line to replace the hardening data due to polychromatic. Based on, the derivative of the function in Eq. (2) at zero is considered to be the slope of the corrected line. That is: K = L (0) = ˛

(3)

Here slope K of the corrected line can be regarded as an equivalent monochromatic projection. Then, the length model for monochromatic based on the polychromatic projection with the same path can be corrected with the Eq. (4) when the hardening curve is obtained. L = K · Pm = L (0) · Pm = ˛ · Pm

(4)

where: • L is the length of the ray crossing object; • a (E, L) is the average absorption coefficient with the length of L; • PP means the polychromatic projection with the Logarithmic transformation; • ˛, ˇ,  is the fitting parameters; • Pm is the monochromatic projection.

4

F. Yang, D. Zhang, H. Zhang et al. / Biomedical Signal Processing and Control 57 (2020) 101823

Fig. 3. The relationship between discrete voxels and rays. (a) is the discrete voxels and rays, (b) shows the matrix coefficient in a 2D image.

Fig. 4. The ray path knowledge. (a) is the preliminary image, (b) shows the binarized slice, (c) shows the length that rays interact with the image.

2.2.2. Weighted-compensation hardening correction (WCHC) algorithm The purpose of the beam hardening correction is to estimate the monochromatic projection which varies with the length of the ray-based on polychromatic projection. To better describe the function of the ray length with projection, an L-based projection model is established. Enlightened by the line integral of path length and attenuation coefficient in Eq. (2), we construct the length model to reflect the relationship based on the polychromatic projection and the attenuation coefficient, such as. L(Pp ) = Pp /a (E, L)

(5)

According to the Taylor formula, the L(Pp ) can be expressed approximately by a derivative near the value Pp0 based on differentiable function, the length model can be rewritten in the following: L(Pp ) = t(Pp0 ) · Pp + t  (Pp0 ) · (Pp − Pp0 ) + o((Pp − Pp0 )n )

(6)

To simplify, the second and subsequent polynomial residual in Eq. (6) are designed as a nonlinear model, approximately. Then a beam hardening curve model can be constructed as: L(Pp ) ≈ c · Pp + a(Pp )b

(7)

Because there is the same path length for the polychromatic and monochromatic, the corrected projection can be addressed from the model that is obtained above: Pm = L(Pp )/L (0) +  = (c · Pp + a(Pp )b )/c + 

Then the final projection (the monochromatic projection) can be calculated as: y = (axb + cx)/c +  ∗ x + 

(9)

At this point, the correction which hauls the polychromatic projection back to monochromatic one has been completed, that is, weighted-compensation hardening correction (WCHC), where, • Pp0 is the mean of polychromatic projection; • t(Pp0 ) refers to the reciprocal of the linear attenuation coefficient, equivalently; • o((Pp − Pp0 )n ) is the polynomial residual; • a, b and c are model coefficients, respectively, and a > 0, b > 0, c > 0; •  is an estimate of noise, which is obtained by Pp − Pp = , and 

Pp is the filtered projection; • x and y state the polychromatic and monochromatic projection, respectively. It should be emphasized that the WCHC algorithm can make full use of the parameters in the original hardened projection model. For example, the LHC algorithm in Eq. (4) uses the parameter A, but the proposed algorithm model in Eq. (9) combines the parameters a, b and c together, and introduces the parameter ␭ as the compensation correction term. Therefore, the WCHC can reflect the projection correction relationship more accurately. 2.3. Performance evaluation

(8)

Moreover, to reduce the influence that the slope of the initial line on the entire correction result, we improve the model and introduce a tradeoff factor  as a compromise to the corrected projection.

To investigate the image characteristics from the corrected algorithm WCHC, we evaluated the image qualities utilizing the root mean square error (RMSE), peak signal to noise ratio (PSNR) [45], and universal quality index (UQI) [46]. RMSE makes an excellent general purpose error metric for numerical predictions. PSNR

F. Yang, D. Zhang, H. Zhang et al. / Biomedical Signal Processing and Control 57 (2020) 101823

describes the difference between the maximum possible signal and the corrupting noise that affect the image. UQI measures both the contrast distortion and the similarity between the recovered and the reference images. They are given by:



RMSE =

 

1 M×N

 PSNR = 20log10

UQI =

(Xij − Xij∗ )2

0≤i≤N 0≤j≤M

 ∗) MAX(X   ∗) RMSE(X, X

 (11)

∗ ¯ × X ¯ 4X X ∗ X

(12)

2 ∗ ¯ + (X ¯ ) ) ( 2 +  2 ∗ )((X) 2

X

(10)

X

where: • Xij and X ∗ denote pixel values of the original image and reconij structed image, respectively; • M and N denote the total number of columns and rows in the image; •   and   ∗ are the variances of Xij and X ∗ ; X X ij •    ∗ is the covariance of Xij and X ∗ ; XX ij

∗ ¯ is the product of the original image and the reconstruction ¯ × X • X image.

For error estimation to the details of the region of interest (ROI) in the selected region, the average gradient (AG) [47], the signalto-noise ratio (SNR) [48], and contrast to noise ratio (CNR) [49] are used. SNR reflects the ratio of signal intensity to noise intensity in ROI. The higher the value, the better the image quality. CNR reflects the contrast of the image in ROI and background. The higher the value, the clearer the image. AG is used to characterize the contour clarity of the image. The larger the AG value is, the more the details are. They are defined as the following: SNR = 10log10



ROI

(14)

S 2 +S 2 ROI BG 2

j

i

j

2 y=1 x=1

M×N

=

(15)

where: • ROI and SROI represent the Mean and Standard Deviation in ROI; • background and SBG represent the Mean and Standard Deviation in the background; • Ix (xi , yj ) = I(xi+1 , yj ) − I(xi , yj ) represents the first derivative of the pixel (xi , yj ) along the x-direction; • Iy (xi , yj ) = I(xi , yj+1 ) − I(xi , yj ) represents the first derivative of the pixel (xi , yj ) along the y-direction. 3. Results and discussion 3.1. Experimental results To avoid the issue that the ray cannot penetrate the objects especially in industry application (It has been a hot topic in CT imaging), we design a #brushing in a diameter of 40 with the material of iron to explain the feasibility and correctness of the proposed modeling based on the cylindrical structure. Furthermore, we design

M 

(Pm )i

i=1 M  i=1

N M   Ix (x ,y )2 +Iy (x ,y )2

AG =

M 

(13)

  ROI − background  CNR =

i

the #wheel in a diameter of 60 with the material of titanium and adding some auxiliary rib structure to change the path of the ray’s attenuation. As has been shown in Fig. 5, the categories of projections (the material of iron named #Bushing, and the material of titanium named #Wheel) are utilized to validate the performance of the proposed method. The distance from the source to the center of rotation is 787 mm and the distance from the source to the detector is 1225 mm. The largest diameter of the #Bushing and #Wheel are 60 mm and 90 mm, respectively. Projections are obtained form 0◦ ∼360◦ scanning using the FPD as XRD 1621 AN15 ES from PerkinElmer. The x-ray tube of the CBCT system is MXR-451HP/11 from Comet. We collect 360 projections in a single scanning with resolution of 512 × 512 and pixel size of 0.2 mm × 0.2 mm. Table 1 shows the scanning parameters and the related value based on the experience of penetration performance, respectively. According to the knowledge and the initial position of the objects, the sampling angle is selected within the range of 70◦ ∼120◦ at the step of 10◦ . Fig. 6 shows the process of curve fitting with the obtained discrete data, where Fig. 6(a) denotes the X-ray sequence which interacts with the binary image of the object at different angles, Fig. 6(b) gives the fitting curve of the sampled discrete data by using the least-squares method. Based on the knowledge of the curve fitting, the model in Eq. (7) is analyzed, where the first and second derivative of the function is taken, and the derivative L (Pp ) > 0, L (Pp ) > 0 is established. It can be concluded that the model L(Pp ) increases monotonically and is a convex where there is no stagnation point and inflection point within the range [0, +∞). Obviously, this function satisfies the relevant requirements for the beam hardening characteristic and the fitting curve, and it is simple and convenient for computing. The model in Eq. (9) has been used to acquire the monochromatic projection, where the tradeoff factor  as a compromise to the model is defined in the following. Table 2 shows the  results and sampling parameter of the curve fitting.



SROI

5

= (Pp )i

(axib + cxi )/c

i=1 M 

(16) xi

i=1

Figs. 7 and 8 show the reconstruction results in different sequences (such as the 140th and the 260th slice) from the original projection (the first column), LHC projection (the second column) and WCHC projection (the third column), respectively. It is found that the region nearby the edge in ROI 1 and ROI 2 are difficult to distinguish from the profile of the uncorrected projection. After the method of LHC, the sharpness of the projection is improved, but the outer boundaries are weakened. In contrast, the WCHC method reduces the deficiencies mentioned above and improves the result greatly. There is no doubt that the cupping artifacts are serious in the slice reconstructed from the original projection, which has made the edge bright and inner structure dark. Although the LHC algorithm has cut down the artifacts, the method makes the grayscale inconsistency, which has resulted in the edge dark but the inner bright. Comparatively, the WCHC algorithm shows a better result, the grayscale is relatively consistent, and the hardening artifacts are suppressed at the same time, as is showed in Figs. 7(d) ∼ (f) and 8 (d) ∼ (f). To quantitatively demonstrate the benefit of the proposed beam hardening correction method, the grayscale with different methods on the reconstruction of ROIs is compared, which has been labeled with l1 , l2 line above. As shown in Fig. 9, the proposed method

6

F. Yang, D. Zhang, H. Zhang et al. / Biomedical Signal Processing and Control 57 (2020) 101823

Fig. 5. The objects and the projections of #Bushing and #Wheel. (a) and (c) are the pictures of the objects; (b) and (d) show the projection, respectively.

Table 1 The sampling parameters for different objects in experimental. Experimental

Material/density(g/cm3 )

Voltage (kVp)

Current (mA)

Projection Number

#Bushing #Wheel

iron /7.86 titanium /4.5

440 420

0.22 0.18

360 360

Fig. 6. Curve fitting with discrete data. (a) shows the interaction of the ray with the binary image, and (b) shows the curve fitting the least-squares method.

Table 2 The related value of the curve fitting. Experimental

#Bushing #Wheel

Angle range/ step

Ray sequence/ step

Modeling parameter a

b

c

70◦ ∼120◦ /10 70◦ ∼120◦ /10

130∼180/10 80∼120/10

0.00348 0.00041

7.6441 7.2152

0.8251 2.9202

effectively reduces the low contrast and improves the distribution of grayscale. Obviously, the edge profile of the section using the WCHC algorithm is clear. And the grayscale of the whole image has a uniform distribution compared to the original and LHC algorithm. That is, the results with the weighted model are better than that of the linear model, where the cupping artifacts are well suppressed. Moreover, Table 3 shows the numerical comparison of the results (the quantitative results are conducted on the whole image). It can be seen from the table that the RMSE of the proposed method is dramatically reduced by 31.6 % and 30.6 % of that without beam hardening correction, respectively. PSNR and UQI also indicate better uniformity of the results with the WCHC method, where the PSNR is increased by 16.8 % and 19.4 %, UQI is increased by 12.8 % and 15.9 % for experimental #Bushing and #Wheel, respectively. To further illustrate the effectiveness of the algorithm, Fig. 10 shows the details of the ROI1 and ROI2 of different algorithms. It is clear that the cupping artifacts of the selected region have

(curve fitting)

Tradeoff factor  4.67519 2.12011

been reduced both in the LHC algorithm and the WCHC algorithm. Although there are still tiny artifacts by using the WCHC algorithm, the grayscale distribution in the edge area has been improved, compared to the uncorrected and LHC algorithm. Table 4 shows the evaluation of local details for experimental #Bushing in ROI1 and ROI2 marked in Fig. 10. Here, the SNR in ROI A is computed, CNR and AG in ROI B are performed. From the numerical comparison of Table 4, the CNR and AG by using LHC and WCHC algorithms are both promoted, but they present a higher numerical value in the CNR and AG indicators within a small window, because the LHC algorithm possesses the characteristic that the outer edge is dark and the inner boundary is bright for the corrected image, which leads to the contrast and the contour clarity of the image in ROI to be slightly higher than the WCHC algorithm. While the SNR by using the WCHC algorithm is going to be higher than the LHC, which reflects the ratio of signal intensity. Although there are still tiny artifacts by using the WCHC algorithm, the grayscale profile which is plotted in Fig. 9 has been greatly improved. Moreover,

F. Yang, D. Zhang, H. Zhang et al. / Biomedical Signal Processing and Control 57 (2020) 101823

7

Fig. 7. The results of #Bushing with different processes and corrections. (a), (b) and (c) are the 140th slice from the original projection, LHC and WCHC, respectively. display window is: [−0.02, 1.10], (d), (e) and (f) are the 260th slice, correspondingly, display window is: [−0.005, 0.86].

Fig. 8. The results of #Wheel with different processes and corrections. (a), (b) and (c) are the cross-view from the original projection, LHC and WCHC, respectively. (d), (e) and (f) are the sagittal-view, correspondingly, the display window is: [−0.005, 0.86].

Fig. 9. The grayscale of the image with a different projection. (a) plots the grayscale along the line l1 in Fig. 7, (b) plots the grayscale along the line l2 in Fig. 8.

8

F. Yang, D. Zhang, H. Zhang et al. / Biomedical Signal Processing and Control 57 (2020) 101823

Table 3 Numerical comparison of image quality. Experimental

#Bushing

#Wheel

Method Original LHC WCHC Original LHC WCHC

PSNR

RMSE

UQI

value

improved/%

value

improved/%

value

improved/%

0.1051 0.0958 0.0719 0.1523 0.1260 0.1056

/ 8.85 31.59 / 17.27 30.66

19.57 20.37 22.86 16.34 17.99 19.52

/ 4.10 16.81 / 10.06 19.44

0.4827 0.5046 0.5448 0.2422 0.2691 0.2809

/ 4.53 12.86 / 11.10 15.97

Fig. 10. The details of the selected region of experimental #Bushing and #Wheel. The first column shows the region of slices, and the columns from the second to the fourth are the details of the region within a different algorithm. The display window is: [−0.005, 0.86].

Table 4 Numerical comparison of image quality.

Experimental

Region 1 Region 2

Reconstructed projection

Original LHC WCHC Original LHC WCHC

Evaluation SNR

CNR

AG

value

improved/%

value

improved/%

value

improved/%

13.7315 15.8199 16.5473 20.5953 22.9782 23.2581

/ 15.21 20.50 / 11.57 12.93

5.1802 6.6017 6.3572 4.8341 5.9673 5.7456

/ 27.44 22.72 / 23.44 19.72

0.0227 0.0509 0.0480 0.0266 0.0459 0.0417

/ 124.49 111.56 / 72.58 56.62

the RMSE, PSNR and UQI indicators of the whole image in Table 3 denotes that the WCHC algorithm shows better characteristics than the other two correction results. Therefore, the WCHC algorithm works better than the LHC algorithm.

3.2. Tradeoff study As has been mention above, the tradeoff factor  is a compromise to the model in Eq. (9). The tradeoff curve is an effective tool to evaluate the presented implementation. In this study, we employ the experimental #Bushing and #Wheel with different materials, including aluminum, titanium, and iron to analysis the performance of the tradeoff factor.

Table 5 gives the tradeoff  which varies with the density of different experimental when using the model in Eq. (16). Obviously, the higher the density, the higher the tradeoff weight. Fig. 11 shows the parameter behavior of  which varies with density. Fig. 11(a) plots the factor  in the model in Eq. (9) by using the proposed method in Eq. (16) on the phantom with different materials both in experimental # Bushing (#A) and experimental # Wheel (#B). Where the density of the material is 2.78 g/cm3 of aluminum, 4.51 g/cm3 of titanium and 7.87 g/cm3 of iron. According to the tendency of the curve in Fig. 11(a), a further study of  is predicted and conducted. Based on the mean value of the experimental data in Table 5, the least-squares fitting is performed by using the power function to obtain the parameter behavior of . Here, a simple power function is selected, such as: () = m

(17)

F. Yang, D. Zhang, H. Zhang et al. / Biomedical Signal Processing and Control 57 (2020) 101823

9

Table 5 The tradeoff factor of different experiments with different materials.

Fig. 11. Parameter behavior of  which varies with density. (a) shows the tendency of  with different materials, (b) is the fitting curve of the .

where the  is the density of the material, m is the model knowledge as the tradeoff estimation. As has been shown in Fig. 11(b), the power function in Eq. (17) represents the behavior of the tradeoff weight, appropriately, so the value is used as the knowledge information to estimate the tradeoff . In order to illustrate the validity of the weight parameter estimation in Eq. (17), the experimental for the Cube named #Slot with different densities are selected for verification. We stack objects of different materials together and obtain the projection of three objects in one scanning. The distance from the source to the center of rotation is 811 mm and the distance from the source to the detector is 1235 mm. The sides of a Cube are 30 mm. Table 6 shows the parameters of the beam hardening (cupping artifacts) correction model in Eq. (9). And the tradeoff  with different densities are computed in Eq. (17). The projection of #Slot is corrected according to the parameters in Table 6, and the results of the CT measurements to investigate the cupping effect as a measure of the amount of beam hardening are shown in Fig. 12. The columns show the entity, original projection and corrected projection by using the WCHC algorithm, respectively. The rows present the corrected results with the tradeoff estimation of  of different densities. Compared with the raw projection, the contrast of the section reconstructed by the WCHC algorithm has greatly improved. Obviously, the gray distribution of the image is more uniform, where the brighter edges and darker interior (that is cupping artifacts) are well suppressed. The results of the chart above show that the tendency of tradeoff estimation by ␭() with different densities is effective and feasible. It demonstrates that the proposed algorithm indeed reduces the cupping artifacts, and the performance of the WCHC model is better than that of the LHC algorithm. Hence, it should be certain that the heuristic assumption made in the derivation is reasonable, and it is induced almost no extra inaccuracy in image reconstruction. There is much room for improvement of the model, and further research is necessary to deal with various metal and bone-related artifacts. The proposed algorithm does not consider metal artifacts that arise from scattering, nonlinear partial volume effects or electric noise. The performance of the proposed sinogram correction method could be improved by enhancing the forward

model, which accurately represents various realistic artifacts. The following research will include further investigation into learning approaches for sinogram correction and clinical research with patients. 4. Conclusion In this paper, an effective cupping artifacts correction method has been described. A method of obtaining the traversing lengths between the ray and the binary image is proposed to obtain discrete beam hardening data. Based on the characteristics of hardening behavior and knowledge, a new weighted and compensated model is constructed. The performance of the proposed method has been verified on the CBCT system. The preliminary data obtained in the experimental investigation shows that the proposed algorithm can at least perform better than the linear corrected algorithm for polychromatic projection. Hence, the proposed algorithm is anticipated to find its utility in nondestructive testing applications wherein the image reconstruction for the complex geometry and high energy X-ray imaging. The key issues discussed in this manuscript concern more about industrial nondestructive imaging. And the phantoms used in the experimental studies are mono-materials. The proposed model is also effective for complicated objects for the adaptive estimation of the tradeoff factor ␭. While there are still deficiencies for estimating the tradeoff factor when rays transmitted through multi-materials. Therefore, the beam hardening correction method for CT based on multiple materials is the key methodology of our next research. Acknowledgments This work was supported in part by the National Natural Science Foundation of China (51675437 and 51605389); Fund of Ministry of Industry and Information Technology of China (MJ2017-F-05); The Fundamental Research Funds for the Central Universities (31020190504006); China Postdoctoral Science Foundation (2019M653749). The authors are grateful to the anonymous reviewers for their valuable comments.

10

F. Yang, D. Zhang, H. Zhang et al. / Biomedical Signal Processing and Control 57 (2020) 101823

Table 6 The related value and tradeoff estimation of test pieces #Slot.

Fig. 12. Verification and comparison of artifacts corrected algorithm. (a), (b) and (c) are the uncorrected slices, (d), (e) and (f) are the corrected results using the WCHC algorithm. The display window is: [−0.05, 1.08].

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References [1] F.Q. Yang, D.H. Zhang, K.D. Huang, Z.Z. Gao, Y.F. Yang, Incomplete projection reconstruction of computed tomography based on modified discrete algebraic reconstruction technique, Meas. Sci. Technol. 29 (2) (2018) 025405. [2] P. Jin, C.A. Bouman, K.D. Sauer, A model-based image reconstruction algorithm with simultaneous beam hardening correction for X-Ray CT, IEEE Trans. Comput. Imaging 1 (3) (2015) 200–216.

[3] O.M.H. Ahmed, Y.S. Song, A review of common beam hardening correction methods for industrial X-ray computed tomography, Sains Malays. 47 (8) (2018) 1883–1890. [4] A. Bevilacqua, D. Barone, S. Malavasi, G. Gavelli, Automatic detection of misleading blood flow values in CT perfusion studies of lung cancer, Biomed. Signal Process. Control 26 (2016) 109–116. [5] S.J. Tang, K.D. Huang, Y.Y. Cheng, X.Q. Mou, X.Y. Tang, Optimization based beam-hardening correction in CT under data integral invariant constraint, Phys. Med. Biol. 63 (13) (2018) 135015. [6] W.C. Cao, T. Sun, G. Kerckhofs, G. Fardell, B. Price, W. Dewulf, A simulation-based study on the influence of the x-ray spectrum on the performance of multi-material beam hardening correction algorithms, Meas. Sci. Technol. 29 (9) (2018) 095002. [7] R.A. Ketcham, R.D. Hanna, Beam hardening correction for X-ray computed tomography of heterogeneous natural materials, Comput. Geosci. 67 (4) (2014) 49–61. [8] K.D. Huang, D.H. Zhang, M.J. Li, Noise suppression methods in beam hardening correction for X-ray computed tomography, in: International Congress on Image and Signal Processing, IEEE, 2009, pp. 1–5.

F. Yang, D. Zhang, H. Zhang et al. / Biomedical Signal Processing and Control 57 (2020) 101823 [9] L.Y. Meng, X.P. Guo, H.G. Li, MRI/CT fusion based on latent low rank representation and gradient transfer, Biomed. Signal Process. Control 53 (2019) 101536. [10] R.B. Ger, D.F. Craft, D.S. Mackin, et al., Practical guidelines for handling head and neck computed tomography artifacts for quantitative image analysis, Comput. Med. Imaging Graph. 69 (2018) 134–139. [11] Y.I. Nesterets, T.E. Gureyev, Noise propagation in x-ray phase-contrast imaging and computed tomography, J. Phys. D Appl. Phys. 47 (10) (2014) 289–290. [12] H.S. Park, D. Hwang, J.K. Seo, Metal artifact reduction for polychromatic X-ray CT based on a beam-hardening corrector, IEEE Trans. Med. Imaging 35 (2) (2016) 480–487. [13] A. Shiras, F. Robert, B. Richard, M. Steffen, B. Oliver, R. Georg, Beam hardening correction using cone beam consistency conditions, IEEE Trans. Med. Imaging 37 (10) (2018) 2266–2277. [14] Z. Christoph, R. Simon, K. Christopher, V. Gloria, K. Florian, D. George, B. Claus, P. Katia, L. Guillaume, Decomposing a prior-CT-based cone-beam CT projection correction algorithm into scatter and beam hardening components, Phys. Imaging Radiat. Oncol. 3 (2017) 49–52. [15] L. Yu, S. Leng, C.H. Mccollough, Dual-energy CT-based monochromatic imaging, AJR Am. J. Roentgenol. 199 (5) (2012) S9. [16] Y.B. Chang, D. Xu, A.A. Zamyatin, Metal artifact reduction algorithm for single energy and dual energy CT scans, in: Nuclear Science Symposium and Medical Imaging Conference, IEEE, 2012, pp. 3426–3429. [17] H. Shi, Z. Yang, S. Luo, Reduce beam hardening artifacts of polychromatic X-ray computed tomography by an iterative approximation approach, J. Xray Sci. Technol. 25 (3) (2017) 417. [18] G.V. Gompel, K.V. Slambrouck, M. Defrise, K.J. Batenburg, J.D. Mey, J. Sijbers, Iterative correction of beam hardening artifacts in CT, Med. Phys. 38 (S1) (2011). [19] S.J. Tang, K.D. Huang, Y.Y. Cheng, X.Q. Mou, X.Y. Tang, Optimization based beam-hardening correction in CT under data integral invariant constraint, Phys. Med. Biol. 63 (13) (2018). [20] H. Wang, Y. Xu, H. Shi, A new approach for reducing beam hardening artifacts in polychromatic X-ray computed tomography using more accurate prior image, J. Xray Sci. Technol. 26 (4) (2018) 593–602. [21] R. Tovey, M. Benning, C. Brune, M.J. Lagerwerf, S.M. Collins, R.K. Leary, P.A. Midgley, C.B. Schonlieb, Directional sinogram inpainting for limited angle tomography, Inverse Probl. 35 (2) (2018) 024004. [22] T. Kano, M. Koseki, A new metal artifact reduction algorithm based on a deteriorated CT image, J. Xray Sci. Technol. 24 (6) (2016) 901–912. [23] H. Zhang, L. Li, L. Wang, Y. Sun, B. Yan, A. Cai, G. Hu, Computed tomography sinogram inpainting with compound prior modelling both sinogram and image sparsity, IEEE Trans. Nucl. Sci. 63 (5) (2016) 2567–2576. [24] X. Zhang, J.W.L. Wan, Image restoration of medical images with streaking artifacts by Euler’s elastica inpainting, in: IEEE International Symposium on Biomedical Imaging, IEEE, Melbourne, VIC, Australia, 2017. [25] C. Yang, P. Wu, S. Gong, S. Gong, J. Wang, Q. Lyu, X. Tang, T. Niu, Shading correction assisted iterative cone-beam CT reconstruction, Phys. Med. Biol. 62 (22) (2017) 8495–8520. [26] P. Wu, X. Sun, H. Hu, T. Mao, W. Zhao, K. Sheng, A. Cheung, T. Niu, Iterative CT shading correction with no prior information, Phys. Med. Biol. 60 (21) (2015) 8437–8455. [27] L. Brabant, E. Pauwels, M. Dierick, D. Van Loo, L. Van Hoorebeke, A novel beam hardening correction method requiring no prior knowledge, incorporated in an iterative reconstruction algorithm, NDT E Int. 51 (2012) 68–73. [28] W. Cao, S. Hawker, G. Fardell, B. Price, W. Dewulf, An improved segmentation method for multi-material beam hardening correction in industrial x-ray computed tomography, Meas. Sci. Technol. 30 (12) (2019) 125403. [29] C. Romano, J.M. Minto, Z.K. Shipton, R.J. Lunn, Automated high accuracy, rapid beam hardening correction in X-Ray Computed Tomography of multi-mineral, heterogeneous core samples, Comput. Geosci. 131 (2019) 144–157. [30] W. Zhao, D.W. Li, K. Niu, W.J. Qin, H. Peng, T.Y. Niu, Robust beam hardening artifacts reduction for computed tomography using spectrum modeling, IEEE Trans. Comput. Imaging 5 (2) (2018) 333–342. [31] S. Sarkar, P. Wahi, P. Munshi, An empirical correction method for beam-hardening artifact in Computerized Tomography (CT) images, NDT E Int. 102 (2019) 104–113. [32] H.S. Park, S.M. Lee, H.P. Kim, J.K. Seo, Y.E. Chung, CT sinogram-consistency learning for metal-induced beam hardening correction, Med. Phys. 45 (12) (2018) 5376–5384.

11

[33] H.M. Zhang, L.Y. Wang, L. Li, A.L. Cai, G.E. Hu, B. Yan, Iterative metal artifact reduction for x-ray computed tomography using unmatched projector/back projector pairs, Med. Phys. 43 (6) (2016) 3019–3033. [34] M. Bazalova, L. Beaulieu, S. Palefsky, F. Verhaegena, Correction of CT artifacts and its influence on Monte Carlo dose calculations, Med. Phys. 34 (6) (2007) 2119–2132. [35] H.S. Park, J.K. Choi, K.R. Park, K.S. Kim, S.H. Lee, J.C. Ye, J.K. Seo, Metal Artifact Reduction in CT by identifying missing data hidden in metals, J. Xray Sci. Technol. 21 (3) (2013) 357–372. [36] A. Mehranian, M.R. Ay, A. Rahmim, H. Zaidi, X-ray CT metal artifact reduction using wavelet domain sparse regularization, IEEE Trans. Med. Imaging 32 (9) (2013) 1707–1722. [37] X. Duan, L. Zhang, Y. Xiao, J. Cheng, Z. Chen, Y. Xing, Metal artifact reduction in CT images by sinogram TV inpainting, in: Nuclear Science Symposium Conference Record, Dresden, Germany, Oct 19, 2008, pp. 4175–4177. [38] Y. Kyriakou, E. Meyer, D. Prell, M. Kachelriess, Empirical beam hardening correction (EBHC) for CT, Med. Phys. 37 (10) (2010) 5179–5187. [39] X.B. Zou, Z.J. Li, TV-based correction for beam hardening in computed tomography, J. Med. Imaging Health Inform. 6 (7) (2016), 1707–1707. [40] A.C. Kak, M. Slaney, G. Wang, Principles of computerized tomographic imaging, Med. Phys. 29 (1) (2002) 107. [41] F.Q. Yang, D.H. Zhang, K.D. Huang, W.L. Shi, C.X. Zhang, Z.Z. Gao, Projection matrix acquisition for cone-beam computed tomography iterative reconstruction, in: Second International Conference on Photonics and Optical Engineering, Xi’an, China, 2017, p. 102560J. [42] W.V. Aarle, K.J. Batenburg, J. Sijbers, Optimal threshold selection for segmentation of dense homogeneous objects in tomographic reconstructions, IEEE Trans. Med. Imaging 30 (4) (2011) 980. [43] K.D. Huang, D.H. Zhang, K. Wang, Beam hardening correction method for cone-beam computed tomography based on registered model simulation, J. Syst. Simul. 21 (4) (2009) 1164–1167. [44] K.D. Huang, D.H. Zhang, Y.M. Kong, K. Wang, Beam hardening correction method for Cone-Beam computed tomography based on reprojection of slice contours, Chin. J. Sci. Instr. 29 (9) (2008) 1873–1877. [45] Q. Huynh-Thu, M. Ghanbari, Scope of validity of PSNR in image/video quality assessment, Electron. Lett 44 (13) (2008) 800. [46] C.J. Hsieh, T.K. Huang, T.H. Hsieh, G.H. Chen, K.L. Shih, Z.Y. Chen, J.C. Chen, W.C. Chu, Compressed sensing based CT reconstruction algorithm combined with modified Canny edge detection, Phys. Med. Biol. 63 (2018) 155011. [47] K.D. Huang, D.H. Zhang, M.J. Li, H. Zhang, Image lag modeling and correction method for flat panel detector in cone-beam CT, Acta Phys. Sin. 62 (21) (2013) 210702. [48] F. Yang, D. Zhang, K. Huang, W. Shi, X. Wang, Scattering estimation for cone-beam CT using local measurement based on compressed sensing, IEEE Trans. Nucl. Sci. 65 (3) (2018) 941–949. [49] H. Zhang, K.D. Huang, Y.K. Shi, M.J. Li, Ring artifact correction method based on air scan for cone-beam CT, Comput. Tomogr. Theory Appl. 21 (2) (2012) 247–254. Fuqiang Yang is with the Key Laboratory of High Performance Manufacturing for Aero Engine (Northwestern Polytechnical University), Ministry of Industry and Information Technology; Engineering Research Center of Advanced Manufacturing Technology for Aero engine (Northwestern Polytechnical University), Ministry of Education, Xi’an 710072, China. Dinghua Zhang is with the Key Laboratory of High Performance Manufacturing for Aero Engine (Northwestern Polytechnical University), Ministry of Industry and Information Technology; Engineering Research Center of Advanced Manufacturing Technology for Aero engine (Northwestern Polytechnical University), Ministry of Education, Xi’an 710072, China. Hua Zhang is with the School of Mechanical Engineering, Northwestern Polytechnical University, Xi’an 710072, China. Kuidong Huang is with the Key Laboratory of High Performance Manufacturing for Aero Engine (Northwestern Polytechnical University), Ministry of Industry and Information Technology; Engineering Research Center of Advanced Manufacturing Technology for Aero engine (Northwestern Polytechnical University), Ministry of Education, Xi’an 710072, China.