Virtual scatter modulation for X-ray CT scatter correction using primary modulator

Abstract

A new scatter estimation algorithm with a concept of virtual scatter modulation for X-ray scatter correction using primary modulator is proposed to reduce the aliasing errors in the estimated scatter. Virtual scatter modulation can be realized through dividing the measured primary-modulated image by the measured modulation function. After the division, the aggravation of the aliasing of primary due to the non-uniformity of the modulation function is largely transferred to that of scatter. Since scatter in general has less high frequencies than primary does, the aggravation of its aliasing is expected to be weaker, and therefore the overall aliasing can be reduced. A CatPhan©600 phantom and an anthropomorphic thorax phantom are scanned on a tabletop X-ray cone-beam computed tomography system to validate our proposed algorithm. On the Catphan phantom, the oscillations that are clearly observed in the central region of the Catphan scatter profile estimated using the original primary-modulation algorithm, are mostly eliminated with the proposed scatter modulation algorithm, leading to less residual artifacts and better CT number uniformity in the reconstructed image. Compared with 38.9 HU of CT nonuniformity in a selected uniform region when the primary-modulation algorithm is used, the new algorithm significantly reduces it to 4.5 HU, reaching the same level of uniformity as the ground truth reference. On the thorax phantom, overall better CT number uniformity is also achieved.

Keywords

Primary modulator scatter modulation scatter correction cone-beam CT

1 Introduction

Scatter is a critical issue in X-ray cone-beam computed tomography (CBCT). It may cause cupping, shading and streaks artifacts, and may also result in CT number shift, as well as image contrast degradation [1, 2]. Scatter correction is an active and challenging topic in research, especially in X-ray imaging using large area detectors [3–10]. Recently, a promising scatter correction method that uses the concept of primary modulation has been developed [11, 12]. In the primary modulation method, a checkerboard pattern consisting of attenuating blockers is inserted between the X-ray source and the object as the primary modulator. With such a primary modulator in place, part of the primary therefore is modulated to high frequency, while scatter still has dominant low-frequency components, making the modulated primary and the scatter highly separable in the frequency domain. Consequently, an accurate scatter estimation can be obtained using filtering techniques. This method is able to achieve effective scatter correction using one single scan without loss of real-time imaging capabilities or increase in patient dose. Validation and evaluation of the proposed method were previously done using Monte Carlo (MC) simulation and phantom studies in Refs. [12–14]. In order to optimize the scatter correction performance, we also investigated considerations of primary modulator design [15, 16].

In practical implementations, aliasings (i.e., frequency overlapping) between the modulated primary and the unmodulated primary and scatter are sources of errors in scatter estimation. The aliasings are rooted in the limited modulation frequency and aggravated by the nonuniformity of the modulation function. Since the primary modulator is inserted between the source and the object, the finite size of focal spot and the geometric magnification remarkably confine the minimum spacing of the modulator blocker. It causes that the achievable modulation frequency is always lower than the sampling frequency of the detector. As a result, scatter estimation is actually carried out on the downsampled projection images [12]: the projection images are downsampled at the centers of the modulator blocker shadows and non-blocker shadows. In the downsampled images, the modulated primary shares the same frequency band with the unmodulated primary and scatter, making the aliasings occur very likely. If a spatially nonuniform modulation function is used, which is true in practice due to the heel effect of the X-ray tube and the beam hardening effect of the primary modulator, the aliasings would be aggravated because the nonuniformity of the modulation function virtually adds more high-frequency components into primary.

The aggravation of aliasing of primary caused by the nonuniformity of the modulation function is not considered in the original algorithm of scatter estimation derived from the principle and hypothesis of primary modulation [12], in which, to estimate a scatter distribution, the modulated primary is first demodulated back, and then scaled and subtracted from the unmodulated primary and scatter. Further improvements have been explored in recent years [17–19], aiming at replacing the original Fourier space filtering technique by non-Fourier space and regularized optimization ones. A clinical study on angiography C-arm CT system has also been reported [20, 21].

In this paper, we still focus on Fourier space filtering-based scatter estimation as it is a unique property of the primary modulator method, computationally efficient and easy to implement. With a concept of virtual scatter modulation, a new algorithm of scatter estimation is derived: the measured X-ray image is first divided by the measured nonuniform modulation function so that a scatter modulation is virtually realized; and then the modulated scatter is demodulated back and scaled to obtain a scatter estimate. After the division, the aggravation of aliasing of primary is transferred to that of scatter. Since scatter in general has less high frequencies than primary does, the aggravation of its aliasing is expected to be weaker, and therefore the overall aliasing can be reduced.

This paper is organized as follows. Section 2 reviews the primary modulation method, where aliasing errors in the original algorithm of scatter estimation are analyzed. The new algorithm using the concept of scatter modulation is proposed in Section 3. Section 4 validates the proposed algorithm on a tabletop CBCT system using a Catphan©600 phantom and an anthropomorphic thorax phantom. Finally, we summarize this paper in Section 5 discussion with a brief discussion.

2 Aliasing errors in primary modulator method

2.1 Primary modulation using primary modulator

In the primary modulation method, a high-spatial-frequency pattern board consisting of attenuating blockers is placed between the X-ray source and the scanning object as a primary modulator [12]. For X-rays passing through the blockers, their intensities are attenuated by the blocker, with a transmission ratio, α, predetermined by the modulator material and thickness. Let p and s be primary and scatter in the spatial domain, respectively, then the measured primary-modulated X-ray image with such a primary modulator in place can be expressed as, $p^{'} (i, j) = p (i, j) m (i, j) + s (i, j),$ (1) where, m (i, j) is the modulation function of the primary modulator. Parameters i and j are the horizontal and vertical indices in the spatial domain. For an ideal modulator, m (i, j) = α, for (i, j) locating within the blocker on the detector; otherwise m (i, j) =1.

According to the derivations in Refs. [12, 14], a scatter distribution in the Fourier domain can be estimated from the primary-modulated image using filtering techniques as $S_{PM} (ω) = P^{'} (ω) H (ω) - \frac{1 + α}{1 - α} P^{'} (ω - π) H (ω),$ (2) where, P′, P and S are the Fourier transforms of p′, p and s, respectively, and ω ∈ [- π, π] × [- π, π] denotes the 2D coordinates of (ω_x, ω_y) in the Fourier domain; while H (ω) is a low-pass filter.

2.2 Aliasing in the primary-modulated image and the related errors

In the primary modulation method, aliasings in the primary-modulated image are sources of errors in scatter estimation. Figure 1 illustrates the assumed ideal situation of primary modulation where primary and scatter are band-limited, and the realistic situation with aliasings occurring between the modulated-primary and the unmodulated primary and scatter, respectively. The aliasings are rooted in the limited modulation frequency and aggravated by the nonuniformity of the modulation function in practical implementations. Since the minimum spacing of the modulator blocker is limited by the finite size of focal spot and the geometric magnification, the achievable modulation frequency of the primary modulator is always lower than the sampling frequency of the detector. Consequently, scatter is initially estimated from the downsampled primary-modulated images. In the downsampled images, the modulated primary shares the same frequency band with the unmodulated primary and scatter, making the aliasings occurs very likely.

Fig.1

An illustration of primary-modulated image in the Fourier domain. (a) The ideal situation; (b) the realistic situation.

If a spatially nonuniform modulation function is used, which is true in practice due to the heel effect of the X-ray tube and the beam hardening effect of the primary modulator, the aliasings would be aggravated. Let n (i, j) be the nonuniformity of the modulation function in the spatial domain, then the realistic modulation function should be m (i, j) = m₀ (i, j) n (i, j). Accordingly, the estimated scatter in Eq. (2) actually is [14], ∥ $\begin{matrix} S_{PM} (ω) & = & S (ω) H (ω) - \frac{1 + α}{1 - α} S (ω - π) H (ω) - \frac{2 α}{1 - α} (P * N) (ω - π) H (ω) . \end{matrix}$ (3) where, N (ω) is the Fourier transform of n (i, j), “∗” denotes the convolution operator. The last two terms in Eq. (3) are the contributions of aliasings, $\begin{matrix} {AE}_{PM} (ω) & = & \frac{1 + α}{1 - α} S (ω - π) H (ω) + \frac{2 α}{1 - α} (P * N) (ω - π) H (ω) . \end{matrix}$ (4)

Since the effective primary, P ∗ N, in scatter estimation generally has more high frequencies than the original primary, P, does, the estimated scatter in Eq. (3) contains more aliasing errors when a nonuniform modulation function is applied.

3 Scatter estimation using virtual scatter modulation

3.1 Virtual scatter modulation

Let p″ (i, j) be the result of the measured transmission image p′ (i, j) divided by the measured modulation function m (i, j) = m₀ (i, j) n (i, j), then from Eq. (1), we have $p^{″} (i, j) = \frac{p^{'} (i, j)}{m (i, j)} = p (i, j) + \frac{s^{'} (i, j)}{m_{0} (i, j)},$ (5) where, s′ (i, j) denotes the effective scatter distribution, s (i, j)/n (i, j). Equation (5) is equivalent to virtually modulating the effective scatter distribution s′ (i, j), with a new ideal modulation function 1/m₀ (i, j).

Following the filtering techniques used in the primary-modulation algorithm above, we can obtain an estimate of the effective scatter distribution as $S_{SM}^{'} (ω) = - \frac{2 α}{1 - α} P^{″} (ω - π) H (ω) .$ (6) where, P″ and S′ are the Fourier transforms of p″ and s′, respectively. Here, it is assumed that S′ (ω) still has dominant low-frequency components.

As s (i, j) = s′ (i, j) n (i, j) [i.e., S (ω) = S′ (ω) ∗ N (ω) in the Fourier space]. an estimate of the original scatter can be obtained as $S_{SM} (ω) = - \frac{2 α}{1 - α} [P^{″} (ω - π) H (ω)] * \tilde{N} (ω),$ (7) where, $\tilde{N} (ω)$ is the Fourier transform of the smoothed n (i, j). The smoothing of n (i, j) is to make sure that the convolution does not introduce significant high-frequency components.

Figure 2 illustrates the main procedures in scatter estimation using the primary-modulation algorithm and the scatter-modulation algorithm, respectively. For the primary-modulation algorithm, scatter is estimated by subtracting the demodulated primary in the primary-modulated image from the unmodulated primary and scatter; while for the scatter-modulation algorithm, scatter is estimated by demodulating scatter in the scatter-modulated image obtained from dividing the measured primary-modulated image by the measured modulation function. They also use different weighting factors.

Fig.2

Major procedures of scatter estimation for the primary-modulation algorithm (a) and for the scatter-modulation algorithm (b), respectively.

Please be noted that for the division in Eq. (5), no geometrical change between the scan of object and the scan of modulator is considered. In case of significant system mechanical vibrations, a more sophisticated calibration on gantry rotation may be required, same as other scatter estimation algorithms. Recently, a modulator-database based registration and interpolation approach has been developed to resolve this problem [20]. Also, a gantry-rotation calibration has been commonly deployed on diagnostic CT [2].

3.2 Aliasing errors in scatter-modulation algorithm

Like primary-modulation algorithm in Eq. (3), the estimated scatter from the scatter-modulation algorithm indeed is

\begin{matrix} S_{SM} (ω) & = & [S^{'} (ω) H (ω)] * \tilde{N} (ω) - \frac{1 + α}{1 - α} [S^{'} (ω - π) H (ω)] * \tilde{N} (ω) \\ - \frac{2 α}{1 - α} [P (ω - π) H (ω)] * \tilde{N} (ω) . \end{matrix}

(8)

The last two terms in Eq. (8) are the contributions of aliasings, $\begin{matrix} {AE}_{SM} (ω) & = & \frac{1 + α}{1 - α} [S^{'} (ω - π) H (ω)] * \tilde{N} (ω) + \frac{2 α}{1 - α} [P (ω - π) H (ω)] * \tilde{N} (ω) . \end{matrix}$ (9)

Unfortunately, the scatter-modulation algorithm cannot bypass the inevitable aliasing errors either. However, compared with the aliasing errors of the primary-modulation algorithm in Eq. (4), the impact of the nonuniformity of the modulation function is mainly transferred from primary to scatter when the scatter-modulation algorithm is used (recall that S′ (ω) is the Fourier transform of the effective scatter, s (i, j)/n (i, j)). Since scatter in general has less high-frequency components than primary does, the increase in high frequencies due to the nonuniformity of the modulation function is generally smaller for scatter than for primary, as illustrated in Fig. 3. As a result, the overall aliasings can be reduced with the proposed scatter-modulation algorithm.

Fig.3

An illustration of the high-frequency increase caused by the nonuniformity of the modulation function. Since scatter in general has less high-frequency components than primary does, the increase in the high-frequency region due to the nonuniformity of the modulation function is expected to be smaller than that of primary.

It is worth noting that, with a uniform (ideal) modulation function, (i.e., n (i, j) =1, N (ω) = δ (ω) and $\tilde{N} (ω) = δ (ω)$ ), S_PM (ω) in Eq. (3) and S_SM (ω) in Eq. (8) are identical. In such cases, both the primary modulation algorithm and the scatter modulation algorithm provide the same scatter estimation. In our analysis above, we do not consider the beam hardening of the scanning object, whose impact to modulation function highly depends on the modulator material [16, 22].

4 Evaluations and results

4.1 The tabletop cbct system

A tabletop CBCT system was used with the system parameters summarized in Table 1. The system consisted of a CPI Indico 100 100 kW programmable X-ray generator (CPI Communication & Medical Products Division, Georgetown, Ontario, Canada), a Varian G-1590SP X-ray tube [Varian X-ray Products (Now Varex), Salt Lake City, UT], a rotation stage, a Huestis collimator (Huestis Medical, Bristol, RI), a Varian PaxScan 4030CB flat panel a-Si large area X-ray detector, and a workstation. The X-ray tube had a tungsten target with a 12 deg of target tilt angle, and an inherent filtration of 1.0 mm Al plus 3 mm of equivalent aluminum due to the collimator. No anti-scatter grid was mounted in front of the detector.

Table 1
Imaging parameters of the physical experiments

X-ray source 0.6 mm nominal focal spot size

120 kVp, 35 mA, 13 ms pulse width

Detector size 397.3×298.0 mm²

1024×768 pixels

Anti-scatter grid none

Rotation circular, 360 deg

Phantom scanned Catphan Thorax

Source-to-detector distance (mm) 1423 1202

Source-to-axis distance (mm) 860 786

Source-to-modulator distance (mm) 231 231

Number of views 800 625

4.2 The primary modulator

During the experiments of the primary modulation method, a one-dimensional (1-D) erbium primary modulator and a two-dimensional (2-D) copper modulator (shown in Fig. 4) were mounted in front of the collimator and used for the Catphan© 600 phantom and the thorax phantom, respectively, with a nominal distance to the X-ray focal spot of 231 mm.

Fig.4

Two primary modulators used in this paper. (a) The 1-D erbium modulator (used for the Catphan phantom): erbium foils (25.4 μm of thickness) were cut into strips and glued upon a uniform aluminum oxide substrate (0.5 mm of thickness), with strip width and spacing being 1.588 mm; (b) the 2-D copper modulator (used for the thorax phantom): checkerboard copper patterns (210 μm of thickness) were manufactured by printed circuit board (PCB) method, with blocker size and spacing being 0.889 mm.

The man-made 1-D erbium modulator consisted of 25.4 μm thick equi-distantly spaced strips. The width and spacing of the strips are both 1.6 mm. The resulting modulation period on the detector due to the geometrical magnification was approximately 12.4 mm. The measured transmission factor of the erbium strip was 81% at 120 kVp. As described in Ref. [23], a 1-D primary modulator aligned parallel to the horizontal direction of the detector, which modulates primary along the vertical direction, is able to avoid ring artifacts in reconstructed images. The material of the 1-D primary, erbium, is specifically selected to minimize beam hardening of the modulator [16]. Compared with a two-dimensional (2D) copper or aluminum primary modulators used in Refs. [12, 14], the 1-D erbium primary modulator is more suitable for the Catphan phantom to demonstrate the performances of the scatter estimation algorithms, because other factors potentially impacting the performance of scatter correction, such as the ring artifacts and the beam hardening effect of the modulator, are already suppressed.

The 2-D copper modulator was built on a circuit board using a conventional etching method. The copper blocker on the modulator has a side length of 889 μm, and a thickness of 210 μm. The measured transmission of the blocker, α, is approximately 75% at 120 kVp. It was used for the thorax phantom due to the relatively more complicated structures and more higher frequency components.

4.3 Experiments and evaluations

Both the primary-modulation algorithm and the scatter-modulation algorithm were evaluated on a CatPhan©600 phantom (The Phantom Laboratory, Salem, NY) with the 1-D erbium primary modulator, and an anthropomorphic thorax phantom with the 2-D copper modulator, respectively. The CatPhan phantom is widely used to evaluate CT number accuracy and uniformity, contrast to noise ratio (CNR), as well as low contrast detectability and spatial resolution, making it the best choice for our study. The thorax phantom, with anatomic structures like human body, can provide helpful information close to real patient scans. The scatter estimations by the two algorithms were implemented using Eq (2) and Eq (7), respectively, where a Hamming window with ω_max = π/4 was chosen as the low-pass filter H (ω). The standard Feldkamp-Davis-Kress (FDK) algorithm [24] was used in image reconstructions. The reconstructed images were presented in Hounsfield units (HU), with a CT number of -1000 HU for air and a CT number of 0 HU for water-equivalent materials.

For the use of comparison, fan-beam scans using a narrow collimator were carried out. Due to the small volume of irradiation, the scatter signals in the narrow-collimator scans were very low, and were used as ground-truth references in the comparisons. The 1-D profile of the measured scatter that was obtained by subtracting the fan-beam image from the cone-beam image in the irradiation area, was also used to exhibit the accuracy of the estimated scatter. The magnitude of the scatter estimation error (SEE) was evaluated as $E_{SEE} = 100 % \frac{1}{N M} \sum_{i = 1}^{M} \sum_{j = 1}^{N} \frac{| s_{est} (i, j) - s_{meas}^{m} (i, j) |}{| s_{meas}^{m} (i, j) |},$ (10) where, s_est (i, j) is the spatial form of the estimated scatter; and $s_{meas}^{m} (i, j)$ is the measured scatter after noise suppressing with a low-pass filter to reduce its impact to the evaluation, respectively.

Besides side-by-side image comparisons, the quality of the reconstructed images is also quantified using the CNR and the error of CT number. The CNR was computed as $CNR (A, B) = | m_{A} - m_{B} | / \sqrt{σ_{A}^{2} + σ_{B}^{2}}$ (11) where, m_A, m_B, σ_A, σ_B are the mean and standard deviation of the CT numbers in the two neighboring regions of interest (ROIs). The error of CT number in the ROIs was calculated as the square root of the mean square error (RMSE), defined as: $E_{RMSE} = \sqrt{mean [(μ_{i} - {\bar{μ}}_{i})^{2}]}$ (12) where i is the index of the ROI, μ_i is the mean reconstructed value in HU inside the ROI, and ${\bar{μ}}_{i}$ is the corresponding value measured in the ground-truth image.

Further more, the CT number nonuniformity in a uniform regions of reconstructed images are quantified using the maximum difference of CT numbers among the selected ROIs, i.e., $Δ_{MAX} = max (μ_{i}) - min (μ_{i}) .$ (13)

4.4 The measured modulation function

The measured modulation function for the 1-D erbium modulator obtained after dark-field and flat-field corrections, is shown in Fig. 5. The mean value of the transmission factor of the erbium strips over the whole downsampled image is 0.81. The nonuniformity of the modulation function is computed by dividing the measured modulation function by a uniform modulation function whose transmission factor is the mean value (here it is 0.81). It has a mean of 0.99, and a standard deviation of 8.34 × 10^-3.

Fig.5

The measured modulation function for the 1-D erbium modulator. Display window: [0.90, 1.03]. (a) The modulation function in the downsampled image (size of 99 × 26 pixels); (b) the nonuniformity of the modulation function. The vertical profiles of the images on the right side are averaged along the horizontal direction.

4.5 Comparisons of the estimated scatter

Figure. 6 shows the profiles of the measured primary, the measured scatter, and the estimated scatter using both the primary-modulation algorithm and the scatter-modulation algorithm, respectively. For the CatPhan©600 phantom, the SPR is up to two. With the primary modulation method, generally speaking both algorithms estimate the scatter distribution quite well. However, technically the scatter-modulation algorithm outperforms the primary-modulation algorithm in terms of estimation accuracy and artifacts. For the primary-modulation algorithm, the scatter estimation error, E_SEE in Eq. (10), is about 3.6%; it is further reduced to about 2.3% when the scatter-modulation algorithm is used. The oscillations that are clearly observed in the central region of the scatter profile estimated using the primary-modulation algorithm, are mostly eliminated with the scatter-modulation algorithm.

Fig.6

The 1-D horizontal profiles of the measured scatter and the estimated scatter for the CatPhan©600 phantom. Scatter estimations using both the primary-modulation (PM) algorithm and the scatter-modulation (SM) algorithm are shown. The measured primary is also provided to imply the corresponding SPR levels.

4.6 Image reconstructions

The scatter correction performances using the two algorithms are more obvious in the comparison of the reconstructed CT images, as shown in Fig. 7. The reconstructions have a size of 512×512 pixels, with 0.4×0.4 mm² for a pixel. The slice thickness in the z direction is 2 mm. The fan-beam reconstruction shown in Fig. 7(b) is used as the reference. Without scatter correction, the reconstructed cone-beam image has severe cupping artifacts. These artifacts are significantly suppressed when scatter is corrected with either the primary-modulation algorithm or the scatter modulation algorithm. However, as also shown in Fig. 6, the scatter estimation using the primary-modulation algorithm has relatively larger errors due to the aggravated aliasing in presence of the nonuniform modulation function. Residual artifacts, such as the doming and ruffling, are present in the image after scatter correction. It is more clearly observed in Fig. 8, where the zoomed-in images in the rectangular region indicated by the dotted lines in Fig. 7(c) are shown.

Fig.7

Image reconstruction of the high contrast module CTP404 of the CatPhan©600 phantom from cone-beam scans with and without scatter correction. Reconstruction size: 512×512 pixels with 0.4×0.4 mm² for a pixel. Display window: [–300, 400] HU. (a) The cone-beam image without scatter correction; (b) the fan-beam image (reference); (c) the cone-beam image with scatter corrected using the primary-modulation algorithm; and (d) the cone-beam image with scatter corrected using the scatter-modulation algorithm.

Fig.8

The zoomed-in images in the rectangular region indicated by the dotted lines in Fig. 7(c). Reconstruction size: 480×330 pixels with 0.16×0.16 mm² for a pixel. Display window: [–67, 120] HU. (a) The fan-beam image (reference); (b) the cone-beam image with scatter corrected using the primary-modulation algorithm; and (c) the cone-beam image with scatter corrected using the scatter-modulation algorithm.

For a quantitative analysis of the reconstructed images, we measured the CT numbers, the noise [standard deviation (STD)] and the CNRs in the selected ROIs labeled in Fig. 7. The results are summarized in Table 2. For the cone-beam image without scatter correction, the RMSE [calculated using Eq. (12) with the fan-beam image as reference] is as high as 312.7 HU, and the averaged CNR is only 12.3 HU. After scatter correction using the primary modulation technique, the primary-modulation algorithm and the scatter-modulation algorithm reduce the RMSE to 19.6 and 19.8 HU, while enhance the averaged CNR to 18.7 and 18.8 HU, respectively. Note that a relatively higher noise level, which is an expected result for all post-processing scatter correction algorithms [25], is seen in the scatter corrected images (the corresponding noise rises from 11.8 HU to 17.1 and 17.0 HU, respectively). Here, a similar level of CT number accuracy and CNR is achieved for high contrast rods, which is expected because the two algorithms differ mainly in the residual aliasing error. The proposed scatter modulation algorithm shows no degradation relative to the original one, reflecting its robustness and effectiveness.

Table 2

Mean, standard deviation (STD) and CNR in the selected ROIs of the CatPhan©600 phantom in Fig. 7

ROIs		Air 1	Air 2	PMP	LDPE	Polysty.	Acrylic	Delrin^TM	Teflon	E _RMSE	Average
Mean (HU)	Cone–beam	–667.1	–676.7	–295.7	–251.6	–228.5	–151.4	–44.9	212.8	312.7	/
	Fan–beam	–929.4	–929.2	–200.9	–117.4	–66.8	82.3	290.8	867.9	/	/
	Primary–Mod.	–931.0	–933.6	–190.6	–115.4	–59.3	87.1	322.1	911.2	19.6	/
	Scatter–Mod.	–931.2	–933.5	–191.0	–116.0	–59.4	87.5	322.7	911.9	19.8	/
STD (HU)	Cone–beam	13.3	10.0	8.6	8.6	8.4	10.6	13.3	18.3	/	11.8
	Fan–beam	12.1	11.8	13.2	13.4	12.6	12.3	12.8	14.3	/	12.8
	Primary–Mod.	14.7	17.0	15.8	16.7	18.0	17.9	19.3	17.0	/	17.1
	Scatter–Mod.	14.6	16.9	15.8	16.6	18.1	17.7	19.2	16.6	/	17.0
CNR	Cone–Beam	24.1	33.2	9.7	6.2	3.6	0.6	6.4	15.0	/	12.3
	Fan–Beam	55.2	57.8	14.2	9.3	6.6	1.5	13.1	43.5	/	25.2
	Primary–Mod.	42.5	40.3	10.5	6.8	5.1	1.3	9.2	34.5	/	18.7
	Scatter–Mod.	42.3	40.3	10.5	6.8	5.1	1.3	9.2	35.0	/	18.8

To quantify the improvement in reduction of aliasing errors, we also measured the CT numbers of a relatively uniform region of the CatPhan phantom at ROIs labeled in Fig. 8. The averaged values in the selected ROIs are summarized in Table. 3, where the maximum differences, Δ_MAX, is calculated using Eq. (13) as the indicator for the CT number nonuniformity in the reconstructed images. Although the overall image quality is much improved after scatter correction using primary-modulation algorithm, the CT nonuniformity in our selected region could still be as high as 38.9 HU, due to the residual aliasing error demonstrated in Fig. 6. Using the proposed virtual scatter modulation algorithm, the CT nonuniformity is significantly reduced to 4.5 HU, at the same level of our ground truth reference (i.e., 3.3 HU).

Table 3

Averaged CT number (HU) in the selected ROIs of a uniform region as shown in Fig. 8

Method	Fan-Beam	Primary-Mod.	Scatter-Mod.
ROI 1	52.2	12.3	44.5
ROI 2	50.1	51.2	49.0
ROI 3	53.4	42.0	48.4
Δ _MAX	3.3	38.9	4.5

To further demonstrate the performance of the virtual scatter modulation, two off-center cross sections (low contrast detectability module CTP515 and high contrast resolution module CTP528) of the CatPhan phantom, along with reformatted images, are also reconstructed and displayed in Fig. 9. Shading artifacts and CT number shift are greatly reduced after scatter correction. The virtual scatter-modulation algorithm consistently outperforms the primary-modulation one with less residual artifacts as highlighted by the “arrows”. Relatively better low contrast detectability is also achieved in the subslice (central) region of the cross section. Specifically, for a contrast target with 7-mm diameter and its background [ROIs shown in Fig. 9(e)], virtual scatter-modulation algorithm has a contrast difference of 8.3 HU, compared with 1.6 HU for the primary-modulation one, which is contaminated by the residual artifact. For the module CTP528, a spatial resolution of 16 lp/cm is achieved without and with scatter correction, suggesting scatter correction barely affect spatial resolution.

Fig.9

Volumetric image reconstructions of the CatPhan©600 phantom from cone-beam scans without and with scatter correction. Reconstruction size: 512×512×320 voxels with 0.4×0.4×0.4 mm³ for a voxel. First row: reformatted images of the Catphan©600 phantom (display window: [–300, 400] HU, averaged over 2-cm thickness), where the dashed-lines in (a) indicate the locations of the cross sections; second row: cross section images of module CTP515 (display window: [–150, 350] HU, averaged over 8-cm thickness), where the low contrast targets are illustrated in (d) according to the manual of the phantom; last row: cross section images of module CTP528 (display window: [–300 2100] HU). First column: the cone-beam images without scatter correction; second column: the cone-beam images with scatter corrected using the primary-modulation algorithm; last column: the cone-beam images with scatter corrected using the scatter-modulation algorithm. Reconstruction kernels: the Ram-Lak for images in the last row and the Shepp-Logan for all the others.

Additionally, a comparative study on an anthropomorphic thorax phantom is also conducted using the 2-D copper modulator, with reconstructed CT images shown in Fig. 10, where a ring correction [26] is applied to suppress the ring/band artifacts caused by the insertion of the 2-D primary modulator. Compared with no scatter correction, scatter corrected transaxial and reformatted (sagittal) images overall have better image quality, with much less shading/streak artifacts, and no obvious CT number shift (taken fan-beam scan as reference). The subtractions of images with primary-modulation algorithm from that of the corresponding scatter-modulation one are also included in Fig. 10, where the differences are dominated by the oscillations that are related to the residual aliasing error in scatter estimation, same observations as on the Catphan phantom study above. Quantitatively, the virtual scatter-modulation algorithm still achieves better CT number uniformity than the original primary-modulation one. As summarized in Table. 4, for the ROIs selected on the transaxial and reformatted images, the CT number non-uniformity measured using Eq. (13) is reduced from 30.4 and 21.4 HU, with the primary-modulation algorithm, to 13.3 and 5.2 HU, with the scatter-modulation algorithm, respectively.

Fig.10

Volumetric image reconstructions of the anthropomorphic thorax phantom from cone-beam scans without and with scatter correction. Reconstruction size: 512×512×256 voxels with 0.4×0.4×0.4 mm³ for a voxel. (a) Fan-beam transaxial image (reference); (b) cone-beam transaxial image without scatter correction; (c) cone-beam reformatted (sagittal) image without scatter correction; (d) cone-beam transaxial images with scatter corrected using the primary-modulation algorithm; (e) cone-beam transaxial images with scatter corrected using the scatter-modulation algorithm; (f) difference image by subtracting (d) from (e); (g) cone-beam reformatted images with scatter corrected using the primary-modulation algorithm; (h) cone-beam reformatted images with scatter corrected using the scatter-modulation algorithm; (i) difference image by subtracting (g) from (h). ROIs in (f) and (i) are in two uniform regions selected for a quantitative analysis in Table. 4. Display windows: [–450 550] HU for all transaxial and reformatted images; [–50 50] HU for the difference images. Reconstruction kernel: the Shepp-Logan.

Table 4

Averaged CT number (HU) in the selected ROIs of the uniform regions as shown in Fig. 10

Method		Fan–Beam	Cone–Beam	Primary–Mod.	Scatter–Mod.
Transaxial	ROI 1	–27.4	–395.3	–51.9	–47.0
	ROI 2	–32.3	–393.1	–30.9	–33.7
	ROI 3	–23.9	–377.6	–21.4	–35.3
	Δ _MAX	8.4	17.7	30.4	13.3
Reformatted	ROI 1	/	–413.7	–47.2	–65.2
	ROI 2	/	–403.0	–68.7	–62.8
	ROI 3	/	–393.3	–55.5	–68.1
	Δ _MAX	/	20.5	21.4	5.2

As demonstrated above, overall both the primary-modulation algorithm and the scatter-modulation one achieve good scatter correction performance, eliminating the cupping artifacts and improving image contrast. The two algorithms have similar levels of improvement in terms of CT number accuracy, CNR and noise for the contrast rods evaluated in Table. 2. However, the scatter-modulation algorithm has less residual artifacts in the estimated scatter and in the reconstructed image as well, which validate our previous argument that the scatter-modulation algorithm is able to reduce the aliasing errors otherwise aggravated by the nonuniformity of the modulation function when the primary-modulation algorithm is used.

5 Conclusion and discussion

In order to obtain good scatter correction performance of the primary modulation method, the primary modulator and the related scatter estimation algorithm needs to be carefully designed and implemented. When designing a primary modulator, the size of the attenuating blocker in the modulator, the transmission factor of the blocker, and the material of the blocker deserve more considerations, because they determine the modulation frequency, the modulation strength and the magnitude of beam hardening of the modulator, respectively. Some relevant investigations can be found in Refs. [15, 16]. Once the modulator is made, these hardware-related parameters usually cannot be easily changed. The scatter estimation algorithm is software-related and can be optimized in data processing. The original algorithm is a natural derivation from the principle and hypothesis of primary modulation. It is straightforward and computationally efficient as the filtering technique is used. However, it may be suboptimal in practical implementations since the principle and hypothesis may not be fully satisfied due to some nonideal issues (e.g., the inaccuracy in manufacture of the modulator, the finite focal spot, and the polychromatic X-ray source). Therefore, it is also important to investigate and improve the scatter estimation algorithm.

In this work, we present a new algorithm to estimate scatter for the primary modulator method, with a concept of virtual scatter modulation, which can significantly reduce the aliasing errors in estimated scatter and therefore improve image quality in reconstructions. In the new algorithm of scatter estimation, the measured X-ray image (i.e., the primary-modulated image) is first divided by the measured modulation function, equivalent to virtually modulating scatter; a scatter distribution is then estimated by demodulating scatter in the scatter-modulated image. Compared with the original algorithm in which scatter is estimated by subtracting the demodulated primary in the primary-modulated image from the unmodulated primary and scatter, the new algorithm possesses better scatter correction performance when a nonuniform modulation function is actually applied in practical implementations.

For the Catphan©600 phantom, the overall scatter estimation error is about 3.6% with the original algorithm. It is further reduced to about 2.3% when the new algorithm is used. Both the two algorithms reduce the CT number error from 312.7 HU to less than 20 HU, and enhance CNR from 12.3 to nearly 19. However, the oscillations that are clearly observed in the central region of the scatter profile estimated using the original algorithm, are mostly eliminated with the new algorithm, leading to less residual artifacts and CT number nonuniformity in the reconstructed image. Compared with 38.9 HU of CT nonuniformity in a uniform region of the Catphan phantom when the primary-modulation algorithm is used, the proposed scatter-modulation algorithm significantly reduces it to 4.5 HU, at the same level of our ground truth reference (3.3 HU). Relatively better low contrast detectability is also achieved in the subslice region of module CTP515, thanks to the less residual artifacts in the center. On the anthropomorphic thorax phantom, overall better CT number uniformity is achieved as well, with the maximum difference among the selected ROIs of uniform regions reduced roughly by a factor of two.

It is worth noting that the virtual scatter modulation shares a common step with a typical beam hardening correction, i.e., the division process using the measured modulation function in Eq. (5). Thus, the virtual scatter modulation can play an important role in compensating for the beam hardening of the modulator and in further reducing the associated scatter estimation bias. A comprehensive spectral compensation with the scatter modulation concept is under development. Without consideration of computational cost, it is also worth investigating iterative algorithms for scatter estimation in the future, where optimization criteria may be utilized, and the nonuniformity of the modulation function may be compensated in each iteration.

Footnotes

Acknowledgments

Research reported in this publication was partially supported by the National Institute of Biomedical Imaging and Bioengineering of the National Institutes of Health under Award Numbers R21EB008186, R21EB019597 and R21EB021545. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health. The project was also partially supported by the Ministry of Science and Technology of China Key Research and Development Project (Grant No. 2016YFC0101400) and by the National Natural Science Foundation of China (Grant No. 81671681). The authors would like to thank N. Robert Bennett for the help in making the 1-D erbium modulator.

References

Kak

A.C.

, Slaney

, Principles of Computerized Tomographic Imaging. New York: IEEE Press, 1988.

Hsieh

, Computed Tomography: Principles, Design, Artifacts,and Recent Advances. Bellingham, WA: SPIE, 2003.

Seibert

and Boone

, X-ray scatter removal bydeconvolution, Med Phys 15(4) (1988), 567–575.

Floyd

C.E.

, Baydush

A.H.

, Lo

J.Y.

, Bowsher

J.E.

and RavinScatter

C.E.

, compensation for digital chest radiography using maximumlikelihood expectation maximization, Invest Radiol 28(5) (1993), 427–433.

Wiegert

, Bertram

, Rose

and Aach

, Model-based scattercorrection for cone-beam computed tomography, Proc of SPIE 5745 (2005), 271–282.

Ning

, Tang

and Conover

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

Boone

, Seibert

, Tang

and Lane

, Grid and slot scanscatter reduction in mammography: Comparison by using Monte Carlotechniques, Radiology 222(2) (2002), 519–527.

Boone

, Scatter correction algorithm for digitally acquiredradiographs: Theory and results, Med Phys 13(3) (1986), 319–328.

Siewerdsen

J.H.

, Daly

M.J.

, Bakhtiar

, Moseley

D.J.

, Richard

, Keller

and Jaffray

D.A.

, A simple, direct method for X-rayscatter estimation and correction in digital radiography andcone-beam CT, Med Phys 33(1) (2006), 187–197.

10.

Kyriakou

, Riedel

and Kalender

W.A.

, Combining deterministicand Monte Carlo calculations for fast estimation of scatterintensities in CT, Phys Med Biol 51 (2006), 4567–4586.

11.

Zhu

and Fahrig

, A new scatter correction method for X-rayand volumetric CT using primary modulation, Proc of RSNA 2005.

12.

Zhu

, Bennett

N.R.

and Fahrig

, Scatter correction method forX-ray CT using primary modulation: Theory and preliminary results, IEEE Trans Med Imag 25(12) (2006), 1573–1587.

13.

Zhu

, Starman

, Bennett

N.R.

and Fahrig

, MTF measurementand a phantom study for scatter correction in CBCT using primarymodulation, IEEE Nuclear Science Symposium Conference Record 3 (2006), 1964–1968.

14.

Gao

, Fahrig

, Bennett

N.R.

, Sun

, Star-Lack

and ZhuScatter

, correction method for x-ray CT using primary modulation:Phantom studies, Med Phys 37(2) (2010), 934–946.

15.

Gao

, Zhu

and Fahrig

, Optimization of system parametersfor modulator design in x-ray scatter correction using primarymodulation, Proc of SPIE 7622 (2010), 76222A.

16.

——, Modulator design for x-ray scatter correctionusing primary modulation: Material selection, Med Phys 37(8) (2010), 4029–4037.

17.

Grimmer

, Fahrig

, Hinshaw

, Gao

and Kachelrie

, Empirical cupping correction for ct scanners with primarymodulation (eccp), Med Phys 39(2) (2012), 825–831.

18.

Zhu

, Local filtration based scatter correction for cone-beam ctusing primary modulation, Med Phys 43(11) (2016), 6199–6209.

19.

Chen

, Song

, Ma

and Zhao

, Optimization-based scatterestimation using primary modulation for computed tomography, Med Phys 43(8) (2016), 4753–4767.

20.

Bier

, Müller

, Berger

, Choi

J.-H.

, Fahrig

and Maier

, Scatter correction for C-arm CT using primary modulation, In Proc 4th Intl Mtg on Image Formation in X-ray CT (2016), pp. 383–6.

21.

Bier

, Berger

, Maier

, Kachelrie

, Ritschl

, Mller

, Choi

J.-H.

and Fahrig

, Scatter correction using a primarymodulator on a clinical angiography c-arm ct system, Med Phys p. In press, 2017.

22.

Gao

, Niu

, Zhu

and Fahrig

, Evaluation of an erbiummodulator in x-ray scatter correction using primary modulation, Proc of SPIE 7961 (2011), 796122.

23.

Zhu

, Gao

, Bennett

N.R.

, Xing

and Fahrig

, Scattercorrection for X-ray conebeam CT using one-dimensional primarymodulation, Proc of SPIE 7258 (2009), 725820.

24.

Feldkamp

L.A.

, Davis

L.C.

and Kress

J.W.

, Practical cone-beamalgorithm, J Opt Soc Am A 1 (1984), 612–619.

25.

Zhu

, Wang

and Xing

, Noise suppression in scattercorrection for cone-beam CT, Med Phys 36(3) (2009), 741–752.

26.

Star-Lack

, Starman

, Munro

, Jeung

, Richters

, Mostafavi

and Pavkovich

, A fast variable-intensity ringsuppression algorithm, Med Phys 33(6) (2006), 1997.

X-ray source	0.6 mm nominal focal spot size
	120 kVp, 35 mA, 13 ms pulse width
Detector size	397.3×298.0 mm²
	1024×768 pixels
Anti-scatter grid	none
Rotation	circular, 360 deg
Phantom scanned	Catphan	Thorax
Source-to-detector distance (mm)	1423	1202
Source-to-axis distance (mm)	860	786
Source-to-modulator distance (mm)	231	231
Number of views	800	625