Scattering measurement and estimation in angular sequence for cone-beam CT based on projection structural tensor and modeling

Abstract

Based on the structural tensor of projection, this study aims to address and test a new improved algorithm applying to the distort projection data to generate a high qualified image by reducing the artifacts and noise from scattering in the cone-beam computed tomography (CBCT). Since the scattering information has a large relationship with the structure of the object, which is reflected by the projection, regional model knowledge for scattering is accomplished by finding the relationship between projection and scattering. As the tensor, the gradient of projection is first calculated in the process for estimating the direction and structural edge of the object. Then, the Determinant and Traces of the tensor map with different characteristics are computed to determine the different regions. By modeling and fitting the regions of scattering distribution, the knowledge of scattering parameters corresponding to a different region is obtained. Based on the similarity of scattering distribution in adjacent angles, the scatterings with angle sequence are completed by interpolating the prior knowledge obtained through the sparse sampling. By performing the studies on polychromatic X-ray to test the performance of the scattering estimation algorithm, the results show a significant improvement in the images that are reconstructed from the corrected projection. The root mean square error (RMSE) of the proposed method is reduced by 21.8% and 39.8%, respectively. Peak signal to noise ratio (PSNR), and universal quality index (UQI) also indicate better uniformity, where the PSNR is increased by 7.4% and 56.7%, UQI is increased by 70.8% and 262.3% for experimental #Wheel and #Cylinder, respectively.

Keywords

Scattering estimation angular sequence structural tensor cone-beam computed tomography corrected projection

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]. It has been a powerful detection technique in applied physics and materials science due to its unique ability to characterize the properties of an object at a high spatial resolution. While the large radiation of CBCT and the signal crosstalk among the detectors, which leads to the detector receiving more scattered photons, resulting in contrast reduction, scattering artifacts and poor CT quantitative accuracy. It has greatly reduced the image quality and limits its further development [3, 4].

As a focused problem to be solved, some strategies for scattering suppression in CBCT has been operated, including knowledgeable selection of projection [3], air gap [5], use of an anti-scattering grid [6, 7], analytical modeling [8], modulation methods [9, 10], Monte Carlo (MC) simulation [11, 12] and blocker-based techniques [13, 14]. The analytical modeling is useful but it fails to heterogeneous and complex objects, Monte Carlo could be accurate but computationally intense, etc. For the projection is the sum of the transmitted rays and the scattered rays, the accurate estimation of the scattering intensity is crucial to the projection correction. Ruhrnschopf et al. [15, 16] divided the scattering estimation into a hybrid method based on measurement, mathematical-physical modeling, and multiple estimations. As involving statistics, this method can simulate complex objects but it takes a long time for computing. Ning et al. [14, 17] showed that the scattering could be suppressed through adding a scattering grating, and the inhomogeneity and linearity of CT values can be greatly improved. Sisniega et al. [18] studied the inhibition effect of geometry and scattering grating to the object, including small animals, hands and feet, head, breast, etc. on CBCT system by using Monte Carlo simulation. It has indicated that the scattering grating could reduce the shaped artifacts, but the correction effect on scattering caused by strip artifacts was limited. Though, increasing the distance between the scanned object and the detector is the best method of scattering suppression [19], the contrast noise ratio (CNR) will lose without changing the intensity of the ray. Ouyang et al. [20] proposed a physical attenuator (i.e., “blocker”) consisting of equal spaced lead strips to simultaneously estimate scatter and reconstruct the complete volume within the field of view (FOV).

To obtain scattering more accurately, the researchers designed the scanning mode of moving grid bars [21, 22] and improved the corresponding reconstruction algorithm, but this process increases the complexity of CT scanning by moving two blockers alternatively. Briggs. et al. [23] proposed the method of beam attenuation grid and obtained the scattering image by formula derivation. This method needs to scan the bright field when there is/no attenuation grid. Shi W. [24] proposed a new scattering correction method based on curve fitting. The method uses the single scattering measurement by the grating to fit the scattering with the projection to estimate all scattering distributions. Although all the scattering distribution is obtained, the estimation error is serious, and the structural impact on scattering is not well interpreted [22, 25]. A vast amount of effort was dedicated in the past few years to improve the scattering estimation of the object; this progress has made it possible to measure the scattering in an irregular object with high density. Thus, the next frontier of scattering correction progress is to improve the signal to noise ratio (SNR) of structural scattering measurements by knowledge sense. In this paper, the regional model knowledge as well as scattering estimation with angle sequence, are explored and accomplished.

To obtain the scattering distribution of different structural, a new regional scattering model based on the projection structural tensor has developed with the guidance of knowledge obtained by sparse sampling, which is interpolated to complete the angular sequence scattering estimation. Finally, the original projection correction is completed by using the estimated scattering distribution.

2 Material and methods

2.1 Method overview

The general scattering estimation in angular sequence for CBCT is shown in Fig. 1. We first estimate few scattering using a partially blocker in the scanning at sparse angle. These scattering are stored as the prior for scattering estimation in angular sequence. To better characterize the effect of structural information on scattering, we use structural tensor (ST) to obtain the region. The projection gradient was first calculated in the process for getting the direction and structural edge of the object. The Determinant (D) and Traces (T) of the tensor map were computed to determine the different region of the projection. Moreover, regional model for scattering is accomplished by finding the relationship between projection and scattering. By fitting different region of prior scattering distribution, the knowledge of scattering parameters corresponding to a different region is obtained. Based on the similarity of scattering distribution in adjacent angles, the scattering estimation with angle sequence is completed. To evaluate the scattering estimation from the premeasured data, we reconstruct the images using the scattering map form different algorithms. These estimated scattering distributions are finally subtracted from the projection data as the correction. Details of these procedures are discussed below.

Above on, the preliminary of the proposed scattering estimation approach contains the following steps:

Fig.1

Work-flow of the scattering estimation using projection structural tensor and modeling from partially prior scattering measurement.

(1)

The scattering measurement in sparse-angle based on the grating/blocker measurement.

(2)

The projection decomposition for the regional division. The trace (T) combined with the determinant (D) of structure tensor is employed.

(3)

A simple and nonlinear projection-based model for scattering estimation is utilized to estimate the scattering in a different region.

2.1.1 The preliminary of scattering measurement

The scattering based on the grating measurement actually is a method that obtains the map through the relationship between the transmission and the projection. As far as we know, the notion of using the aid of projection to reflect the structural information is suitable. Since the scattering information has a large relationship with the structure of the object which has been reflected by the projection. And the validity and correctness of this hypothesis have been confirmed in references mentioned in [26, 27].

The geometry of the blocked CBCT system is shown in Fig. 2, the scattering is obtained by collecting two projection information:

Fig.2

Schematic diagram of scattering from a grating. (a) shows the sampling without blocker;(b) shows the sampling with a blocker.

Obtaining the projection of the object without grating. As shown in Fig. 2(a), the signal that reached to the detector is the transmission information and the object scattering, which is denoted as P + S _object.

Obtaining the projection of the object with the grating in front of the detector. As shown in Fig. 2(b), the signal that reached to the detector is the transmission information and grating additional scattering, which is denoted as P + S _blocker.

Considering the characteristics that low frequency and smoothness of the scattering, the scattering acquisition is based on the following assumptions: (1) The scattering of object blocked by the lead strip is all absorbed; (2) The signal reached to the detector is corresponding to the blocker, where the blocked part only contains the scattering itself. For the transmitted photons can pass through the slit portion, the photons received by the detector is the transmission plus blocker scattering.

The flowchart of the scattering obtained by the blocker is shown in Fig. 3, and the specific operation and post-processing are carried out according to the steps in Ref. [24], which has been published as academic achievement. For the scattering information obtained by grating shows a strip distribution, the scattering map is filtered after interpolating. The scattering obtained by beam blocker can be used as a priori and reference for scattering estimation below.

Fig.3

Flowchart of the scattering obtains by the blocker.

2.1.2 Structural tensor description

Structure Tensor (ST) is an extraction algorithm which aims to estimate the direction and local structure of the image edge, and has been widely used in image processing and computer vision. Its cover corner detection [28], texture analysis [29] and so on.

For a two-dimensional (2D) image I (x, y), the gradient vector of the image is defined as: $\nabla I = {(I_{x}, I_{y})}^{Tr}$ (1) where I _x and I _y indicate the partial derivative of the image along the x and y directions. Tr indicates the transposition operation. Then the ST of the image can be expressed as: $E_{0} = \nabla I \cdot \nabla I^{Tr} = (\begin{matrix} I_{x}^{2} & I_{x} I_{y} \\ I_{y} I_{x} & I_{y}^{2} \end{matrix}) = (\begin{matrix} u_{11} & u_{12} \\ u_{21} & u_{22} \end{matrix})$ (2)

According to the definition of the gradient vector in Equation (1), the vertical gradient vector is defined as: $\nabla I_{⊥} = {(- I_{y}, I_{x})}^{Tr}$ (3)

From which the gradient of ST can be calculated: $E_{0} \cdot \nabla I = (\begin{matrix} I_{x} \\ I_{y} \end{matrix}) \cdot (I_{x}, I_{y}) \cdot (\begin{matrix} I_{x} \\ I_{y} \end{matrix}) = {∥ \nabla I ∥}^{2} \nabla I$ (4) $E_{0} \cdot \nabla I_{⊥} = (\begin{matrix} I_{x}^{2} & I_{x} I_{y} \\ I_{y} I_{x} & I_{y}^{2} \end{matrix}) \cdot (\begin{matrix} - I_{y} \\ I_{x} \end{matrix}) = (\begin{matrix} I_{x} \\ I_{y} \end{matrix}) \cdot (I_{x}, I_{y}) \cdot (\begin{matrix} - I_{y} \\ I_{x} \end{matrix}) = 0$ (5)

Therefore, the tensor information represents the strength of the image grayscale change. The larger the eigenvalue is, the more obvious the gradation change is in the direction of the normal vector; the smaller the eigenvalue is the less change in vector of the edge. Two eigenvalues λ ₁, λ ₂ and eigenvectors v ₁, v ₂ (where, λ ₁ ⩾ λ ₂, v ₁ ⊥ v ₂) are obtained by decomposition the ST matrix, and the corresponding computing are in Equations (6 and 7). $λ_{1, 2} = \frac{1}{2} (u_{11} + u_{22} \pm \sqrt{{(u_{11} - u_{22})}^{2} + 4 u_{12}^{2}})$ (6) $v_{1} = {(2 u_{12}, \sqrt{{(u_{11} - u_{22})}^{2} + 4 u_{12}^{2}} - u_{11} + u_{22})}^{Tr}$ (7)

Although the process of X-ray image generation is usually ill-posed, for example, the projection is affected by noises somewhere, the values on the detector are much larger than zero, so the determinant is finite. Based on a 2D projection image, the trace (T) combined with the determinant (D) of the ST can reflect different regions. Since the eigenvalues of the ST (such as trace T) describe the contrast change along the feature direction, the determinant (D) of the ST reflects the coherence of the surroundings, so, the flat region, edge region, and corner region of the image can be distinguished. The local features such as edges and corner points can be better identified.

To acquire the local structure information quickly through the structure tensor, the D and T are obtained from the obtained ST matrix, and the regional pattern of the image can be distinguished based on the relationship between D and T. Assuming a given threshold τ, the pattern is divided into the following categories: ①The flat region: T < = τ & |D| ≈ 0; ② The edge region: T > τ & |D| ≈ 0; ③The corner region: T > τ & |D| > 0. Figure 4(a) shows the relationship between the D and T of the ST. The area Ω ₁ indicates a flat area, the area Ω ₂ indicates an edge area, and the area Ω ₃ indicates a corner area. Of course, there are no strict distinctions between the regional boundaries, because there is a certain overlap in the connected boundaries.

Fig.4

Decomposition for the projection, (a) shows the relationship between D and T of structure tensor, (b) is the background; (c) ∼ (e) are the flat region, edge region and corner region of projection.

Since there exists a great gap of the grayscale value between the object and the background on the detector, the flat region in the background region can be obtained based on the law T < = τ & D ≈ 0, approximately. In this study, the threshold is selected as 12000, such as Fig. 4(b) and 4(c) show the background and flat region in the object, and Fig. 4(d) and 4(e) are corresponding to the edge region and corner region of projection. Such that the projection-scattering mapping relationship based on the structural information can be found.

2.1.3 Projection-scattering modeling

As far as we know, the projection information is used twice in the full text, first of all, the notion of using the aid of projection to reflect the structural information is suitable, since the scattering information has a large relationship with the structure of the object which has been reflected by the projection. Secondly, the validity and correctness of this hypothesis that scattering can be modeled into the convolution of a kernel with projection have been confirmed in references mentioned in [26, 27]. Generally, the analytical model can be expressed as: $I = P + S = P + ζ (u, v) \otimes P$ (8) where I , P , and S denotes the total intensity, the projection intensity, and the scattering intensity, respectively. Here, ζ (u, v) is the convolution kernel for the estimation of scattering.

However, kernel selection complicates the scattering estimation, and more sophisticated kernels may require additional assumptions which may not be always valid in general situations. To ease the process, a simple, nonlinear, and projection-based model for scattering estimation is proposed [27], which is shown as: $S (i, j) = a * P (i, j)^{b}$ (9) where S (i, j) denotes the scattering distribution at the pixel (i, j), and P (i, j) is the projection. a and b are model parameters which depict the overall magnitude adjustment and local adjustment, respectively.

Based on the structural information and the model above, the projection-scattering (PS) information at a discrete angle position is fitted. Then, scattering can be obtained and estimated by parameters of different regions, and the final scattering field S _T can be represented as the combination of scattering in different regions, such as: $S_{T} = S_{background} + S_{flat} + S_{edge + corner}$ (10)

2.2 Evaluation methods

To verify the validity, the categories of projections (named #Wheel and #Cylinder with the material of titanium) were employed to validate the performance of the proposed method. Projections were obtained from 360° scan using the FPD as XRD 1621 AN15 ES from PerkinElmer. The x-ray tube of the CBCT system was MXR-451HP/11 from Comet, and was operated at 420kVp and 0.18 mA. The projection is of 512×512 pixels and pixel size of 0.2×0.2 mm². Figure 5 shows the photograph of the phantom, respectively.

Fig.5

The objects and the projections of #Wheel and #Cylinder. (a) is the objects of #Wheel and (b) is the objects of # Cylinder.

For quantitative analysis of the image, we use the root mean square error (RMSE), peak signal to noise ratio (PSNR), and universal quality index (UQI) to assess and compare the algorithm results. These are defined below.

$RMSE = \sqrt{\frac{1}{M \times N} \sum_{0 ⩽ i ⩽ N} \sum_{0 ⩽ j ⩽ M} (X_{ij} - X_{ij}^{*})^{2}}$ (11)

$PSNR = 20 {log}_{10} (\frac{MAX ({\vec{X}}^{*})}{RMSE (\vec{X}, {\vec{X}}^{*})})$ (12)

$UQI = \frac{4 σ_{\vec{X} {\vec{X}}^{*}} \bar{\vec{X}} \times {\bar{\vec{X}}}^{*}}{(σ_{\vec{X}}^{2} + σ_{{\vec{X}}^{*}}^{2}) ((\bar{\vec{X}})^{2} + ({\bar{\vec{X}}}^{*})^{2})}$ (13)

where: •

X _ij and $X_{ij}^{*}$ denote pixel values of the original and reconstructed image, respectively;

•

M and N denote the total number of columns and rows in the image;

•

$σ_{\vec{X}}$ , $σ_{{\vec{X}}^{*}}$ and $σ_{\vec{X} {\vec{X}}^{*}}$ are the variances and covariance of X _ij and $X_{ij}^{*}$ ;

•

$\bar{\vec{X}} \times {\bar{\vec{X}}}^{*}$ in Equation (13) means the average value of the original image multiplied by the average value of the reconstructed image [30].

3 Results

To obtain the prior knowledge of the model for scattering, projections and corresponding scatterings in different positions are selected and calculated. Here, 6 sampling angles are selected for experimental in range of the 0∼360° evenly (such as 0°, 60°, 120°, 180°, 240°, 300°). Table 1 shows the model parameters information of scattering in each region at selected locations.

Table 1
Parameter information of different regional

No. Location Background Flat region Boundary region Corner region

1 0° a = 20.87 b = 0.47 a = 16.90 b = 0.51 a = 18.20 b = 0.52 a = 20.91 b = 0.56

2 60° a = 14.74 b = 0.51 a = 14.09 b = 0.52 a = 10.53 b = 0.58 a = 9.65 b = 0.63

3 120° a = 13.72 b = 0.51 a = 16.34 b = 0.51 a = 9.71 b = 0.60 a = 9.89 b = 0.63

4 180° a = 15.40 b = 0.50 a = 13.39 b = 0.53 a = 10.49 b = 0.57 a = 11.04 b = 0.62

5 240° a = 17.75 b = 0.49 a = 17.77 b = 0.50 a = 13.33 b = 0.55 a = 15.18 b = 0.63

6 300° a = 17.71 b = 0.49 a = 19.54 b = 0.49 a = 13.66 b = 0.55 a = 16.05 b = 0.60

No.	Location	Background	Flat region	Boundary region	Corner region
1	0°	a = 20.87	b = 0.47	a = 16.90	b = 0.51	a = 18.20	b = 0.52	a = 20.91	b = 0.56
2	60°	a = 14.74	b = 0.51	a = 14.09	b = 0.52	a = 10.53	b = 0.58	a = 9.65	b = 0.63
3	120°	a = 13.72	b = 0.51	a = 16.34	b = 0.51	a = 9.71	b = 0.60	a = 9.89	b = 0.63
4	180°	a = 15.40	b = 0.50	a = 13.39	b = 0.53	a = 10.49	b = 0.57	a = 11.04	b = 0.62
5	240°	a = 17.75	b = 0.49	a = 17.77	b = 0.50	a = 13.33	b = 0.55	a = 15.18	b = 0.63
6	300°	a = 17.71	b = 0.49	a = 19.54	b = 0.49	a = 13.66	b = 0.55	a = 16.05	b = 0.60

Figure 6 shows the variation tendency of the model parameters at different locations within different structural region. In the sub-graph above, the parameter tendency is basically the same that the variation of a is opposite to that of b . When parameter a increases, parameter b decreases; when parameter a decreases, parameter b increases. To obtain the model information in all angular sequence, the prior coefficient knowledge is interpolated with a spline. Then, the corresponding model knowledge for scattering estimation in the angular sequence is completed. The results are shown in Fig. 7.

Fig.6

Variation of parameters in different structural regions. (a) ∼ (d) correspond to the different regions, including the background region (a), the object flat region (b), the edge region (c), and the intersection region (d).

Fig.7

Model parameters obtained by spline interpolation

According to the model parameter at the angular sequence, the scattering in the different structural region is estimated in Equation (9), including the background, flat, boundary and corner region, respectively. By using S _θbackground, S _θflat, S _θboundary, S _θcorner to represent the scattering estimation, the final scattering distribution S _θ (x) of regional feature is obtained. $S_{θ} (x) = S_{θ background} + S_{θ flat} + S_{θ boundary} + S_{θ corner}$ (14)

Due to the inevitable model estimation error, and to reduce the numerical mutation in different regions, the overall scattering distribution is smoothed as well. The scattering distribution obtained by a grating blocker in all angle range is used as the reference. we compare the proposed scattering estimation with the non-region fitting estimation method, As shown in Fig. 8, the scattering distribution that is obtained by different estimation algorithms is compared (here we randomly select the image at the position of 143°). The first row is corresponding to experiment #Wheel, and the second row is corresponding to experiment # Cylinder. The characteristic information of the scattering estimation is prominent through fitting the undistinguished data as a whole, which is easy to weaken the self-structure of the projection with subtraction.

Fig.8

The scattering distribution obtained by different scattering estimation algorithms. (a) ∼ (c) is the experimental #Wheel, and (d) ∼ (f) is the experimental #Cylinder, respectively. The first column shows the grating scattering; the second column shows the fitting scattering, and the third column shows the proposed structural scattering. The display window is: [1200, 6420].

Then, the estimated scattering is used to correct the projection. Figure 9 shows the slice images reconstructed from the corrected projections of #Wheel (the first row) and #Cylinder (the second row) in cross view. The first column is the reference slice, and the second column shows the uncorrected slice. It is obvious that the image has been contaminated, seriously, where the contour sharpness is not clear enough. The third and fourth column are the corrected image with a different method, where the third column is the results using the curve corrected, and the fourth column is the results using the proposed above, it can be seen that both the scattering correction are useful. And the fourth column has a better performance than the middle two columns.

Moreover, Fig. 10 plots the grayscale of the reconstruction image with different scattering correction, along the line L1 and L2, respectively. The grayscale of the image with the region scattering correction has superior results to the non-region scattering correction, which shows a higher contour clarity.

Fig.9

The reconstructed slice from the scattering corrected projection. (a) ∼ (d) are corresponding to the experimental #Wheel, (e) ∼ (h) are corresponding to the experimental #Cylinder. Display window: [– 0.02, 0.18].

Fig.10

The profile of the reconstructed slice from scattering corrected projection. (a) is the grayscale along the L1 in the figure; (b) is the grayscale along the L2 in the figure.

Table 2 shows the numerical comparison of the results in Fig. 9 (the quantitative results were conducted on the whole image). The RMSE of the proposed method for the slice was dramatically reduced by 21.83% and 39.85% of that without scattering correction, respectively. PSNR and UQI also indicate better uniformity of the results, where the PSNR was increased by 7.42% and 56.77%, UQI was increased by 70.82% and 262.36% for experimental #Wheel and #Cylinder, respectively.

Table 2

Image quality comparison with scattering corrected

Experimental	Method	RMSE		PSNR		UQI
	value	improved/%	value	improved/%	value	improved/%
#Wheel	original	0.1292	/	17.7767	/	0.4412	/
	non-region fitting	0.1149	11.07	18.7945	5.73	0.6032	36.71
	the proposed	0.1010	21.83	19.0963	7.42	0.7537	70.82
#Cylinder	original	0.4085	/	7.7759	/	0.0890	/
	non-region fitting	0.2835	30.60	11.5839	48.97	0.2508	181.80
	the proposed	0.2457	39.85	12.1906	56.77	0.3225	262.36

It can be seen that the quality of the reconstructed image after scattering correction is greatly improved, and the results with the proposed algorithm are higher than that of the non-region fitting algorithm.

4 Discussion

4.1 Effect of interpolation in proposed scattering estimation

This part mainly analyzes the error of scattering estimation of different algorithms. Similarly, taking the grating scattering distribution as a reference, the ratio of the scattering obtained by different methods to the reference was determined. Here, the logarithmic transformation was taken as the error image to analyze the entire scattering estimation. The calculation is defined as follows: $ξ = log (\frac{S_{estimated_scttering}}{S_{grating_scattering}})$ (15) where S _{estimated_scttering} denotes the scattering estimation of the model, S _{grating_scattering} represents the grating scattering distribution, ξ is used as the scattering estimation error.

As is showed in Fig. 11, the scattering distribution at a random angle, such as the position 55° is selected, the results show the approximation or anastomosis of different scattering estimates with the grating scattering distribution. The first row, including Fig. 11(a∼b), plot the error distribution of the scattering estimation of phantom named #Wheel, the second column, including Fig. 11(c∼d), plot the error distribution to the phantom named #Cylinder. Obviously, the scattering estimation error using the non-region fitting is serious, while the error using the proposed algorithm has been relieved, which shows a small amplitude and slower change in error.

Fig.11

Comparison of error images of different scattering estimation methods. (a) ∼ (d) plot the error distribution of Wheel and Cylinder, respectively.

To illustrate the effectiveness and practicability of this algorithm, a far more quantitative description, and comparison are given. In Fig. 12, the Box-whisker Plot is used. Figure 12(a∼b) show the statistics on scattering data of three different methods, which is corresponding to the two different objects in Fig. 5. It is obvious that the dispersion behavior of the region fitting scattering estimation is consistent with the one in grating scattering, but there is a certain difference for the non-region fitting scattering estimation. Meanwhile, Fig. 12(c∼d) plot the profile along the line L 3 and L 4 in Fig. 11. Clearly, the estimated error of structural scattering is much smaller than that of non-region fitting scattering estimation, and it is closer to the scattering distribution with the grating as the reference. Therefore, the structural scattering estimation model can be used to estimate the scattering in angular sequence.

Fig.12

Error profile of different scattering estimation. (a) ∼ (b) plot the statistical characteristics of different scattering data, and they are corresponding to #Wheel and #Cylinder, respectively, (c) plots the relative error profile of #Wheel along the line L3, (d) plots the relative error profile of #Cylinder along the line L4.

4.2 Parameters analysis of modeling

In order to analyze the model coefficient tendency, the experimental of different materials of aluminum (2.7 g/cm³), titanium (4.5 g/cm³), and iron (7.86 g/cm³) were scanned under the same exposure. With the scattering as S (i, j) in Equation (9), the parameters a and b can be solved using linear regression. Figures 13 and 14 show the projection-scattering fitting of #Wheel and #Cylinder in a different region, respectively.

Fig.13

Scattering distribution trends with different densities in different areas of #Wheel. (a) is corresponding to the flat region of background; (b) is corresponding to the flat region of object projection; (c) is corresponding to the edge and corner region of object projection, all of which based on different materials.

Fig.14

Scattering distribution trends with different densities in different areas of #Cylinder. (a) is corresponding to the flat region of background; (b) is corresponding to the flat region of object projection; (c) is corresponding to the edge and corner region of object projection, all of which based on different materials.

Obviously, the scattering of projection corresponding to different regions have a different tendency. The scattering in the background is more serious, the scattering in the flat region increase smoothly, and the scattering in the edge and corner region shows disunity. From the plots in Figs. 13 and 14, it can be seen that the scattering rises with the material density increasing, but different areas are not the same. So, the parameters fitted by modeling are slightly changed.

Form the information showed in Fig. 7, that the modeling parameters is a function of angle. To further explore the behavior and laws of the model, the scattering can be rewritten as $S_{θ} (x) = {ax}_{θ}^{b}$ (16)

Taking logarithmic transformation on both sides in Equation (16), the expression can be expressed as: $ln (S_{θ} (x)) = b ln x_{θ} + ln a$ (17) where x _θ represents the projection value at the position θ, and S _θ (x) represents the scattering value at the position θ. If ln a = A, b = B, then the coefficients A and B have the following relations: $A = ln (S_{θ} (x)) - B ln x_{θ}$ (18)

If the projection value and scattering value at the position θ are known, the linear relation in Equation (18) exists between coefficients A and B. It can be seen that the coefficient A decreases with the increase of B, which is consistent with the tendency of parameters in different regions in Figs. 6 and 7. As a knowledge, scattering parameters can be obtained by interpolation between sampling locations.

5 Conclusion

An effective scattering estimation in the angular sequence has been described in this paper. Based on the estimated scattering completed by interpolating the model parameters of prior scattering obtained by sparse sampling, the performance of the proposed algorithm has been evaluated and verified by using experimental #Wheel and #Cylinder on CBCT system. In the study, the results show a significant improvement in the RMSE, PSNR, and UQI of the images reconstructed from scattering corrected projection. In comparison, we have got the approximation or anastomosis of corrected projection as the full angle grating scattering process, by interpolating the sparse angle scattering parameters. This can reduce the scanning cost and perform dense scattering estimation on the scanning angle of interest, and improve the accurate correction of the projection further.

While scattering is not the only factor that affects image quality. Beam hardening which couples with scattering is also an important factor to affect and distort the projection. So, there are the cupping artifacts on the reconstructed image in Fig. 9. However, the correction of beam hardening is usually after the scattering estimation, which does not affect the estimated results of scattering by the method that the paper proposed.

Footnotes

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 (MJ-2017-F-05). Project funded by China Postdoctoral Science Foundation (2019M653749). The authors are grateful to the anonymous reviewers for their valuable comments.

References

Yang

F.Q.

Zhang

D.H.

Huang

K.D.

Gao

Z.Z.

and Yang

Y.F.

, Incomplete projection reconstruction of computed tomography based on the modified discrete algebraic reconstruction technique, Measurement Science & Technology 29(2) (2018), 025405.

Lan

Liu

Xiao

and Yao

, Internal-illumination photoacoustic computed tomography, Journal of Biomedical Optics 23(3) (2018), 1–4.

Huang

K.D.

Zhang

Shi

Y.K.

Zhang

and Xu

, Scattering correction method for cone-beam CT based on interlacing-slit scan, Chin Phys B 23(9) (2014), 098106.

Yang

F.Q.

Zhang

D.H.

Huang

K.D.

Yang

Y.F.

and Liao

J.M.

, Image artifacts and noise reduction algorithm for cone-beam computed tomography with low-signal projections, J Xray Sci Technol 26(2) (2018), 227–240.

Siewerdsen

J.H.

and Jaffray

D.A.

, Optimization of x-ray imaging geometry (with specific application to flat-panel cone-beam computed tomography, Med Phys 27(8) (2000), 1903–1914.

Zhao

Vernekohl

Zhu

Wang

and Xing

, A model-based scatter artifacts correction for cone beam CT, Med Phys 43(4) (2016), 1736–1753.

Min

Pua

Kim

Han

and Cho

, Analytic image reconstruction from partial data for a single-scan cone-beam CT with scatter correction, Med Phys 42(11) (2015), 6625–6640.

Star-Lack

Sun

Kaestner

Hassanein

Virshup

Berkus

and Oelhafen

, Efficient scattering correction using asymmetric kernels, Proc SPIE 7258 (2009), 1Z.

Maltz

J.S.

Blanz

W.E.

Hristov

and Bani-Hashemi

, Cone beam x-ray scattering removal via image frequency modulation and filtering, Med Phys 2(6) (2005), 1854–1857.

10.

Bae

and Ye

J.C.

, Translational motion correction algorithm for truncated cone-beam CT using opposite projections, J Xray Sci Technol 25(6) (2017), 927–944.

11.

Lee

Xing

Lee

and Fahimian

, Scatter correction in cone-beam CT via a half beam blocker technique allowing simultaneous acquisition of scatter and image information, Med Phys 39(5) (2012), 2386–2395.

12.

Zhang

Jacobs

and Bosmans

, A model-based volume restoration approach for Monte Carlo scatter correction in image reconstruction of cone beam CT with limited field of view, Proc SPIE (2013), p. 86682.

13.

Zhu

Strobel

and Fahrig

, X-ray scattering correction for cone beam CT using moving blocker array, Proc SPIE 5745 (2005), 251.

14.

Ning

Tang

and Conover

, X-ray scattering correction algorithm for cone beam CT imaging, Med Phys 31(5) (2004), 1195–1202.

15.

Ruhrnschopf

E.P.

and Klingenbeck

, A general framework and review of scatter correction methods in X-ray cone-beam computerized tomography. Part 1: Scatter compensation approaches, Med Phys 38(7) (2011), 4296–4311.

16.

Ruhrnschopf

E.P.

and Klingenbeck

, A general framework and review of scatter correction methods in cone beam CT. Part 2: Scatter estimation approaches, Med Phys 38(9) (2011), 5186–5199.

17.

Yan

Mou

and Tang

, et al., Projection correlation based view interpolation for cone beam CT: Primary fluence restoration in scatter measurement with a moving beam stop array, Phys Med Biol 55(21) (2010), 6353–6375.

18.

Sisniega

Zbijewski

Badal

, et al., Monte Carlo study of the effects of system geometry and antiscatter grids on cone-beam CT scatter distributions, Med Phys 40(5) (2013), 051915.

19.

Rong

J.Y.

Liu

W.L.

and Gao

, et al., The review of scatter suppression methods in cone beam CT, CT Theory and Applications 25(2) (2016), 235–250.

20.

Ouyang

Song

Solberg

and Wang

, A moving blocker system for cone-beam computed tomography scatter correction, Med Phys 40(7) (2013), 071903.

21.

Niu

and Zhu

, Scatter correction for full-fan volumetric CT using a stationary beam blocker in a single full scan, Med Phys 38(11) (2011), 6027–6038.

22.

Huang

, Scatter measurement and correction method for cone-beam CT based on single grating scan, Optical Engineering 56(6) (2017), 064106.

23.

Briggs

N.T.

Slade

W.H.

Boss

and Perry

M.J.

, Method for estimating mean particle size from high-frequency fluctuations in beam attenuation or scattering measurements, Applied Optics 52(27) (2013), 6710.

24.

Shi

W.L.

, Research on the scatter measurement and analysis for CBCT using a grating, NWPU (2016), 43–50.

25.

Yang

F.Q.

Zhang

D.H.

Huang

K.D.

Liao

J.M.

and Wang

X.Y.

, Scattering estimation for cone-beam CT using local measurement based on compressed sensing, IEEE Transactions on Nuclear Science 65(3) (2018), 941–949.

26.

Sun

and Star-Lack

J.M.

, Improved scatter correction using adaptive scatter kernel superposition, Phys Med Biol 55 (2010), 6695–6720.

27.

Meng

Lee

Xing

and Fahimian

, Single-scan patient-specific scatter correction in computed tomography using peripheral detection of scatter and compressed sensing scatter retrieval, Med Phys 40(1) (2013), 011907.

28.

Yan

W.J.

Yang

Wang

H.X.

and Zhang

X.H.

, Corner detection algorithm based on B-spline multi-scale boundary structure tensor, Computer Technology & Development 4 (2010), 80–83.

29.

Alizadeh

Conklin

C.J.

Middleton

D.M.

Shah

Saksena

Krisa

Finsterbusch

Faro

S.H.

Mulcahey

M.J.

and Mohamed

F.B.

, Identification of ghost artifact using texture analysis in pediatric spinal cord diffusion tensor images, Magnetic Resonance Imaging 47 (2018), 7–15.

30.

Hsieh

C.J.

Huang

T.K.

Hsieh

T.H.

Chen

G.H.

Shih

K.L.

Chen

Z.Y.

Chen

J.C.

and Chu

W.C.

, Compressed sensing based CT reconstruction algorithm combined with modified Canny edge detection, Phys Med Biol 63 (2018), 155011.