1. Introduction

Cone-beam computed tomography (CBCT) plays an essential role in medical diagnosis with its great advantages such as non-destructive, high-resolution and non-overlapping images [1]. However, photons passing through the phantom collide, and both scattered and non-scattered signals (which are useful) are recorded by the detector. So the scattered signals affect the useful signals, especially when using high-energy X-ray and flat-panel detectors. Photon scattering reduces the contrast resolution of the image and produces streak artifacts and makes the cup-shaped artifacts more serious [2]; this affects the accuracy of CBCT image. Therefore, the elimination of photon scattering is fundamental.

Numerous non-Monte Carlo (MC) methods and MC methods have been proposed to estimate and correct the scattered signal for improving the image quality. The non-MC methods include deterministic approach, machine learning, etc. [3,4]. Acuros CTS is one of the well-known algorithms for deterministic approach, which was applied to the clinic almost immediately for treatment planning for cancer radiotherapy [5–7]. MC method was first proposed in the 1940s for solving numerical problems concerning random neutron diffusion in fissible material in atomic bomb designs [8]. Over these years, it plays an important role in complex physical systems and computational biology because of its extreme flexibility and power. MC method is commonly considered as the most accurate method for solving problems of radiation transport because it can accurately describe the underlying physical interactions between the radiation and matter. However, MC method has a very low convergence rate $O(N^{-1/2})$, where $N$ is the number of simulations. When greater accuracy is required, the number of simulations needed increases. Many efforts have been devoted to improve the efficiency and convergence rate of MC method. A large number of MC tools have been developed, such as PENELOPE [9], fixed forced detection (FFD) [10,11], PENELOPE subroutine: MC-graphics processing unit (GPU) [12], GPU-based MC tool (gMCDRR) [13,14], GPU-based Metropolis MC (gMMC) [15] and the method by smoothing the scattering functions to optimize the calculation [16]. FFD only considers the path to the detector to improve efficiency. PENELOPE subroutine: MC-GPU uses GPU instead of the central processing unit as the main computing device. However, MC-GPU spends lots of time on sampling and calculating the Compton scattering direction on the PENELOPE model. Besides, some realistic features, such as the detector response, are missed in the MC-GPU. The gMCDRR method establishes the discrete relationship between the inverse function of the cumulative density function of the cosine of the scattering angle and the energy and its cumulative density function, and then samples the scattering angle based on this relationship, which is an improvement of MC-GPU. However, a large number of particles do not reach the detector, causing them to get wasted, deteriorating overall computational efficiency. gMMC uses the Metropolis-Hasting algorithm to sample the entire photon path from the X-ray source to the detector. The steps for importance sampling according to the Metropolis-Hasting algorithm are as follows: (i) Given the current state $x_{i-1}$, sample $y$ from the proposal distribution $T(x_{i-1},y)$; (ii) Generate random number $u$, and update $x_{i}=y$, if $u<\textrm {min}\Big \{1, \frac {p(y)T(y,x_{i-1})}{p(x_{i-1})T(x_{i-1},y)}\Big \}$; otherwise, $x_{i}=x_{i-1}$, where $p(x_{i-1})$ is the probability density function (PDF) of the particle path. However, when gMMC uses an independent sampling strategy, it will encounter a path $l_{x_{i-1}}$ with small $T(y,x_{i-1})$ that is close to $0$ and a large value of PDF $p(x_{i-1})$, $\frac {p(y)T(y,x_{i-1})}{p(x_{i-1})T(x_{i-1},y)}$ is close to 0, thus the path $l_{x_{i-1}}$ may be selected a hundred times or even more. In this situation, the path such as $l_{x_{i-1}}$ will cause large noise in the scattered image and slow down the convergence rate. Such event $l_{x_{i-1}}$ usually occurs when the intersection length between path $l_{x_{i-1}}$ and phantom is short and with low energy; for example, path $l_{x_{i-1}}$ is a photon scatter in the surface of phantom and reaches to the detector area outside the phantom projection, while most paths $l_y$ go through phantom with low probability, so the $p(x_{i-1})$ would be far greater than $p(y)$. Besides, the proportion of photon scattering will increase with increasing energy, the gMMC needs to calculate the scatter ratio in each segment energy by other MC tools when simulating with a polyenergetic spectrum, requiring additional simulation time. Finally, the efficiency of the gMMC method is affected by the proposal function. The closer the proposal function is to the target distribution, the more efficient the Metropolis-Hasting algorithm is. However, it is difficult to choose the proposal function close to the PDF of the particle path.

Compared with MC, quasi-Monte Carlo (QMC) uses low-discrepancy sequences [17] instead of random sequences for simulation. In this paper, Sobol’ sequence [18,19] is used for simulation. The convergence rate of MC is $(O(N^{-1/2}))$, which is not affected by the problem dimension, where $N$ is the number of simulations. However, under appropriate conditions, the convergence rate of QMC is close to $(O(N^{-1+\epsilon }))$ [17], which depends on the problem dimension, where $\epsilon >0$. The convergence order of QMC method is asymptotically much faster than that of MC for fixed dimension.

In this paper, we transform the problem of calculating the X-ray scattering signal in CT into a high-dimensional integration problem and give the integral formula for the probability of photons emitted from the X-ray source through the phantom to reach the detector. In addition, we propose a new algorithm called gQMCFRD, which combine GPU-based gQMC method with forced random detection technique (FRD) for simulation. Specifically, according to the integral formula of the photon scattering probability along each path, we first calculate the probability of the photon passing through the phantom without scattering and determine the position of the first scattering point $A_1$ by sampling. Second, sample the scattering type and the direction of photon after scattering by using the weighted statistics and the Rational Inverse Transform with Aliasing (RITA) algorithm [20], then we sample the position of the next-order scattering point $A_{i},i=2,\dots ,n$ and calculate the weight of the photon moving from the current interaction point to the next-order interaction point $A_i$. Finally, force a random selection of a patch $D_i$ of size $k \times k$, corresponding to the interaction point $A_i$, and then accumulating the fluence from $A_i$ to each pixel of $D_i$.

The rest of this article is organized as follows. Section 2 introduces the basic idea and the overall structure of the proposed method. Key gQMCFRD technologies and implementation are described in detail. Section 3 contains results of the scatter calculation in both Shepp-Logan phantom and Head phantom. Finally, the article is summarized and discussed in Section 4.

2. Methods and implementation

2.1 Quasi-Monte Carlo simulation

Many problems in particle transport can be formulated as high-dimensional integrals, and only in rare cases do explicit solutions exist, simulation is often the only method of solution. Because of the extreme flexibility and power, MC or QMC methods are used to approximate high-dimensional integrals. Consider an integral $I_d=\int _{\Omega }g(\boldsymbol X)\mathrm{d}\boldsymbol X.$ Through suitable transformations, such as inverse transform [17], the integral can be transformed into an integration problem over $[0,1]^{d}$,

The MC or QMC estimate of $I_d(f)$ is of the form

For low discrepancy point set, $D^*_{N}$ is in the order $O(N^{-1}[log(N)]^{d-1})$ [21], where

2.2 Overall structure of gQMCFRD

This paper considers the geometric model of the CT imaging system shown in Fig. 1(a). It is composed of X-ray source ${\rm S}$, a phantom ${\rm M}$ and a detector ${\rm D}$. In Fig. 1(a), the photon penetrating phantom ${\rm M}$ undergoes $n$-order scattering before hitting the detector ${\rm D}$, and $A_1, A_2, \dots ,A_n$ represent the scattering points, respectively. After the photon reaches each scattering point $A_i,i=1,2,\dots ,n$, its direction changes and it passes through the phantom ${\rm M}$, reaching the detector ${\rm D}$ with a certain probability. In order to improve photon utilization efficiency, the $(6i-1)$-th and $6i$-th components of the Sobol’ sequence are used to force the random selection of the patch $D_i$, which is the detector pixel $D_ {i_0}$ in the upper left corner and has a size of $k\times k$, corresponding to the interaction point $A_i$. Other detector pixels in patch $D_i$ are recorded as $D_{i_j},j=1,\dots ,k^2-1$. Denote the path of the photon to the detector pixel $D_{i_j},j=0,1,\dots ,k^2-1$ as $l_{D_{i_{j}}}:{\rm S}\rightarrow A_1 \rightarrow \dots \rightarrow A_i \rightarrow B_{i_{j}} \rightarrow D_{i_{j}}$, where $A_0$ is the intersection of the initial incident X-ray and the phantom ${\rm M}$, $B_{i_{j}}$ is the intersection on boundary of ${\rm M}$ along with $\overrightarrow {A_{i}D_{i_{j}}}$ and $C_i,i=0,\dots ,n$, is the intersection on the boundary of the phantom ${\rm M}$ along with $\overrightarrow {A_{i}A_{i+1}}$. When $k$ increases, the calculation will increase. An illustration of the photon paths $l_{D_{i_j}}$ is shown in Fig. 1(b).

 

Fig. 1. (a) Geometry illustration of gQMCFRD simulation; (b) Geometry illustration of photon paths.

Download Full Size | PDF

The main structure of gQMCFRD is shown in Algorithm 1. Before the simulation begins, initialize the relevant parameters with necessary data such as X-ray spectrum and physics data. All of these data are transferred to GPU memory during this step. Specifically, the voxelized phantom data and attenuation coefficients are stored in the texture memory of GPU to achieve quick access because of cache-function of GPU. The differential cross section (DCS) data are stored in GPU global memory and the rest of the variables are stored in GPU constant memory. All image tally counters are set to zero. Next, segment the energy spectrum, simulate the scattering results of each energy spectrum in batches, and then accumulate the scattering results of all batches. In the simulation, Rayleigh scattering, Compton scattering and photoelectric effect are considered [27]. However, because the photon is absorbed and does not reach the detector when the photoelectric effect occurs, the simulation results in this paper only record the probability and energy of photon reaching the detector after Rayleigh scattering and Compton scattering. Besides, if Rayleigh scatter occurs, the energy does not change; if Compton scatter occurs, the energy changes and calculate the energy after each collision according to formula (9).

In the algorithm, we use $6n$-dimensional Sobol’ sequence $\{u_1,u_2,\dots ,u_{6n}\}$ in one simulation, where $n$ is the scattering order. When simulating the trajectory of photons in the phantom after n collisions, the constructive dimension of the gQMCFRD algorithm is $6n$, where $n$ is a positive integer. Specifically, the steps of simulating the trajectory of the photon with single scatter are as follows: use the first component $u_1$ of the pregenerated Sobol’ sequence to sample the initial energy of the incident light from the predetermined energy spectrum by the inverse transform method [17], use $u_2, u_3$ to sample the incident direction of the photon, use $u_4$ to sample the first interaction point and $u_5, u_6$ to force random sampling of the detector pixels corresponding to the first-order interaction point. At this time, the problem dimension is $6$. If the photon undergoes second-order scattering, the simulation of the trajectory of the photon from the X-ray source to the first-order scattering point is the same as when the photon only collides once. Then use $u_7$ to sample the scatter type, $u_8, u_9$ to sample the scattering angle direction of the photon at the first-order scattering point, $u_{10}$ to sample the position of the second-order scattering point, and use $u_{11}, u_{12}$ to force the random selection of the detector pixels corresponding to the second-order interaction point. At this time, the problem dimension is $12$. In the same way, it can be inferred that the problem dimension of the $n$-th order scattering path of a photon is $6n$.

2.3 Derivation of photon path probability

Consider the path $l_{D_{i_{j}}}, i\in \{i,\dots ,n\}, j=0,1,\dots , k^2-1$ of a photon with an initial energy of $E_0$ through a three-dimensional phantom ${\rm M}$ (not necessarily uniform), as shown in Fig. 1(b). The calculation of the probabilities of each path is divided into the following steps. (i) Calculate the probability of the initial energy spectrum and the initial incident direction of the photon and calculate the non-attenuated probability. (ii) Sample the first interaction point $A_1$ in the initial incident direction and calculate the probability $A_0\rightarrow A_1$. (iii) Sample the scattering interaction type and the scattering angle distribution of the photon at the interaction point $A_{i-1},i=2,\dots ,n$, sample the position of the next-order scattering point $A_i$ in the scattering direction of $A_{i-1}$ and calculate the probability $A_{i-1}\rightarrow A_{i}$. (iv) Force the random selection of the patch $D_i$ of size $k \times k$ corresponding to the interaction point $A_i$. (v) Calculate the probability $A_{i}\rightarrow D_{i_j}$.

Let the unit initial incident direction be $ {\mathop{\omega}\limits ^{\rightarrow}}_0$, $C_{0}=A_{0}+c_0 {\mathop{\omega}\limits ^{\rightarrow}}_0$. After the photon enters the phantom ${\rm M}$, the path length of the photon from its current position $A_0$ to the intersection point $C_0$ of the photon travelling along ${\mathop{\omega}\limits ^{\rightarrow}}_0$ on the boundary of the phantom ${\rm M}$ follows the distribution $\mu _{tot}(x,E_0)\textrm {exp}{\big [}{-\int _{0}^{t_1}\mu _{tot}(A_0+l {\mathop{\omega}\limits ^{\rightarrow}}_0,E_0){\textrm {d}}l}{\big ]},x=A_0+t_1{\mathop{\omega}\limits ^{\rightarrow}}_0,0\leq t_1\leq c_0$, where $\mu _{tot}(x,E)$ is the linear attenuation coefficient which is a function of both the spatial coordinates $x$ and the energy $E$ (namely the macroscopic total interaction cross section) and is the sum of Compton scattering coefficient $\mu _{T_0}(x,E)$, Rayleigh scattering coefficient $\mu _{T_1}(x,E)$ and photoelectric coefficient $\mu _{a}(x,E)$. So the probability of the photon escaping the phantom ${\rm M}$ along the unit initial incident direction $ {\mathop{\omega}\limits ^{\rightarrow}}_0$ (namely the non-attenuated probability) is $p_0({\mathop{\omega}\limits ^{\rightarrow}}_{0},E_0)=\textrm {exp}{\big [}{-\int _{0}^{c_0}\mu _{tot}(A_0+l{\mathop{\omega}\limits ^{\rightarrow}}_0,E_0){\textrm dl}}{\big ]}$ and the probability of the photon from the X-ray source ${\rm S}$ to the first interaction point $A_1$ along the unit initial incident direction ${\mathop{\omega}\limits ^{\rightarrow}}_0$ is

After the first collision, we need to sample the scatter type, calculate the remaining energy $E_1^{\delta _1} ,\delta _1=0,1$ and sample the unit scatter direction $ {\mathop{\omega}\limits ^{\rightarrow}}_1$. The unit scatter direction is controlled by the photon interaction type and corresponding DCS at $A_1$. The PDF of the scattering angle is a normalized DCS function. Let $T_0, T_1$ represent the scatter type of Compton and Rayleigh respectively, $p_{T_0},p_{T_1}$ represent the corresponding probability, and $p_{\theta _0},p_{\theta _1}$ are the corresponding PDF. So

Similarly, the probability of the photon ${\rm S}\rightarrow A_1 \rightarrow A_2 \dots \rightarrow A_i$ is

Let the unit direction of the photon from the interaction point $A_i,i\in \{1,\dots ,n\}$ to the detector pixel $D_{i_j}$ $ {\mathop{r}\limits ^{\rightarrow}}_{i_j}$, which $D_{i_j}$ is a detector pixel of the patch $D_i$ corresponding to $A_i$. The probability of the photon from $A_i,i\in \{1,\dots ,n\}$ to the detector pixel $D_{i_j},j=0,\dots ,k^2-1$ ($A_{i}\rightarrow B_{i_j} \rightarrow D_{i_j}$) is

So the whole probability that the photon follows the path $l_{D_{i_j}}: {\rm S} \rightarrow A_1 \rightarrow \dots \rightarrow A_i \rightarrow B_{i_j} \rightarrow D_{i_j}, i\in \{1,\dots , n\},j=0,\dots ,k^2-1$, is

The $A_0,A_1,\dots ,A_i,i=1,\dots ,n$, in Eq. (6) are random variables, which are related to the initial energy and the initial incident direction, $A_i$ is related with $A_{i-1}$ and force the random selection of the patch $D_i$. The specific implementation method of sampling the initial energy spectrum and unit initial incident direction ${\mathop{\omega}\limits ^{\rightarrow}}_0$ of the X-source ${\rm S}$, selecting $A_i,i=1,\dots ,n$, sampling the scattering type $T_0, T_1$ and scattering direction $ {\mathop{\omega}\limits ^{\rightarrow}}_{i},i=1,\dots ,n$, calculating the remaining energy $E_i^{\delta },i=1,\dots ,n,\delta =0,1$ of the photon after $i$- collision, and forcing random selection of the patch $D_i$ corresponding to $A_i$, see Section 2.4 in this article.

We transform the abstract particle transportation path problem into a high-dimensional integration problem. However, it is often difficult to calculate high-dimensional integrals directly. To calculate such integrals, we propose gQMCFRD algorithm. When calculating the probability of path $l_{i_j}$, the algorithm can obtain the probability of path $l_{m_j},m=1,\dots ,i$ at the same time. That is, by calculating the path probability of the $i$-th order scattering of the photon, the information of the $m(\leq i)$-th order scattering of the photon can be clearly known. This helps to clearly know the contribution of each order of scattering to the total scattering and to clarify the structure of the total scattering. For comparison, the traditional Monte Carlo method calculates a specific scattering event, and cannot obtain a scattering information lower than the scattering order at the same time.

2.4 gQMCFRD algorithm and implementation

In this section, we use the gQMCFRD algorithm to calculate the probability in (6). The initial photon energies, directions, interaction position, interaction types, and location of the detector are involved in the simulation. The detector pitch $D_i$ is randomly selected at each simulation. Therefore, to calculate the weight of the photon paths $l_{D_{i_j}}$, we need to implement Algorithm 1.

2.4.1 Initialize photon energy and direction

We consider polyenergetic source in the simulation, $u_1$ is used to sample the photon energy from the given energy spectrum by Walker’s aliasing method [28] (see Fig. 2(a)). The $\Omega _0$ of the direct light to the detector is usually rectangular or circular. The flow intensity in $\Omega _0$ may be complicated, it is difficult to directly sample particles that obey the flow intensity distribution. We will instead sample uniformly in $\Omega _0$ by using $u_2$ and $u_3$, then the photon incident direction ${\mathop {\omega } ^{\rightarrow }}_0$ and $A_0$ can be calculated. Afterwards, we calculate $f_0(A_{0})$ as the initialization weight.

 

Fig. 2. (a) Illustration of photon energy distribution; (b) Angular deflections in single-scattering events.

Download Full Size | PDF

2.4.2 Interaction type and the weight of the corresponding type

Let $p_{T_0}(A_{i-1}), p_{T_1}(A_{i-1}), p_{a}(A_{i-1}), i\in \{2,\dots ,n+1\}$ denote respectively the proportions of Compton scattering, Rayleigh scattering and photoelectric effect in the $(i-1)$-th interaction point $A_{i-1}$, $p_{T_0}(A_{i-1})+p_{T_1}(A_{i-1})+p_{a}(A_{i-1})=1$. In the simulation, three interactions of photons are considered. But when the photoelectric effect occurs, photons are absorbed and does not reach the detector. Therefore, only record the probability and energy of photon reaching the detector after Rayleigh scattering and Compton scattering in the simulation. We assign Rayleigh scatter to an interaction point with a sample probability $\gamma$ and use the $(6(i-1)+1)$-th component of $6n$-dimensional Sobol’ sequence, if $u_{6(i-1)+1}<\gamma$, it is considered that Rayleigh scattering occurs; otherwise, Compton scattering occurs. $\gamma =0.5$ is used in numerical experiments. Because importance sampling is used, the proportions of Compton scattering $p'_{T_0}(A_{i-1})=p_{T_0}(A_{i-1})/(1-\gamma )$, the proportions of Rayleigh scattering $p'_{T_1}(A_{i-1})=p_{T_1}(A_{i-1})/\gamma$.

2.4.3 Distribution of scattering angle

In this work, PENELOPE physical database used in GEANT4 (GEometry ANd Tracking) [29] is employed. After $(i-1)$-th interaction, the polar scattering angle PDF of Rayleigh scattering is given by the following formula [30],

In addition, the Klein-Nishina DCS is used to characterize the polar scattering angle PDF of Compton scattering [31],

In case of Rayleigh scatter photon energy remains the same $E_i=E_{i-1}$, in case of Compton scatter photon energy changes after each interaction,

After $(i-1)$-th interaction, the polar scattering angle $\theta _{i-1}$, in fact $\textrm {cos}(\theta _{i-1})$, is sampled by RITA algorithm [20] for both Rayleigh and Compton scatter which uses the $(6(i-1)+2)$-th component of $6n$-dimensional Sobol’ sequence. The azimuthal scattering angle $\phi _{i-1}$ is generated according to the uniform distribution in $(0,2\pi )$, i.e. $\phi _{i-1}=2\pi u_{6(i-1)+3}$, as show in Fig. 2(b). We will give the experimental results of using the RITA algorithm to sample the scattering polar angle in the subsequent Section 3.3.

2.4.4 Interaction position

Since the photon from $A_{i-1} \rightarrow A_{i}, i=1,\dots ,n$, obeys the distribution:

2.4.5 Location on detector

We force the random selection of the detector pitch $D_{i},i=1,\dots ,n$, corresponding $A_{i}$ by $u_{6i-1}$ and $u_{6i}$. Then, from Eq. (5) the scatter angle probability density function $p_{\theta _{\delta _i}}(A_i,E_{i}^{\delta _{i}},{\mathop{r}\limits ^{\rightarrow}}_{i_j})$ and the solid angle $\Omega _{A_{i},D_{i_j} }$, and through probability $\hat {p}_i({{\mathop{r}\limits ^{\rightarrow}}_{i_j}},E_{i}^{\delta _{i}})=\textrm {exp}{\big [}{-\int _{0}^{b_{i_j}}\mu _{tot}(A_i+l{\mathop{r}\limits ^{\rightarrow}}_{i_j},E_i^{\delta _i})\mathrm{d}l}{\big ]}$ are calculated to get the corresponded weight $P({l_{D_{i_j}}})$ which can be expressed as a function $f_{i}(u_1,\dots ,u_4,\dots ,u_{6(i-1)+1},\dots ,u_{6(i-1)+4})g_{i_j}(u_{6i-1},u_{6i}), i\in \{1,\dots ,n\},j=0,\dots ,k^2-1$. The required probability $P(l_{D_{i_j}})$ is the mathematical expectation $P(l_{D_{i_j}}) =E(f_{i}g_{i_j})$.

The QMC method uses the average of the function values to estimate the integral $E(f_ig_{i_j})$:

Let $\boldsymbol {u^{m}}=\{u^m_1,\dots ,u^m_{6i}\}$. $\boldsymbol {u^{1}}, \dots , \boldsymbol {u^{N}}$ are a Sobol’ sequence on $\boldsymbol {[0,1]}^d$, $d=6i, i\in \{1,\dots ,n\}$, $n$ is the highest scattering order.

Let $P(l^m_{D_{i_j}}) =f_{i}(u^m_1,\dots ,u^m_4,\dots ,u^m_{6(i-1)+1},\dots ,u^m_{6(i-1)+4})g_{i_j}(u^m_{6i-1},u^m_{6i})$ and

3. Numerical results

We perform scattering calculations on the Shepp-Logan phantom [32] and Head phantom to demonstrate the effectiveness of the proposed gQMCFRD algorithm. The setting of in-field of view (IFOV) is simulated and the X-ray is point source with energy spectrums: $5-120$keV. For the IFOV case, the whole phantom is under the X-ray illumination and the X-ray source S-to-phantoms distance is $500$ mm and phantoms-to-detector distance is $500$ mm. The detector resolution is $512 \times 512$ pixels with a pixel size of $0.8 \times 0.8 {\textrm {mm}}^2$. With regard to computational hardware, we use a GeForce RTX 2080 Ti graphics card that is equipped with 4352 processors and 11.0 GB global memory. In the numerical experiment of this article, the detector patch size forcibly selected in the gQMCFRD algorithm is $4\times 4$. We calculated the efficiency of the gQMCFRD algorithm with the patch size $1\times 1$, $2\times 2$, $4\times 4$ and $8\times 8$, respectively, and found that the efficiency of the gQMCFRD algorithm is more suitable with the patch size is $4\times 4$. As the patch size increases, the computational complexity of gQMCFRD will increase.

gQMCFRD, gMCDRR, gMMC and gMCFRD are benchmarked against MC-GPU on the same GPU card. MC-GPU code is publicly available via $http://code.google.com/p/mcgpu/$ and is described in detail in [12]. The gMCDRR and gMMC are implemented according to [13,15]. The difference between gMCFRD and gQMCFDR is to replace Sobol’ points with MC points. To quantify the speed performances of gQMCFRD, we use the so-called figure of merit (FOM) [33]

3.1 Shepp-Logan phantom

In this section, we perform the experiment to compare the practical performance of the proposed method in Shepp-Logan phantom, which is composed of bone, water and air. The phantom dimension is $320\times 400\times 360$ voxels with a voxel size of $0.5\times 0.5\times 0.5\textrm {mm}^3$. Figures 3(a), 3(a1), 3(a2), 3(a3) and 3(c) are the geometric illustration. gQMCFRD evaluate the path integrals in the exponential term in Eq. (6) by using Siddon’s ray-tracing algorithm [34].

 

Fig. 3. (a) Shepp-Logan phantom; (a1),(a2) and (a3) the transversal, coronal and sagittal of the Shepp-Logan phantom respectively; (c) Shepp-Logan phantom under under X-ray point source; (b) Head phantom; (b1),(b2) and (b3) the transversal, coronal and sagittal of the head phantom respectively; (d) Head phantom under under X-ray point source.

Download Full Size | PDF

The scatter signals computed by MC-GPU, gMCDRR, gMMC, gMCFRD and gQMCFDR are presented in Fig. 4. These images visually demonstrate good agreements among the gQMCFRD results, the MC-GPU results, the gMCDRR results, the gMMC results and the gMCFRD results, indicating the accuracy of gQMCFRD code. And it can be seen that FOM is $371$ when $2^{38}$ photons transported in MC-GPU, the relative differences of the total scatter signals are $1.21\%$, FOM is $9961$ and EIF is $26.85$ when $2^{29}$ photons transported in gQMCFRD, $1.13\%$, $404$ and $1.09$ when $2^{38}$ photons transported in gMCDRR, $2.48\%$, $45$ and $0.12$ when $2^{34}$ photons transported in gMMC, $1.68\%$, $1343$ and $3.62$ when $2^{29}$ photons transported in gMCFRD from the Table 1. Thus gQMCFRD is $27$, $25$, $221$ and $7$ times faster than MC-GPU, gMCDRR, gMMC and gMCFRD, respectively, for the Shepp-Logan phantom. Furthermore, it indicates that QMC method is more effective than MC method for Shepp-Logan phantom.

 

Fig. 4. The primary signals (namely non-scattered singals) and total scattering signals of the Shepp-Logan phantom. (a1)-(a4) are the primary signals, total scatter signals (sum of first- and second-order scatter), first-order Rayleigh scatter, first-order Compton scatter of MC-GPU; (b1)-(b4) are corresponding results of gMCDRR; (c2)-(c4) are corresponding results of gMMC; (d1)-(d4) are corresponding results of gMCFRD; (e1)-(e4) are corresponding results of gQMCFRD; (f2)-(f4) are the signal profiles of Shepp-Logan phantom through the center of the corresponding detector which parallel to the x-axis.

Download Full Size | PDF

Table 1. The total scattering signal results of gMCDRR, gMMC, gMCFRD and gQMCFRD are compared with the result of MC-GPU for Shepp-Logan phantom.^a

View Table | View all tables in this article

3.2 Head phantom

Head phantom is composed of tissue, bone and air. The phantom dimension is $270\times 400\times 430$ voxels with a voxel size of $0.5\times 0.5\times 0.5\textrm {mm}^3$. Figures 3(b), 3(b1), (b2), 3(b3) and 3(d) are the geometric illustration.

Computational results of Head phantom are presented in Fig. 5. And it can be seen that gQMCFRD is $37$, $32$, $253$ and $10$ times faster than MC-GPU, gMCDRR, gMMC and gMCFRD, respectively, and the relative differences are less than $2\%$ for the Head phantom from the Table 2.

 

Fig. 5. The primary signals (namely non-scattered singals) and total scattering signals of the Head phantom. (a1)-(a4) are the primary signals, total scatter signals (sum of first- and second-order scatter), first-order Rayleigh scatter, first-order Compton scatter of MC-GPU; (b1)-(b4) are corresponding results of gMCDRR; (c2)-(c4) are corresponding results of gMMC; (d1)-(d4) are corresponding results of gMCFRD; (e1)-(e4) are corresponding results of gQMCFRD; (f2)-(f4) are the signal profiles of Head phantom through the center of the corresponding detector which parallel to the x-axis.

Download Full Size | PDF

Table 2. The total scattering signal results of gMCDRR, gMMC, gMCFRD and gQMCFRD are compared with the result of MC-GPU for Head phantom.^a

View Table | View all tables in this article

The number of $i(i=0,1,2)$- scattered photons in Shepp-Logan phantom and Head phantom is shown in Table 3. It indicates that the proportion of higher order scatter will decrease.

Table 3. The number of $i$-order scattered photons in Shepp-Logan phantom and Head phantom when $2^{29}$ photons are uesd in one simulation.

View Table | View all tables in this article

3.3 RITA vs GITM

In this section, we compare the accuracy of scattering polar angles sampled by RITA [20] and a generalized inverse transform method (GITM) [13] in bone materials. In gMCDRR, using GITM which uses linear interpolation to sample the polar scattering angle $\theta _{i},i=1,\dots ,n$ [13]. In our implementation, the data used in RITA are precalculated by DCS. At each energy, $64$ suitable grid points are used to minimize the interpolation error between the interpolating PDF and the original DCS. The DCS of Rayleigh scatter can be directly calculated by the Eq. (7). While the DCS of Compton scatter cannot be directly calculated by the Eq. (8), we use Geant4 to generate Compton scattering polar angle to get a histogram at each energy, and then the DCS can be calculated.

To illustrate the accuracy of the sampling method, we introduce the interpolation error in the $s$-th grid $x_s$, which is defined as $\epsilon _s=p(x_s)-\hat {p}(x_s)$. For Rayleigh scattering polar angle, it can be seen that the mean error is $1.15\times 10^{-3}$ at 30keV when $4096$ grid points is used by GITM, and $3.18\times 10^{-4}$ when $64$ grid points is used by RITA from Table 4. This shows that RITA has higher accuracy than GITM. Besides, the interpolation mean error and maximum error increase with increasing energy and when single-precision floating-point numbers are used, GITM works badly, the interpolation error increases as the grid points increase. When sampling the Rayleigh scattering polar angle distribution by GITM in this section, double-precision floating-point numbers are used. Rayleigh DCS with the energy range of $5-150$keV for bone is shown in Fig. 6(a). Figures 6(b) and 6(c) demonstrate the interpolation error between the interpolating PDF and the original Rayleigh DCS which sampled from the Rayleigh DCS at 30 keV and 60 keV, respectively and the number of interpolation grid points is $2048$ by GITM. Those indicate the effectiveness of RITA.

 

Fig. 6. (a) Rayleigh DCS with the energy range of $5-150$keV for bone; (b) and (c) the interpolation errors of Rayleigh DCS at 30 keV and 60 keV, respectively.

Download Full Size | PDF

Table 4. Interpolation errors of sampling Rayleigh scattering polar angle distribution for bone^a

View Table | View all tables in this article

As for Compton scattering polar angle, it can be seen that the mean error is $1.08\times 10^{-4}$ at 60keV when $4096$ grid points is used by GITM, and $1.57\times 10^{-4}$ when $64$ grid points is used by RITA from Table 5. Compton DCS with the energy range of $5-150$keV for bone is shown in Fig. 7(a). Figures 7(b) and 7(c) demonstrate the interpolation errors between the interpolating PDF and the original Compton DCS which sampled from the Compton DCS at 30 keV and 60 keV, respectively and the number of interpolation grid points is $2048$ by GITM. Those indicate the effectiveness of RITA.

 

Fig. 7. (a) Compton DCS with the energy range of $5-150$keV for bone; (b) and (c) the interpolation error of Compton DCS at 30 keV and 60 keV, respectively.

Download Full Size | PDF

Table 5. Interpolation errors of sampling Compton polar scattering angle distribution for bone^a

View Table | View all tables in this article

4. Discussion and conclusion

In all the cases studied, gQMCFRD only needs to simulate $2^{29}$ photon paths to achieve good results in which the relative difference of the total scatter signals is less than $2\%$ and the statistical uncertainty is of the order of $10^{-3}$, whereas the gMCDRR, gMMC and gMCFRD methods need to simulate $2^{38}$, $2^{34}$ and $2^{29}$ photons, respectively, to achieve the relative difference is about $2\%$ and the statistical uncertainty is of the order of $10^{-2}$. The photon paths in gMCDRR and gMMC are $512$ and $32$ times, respectively, the number paths in gQMCFRD and the Efficiency Improvement Factors (EIF) range $27\sim 37$ (or $1.09\sim 1.16$, or $0.12\sim 0.15$) by gQMCFRD (or gMCDRR, or gMMC) with comparison to MC-GPU. These indicate that gQMCFRD has a faster convergence rate and can reach a high level of accuracy using less photons. However, the computations in gQMCFRD are more complex than gMCDRR. In particular, gMCDRR eliminates the need for complex voxel crossing computations in photons transports by using the Woodcock scheme [35]. In contrast, gQMCFRD needs to perform ray-tracing operations to go through phantoms and calculate the next scatter point by a specified integral value, as shown in Eqs. (6) and (11), which is the major burden. In the implementation, we have optimized the algorithm by reducing the number of calculations. Two parts of the integral need to be calculated in Eq. (11) but the integral path $\overrightarrow {A_iA_{i+1}}$ is part of the integral path $\overrightarrow {A_iC_i}$. The calculation of integral on path $\overrightarrow {A_iA_{i+1}}$ can be reduced by a look-up table that is the precalculated integrated result at a different position, which is usually equally divided distance on the path $\overrightarrow {A_iC_i}$. The main computation of gMMC is the ray-tracing operations on a photon path probability evaluation and the corresponded accepted probability and the efficiency is affected by the selected proposal function. The EIF range $3.62\sim 3.70$ by gMCFRD with comparison to MC-GPU for these phantoms. The implementation of gMCFRD is the same as gQMCFRD, the only difference is that random points are used instead of Sobol’ points. That indicates that gQMCFRD is more effective for these phantoms. Besides, gMCFRD is faster than MC-GPU, gMCDRR and gMMC for these all cases.

In this study, we transform the photon transport simulation into high-dimensional integration problem and propose a new algorithm called gQMCFRD for simulation which combines GPU-based QMC method with forced random detection technique. We demonstrate the effectiveness of this method in the context of computing scattered photon signals in CT. In all the cases studied, the EIF range $27\sim 37$ by gQMCFRD with comparison to MC-GPU, while the relative difference is less than $2\%$ and the statiscal uncertainty is the order of $10^{-3}$. The EIF range $1.09\sim 1.16$ (or $0.12\sim 0.15$, or $3.62\sim 3.70$) by gMCDRR (or gMMC, or gMCFRD) with comparison to MC-GPU, the relative difference is about $2\%$ and the statiscal uncertainty is the order of $10^{-2}$.

For efficiency considerations, this paper only calculates the probability of first- and second-order scattering photons reaching the detector. Higher scattering order means that we need to simulate the photon path in a higher-dimensional space, because simulating the $n$-th order scatter is equivalent in calculating an integral to $d=6n$ dimensional space. As a result, the convergence rate of gQMCFRD and the effective dimension of the problem, in theory, are needed for future study. Also, using this method in a real CT system for scatter correction will be our future work.