Matrixed-Spectrum Decomposition Accelerated Linear Boltzmann Transport Equation Solver for Fast Scatter Correction in Multi-Spectral CT

TL;DR

This paper introduces a novel matrixed-spectrum LBTE solver that efficiently computes multi-spectral X-ray scatter in CT, significantly reducing computational cost while maintaining high accuracy for scatter correction.

Contribution

The paper proposes a unified, spectrum-agnostic LBTE solver using matrixed-spectrum decomposition and a spectrum basis, enabling simultaneous scatter estimation across multiple spectra with minimal additional computation.

Findings

Reduces computational time for multi-spectral scatter estimation

Maintains high accuracy comparable to Monte Carlo methods

Effectively reduces scatter artifacts in reconstructed CT images

Abstract

X-ray scatter has been a serious concern in computed tomography (CT), leading to image artifacts and distortion of CT values. The linear Boltzmann transport equation (LBTE) is recognized as a fast and accurate approach for scatter estimation. However, for multi-spectral CT, it is cumbersome to compute multiple scattering components for different spectra separately when applying LBTE-based scatter correction. In this work, we propose a Matrixed-Spectrum Decomposition accelerated LBTE solver (MSD-LBTE) that can be used to compute X-ray scatter distributions from CT acquisitions at two or more different spectra simultaneously, in a unified framework with no sacrifice in accuracy and nearly no increase in computation in theory. First, a matrixed-spectrum solver of LBTE is obtained by introducing an additional label dimension to expand the phase space. Then, we propose a ``spectrum basis'' for LBTE and a principle of selection of basis using the QR decomposition, along with the above solver to construct the MSD-LBTE. Based on MSD-LBTE, a unified scatter correction method can be established for multi-spectral CT. We validate the effectiveness and accuracy of our method by comparing it with the Monte Carlo method, including the computational time. We also evaluate the scatter correction performance using two different phantoms for fast-kV switching based dual-energy CT, and using an elliptical phantom in a numerical simulation for kV-modulation enabled CT scans, validating that our proposed method can significantly reduce the computational cost at multiple spectra and effectively reduce scatter artifact in reconstructed CT images.

Tables

Table 1. Table 1 : Experiment Setup

Experiment	kV-modulation	kV-switching
phantom	water Phantom	Multi-Energy Phantom	Head Phantom
SID	750mm	750mm	755mm
SDD	1184mm	1184mm	1155mm
Pixel size	$0.45$mm	$0.388$mm	$0.388$mm
Detector	$960\times 960$	$1024\times 768$	$1024\times 768$
Voltage	70 to 140kVp	80/140kVp	80/120kVp
Filter	2mm Al	2mm Al	0.3mm Cu

Table 2. Table 2 : Accuracy of scatter distributions between Geant4 and MSD-LBTE for different spectra in Figs. 6 , 7

Spectrum	Cylinder		Head
	MRE(%)	SSIM	MRE(%)	SSIM
Spec1	2.39	0.9741	2.94	0.9658
Spec2	1.13	0.9847	1.03	0.9769
Spec3	1.81	0.9796	1.79	0.9723
Spec4	1.32	0.9825	1.34	0.9739
Spec5	1.99	0.9789	2.40	0.9718
Spec6	2.38	0.9668	2.90	0.9546
Spec7	2.28	0.9731	2.33	0.9617

Table 3. Table 3 : Averaged density (mg/ml), R ​ M ​ S ​ E ¯ \overline{RMSE} and P ​ S ​ N ​ R ¯ \overline{PSNR} for ROIs of water and iodine density images in Fig. 8 .

ROI	Iodine Density (mg/ml)
	Ground truth	Uncorrected	Kernel-based	MSD-LBTE	Fan beam
1	0.0	-0.3	0.7	0.1	0.1
2	2.0	0.4	2.5	2.0	2.0
3	5.0	1.7	5.7	5.2	5.1
4	10.0	3.8	10.6	9.7	10.0
5	15.0	5.6	15.7	14.9	14.9
6	0.0	-1.1	2.9	0.5	0.6
7	0.0	-1.2	2.4	0.6	0.4
$\overline{RMSE}$	0.0	3.30	1.33	0.47	0.36
$\overline{PSNR}$	-	0.98	13.73	18.76	20.20
ROI	Water Density (mg/ml)
	Ground truth	Uncorrected	Kernel-based	MSD-LBTE	Fan beam
1	1000.0	919.8	994.9	1012.1	1010.7
2	1000.0	936.3	997.4	1018.3	1017.0
3	1000.0	953.4	995.6	1019.0	1023.0
4	1000.0	986.4	997.5	1036.6	1031.0
5	1000.0	1019.6	993.1	1037.3	1036.7
6	1000.0	802.8	921.7	1004.7	991.4
7	1000.0	834.8	943.4	999.6	995.7
$\overline{RMSE}$	0.0	87.67	31.07	23.09	21.65
$\overline{PSNR}$	-	23.08	32.71	34.03	34.60

Table 4. Table 4 : Averaged CT value (HU), E R ​ M ​ E ​ S E_{RMES} adn P ​ S ​ N ​ R ¯ \overline{PSNR} for ROIs of VMIs (70 keV) in Fig. 9 .

ROI	VMI (HU)
	Uncorrected	Kernel-based	CST-LBTE	MSD-LBTE	Fan beam
1	1566.6	1706.8	2002.8	1975.0	1957.0
2	1252.1	1435.6	1765.1	1746.5	1699.2
3	1615.3	1731.6	1903.0	1871.0	1925.8
4	917.0	975.6	1038.4	1034.1	1024.9
5	863.9	952.2	1048.0	1048.2	1024.1
6	783.2	952.1	1062.8	1073.8	1025.3
7	772.3	915.0	1022.5	1026.3	987.1
$\overline{RMSE}$	268.72	140.61	39.72	37.07	0
$\overline{PSNR}$	13.15	20.72	31.80	32.53	-

Equations

$$\overrightarrow{\Omega }\cdot \nabla \varphi (\overrightarrow{r},E,\overrightarrow{\Omega })+{\mu }_t(\overrightarrow{r},E)\varphi (\overrightarrow{r},E,\overrightarrow{\Omega })=S(\overrightarrow{r},E,\overrightarrow{\Omega })$$

$$+$$

$$S^{(k)}(\overrightarrow{r},E,\overrightarrow{\Omega })=$$

$$S^{(k+1)}(\overrightarrow{r},E,\overrightarrow{\Omega })=$$

$${\varphi }^{(k)}(\overrightarrow{r},E,\overrightarrow{\Omega })$$

$$I(\overrightarrow{r})={\int }_0^{E_{max}}dE{\int }_{4\pi }\varphi (\overrightarrow{r},E,\overrightarrow{\Omega })D^{res}(\overrightarrow{r},E)\cdot (\overrightarrow{\Omega }\cdot \overrightarrow{n_d})d\overrightarrow{\Omega },$$

$${\theta }_i$$

$${\phi }_j$$

$${\varphi }_{n,g,m_{n^{\prime }\to n}}^{(k)}$$

$$S_{n,g,m}^{(k+1)}$$

$${\varphi }_{n,g,m_{\overrightarrow{r_s}\to \overrightarrow{r_n}}}^{(0)}$$

$$I_d^S$$

$$I_d^P$$

$$\overrightarrow{{\varphi }^{(k)}}$$

$$\overrightarrow{{\varphi }^{(0)}}$$

$$\overrightarrow{I^S}$$

$$\overrightarrow{I^P}$$

$$\overrightarrow{I^S}$$

$$\overrightarrow{I_l^S}$$

$$\overrightarrow{I_l^P}$$

$${\mathbf{I}}^{\mathbf{p}}$$

$${\mathbf{I}}^{\mathbf{s}}$$

$$\overrightarrow{S}=c_0^{basis}\overrightarrow{S_0^{basis}}+c_1^{basis}\overrightarrow{S_1^{basis}}+.+c_{N_E-1}^{basis}\overrightarrow{S_{N_E-1}^{basis}}$$

$$\overrightarrow{I}$$

$$=f(c_0^{basis}\overrightarrow{S_0^{basis}}+c_1^{basis}\overrightarrow{S_1^{basis}}+.+c_{N_E-1}^{basis}\overrightarrow{S_{N_E-1}^{basis}})$$

$$=c_0^{basis}f(\overrightarrow{S_0^{basis}})+c_1^{basis}f(\overrightarrow{S_1^{basis}})+..+c_{N_E-1}^{basis}f(\overrightarrow{S_{N_E-1}^{basis}})$$

$$=c_0^{basis}\overrightarrow{I_0^{basis}}+c_1^{basis}\overrightarrow{I_1^{basis}}+.+c_{N_E-1}^{basis}\overrightarrow{I_{N_E-1}^{basis}}$$

$$\overrightarrow{I_i^{basis}}=f(\overrightarrow{S_i^{basis}}),i=0,1,.,N_E-1$$

$$\mathbf{S}$$

$${\mathbf{I}}^{\mathbf{basis}}$$

$$\mathbf{I}$$

$${\mathbf{S}}^{\mathbf{T}}={\mathbf{Q}}_{\mathbf{s}}{\mathbf{R}}_{\mathbf{s}}$$

$$\mathbf{S}={\mathbf{R}}_{\mathbf{s}}^{\mathbf{T}}{\mathbf{Q}}_{\mathbf{s}}^{\mathbf{T}}$$

$$MRE=\frac{1}{N_D}\sum _{i=1}^{N_D}|\frac{I_i^{ref}-I_i}{I_i^{ref}}|,$$

$$SSIM=\frac{(2{\mu }_x{\mu }_y+C_1)(2{\sigma }_{xy}+C_2)}{({\mu }_x^2+{\mu }_y^2+C_1)({\sigma }_x^2+{\sigma }_y^2+C_2)}$$

$$\overline{RMSE}=\frac{1}{M}\sum _{i=1}^M\sqrt{\frac{1}{N_{RO{I}_i}}\sum _{j=1}^{N_{RO{I}_i}}(I_{i:j}^{ref}-I_{i;j}{)}^2}$$

$$\overline{PSNR}=\frac{1}{M}\sum _{i=1}^M\left\{20\cdot \lg (\frac{{\max }_{j\in \mathrm{RO}{\mathrm{I}}_{\mathrm{i}}}(I_{i;j})}{RMS{E}_i})\right\}$$

Full text

Linear Boltzmann Transport Equation, Scatter Correction, Multi-Spectral CT

I Introduction

X-ray scatter has been a serious concern since the invention of computed tomography (CT), as it could cause a variety of artifacts in reconstruction and significantly degrade CT image quality if no effective scatter correction or rejection is applied. In the literature, numerous scatter correction or rejection methods have been developed that can be mainly divided into two categories: hardware techniques of scatter rejection and software algorithms of scatter correction. Hardware techniques typically require some sort of interventions in CT system design, including bowtie filter, air-gap and anti-scatter grid[1], which may be less convenient to deploy in clinic.

For software algorithms, statistical methods constitute a category, among which the Monte Carlo (MC) simulation is a classical scattering estimation technique and is considered a gold standard for particle transport simulation[2, 3]. However, substantial computational requirements and random noise make it quite slowly in order to reach signal levels with low statistic errors, although some acceleration strategies have been developed [4, 5, 6, 7]. Recent advancements in deep learning have facilitated the development of learning-based statistical scatter correction methods, and several networks based on supervised and unsupervised learning have been developed [8, 9, 10, 11, 12] which perform quite well in scatter correction, but the limited physical interpretability and the needs of training dataset may restrict their applicability in many aspects. Another category of scatter correction involves analytic methods. The scatter kernel superposition (SKS) method has been developed in the literature with acceptable accuracy and rapid computational speed[13, 14, 15]. However, the design of scatter kernels that can both achieve generality and accuracy remains challenging, and its efficacy may be limited when applied to complex objects.

It is well known that the linear Boltzmann transport equation (LBTE) can be used to describe the statistical behaviors of photons/particles as they transport in space and interact with matter. In contrast to the MC method, which focuses on individual photon behavior at the microscopic level, the LBTE method tries to solve macroscopic photon distribution throughout the entire spatial domain. Scatter correction method based on it has characteristics of high accuracy, good generality, and excellent physical interpretability. Although computationally intensive, the superior parallelism of the equations makes it possible to use GPU and parallel technologies to accelerate the computation process, which has been demonstrated by previous investigations[16, 17].

However, for the conventional LBTE based scatter estimation and correction, it is cumbersome if not impossible to calculate the scatter signals at multiple spectra one by one for multi-spectral CT scans such as KV-switching CT and KV modulation CT. As a result, it is highly desired to develop a method that can unitedly compute the scatter signals at multiple spectra in a single calculation with an improved computational efficiency.

In this work, we proposed a novel we propose a Matrixed-Spectrum Decomposition accelerated LBTE solver (MSD-LBTE). First, we derive a semi-analytic solution of the LBTE in a form based on series expansion according to the times of multiple scattering and phase space discretization scheme. Then, we transform the implicit mapping hidden in LBTE into an explicit matrix. After that, we expand the one-dimensional energy group space into two dimensions and obtain a spectrum-matrixed solver of LBTE. Based upon these findings, we further bring up a new concept of spectrum basis, and a principle of basis selection using QR decomposition, which makes spectrum basis matrix a lower triangular one and significantly reduces computational redundancy.

In our previous conference paper[18], only a preliminary solution was proposed with lack of experimental validations. while in this study, we further optimized the computational efficiency and provided a rigorous mathematical solution formulation. In addition, we have conducted several validations of the computational efficiency and accuracy, and evaluated the feasibility and effectiveness of our proposed method in kV-switching based dual-energy spectral CT and kV-modulation CT scans.

II Methods

II.1 Linear Boltzmann Transport Equation and Its Numerical Solution

Considering the essential characteristics of photons and the specific requirements of a CT scan, the state of photons can be distinctly defined by their spatial position ($\overrightarrow{r}=(x,y,z)$), flight direction ($\overrightarrow{\Omega}=(\theta,\varphi)$), and energy ($E$), collectively forming the phase space $(\overrightarrow{r},\overrightarrow{\Omega},E)$. The behavior of photons in the phase space during CT scanning can be described by the linear Boltzmann transport equation (LBTE) as follow[16]:

[TABLE]

Here, $\phi$ represents the photon fluence, denoting the number of photons passing through a unit section at position $\overrightarrow{r}$ with energy $E$ and flight direction $\overrightarrow{\Omega}$, measured in $m^{-2}sr^{-1}J^{-1}$; $S$ signifies the source term, indicating the number of photons newly generated with energy $E$ and flight direction $\overrightarrow{\Omega}$ in a unit space at position $\overrightarrow{r}$, measured in $m^{-3}sr^{-1}J^{-1}$; $\mu_{t}$ stands for the linear attenuation coefficient with units of $m^{-1}$; $\mu_{s}$ denotes the linear scatter coefficient, representing the fraction of photons with energy $E^{\prime}$ and flight direction $\overrightarrow{\Omega^{\prime}}$ scattered as photons with energy $E$ and flight direction $\overrightarrow{\Omega}$ in a unit space at position $\overrightarrow{r}$; $E_{max}$ corresponds to the maximum energy of photons. In the field of medical imaging, X-ray energies typically do not exceed 200 keV, well below the pair-production effect whose threshold is 1.022 MeV. Within this energy range, photon-matter interactions are primarily characterized by the photoelectric effect, Compton scattering, and Rayleigh scattering.

II.1.1 Semi-Analytic Solution Based on Series Expansion

As an integral-differential equation, it is often difficult to find an analytic solution for LBTE. To solve such as equation, one way is to classify the photons according to the incident times of multiple scattering, so that source distribution and fluence distribution can be unrolled into series form: $\phi(\overrightarrow{r},E,\overrightarrow{\Omega})=\sum_{k=0}^{\infty}\phi^{(k)}(\overrightarrow{r},E,\overrightarrow{\Omega})$ and $S(\overrightarrow{r},E,\overrightarrow{\Omega})=\sum_{k=0}^{\infty}S^{(k)}(\overrightarrow{r},E,\overrightarrow{\Omega})$[16]. When $k=0$, $S^{(0)}$ represents the X-ray source in a CT scan. After that, the integral-difference equation can be divided into an integral equation and a difference equation, separately.

[TABLE]

Here, Eq.(2b) shows the explicit transformation from $\phi^{(k)}(\overrightarrow{r},E,\overrightarrow{\Omega})$ to $S^{(k+1)}(\overrightarrow{r},E,\overrightarrow{\Omega})$. It is well know that Eq.(2a) has an analytic solution as:

[TABLE]

where, both $u$ and $v$ represent integral variables. Since $S^{(0)}$ is determined by the geometry of a CT scan and the spectrum of the X-ray source, we can calculate $\phi^{(0)}$ and then $S^{(1)}$ using Eqs.(3) and (2b), iteratively. Subsequently, by repeatedly applying the same process, we can determine $\phi^{(1)}$ and $S^{(2)}$ based on the previously calculated $S^{(1)}$. This procedure allows us to obtain $S^{(i)}$ and $\phi^{(i)}$ for all $i\in\mathbf{N}$. The detected signals $I(\overrightarrow{r})$ can be computed as indicated in Eq.(4):

[TABLE]

where, $\overrightarrow{n_{d}}$ represents the normal vector of the detector plane, and $D^{res}(\overrightarrow{r},E)$ represents the detector response function. So far, we have derived a series expansion based semi-analytic solution of LBTE, along with an expression for the detector received primary signals and its scatter counterparts.

II.1.2 Discretization of Solution

The LBTE defined in the continuous phase space cannot be directly solved in a digital computing system. Instead, we must discretize the continuous phase space as illustrated in Fig. 1. In reality, the pixelated detectors are naturally discrete. We denote the total number of pixels as $N_{D}$, and use the index $d$ to represent $d$-th pixel. In the spatial position domain, we focus mainly on the scanned object and ignore other external factors. One can construct a partition grid within the scanned object, divide it into $N_{V}$ congruent cuboid voxels, and assume that the material composition and photon distribution in a voxel are uniform. Here, index $n$ denotes $n$-th voxel. For energy domains, the continuous spectrum is discretized into a finite number $N_{E}$ of separate energy groups. Suppose that for the $g$-th energy group, the lower bound is $E_{g}$ and the upper bound is $E_{g+1}$, then the relative intensity $S_{g}=(\int_{E_{g}}^{E_{g+1}}S(E)dE)/(E_{g+1}-E_{g})$. For flight direction domains, we discretized the $4\pi$ solid angle into $M$ discrete angles based on the distribution of longitude and latitude lines. The $(\theta_{i},\phi_{j})$ coordinates of discrete angle $\overrightarrow{\Omega_{m}}$ can be written as:

[TABLE]

When $i=0$ or $N_{\Omega}$, the discrete angle is equal to the "pole" in longitude and latitude, and for different $j$ it represents the same direction. Here, the values of $N_{A}$ and $N_{\Omega}$ are associated as $N_{A}=2N_{\Omega}(N_{\Omega}-1)+2$. Based on the above discretization scheme, we divide the continuous phase space $(\overrightarrow{r},E,\overrightarrow{\Omega})$ into the discrete phase space $(\overrightarrow{r_{n}},E_{g},\overrightarrow{\Omega_{m}})$.

In the discrete phase space, the discrete forms of Eqs.(3) and (2b) can be rewritten as:

[TABLE]

where, $w^{\prime}$ represents the path-length in voxel $i$, $\mu_{t:i,g}$ represents the average linear attenuation coefficient in voxel $i$ in energy group $g$. For voxel $i$ containing more than one material, $\mu_{t:i,g}=\sum_{n=1}^{N}\mu_{m:n,g}\rho_{n}$, where $\mu_{m:n,g}$ represents the average mass attenuation coefficient of material $n$ in energy group $g$ and $\rho_{n}$ represents the density of material $n$ in voxel $i$.

In particular, $S^{(0)}(\overrightarrow{r},E,\overrightarrow{\Omega})$ differs from $S^{(k)}(\overrightarrow{r},E,\overrightarrow{\Omega}),k\geq 1$, because $S^{(0)}$ represents the zero-order X-ray source in a CT scan rather than the higher-order X-ray source in the voxel. Therefore, the spatial position of $S^{(0)}$ cannot be represented by discrete voxel positions because $S^{(0)}$ is outside the scanned object. In continuous phase space, the expression of $S^{(0)}(\overrightarrow{r},E,\overrightarrow{\Omega})$ is given by $S^{(0)}(\overrightarrow{r},E,\overrightarrow{\Omega})=O(\overrightarrow{\Omega})S(E)\delta(\overrightarrow{r}-\overrightarrow{r_{s}})$, where $\overrightarrow{r_{s}}$ represents the spatial location of the X-ray source, $S(E)$ represents the spectrum of X-ray source and $O(\overrightarrow{\Omega})$ represents the focus range of the X-ray source. After discretization, it is seen that $S^{(0)}_{g,m}=O_{m}\cdot S_{g}$. In the case of $k=0$, Eq.(6a) is simplified as:

[TABLE]

Finally, the discrete forms of the primary and scatter signals received by the detector can be written as:

[TABLE]

where, $\overrightarrow{r_{d}}$ represents the spatial position of pixel $d$, $a_{d}$ represents the area of pixel $d$, $D^{res}_{d,g}$ represents the discrete form of $D^{res}(\overrightarrow{r},E)$, and $N_{S}$ represents the upper limit of scattering. In practice, the number of photons in a higher-order of multiple scattering is smaller and their contributions to the overall scattered photons is mostly negligible, so one can truncate at a finite order of multiple scattering. So far, we have turned formulas defined in continuous phase space into discrete forms, which can be solved in a computer.

II.2 Spectrum-Matrixed Solver for LBTE and a Concept of Spectrum Basis

For ease of narration, we convert the above discrete expressions into matrix form. Consider $\phi^{(k)}_{n,g,m},k\geq 0$ and $S^{(k)}_{n,g,m},k\geq 1$ as ${N_{V}\times N_{E}\times N_{A}}$-dimension vectors, respectively, denoted by $\overrightarrow{\phi^{(k)}}$ and $\overrightarrow{S^{(k)}}$, where $\overrightarrow{\phi^{(k)}},\overrightarrow{\phi^{(k)}}\in\mathbb{R}^{N_{V}N_{E}N_{A}}$. Similarly, $I^{S}_{d}$, $I^{P}_{d}$ are represented by vectors $\overrightarrow{I^{S}}$, $\overrightarrow{I^{P}}$ $\in\mathbb{R}^{N_{D}}$. and $S_{g}$ is represented by the vector $\overrightarrow{S}\in\mathbb{R}^{N_{E}}$. Then the Eqs.(6)-(8) can be rewritten to the following matrix forms:

[TABLE]

combining Eqs.((9a), (9b), and (9c), one gets:

[TABLE]

Here,

$\mathbf{T_{0}}\in\mathbb{R}^{N_{V}N_{E}N_{A}\times N_{E}N_{A}}$ is the primary transport matrix, and $\mathbf{R_{0}}\in\mathbb{R}^{N_{D}\times N_{E}N_{A}}$ is primary detector matrix. Both of them are sparse matrices determined by the source position and material density distribution. $\mathbf{R_{0}}$ also depends on the detector position and response function.

$\mathbf{O}\in\mathbb{R}^{N_{E}N_{A}\times N_{E}}$ is primary X-ray source focus range matrix, which determined by X-ray source position and its field of view (FOV).

$\mathbf{T}\in\mathbb{R}^{N_{V}N_{E}N_{A}\times N_{V}N_{E}N_{A}}$ is the secondary transport matrix, $\mathbf{C}\in\mathbb{R}^{N_{V}N_{E}N_{A}\times N_{V}N_{E}N_{A}}$ is secondary scatter matrix, and $\mathbf{R}\in\mathbb{R}^{N_{D}\times N_{V}N_{E}N_{A}}$ is secondary detector matrix. All of them are sparse matrices determined by material density distribution. $\mathbf{R}$ also depends on the detector position and response function.

II.2.1 Phase Space Expansion and Spectrum-Matrixing

In order to fully utilize the parallelism of LBTE, we introduce a label dimension ($l$) to the phase space, which assigns an additional information to the photon emitted from the X-ray source. As a result, the phase space expands from $(\overrightarrow{r_{k}},E_{g},\overrightarrow{\Omega_{m}})$ to $(\overrightarrow{r_{n}},E_{g},\overrightarrow{\Omega_{m}},l)$, $l\in\{0,1,...,N_{L}-1\}$. Then, $S^{(k)},\phi^{(k)}$ in Eq.(9) are expanded as $\overrightarrow{S^{(k)}_{l}},\overrightarrow{\phi^{(k)}_{l}},l=0,1,2,...,N_{L}-1$, $\overrightarrow{I^{S}}$, $\overrightarrow{I^{P}}$ are expanded as $\overrightarrow{I^{S}_{l}},\overrightarrow{I^{P}_{l}},l=0,1,2,...,N_{L}-1$, $\overrightarrow{S}$ is expanded as $\overrightarrow{S_{l}},l=0,1,2,...,N_{L}-1$, and $\mathbf{T},\mathbf{T_{0}},\mathbf{R},\mathbf{R_{0}},\mathbf{C}$ and $\mathbf{O}$ will not be expanded because they are independent of the information of photons and instead only dependent of geometry, detector response, X-ray source FOV or material density distribution. After such as phase space expansion, Eqs.(9d) and (10) can be rewritten as:

[TABLE]

and a series of spectrum vectors $\overrightarrow{S_{l}},l=0,1,2,..,N_{L}-1$ can be written as a single spectrum matrix $\mathbf{S}=[\overrightarrow{S_{0}},\overrightarrow{S_{1}},...,\overrightarrow{S_{N_{L}-1}}]$, where each column vector of spectrum matrix $\mathbf{S}$ represents a spectrum vector $\overrightarrow{S}$ under different $l$. As such a process means that the spectrum space changes from $\mathbb{R}^{N_{E}}$ to $\mathbb{R}^{N_{E}\times N_{L}}$, it is called “spectra-matrixing” in this paper.

Similarly, $\overrightarrow{I_{l}},\overrightarrow{I_{s}},l=0,1,2,..,N_{L}-1$ can also be written as matrix $\mathbf{I^{S}}=[\overrightarrow{I^{S}_{0}},\overrightarrow{I^{S}_{1}},...,\overrightarrow{I^{S}_{N_{L}-1}}],\mathbf{I^{P}}=[\overrightarrow{I^{P}_{0}},\overrightarrow{I^{P}_{1}},...,\overrightarrow{I^{P}_{N_{L}-1}}]$ and the detected signal space changes from $\mathbb{R}^{N_{D}}$ to $\mathbb{R}^{N_{D}\times N_{L}}$. The same applies to $\overrightarrow{\phi^{(k)}}$ and $\overrightarrow{S^{(k)}}$, as they will be expanded into matrices $\mathbf{\Phi^{(k)}}$ and $\mathbf{S^{(k)}}$ with $\mathbb{R}^{N_{V}N_{E}N_{A}}$ expanded into $\mathbb{R}^{N_{V}N_{E}N_{A}\times N_{L}}$.

If we set $\mathbf{A^{s}}=\mathbf{R}\left(\sum_{k=1}^{N_{S}}\mathbf{(TC)}^{k-1}\right)\mathbf{C}\mathbf{T_{0}}\mathbf{O}$ and $\mathbf{A^{p}}=\mathbf{R_{0}O}$, then Eqs.(9d) and (10) become:

[TABLE]

In such a way above, the mappings from the spectrum to primary and scatter signals ,implied in LBTE ,are represented explicitly as matrices $\mathbf{A^{p}}$ and $\mathbf{A^{s}}$. In order to facilitate the narration, we use $\mathbf{A}\in\mathbb{R}^{N_{D}\times N_{E}}$ to represent these two matrices, and use $\mathbf{I}\in\mathbb{R}^{N_{D}\times N_{L}}$ to represent $\mathbf{I^{P}}$ and $\mathbf{I^{S}}$, where $\mathbf{S}\in\mathbb{R}^{N_{E}\times N_{L}}$.

Here, it is worth noting that the introduction of label dimension adds just a small computational overhead. When $N_{L}=1$, let the computation cost be 1 unit. When $N_{L}=2$, the computation cost does not increase to 2 unit but to $1+t$ unit with $t<1$ (for $\forall N_{L}$, the computation cost equals $1+(N_{L}-1)t$ unit). Because $\mathbf{A}$ are independent of dimension $l$, we only need to calculate them once. Adding $N_{L}$ labels only increases about $2N_{S}(N_{L}-1)$ times of sparse matrix multiplication as $(N_{L}-1)t$ unit computation cost (refer to Eq.(9)), which is relatively small when compared with the cost of computing these sparse matrices. The subsequent experiments will validate this assumption. So far, we have established a matrixed-spectrum solver for LBTE (MS-LBTE).

II.2.2 A Novel Concept with Spectrum Basis Vectors

In order to obtain different scatter signals for multiple spectra in a single computation, we need to optimize the process of LBTE solving and reduce the computational redundancy. Fortunately, we notice that process of solving LBTE is actually to find a linear implicit mapping from the spectrum $S(E)$ to the detected signal $I(\vec{r})$ as $f:S(E)\rightarrow I(\overrightarrow{r})$. After discretization as shown in Sec. II.1, the continuous variable $S$ becomes an $N_{E}$-dimensional vector $\overrightarrow{S}$, and $I$ becomes $N_{D}$-dimensional vectors $\overrightarrow{I}$. The following relationship still holds: $\overrightarrow{I}=f(\overrightarrow{S})$.

Since $\overrightarrow{S}$ is in $\mathbb{R}^{G}$, it is easy to pick a set of linearly independent $\overrightarrow{S^{basis}_{0}},\overrightarrow{S^{basis}_{1}},...\overrightarrow{S^{basis}_{G-1}}$ as the spectrum basis vectors in $\mathbb{R}^{G}$. Then, for $\forall\overrightarrow{S}\in\mathbb{R}^{G}$, we can get a set of coefficients $c^{basis}_{0},c^{basis}_{1},...,c^{basis}_{N_{E}-1}$ that satisfy:

[TABLE]

Given that $f$ is a linear mapping, we have:

[TABLE]

where,

[TABLE]

with $\overrightarrow{I^{basis}_{i}},i=0,1,...,N_{E}-1$ being detected signal basis vectors.

By using such spectrum bases, for any spectrum $\overrightarrow{S}$, we do not need to calculate the corresponding $\overrightarrow{I}$ one by one. Instead, we can first select a set of spectrum basis vectors and obtain a set of coefficients that $\overrightarrow{S}$ resolves under that bases as in Eq.(13). Then we calculate a set of detected signals basis vectors based on spectrum basis vectors as in Eq.(15). Finally, we can get the detected signal $\overrightarrow{I}$ corresponding to $\overrightarrow{S}$, which is equal to a linear combination of detected signal bases as in Eq.(14).

An obvious advantage of such a scheme above is that no matter how many spectra ($N_{spec}$) we used (even $N_{spec}>N_{E}$), we only need to compute $N_{E}$ detected signal vectors $\overrightarrow{I^{basis}_{i}},i=0,1,...,N_{E}-1$ based on G selected spectrum basis vectors $\overrightarrow{S^{basis}_{i}},i=0,1,...,N_{E}-1$ as in Eq.(15). In addition, if we only use $N_{spec}$ different spectra $\overrightarrow{S_{i}},i=0,1,...,N_{spec}-1$ and $N_{spec}<N_{E}$, these vectors cannot fill the whole $\mathbb{R}^{N_{E}}$ space. In fact, the ${N_{spec}}$ spectra vectors span only the lower dimensional $\mathbb{R}^{N_{spec}}$ space ($\mathbb{R}^{N_{spec}}\subset\mathbb{R}^{N_{E}}$) in fact. As a result, we can select fewer spectrum basis vectors $\overrightarrow{S^{basis}_{i}},i=0,1,...,{N_{spec}}-1$, leading to fewer computations of detected signal basis vectors $\overrightarrow{I^{basis}_{i}},i=0,1,...,{N_{spec}}-1$.

In MS-LBTE, the implicit mapping $f$ is transformed into an explicit matrix $\mathbf{A}\in\mathbb{R}^{N_{D}\times N_{E}}$, and spectrum vector $\overrightarrow{S}\in\mathbb{R}^{N_{E}}$ and detected signal vector $\overrightarrow{I}\in\mathbb{R}^{N_{D}}$ are expanded into $\mathbb{R}^{N_{E}\times N_{L}}$ and $\mathbb{R}^{N_{E}\times N_{L}}$. With no spectrum basis strategy above, if there are $N_{spec}$ different spectra, it is considered that $l$ marks different spectrum vectors used so that $N_{L}=N_{spec}$, which can be written as $\mathbf{I}=\mathbf{A}\cdot\mathbf{S}$, with $\mathbf{I}=[\overrightarrow{I_{0}},\overrightarrow{I_{1}},...,\overrightarrow{I_{N_{spec}-1}}]\in\mathbb{R}^{N_{D}\times N_{spec}}$ and $\mathbf{S}=[\overrightarrow{S_{0}},\overrightarrow{S_{1}},...,\overrightarrow{S_{N_{spec}-1}}]\in\mathbb{R}^{N_{E}\times N_{spec}}$.

Using the spectrum basis strategy, we can write multiple spectrum basis vectors as a matrix $\mathbf{S^{basis}}=[\overrightarrow{S^{basis}_{0}},\overrightarrow{S^{basis}_{1}},\overrightarrow{S^{basis}_{2}},...]$. Similarly we have $\mathbf{I^{basis}}=[\overrightarrow{I^{basis}_{0}},\overrightarrow{I^{basis}_{1}},\overrightarrow{I^{basis}_{2}},...]$ and $\overrightarrow{C^{basis}}=[c^{basis}_{0};c^{basis}_{1};c^{basis}_{2};...]$, then we can obtain the basis decomposition of $\overrightarrow{S}$ as $\overrightarrow{S}=\mathbf{S^{basis}}\overrightarrow{C^{basis}}$ corresponding to Equation (13) and the linear combination of $\overrightarrow{I}$ as $\overrightarrow{I}=\mathbf{I^{basis}}\overrightarrow{C^{basis}}$ corresponding to Eq.(14).

If there are $N_{spec}$ different spectra, now $\overrightarrow{I_{i}}=\mathbf{A}\overrightarrow{S_{i}},i=0,1,...,N_{spec}-1$. Let’s $\mathbf{I}=[\overrightarrow{I_{0}},\overrightarrow{I_{1}},...,\overrightarrow{I_{N_{spec}-1}}]$ ($\mathbf{I}\in\mathbb{R}^{N_{D}\times N_{spec}}$) and $\mathbf{S}=[\overrightarrow{S_{0}},\overrightarrow{S_{1}},...,\overrightarrow{S_{N_{spec}-1}}]$ ($\mathbf{S}\in\mathbb{R}^{N_{E}\times N_{spec}}$), we can obtain $\mathbf{S}=\mathbf{S^{basis}C^{basis}}$ and $\mathbf{I}=\mathbf{I^{basis}C^{basis}}$, with $\mathbf{C^{basis}}=[\overrightarrow{C^{basis}_{0}},\overrightarrow{C^{basis}_{1}},\overrightarrow{C^{basis}_{2}},...]$ and $\overrightarrow{S_{i}}=\mathbf{S}\overrightarrow{C^{basis}_{i}}$. Then we can finally obtain the following expressions,

[TABLE]

Such an introduction of spectrum basis matrix serves as a transitional position, which enable us to freely select the spectrum matrix used in LBTE solving without being restricted by the actual spectrum matrix determined by X-ray source as in the following Sec.II.3.

II.3 Spectrum Basis Optimization with QR Decomposition

When $N_{spec}\geq N_{E}$, we can set $N_{L}=N_{E}$. The spectrum matrix space is $\mathbb{R}^{N_{E}\times N_{E}}$. We select a set of unit orthogonal spectrum basis vectors $\overrightarrow{e_{0}},\overrightarrow{e_{1}},...,\overrightarrow{e_{N_{E}-1}}$ in $\mathbb{R}^{G}$, where $\overrightarrow{e_{0}}=[1,0,...,0],\overrightarrow{e_{1}}=[0,1,...,0],...,\overrightarrow{e_{N_{E}-1}}=[0,0,..,1]$. Then,$\mathbf{S^{basis}}=[\overrightarrow{e_{0}},\overrightarrow{e_{1}},...,\overrightarrow{e_{N_{E}-1}}]=\mathbf{E}$, which is an identity matrix in $\mathbb{R}^{N_{E}\times N_{E}}$. In theory, the computational cost is $1+(N_{E}-1)t$ in this case. Under such circumstances, as the photoelectric effect, Rayleigh scatter and Compton scattering effect are considered in the equation, the energy of photons will merely decrease rather gradually than increase. Let energy group [math] be the highest energy group and energy group $G-1$ be the lowest energy group. Then, for the spectrum basis vectors $\overrightarrow{e_{l}}$ (i.e., the energy spectrum where only energy group $l$ has a non-zero value), only the behavior of photons in the lower energy groups $l,l+1,l+2,...,N_{E}-1$ needs to be taken into account, while the behavior of the higher energy groups $1,2,...,l-1$ need no consideration as shown in Fig. 3(a). This enables the computational overhead to be further reduced, from $1+(N_{E}-1)t$ unit to $1+(N_{E}-1)t\cdot\frac{N_{E}+1}{2N_{E}})$ unit instead.

When $N_{spec}<N_{E}$, we set $N_{L}=N_{spec}$, the spectrum matrix space is $\mathbb{R}^{N_{E}\times N_{spec}}$. We can also select a set of spectrum basis vectors $\overrightarrow{b_{0}},\overrightarrow{b_{1}},...,\overrightarrow{b_{N_{E}-1}}$ which form a lower triangular matrix $\mathbf{B}$ as shown in Fig. 3(b).

Luckily, the validity of such a selection can be guaranteed by the QR decomposition of a matrix. In other words, We can obtain $\mathbf{S^{basis}}$ and $\mathbf{C^{basis}}$ by finding the QR decomposition of the transpose of the matrix $\mathbf{S}$ as,

[TABLE]

where $\mathbf{Q_{s}}$ is an orthogonal of size $N_{spec}\times N_{spec}$ and $\mathbf{R_{s}}$ is an upper triangular matrix. Taking the transpose of both sides, we obtain,

[TABLE]

Then, compared with Eq.(16a), we can get $\mathbf{S^{basis}}=\mathbf{R_{s}^{T}},\mathbf{C^{basis}}=\mathbf{Q_{s}^{T}}$.

Now, $\mathbf{S^{basis}}$ is a lower-triangular matrix. Considering the differences between different spectra, we can think of $\mathbf{S}$ as a column full rank matrix, which avoids column pivoting in the QR decomposition. It is worth noting that the $\mathbf{S^{basis}}$ obtained by the QR decomposition may contain negative values which has no mean in physics, but will have no negative impact on final outcome in computation as our LBTE solving process is linear with respect to the number of photon at a specific energy level. In such a case, the computational cost is no longer $1+(N_{spec}-1)t$ unit. It becomes $1+(N_{spec}-1)t\cdot(1-\frac{N_{spec}-1}{2N_{E}})$ unit instead, leading to an LBTE solver with better computational efficiency. In this paper, such a matrixed-spectrum decomposition accelerated solver of LBTE is called MSD-LBTE.

II.4 Unified and Accelerated Scatter Correction in Multi-Spectral CT

Based on the above deviations and findings, a unified scatter correction for multi-spectral CT can be straightforwardly implemented, as illustrated in Fig. 4. Multiple steps are usually involved. First, it is essential to preprocess the projection data to generate the initial CT images. Subsequently, the distributions of CT values need to be converted into material densities, serving as inputs of scanned objects for MSD-LBTE to. After that, the spectrum matrix, consisting of multiple spectrum vectors, is decomposed into a product of the spectrum basis matrix and a coefficient matrix. The spectrum basis matrix, along with the multi-material densities of the object, is fed into MSD-LBTE to generate the scatter signal bases. The bases are then multiplied by the coefficient matrix to derive the actual scatter signals needed for each specific spectrum. Ultimately, a set of scatter signals are obtained, usually with an appropriate scaling factor for scatter correction determined in advance practically. In the end, the scatter signals are then subtracted from the total transmission data to generate scatter-corrected projections from which final reconstructions can be done. In comparison with the conventional LBTE based ones, our method offers a clear advantage of simultaneously scatter computation for multiple spectra, leading to a significantly improved computational efficiency, with no sacrifice in computational accuracy.

II.5 Experimental Setups

In order to validate our proposed MSD-LBTE method, we measured the computational time of using conventional LBTE methods to estimate scatter distributions, along with convetional LBTE and MS-LBTE methods, with different numbers of spectra in CT scans in the same hardware environment (Nvidia RTX A5000, $N_{V}=16^{3}$, $N_{E}=10$ and $N_{A}=114$). In terms of computation accuracy, we used GEANT4 to conduct MC simulations to generate primary and scatter signals at different spectra as references. The detector consisted of $80\times 80$ pixels with a pixel size of $5.4$ mm; seven different spectra were used in MC simulations; the angular flux of X-ray source was set to $4.37142\times 10^{9}$ $\mathrm{sr^{-1}}$. The source-to-isocenter distance (SID) and the source-to-detector distance (SDD) were $750.2$ $\mathrm{mm}$ and $1184.04$ $\mathrm{mm}$, respectively, with 7.4∘ cone angle of the X-ray source. In MSD-LBTE, $N_{E}=9$, $N_{V}=24^{3}$ and $N_{A}=114$. The phantoms used were a water cylinder with a diameter in $200$ $\mathrm{mm}$ and a height in $200$ $\mathrm{mm}$, and a voxelized head phantom consisting of water and bone with different densities.

To further validate efficiency of the unified scatter correction method based on our proposed MSD-LBTE, we conducted a set of physical experiments, on an experimental CBCT tabletop system with a Varian 4030 DX flat-panel detector. The GammexTM multi-energy CT phantom and the Kyoto head phantom were used in our experiments. The kernel-based method[15], and the CBCT Software Tools (CST) [16, 19] with embedded LBTE-based scatter correction were also conducted as comparison. Some key parameters are listed in TABLE 1, with kV-switching mimicked by dual-kVp scans at 80/140 kVp for the Multi-Energy phantom and at 80/120 kVp for the head phantom, respectively.

In addition, we also conducted a numerical simulation to explore potential of applying MSD-LBTE in kV-modulation enabled CT [20, 21]. Here, a virtual kV-modulation CBCT scan was carried out, using an elliptical cylinder whose long axis is set as $24$ cm, short axis as $12$ cm, and height as $20$ cm. Projection data with and without scatter signals were generated with a kV-modulation designed to optimize the SNR over dose [22]. Some key parameters are also listed in TABLE 1, with X-ray source varying in the range of 70 to 140 Kilo-voltage peak (kVp) with an interval of 0.5 keV. Meanwhile, in order for comparison, we also approximate the continuously varying tube voltages to one fixed energy level of 105 kVp, and two fixed energy levels of 87.5 and 122.5 kVp, respectively, for the conventional LBTE-based scatter estimation and correction.

II.6 Evaluation metrics

In the comparison of scatter estimation accuracy, the mean relative error (MRE) and the structural similarity index (SSIM) are calculated in the images of scatter signals. The MRE is defined as,

[TABLE]

where $I^{ref}_{i}$ is the $i$-th pixel value of the reference image, and $I_{i}$ is the $i$-th pixel value of the obtained image. SSIM values were calculated by,

[TABLE]

where, $\mu_{x}$, $\mu_{y}$, $\sigma_{x}$, $\sigma_{y}$ and $\sigma_{xy}$ represent the local mean, standard deviation and cross-covariance of images x and y. In addition, in order to avoid errors due to the dynamic range difference,, we normalized the scatter signals to [0,1] when calculating SSIM.

For the numerical simulations and physical experiments, we use an average error of water or iodine density among the regions of interest (ROIs) ($\overline{RMSE}$) defined as,

[TABLE]

Where, $I^{ref}_{i;j}$ is the $j$-th pixel value of the $i$-th ROI of the reference image, $I_{i;j}$ is the $j$-th pixel value of the $i$-th ROI of the obtain image. Likewise, we use an average peak signal-to-noise ratio of water or iodine density in ROIs ($\overline{PSNR}$) defined as,

[TABLE]

where $RMSE_{i}$ is the root mean square error of the $i$-th ROI.

III Results

III.1 Computational Efficiency

The actual computational costs for different number of spectra are shown in Fig. 5 with energy group $N_{E}=9$. For the conventional LBTE-based scatter correction (convLBTE), the computation time increases linearly as the number of spectra increases, as expected, with about $0.24$ s per spectrum. For MS-LBTE, it still increases linearly but computation time per spectrum is decreased by a factor of two-third, thanks to the sharing of common calculations across different spectra being considered.

Using our proposed MSD-LBTE, the computational cost further reduced significantly even compared with MS-LBTE. It is observed that as the number of spectra increases, the computation time will reaches a saturation point of maximum limit. In the case of two spectra, the computation time can be reduced to 66.7% of that of the conventional LBTE. The number drops to 52.8% for three spectra. When 9 spectra used, the computation time of MSD-LBTE saturates, which is equivalent to 26.5% of that of the conventional LBTE.

III.2 Computational Accuracy

In Figs. 6 and 7, we compare the simulated scatter signals of the MC simulation with that of MSD-LBTE, using the water cylinder phantom and the voxelized head phantom, respectively, where the horizontal and vertical line profiles across the center of the primary and scatter signal images are also plotted. It is seen that overall results of MSD-LBTE agrees well with the MC simulation. The MRE and SSIM of the scatter signals are summarized in TABLE 2. Quantitatively, for the water cylinder phantom, MRE is less than $2.39\%$ and SSIM is greater than $96.68\%$ for every spectrum, respectively; while for the voxelized head phantom, they are less than $2.94\%$ and greater than $95.46\%$ across all spectra, respectively.

It is worth noting that although our proposed MSD-LBTE achieved a comparable accuracy to the MC method, in scatter computation, with much faster speed, some noticeable discrepancies between MC and MSD-LBTE appear mainly density in the peripheries. This is because that the sampling (discretization) is not dense enough in those regions, where the primary/scatter signal changes sharply. Increasing the discretization density (such as from $N=24^{3}$ to $N=32^{3}$) can reduce those discrepancies. Fortunately, however, they will not affect scatter correction much in the end, as the primary signals are relatively much bigger (resulting in a negligible residual error in SPR).

III.3 Scatter Correction in KV-Switching Spectral CT scan

III.3.1 Multi-Energy phantom study

The reconstructed CT images from low- and high-energy projections, as well as basis material projections after material decomposition after material decomposition are shown in Fig. 8. It is seen that, with scatter correction there are serious artifacts in the CT images as expected, and the water and iodine densities are way off their nominal values. Our proposed MSD-LBTE method significantly removes the artifacts and recoveries water and iodine densities. The plotted profiles show that after scatter correction, the CT values from cone-beam scans, along with their water and iodine densities after material decomposition are in good agreement with the corresponding fan-beam scan (reference).

Quantitatively, TABLE 3 lists the averaged density, $\overline{RMSE}$ and $\overline{PSNR}$ in the selected ROIs, where the groundtruth of iodine and water density are referred to the user manual of the Gammex Multi-Energy CT phantom. It is seen that the results from uncorrected cone-beam scan have significant bias. After scatter correction by our proposed method, in iodine density, the $\overline{RMSE}$ is reduced from $3.30$ mg/ml to $0.47$ mg/ml and the $\overline{PSNR}$ is increased from $0.98$ dB to $18.76$ dB; in water density, the $\overline{RMSE}$ is reduced from $87.67$ mg/ml to $23.09$ mg/ml and the $\overline{PSNR}$ is increased from $23.08$ dB to $34.03$ dB. Besides, compared with the conventional LBTE method of separately conducting scatter correction for high and low energy in dual-energy kV-switching CT, MSD-LBTE reduces the computation time by $1/3$.

III.3.2 Head phantom study

The low- and high-energy reconstruction image along with material decomposition results using different scatter correction methods and their profiles are shown in Fig. 9. It also shows the difference images between cone-beam VMIs with or without scatter correction with respect to the fan-beam VMI. Preliminary results show that the CT reconstructions after scatter correction by using either CST or our proposed method are in good agreement with that from the fan-beam scans. Quantitatively, TABLE 3 lists the averaged CT values, $\overline{RMSE}$ ,and $\overline{PSNR}$ of the selected ROIs of VMIs, where the fan-beam scan results are utilized as a reference. Both CST and our proposed method show a good performance in scatter correction, demonstrating the effectiveness and accuracy of the LBTE-type method.

III.4 Potential in KV-Modulation CT scan

As shown in Fig. 10, before scatter correction, there are serious scatter artifacts in CT images reconstructed from the kV-modulated projections with $RMSE=187.6$, significantly deviating from the reference, which emphasizes the need of a proper scatter correction. After single-kVp-approximated LBTE correction which costs 0.24s (each view), the $RMSE$ is reduced to $11.41$ but residual scatter artifacts are quite obvious. Double-kVp-approximated LBTE correction which costs 0.48s further improves the image quality, with the $RMSE$ decreased to $5.89$ while some residual scatter artifacts still noticeable. Using our proposed MSD-LBTE which costs 0.57s, the reconstructed CT images have the best uniformity visually and the lowest $RMSE=4.78$, with no obvious scatter artifacts. Compared with double-kVp-approximation LBTE correction, the computational time of using MSD-LBTE, which obtains scatter signals for all kVps without approximation, just increases 18% overall. It is seen that, with our proposed method, scatter correction will no longer be a bottle-neck for CT scan with kV-modulation, in terms of computation at least.

IV Discussions and Conclusion

In this work, we proposed a matrixed-spectrum decomposition accelerated linear Boltzmann transport equation solver, and established a unified and enhanced scatter correction basing on it. Our proposed method significantly improves computational efficiency at multiple spectra while preserving high levels of accuracy.

For our implementation of the LBTE method, while the computational load introduced by the phase space expansion is substantially small portion of the entire runtime, achieving the desired optimization is hindered by the time-consuming IO read and write conflicts in multi-threaded computer programming. Further balancing the IO read and write operations to mitigate access conflicts in multi-thread process could further improve the acceleration in computation speed. Moreover, the MSD-LBTE can be applied not only in multi-spectral CT but also in multi-source CT configurations such as stationary CT, where $l$ can be employed to distinguish between different positions of sources with different spectra potentially. Also, our work could play a significant role for learning-based tasks where efficiently generating CT scans from publicly available CT datasets for a range of different kVp’s, with realistic scatter signals, becomes more and more important.

There exists some limitations for the LBTE-type method. First, a challenge lies in performing scatter correction on truncated projection data, regardless of whether the truncation occurs in the $\mathrm{x-y}$ direction or the $\mathrm{z}$ direction. Portions of the scanned objects that are outside the feild of view (FOV) still contribute to the scatter signals. Second, it is difficult to conduct scatter correction for objects with intricate structure composed of various materials with significant differences in attenuation and scattering characteristics. In such a case, the accuracy of material conversion could be notably compromised. These two issues exists for most post-reconstruction methods.

Future endeavors includes further optimization of the implementation to enhance both computation speed and accuracy, exploration of new application, and more in-depth assessments and analyses in a variety of multi-spectral CT scenarios beyond fast-kV switching based spectral CT and kV-modulation enabled CT scans.

V Acknowledgment

This project was supported in part by Grants from the National Key R&D Program of China (No. 2022YFE0131100 and 2024YFC2417500), and in part by Beijing Natural Science Foundation (No. L252021).