MB-DECTNet: a model-based unrolling network for accurate 3D dual-energy CT reconstruction from clinically acquired helical scans

In this study, we modeled the linear attenuation coefficient (LAC) as a sum of products of spatial functions times energy functions as

$$\begin{eqnarray}&&\mu (x,E)=\displaystyle \sum _{i}{\mu }_{i}(E){c}_{i}(x),\end{eqnarray} \tag{ 1 }$$

where $x\in {{\mathbb{R}}}^{3}$ denotes the spatial location in the image domain, E denotes the energy index, μ(x, E) denotes the LAC value at location x and energy E, i indexes the basis component, μ_i (E) denotes the LAC of the selected material, and c_i (x) denotes the basis component weights that need to be reconstructed.

In DECT statistical iterative reconstruction (SIR), the stationary point is reached by minimizing the objective function with respect to the target image c = {c₁, c₂} as

$$\begin{eqnarray}&&\hat{{\boldsymbol{c}}}=\arg \,{\min }_{{\boldsymbol{c}}}{ \mathcal D }({\boldsymbol{c}},{\boldsymbol{d}})+{ \mathcal R }({\boldsymbol{c}}),\end{eqnarray} \tag{ 2 }$$

where ${ \mathcal D }$ denotes the data-fidelity term, ${ \mathcal R }$ denotes the penalty term, d = {d_L , d_H } denotes the set of measurements acquired at low- and high-kVp spectra.

The choice of data fidelity term depends on the assumed nature of the measured noise. If Poisson noise is assumed, I-divergence may be employed (O'Sullivan and Benac 2007, Long and Fessler 2014, Chen et al 2016), whereas a weighted least square is utilized if Gaussian noise is assumed (Huh and Fessler 2009, Xu et al 2009, Zhang et al 2014). Consequently, the target image c undergoes iterative updates based on the penalized selected fidelity function.

Compared to image-domain techniques, SIR reconstructs more accurate images with higher resolution and less uncertainty (Zhang et al 2018, Medrano et al 2022). However, DECT SIRs are usually computationally intensive due to the large system operator and the low convergence rate caused by the surrogate function that addresses the nonlinearity of the forward model.

The framework of the traditional DECT SIR is shown in the upper part of figure 2, where (2) is solved by iteratively updating the target image based on the update direction. The update step for such algorithms can be written as

$$\begin{eqnarray}&&{{\boldsymbol{c}}}^{k}={{\boldsymbol{c}}}^{k-1}+{\gamma }^{k-1}{\rm{\Delta }}{{\boldsymbol{c}}}^{k-1},\end{eqnarray} \tag{ 3 }$$

where k denotes the iteration number, Δ c ^k−1 denotes the update direction, and γ^k−1 is the step size. This kind of algorithm usually requires thousands of iterations to converge.

Zoom In Zoom Out Reset image size

Our model-based deep unrolling network (MB-DECTNet) replaces the update step with a neural network, and Δc with the data-consistency (DC) layer, which ensures that the unrolled network, (lower picture of figure 2) is model-based. Here Ψ_k denotes the kth deep image updater. MB-DECTNet consists of several update blocks that are analogous to iterations in SIR. Intuitively, these networks learn a more aggressive, high-dimensional step size, as well as a denoiser from the training set. With this model-based unrolled network, we are trying to get the stationary point in just a few iterations.

The structure of the update block is shown in figure 3. The DC layer takes the image from the previous block and the transmission data as the inputs to derive the update direction, which quantifies the inconsistency between the current estimate and the measurement. The update direction can vary with respect to the algorithm type. For gradient descent optimizers, it is defined as the simple gradient of the objective function with respect to the previous estimate c ^k−1,

$$\begin{eqnarray}&&{\mathscr{DC}}({{\boldsymbol{c}}}^{k-1},d)={\left.\displaystyle \frac{\partial { \mathcal D }({\boldsymbol{c}},{\boldsymbol{d}})}{\partial {\boldsymbol{c}}}\right|}_{{\boldsymbol{c}}={{\boldsymbol{c}}}^{k-1}},\end{eqnarray} \tag{ 4 }$$

or as the stationary point of the surrogate function minus the previous estimate c ^k−1 when SIR minimizes a surrogate of the objective function

$$\begin{eqnarray}&&{\mathscr{DC}}({{\boldsymbol{c}}}^{k-1},d)=\left\{\arg \mathop{\min }\limits_{{\boldsymbol{c}}}\tilde{{ \mathcal D }}\{{\boldsymbol{c}},{\boldsymbol{d}},{{\boldsymbol{c}}}^{k-1}\}\right\}-{{\boldsymbol{c}}}^{k-1},\end{eqnarray} \tag{ 5 }$$

where k denotes the index of the current block, and $\tilde{{ \mathcal D }}$ denotes the surrogate function of the data fidelity term at c ^k−1. Since the update block can be regarded as a denoiser, an explicit penalty term ${ \mathcal R }(c)$ in (2) is not necessary.

Zoom In Zoom Out Reset image size

The DECT image and DC output are stacked along the channel dimension as the input to the deep learning-based image updater Ψ_k . The subsequent estimate is evaluated by the pointwise summation of the previous estimate and the network output. With this residual learning strategy, MB-DECTNet is able to focus on the inconsistency of the estimate against the training target (He et al 2015, Zhang et al 2016), making training more efficient. To sum up, in the proposed model-based unrolling network, (3) in DECT SIR is rewritten as

$$\begin{eqnarray}&&{{\boldsymbol{c}}}^{k}={{\boldsymbol{c}}}^{k-1}+{{\rm{\Psi }}}_{k}\left({\mathtt{stack}}({\mathscr{DC}}({{\boldsymbol{c}}}^{k-1},d),{{\boldsymbol{c}}}^{k-1})\right).\end{eqnarray} \tag{ 6 }$$

In this work, we adapted a shrunken U-Net as the deep learning-based image updater. The channel number of the first layer N_chl varies according to the memory footprint of each update block. With the increase in the depth, the number of layers in the encoder part increased to 2N_chl, 4N_chl, and 8N_chl, and then decreased by a factor of 2 after each concatenation operation in the decoder part.

Attention gates (Oktay et al 2018) were deployed before the concatenation operations in the shrunk UNet. This implementation enables the image updater to capture salient features and learn to focus on image structures exhibiting a higher level of inconsistency with the stationary point. In the attention gate, the input from the skip connection is scaled by an attention coefficient. The attention coefficient is evaluated by combining the input itself with the gating signal from a coarser scale, followed by a ReLU layer, a linear layer, and a soft activation function.

In order to stabilize and speed up the training process, our network is trained under the supervision of the weighted sum of a set of intermediate loss functions. For the sake of simplicity, we rewrite RHS in (6) as ${{\rm{\Gamma }}}_{{\theta }_{k}}({{\boldsymbol{c}}}^{k-1})$. The final loss function is

$$\begin{eqnarray}&&\arg \,{\min }_{\{{\theta }_{0}\ldots {\theta }_{n}\}}\displaystyle \sum _{n=0}^{N-1}{{\rm{w}}}_{t,n}\cdot { \mathcal L }\left({{\boldsymbol{c}}}_{\mathrm{truth}},{{\rm{\Gamma }}}_{{\theta }_{n}}\circ {{\rm{\Gamma }}}_{{\theta }_{n-1}}...\circ {{\rm{\Gamma }}}_{{\theta }_{0}}({{\boldsymbol{c}}}_{\mathrm{init}})\right),\end{eqnarray} \tag{ 7 }$$

where w_t,n ∈ [0, 1] denotes the scalar that controls the back-propagation weight, N denotes the number of update blocks, t denotes the index of the current training epoch, and {θ₀, θ₁} denotes the set of trainable parameters in all update blocks. ◦denotes the function composition, and ${\mathfrak{L}}$ denotes the customized loss function for the dual-energy BVM model

$$\begin{eqnarray}&&{ \mathcal L }({{\boldsymbol{c}}}_{\mathrm{truth}},{{\boldsymbol{c}}}_{\mathrm{est}})=\displaystyle \sum _{i}\parallel {c}_{i,\mathrm{truth}}-{c}_{i,\mathrm{est}}{\parallel }_{2}^{2}+\displaystyle \sum _{E}{\psi }_{L}(E)\parallel \displaystyle \sum _{i}{\mu }_{i}(E){c}_{i,\mathrm{truth}}-\displaystyle \sum _{i}{\mu }_{i}(E){c}_{i,\mathrm{est}}{\parallel }_{2}^{2},\end{eqnarray} \tag{ 8 }$$

where ∥ · ∥₂ is the ℓ₂ norm, and ψ_L (E) denotes the lower-energy spectrum at energy E. The second physics-based loss term is introduced to account for the correlation between c₁ and c₂.

In this study, we used dual-energy alternating minimization (DEAM) to demonstrate the performance of MB-DECTNet. DEAM is a jointly statistical image reconstruction technique that maximizes the log-likelihood of Poisson random variables with realizations d and estimated mean polychromatic sinogram g ( c ) = {g_L ( c ), g_H ( c )} (O'Sullivan and Benac 2007).

Maximizing the Poisson likelihood is equivalent to minimizing the I-divergence

$$\begin{eqnarray}&&{ \mathcal I }({\boldsymbol{d}}\parallel {\boldsymbol{g}}({\boldsymbol{c}}))=\displaystyle \sum _{j}\displaystyle \sum _{y}{d}_{j}(y)\mathrm{log}\displaystyle \frac{{d}_{j}(y)}{{g}_{j}(y:{\boldsymbol{c}})}-{d}_{j}(y)+{g}_{j}(y:{\boldsymbol{c}}),\end{eqnarray} \tag{ 9 }$$

where j = {L, H}, and g_j (y: c ) is given by the polychromatic forward model as

$$\begin{eqnarray}&&{g}_{j}(y:{\boldsymbol{c}})=\int {I}_{0,j}(E)\exp \left(-\int h(x,y)\displaystyle \sum _{i}{c}_{i}(x){\mu }_{i}(E){dx}\right){dE}+{\xi }_{j}(y),\end{eqnarray} \tag{ 10 }$$

where $y\in {{\mathbb{R}}}^{3}$ denotes the location in the measurement space, h(x, y) denotes the system operator, and I₀(y, E) denotes the mean photon count in the absence of the object, including the spectrum and bowtie filter. In the presence of substantial scatter, a scatter term ξ_j (y) may be introduced. Derivation of an ideal scatter model remains an active research topic (Liu et al 2021, Maier et al 2019, Medrano et al 2022).

Decoupling the I-divergence in energy, spatial, and basis domain yields an analytical solution to the surrogate function (refer to the appendix for derivation and step size Z_i (x)):

$$\begin{eqnarray}&&{c}_{i}^{k+1}(x)={c}_{i}^{k}(x)-\displaystyle \frac{1}{{Z}_{i}(x)}\mathrm{log}\displaystyle \frac{{\sum }_{j}{\sum }_{y}{\sum }_{E}{\mu }_{i}(E)h(x,y){p}_{j}^{k}(y,E)}{{\sum }_{j}{\sum }_{y}{\sum }_{E}{\mu }_{i}(E)h(x,y){q}_{j}^{k}(y,E)},\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&{q}_{j}^{k}(y,E)={I}_{0,j}\exp \left(-\displaystyle \sum _{x}h(x,y)\displaystyle \sum _{i}{\mu }_{i}(E){c}_{i}^{k-1}(x)\right)\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{p}_{j}^{k}(y,E)={q}_{j}^{k}(y,E)\displaystyle \frac{{d}_{j}(y)}{{\sum }_{E^{\prime} }{q}_{j}^{k}(y,E^{\prime} )}.\end{eqnarray} \tag{ 13 }$$

Therefore, for DEAM, the DC layer in MB-DECTNet is given by the log ratio between the forward projection of the attenuation-weighted measured sinogram and the forward projection of the attenuation-weighted mean sinogram that are summed over measurement index j and energy E, and weighted by LACs of the selected materials, formulated as

$$\begin{eqnarray}&&{\mathscr{DC}}({{\boldsymbol{c}}}^{k-1},d)=\left[\begin{array}{c}-\displaystyle \frac{1}{{Z}_{1}(x)}\mathrm{log}\displaystyle \frac{{\sum }_{j}{\sum }_{y}{\sum }_{E}{\mu }_{1}(E)h(x,y){p}_{j}^{k}(y,E)}{{\sum }_{j}{\sum }_{y}{\sum }_{E}{\mu }_{1}(E)h(x,y){q}_{j}^{k}(y,E)}\\ -\displaystyle \frac{1}{{Z}_{2}(x)}\mathrm{log}\displaystyle \frac{{\sum }_{j}{\sum }_{y}{\sum }_{E}{\mu }_{2}(E)h(x,y){p}_{j}^{k}(y,E)}{{\sum }_{j}{\sum }_{y}{\sum }_{E}{\mu }_{2}(E)h(x,y){q}_{j}^{k}(y,E)}\end{array}\right].\end{eqnarray} \tag{ 14 }$$

The GPU memory footprint is a critical design issue for 3D unrolled networks due to the large size of the image and multiple unrolling blocks. Besides the shrunken UNet mentioned in the previous subsection, we partition the image volume into multiple small stacks consisting of 8 slices to reduce memory consumption. The transmission data corresponding to each image stack is also truncated accordingly. For the geometry used in this work, 8 slices correspond exactly to a gantry rotation of the scanner. While the term 'slice' typically has a different meaning in CT reconstruction, in the context of this paper, it refers to a collection of voxels that share the same z-position.

We used 4×NVIDIA Tesla V100 to train the model. Each V100 has 32 GB memory, which can only accommodate a single instance of the training data. Therefore, we employed the gradient accumulation technique to compensate for the small batch size used in our training process. Instead of updating the model parameters after each singleton batch, we accumulated the gradients across multiple batches and updated the model every four samples to simulate a training process with a batch size of 4. Moreover, group normalization (Wu and He 2018) is utilized rather than the widely-used batch normalization in order to increase the accuracy of statistics estimation. These two techniques help reduce the memory requirements during training and improve the speed and stability of training, especially when GPU memory is limited.

In helical CT reconstruction, it is crucial to account for the margin effect, especially when a thin slice is reconstructed. As depicted in figure 4, the image slices corresponding to a single gantry rotation (highlighted in orange) are determined by the central row of the first and last views, where ${{ \mathcal Z }}_{d}$ represents the collimation width, and ${{ \mathcal Z }}_{\mathrm{feed}}$ denotes the bed travel distance along the z-direction per rotation. However, light orange segments outside the ROI will still influence the data acquired, and vice versa. Therefore, if the DC layer only evaluates the update direction based on the brighter orange area, the update direction will be much larger than expected to compensate for the absence of the light orange area. Thus, it is essential to pad the input image volume along the z-direction to ensure accurate results from the DC layer. Note that the diagram assumes parallel x-ray beams due to the large ratio between the source-to-detector distance and ${{ \mathcal Z }}_{d}$, which does not accurately represent the actual geometry in our implementations.

Zoom In Zoom Out Reset image size

Furthermore, each slice is not reconstructed from the same number of views. Let z be the distance between the centers of the first and an arbitrary slice. The number of gantry steps contributing to a slice can be calculated as

$$\begin{eqnarray}{n}_{\mathrm{vps}}(z)=\left\{\begin{array}{cc}{n}_{\mathrm{vpr}}\left(\frac{{{ \mathcal Z }}_{d}}{2{{ \mathcal Z }}_{\mathrm{feed}}}+\frac{z}{{{ \mathcal Z }}_{\mathrm{feed}}}\right) & z\lt {{ \mathcal Z }}_{\mathrm{feed}}-\frac{{{ \mathcal Z }}_{d}}{2}\\ {n}_{\mathrm{vpr}}\left(\frac{{{ \mathcal Z }}_{d}}{2{{ \mathcal Z }}_{\mathrm{feed}}}+\frac{{{ \mathcal Z }}_{\mathrm{feed}}-1-z}{{{ \mathcal Z }}_{\mathrm{feed}}}\right) & z\gt \frac{{{ \mathcal Z }}_{d}}{2}-1\\ {n}_{\mathrm{vpr}} & \mathrm{else},\end{array}\right.\end{eqnarray} \tag{ 15 }$$

where n_vpr denotes the number of views per rotation.

The reference gradient for peripheral slices will be updated based on fewer views compared to the gradient for central slices. This implies that the peripheral slices of the reconstructed image may exhibit higher error margins than the central slices. To mitigate this issue, only the central slices of the MB-DECTNet output were utilized to evaluate the performance of the proposed method, and the contribution of peripheral slices to the loss function was down-weighted. Ideally, we want to extract the slices that correspond to as many as possible views. In the clinical setup, ${{ \mathcal Z }}_{\mathrm{feed}}$ is 8.304 mm and ${{ \mathcal Z }}_{d}$ is 12 mm. Therefore, the slice index must satisfy $2.2\leqslant \tilde{z}\leqslant 4.8$ to ensure the image slice is reconstructed by a full rotation of views.

Figure 5 shows the plots of estimation error by slices along the z-direction. A periodic pattern can be seen in the MSE plot of output with the margin effect. The model-based network with the proposed strategy significantly outperforms the same network that is trained and tested based on all the slices of the samples.

Zoom In Zoom Out Reset image size

During the training process, the weight of the neural network is updated based on the gradient of the loss function with respect to the corresponding parameter. Following the chain rule, this gradient can be propagated from the loss function to the foremost layer, a technique commonly referred to as back-propagation.

The Jacobian of the DC layer can be expressed as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial {c}_{i^{\prime} }(x^{\prime} )}{\mathscr{DC}}{\left({\boldsymbol{c}},d\right)}_{(i,x)}={\left[\displaystyle \sum _{j}\displaystyle \sum _{y}\displaystyle \sum _{E}\displaystyle \frac{\partial {\mathscr{DC}}({\boldsymbol{c}},d)}{\partial {q}_{j}(y,E)}\displaystyle \frac{\partial {q}_{j}(y,E)}{\partial {\boldsymbol{c}}}\right]}_{(x,i),(x^{\prime} ,i^{\prime} )}\\ \quad =-\displaystyle \frac{1}{{Z}_{i}(x)}\displaystyle \sum _{j^{\prime} }\displaystyle \sum _{y^{\prime} }\displaystyle \sum _{E^{\prime} }\displaystyle \frac{\partial }{\partial {q}_{j^{\prime} }(y^{\prime} ,E^{\prime} )}\mathrm{log}\\ \quad \cdot \displaystyle \frac{{\sum }_{j}{\sum }_{y}{\sum }_{E}{\mu }_{i}(E)h(x,y){q}_{j}(y,E)\tfrac{{d}_{j}(y)}{{\sum }_{E^{\prime\prime} }q(y,E^{\prime\prime} )}}{{\sum }_{j}{\sum }_{y}{\sum }_{E}{\mu }_{i}(E)h(x,y){q}_{j}(y,E)}\displaystyle \frac{\partial {q}_{j^{\prime} }(y^{\prime} ,E^{\prime} )}{\partial {c}_{i^{\prime} }(x^{\prime} )}\\ \quad =\displaystyle \frac{1}{{Z}_{i}(x){\left({\hat{b}}_{i}(x)\right)}^{2}}\displaystyle \sum _{j}\displaystyle \sum _{y}\displaystyle \sum _{E}{\mu }_{i}(E)h(x,y)\left({A}_{j}(y,E){\hat{b}}_{i}(x)-{\tilde{b}}_{i}(x)\right){q}_{j}(y,E){\mu }_{i^{\prime} }(E)h(x^{\prime} ,y)\\ \quad \approx -\displaystyle \frac{1}{{Z}_{i}(x){\hat{b}}_{i}(x)}\displaystyle \sum _{y}h(x,y)\left(\displaystyle \sum _{j}\displaystyle \sum _{E}{\mu }_{i}(E)\displaystyle \frac{{q}_{j}{\left(y,E\right)}^{2}}{{d}_{j}(y)}{\mu }_{i^{\prime} }(E)\right)h(x^{\prime} ,y)\end{array}\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&{\hat{b}}_{i}(x)=\displaystyle \sum _{j}\displaystyle \sum _{y}\displaystyle \sum _{E}{\mu }_{i}(E)h(x,y){q}_{j}(y,E)\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\tilde{b}}_{i}(x)=\displaystyle \sum _{j}\displaystyle \sum _{y}\displaystyle \sum _{E}{\mu }_{i}(E)h(x,y){p}_{j}(y,E)\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{A}_{j}(y,E)=\displaystyle \frac{{d}_{j}(y)\left({\sum }_{E^{\prime} }{q}_{j}(y,E^{\prime} )-{q}_{j}(y,E)\right)}{{\left({\sum }_{E^{\prime} }{q}_{j}(y,E^{\prime} )\right)}^{2}}.\end{eqnarray} \tag{ 19 }$$

The last expression in (16) follows by assuming g_j (y) ≈ d_j (y) and q_j (y, E) ≈ p_j (y, E). Note that q_j (y, E) and ${\hat{b}}_{i}(x)$ have been computed during forward propagation, so only 2 forward-projections and 2 back-projections are required for an update block per iteration. Figure 6 shows an example of the gradients before and after the DC layer.

Zoom In Zoom Out Reset image size

During our experimentation, we observed that MB-DECTNet utilizing the exact DC gradient demonstrated a higher convergence rate than the MB-DECTNet implementation that relied on the zero DC gradient. However, it is noteworthy that both networks reached comparable levels of performance upon convergence.

In our study, we have performed a quantitative comparison between three MB-DECTNet variants and two benchmark methods: the prevalent image-domain decomposition (IDD) (Yan et al 2000) and a standalone deep learning technique (Pure Net). To comprehend the interplay between the model performance and the physical information incorporated, we introduced three MB-DECTNet variants with 1, 2, and 4 update blocks, referred to as MB-DECTNet-1, MB-DECTNet-2, and MB-DECTNet-4. Each variant comprises 1, 2, or 4 DC layers, representing increasing levels of synergy with the physical model. To maintain a balanced comparison, the channel count, N_chl, in the MB-DECTNets was fine-tuned, ensuring roughly equivalent learnable parameters across models. For the Pure Net method, we employed the widely used U-Net structure (Ronneberger et al 2015) with both input and output channel numbers set to 2, representing a model devoid of physical knowledge.

To quantitatively compare the accuracy of various models, we plotted histograms of the estimated LAC errors at 30, 60, and 150 keV, where the error was calculated as

$$\begin{eqnarray}&&\displaystyle \frac{{\mu }_{{est}}(x,E)-{\mu }_{{GT}}(x,E)}{{\mu }_{{GT}}(x,E)}\times 100 \% .\end{eqnarray} \tag{ 20 }$$

To enhance training efficiency and promote training stability, we employed a pre-training strategy in this study. Specifically, we first trained the first block of the MB-DECTNet independently for 100 epochs using the ADAM optimizer with a starting learning rate of 10⁻⁴ and a decay rate of 0.8 every 10 epochs. During this stage, the input remained constant throughout epochs, and the reference update from the DC layer was reusable, resulting in a relatively quick training process. We then propagated the pre-trained weights to other blocks to initialize the complete unrolling network. Finally, we trained the entire network end-to-end for an additional 50 epochs, minimizing the loss function in equation (7) with a learning rate of 10⁻⁴ and a decay rate of 0.8 every 5 epochs. The loss weights w_t,n employed for MB-DECTNet with varying numbers of blocks are shown in table 1.

In this study, we utilized images reconstructed by DEAM as the ground truth. The inputs are the image domain decomposition method with an iterative scheme (Yan et al 2000). The dual-energy measurements were obtained at 90 and 140 kVp using a Philips Brilliance Big Bore CT scanner located at the Alvin J Siteman Cancer Center in Washington University School of Medicine, following a helical protocol. In order to obtain the raw data necessary for our approach, beam hardening corrections in these measurements have been removed, which requires special agreement with the vendor. The number of patients available within the human subject protocol was limited. The detector array contains 16 rows and 816 channels, with a row spacing of 0.75 mm and channel spacing of 0.001 189 71 radians. Each rotation contains 1320 gantry steps with a z-direction feed of 8.304 mm.

These measurements were subsequently reconstructed through the application of a penalized DEAM within a motion-compensation framework (Ge et al 2023), with a sufficient number of iterations. Basis materials are 23% aqueous solution of CaCl₂ and polystyrene (Williamson et al 2006). The size of the reconstructed head-neck image is 610 × 610 × 340, and the image resolution is 1 × 1 × 1.034 mm³. The penalty term is written as

$$\begin{eqnarray}&&{ \mathcal R }({\boldsymbol{c}})=\displaystyle \sum _{x}\displaystyle \sum _{i}\displaystyle \sum _{x^{\prime} \in { \mathcal N }(x)}{w}_{R}(x,x^{\prime} )\psi \left({c}_{i}(x)-{c}_{i}(x^{\prime} )\right),\end{eqnarray} \tag{ 21 }$$

where ${ \mathcal N }(x)$ denotes the neighborhood of x, w_R denotes the weight addressing the distance between x and $x^{\prime}$. The derivation of the image update with the penalty of the form (21) is shown in the appendix. The potential function ψ is an element-wise convex symmetric smooth function:

$$\begin{eqnarray}&&\psi (t)={\delta }^{2}\left(\displaystyle \frac{| t| }{\delta }-\mathrm{ln}\left(1+\displaystyle \frac{| t| }{\delta }\right)\right),\end{eqnarray} \tag{ 22 }$$

where δ controls the transition between approximate ℓ₂-norm and approximate ℓ₁-norm regions.

Each DEAM ground-truth image set took 20 h to converge on 4×NVIDIA V100 32 GB GPUs (400 iterations with 33 ordered subsets, followed by an additional 1000 iterations without ordered subsets).

The training set contains the images and measurements for 4 head-neck patients, one lung patient, and a customized head-sized phantom. The number of epochs was determined through a validation set containing out-of-sample lung patient data. Our test set contains 320 samples extracted from an out-of-sample head-neck patient scan.

Scans of patients treated with proton or photon therapy are acquired under IRB study NCT03403361. The scanned region of the head-neck patient extended from the top of the head to the base of the skull to model a typical central-nervous-system scan, and the scanned region of the lung patient extended from the base of the skull to the top of the pelvis to represent a typical lung scan.

The customized phantom has a cylindrical acrylic shell with a diameter of 215 mm and length of 20 mm, which is filled with water. Eight inserts are plastic bottles with a diameter of 31 mm and filled with liquid samples, including ethanol, propanol, butanol, and K₂HPO₄ solutions with different concentrations.

Each sample is the truncated portion of the spiral sinogram with dimensions 16 × 816 × 1320, which corresponds to a rotation of the scanner gantry. The image size of a sample was 610 × 610 × 8, with an image resolution of 1 × 1 × 1.034 mm³. In total, the training set contained 1280 samples.

Figure 7 illustrates the estimated component weights and virtual monoenergetic images (VMIs) of the selected slice in the test set. The IDD result has the highest noise level. The pure deep-learning model offers some noise reduction but exhibits an underestimation of the soft tissue in VMIs at 60 and 150 keV. MB-DECTNet with 1 update block (MB-DECTNet-1) offers better denoising capabilities compared to the pure network, but it tends to underestimate the LAC of bony tissues, particularly near the boundary between the skull and brain tissue. MB-DECTNet with 2 and 4 update blocks (MB-DECTNet-2 and MB-DECTNet-4) deliver the best visual performance and demonstrate the potential for reducing noise while preserving structural information. In inference mode, MB-DECTNet-2 and MB-DECTNet-4 generate a 610 × 610 × 340 image in roughly 352 and 700 s, which are 205 and 103 times shorter than the elapsed time of DEAM, respectively.

Zoom In Zoom Out Reset image size

The histograms of percentage errors are depicted in figure 8. The region of interest was limited to all soft and bony tissues (LAC threshold range [0.01, 0.05] at 60 keV) to avoid extremely small values being used as the denominator. The IDD estimation errors at 30 keV and 150 keV span a wide range, indicating a large variance in the corresponding VMIs. Additionally, the peak at −48% in the 30 keV plot suggests a significant estimation bias in IDD. All deep-learning models demonstrate the ability to correct bias and reduce noise in the initial image. As the pure network does not utilize data-consistency information, its performance is inferior to that of model-based networks. Among the model-based learning methods, MB-DECTNet-4 exhibits the most near-zero errors at the three selected energies, indicating that the proposed network with the most update blocks achieves the best performance.

Zoom In Zoom Out Reset image size

Table 2 also reflects the same trend, where MB-DECTNet-4 emerges as the top performer across most evaluated metrics. MB-DECTNet-4 attains the lowest MAE and highest PSNR at 60 and 150 keV. The percentage bias of MB-DECTNet-2 at 150 keV is lower than that of MB-DECTNet-4, and their biases are comparable at 60 keV. Since the component weight estimates the LAC over a broad energy range, metrics at 60 and 150 keV may be insufficient to fully demonstrate the performance of the models. Hence, we utilized the original objective function of DECT SIR as an alternate performance measure that projects the estimates to the sinogram domain and directly compares them with the measurement. Among the deep-learning models, MB-DECTNet-4 achieves the lowest sinogram-domain objective function value.

A noticeable disparity can be observed between the objective function of the reference image and the MB-DECTNet-4 estimate. However, this discrepancy does not necessarily imply suboptimal performance by MB-DECTNet-4. As commonly acknowledged, achieving a neural network output with precisely zero background is challenging. In the calculation of the objective function, the nonzero background values for air will be projected and accumulated in the sinogram domain, and undergo nonlinear amplification through the exponentiation term in the forward model. These background errors can be mitigated by post-processing.

The previous subsection demonstrated the potential of MB-DECTNet to substitute for model-based statistical algorithms for DECT reconstruction. However, it can be seen that due to the limited number of accessible unpreprocessed clinical measurements in the training set, the estimated image by MB-DECTNet did not meet the situation where highly accurate reconstruction results with the subpercentage error were desired.

In order to further improve the result from the neural network, we initialized the traditional DECT SIR using the deep learning-based model and compared the performance with the traditional IDD method on the same task. Specifically, we used IDD, Pure Net, and MB-DECTNet-4 as the initial condition of DEAM with 33 ordered subsets, and plotted the objective function values of DEAM versus iterations. The plot is shown in figure 9. As it is shown, DEAM initialized by MB-DECTNet-4 is close to convergence after 80 iterations, while DEAM initialized by IDD takes 100 iterations to reach the objective function value of DEAM initialized by MB-DECTNet-4 at the 20th iteration.

Zoom In Zoom Out Reset image size

Figure 10 shows the percentage MAE of the estimation against the reference image versus iteration, which provides another view of the convergence analysis. As we can see, MAEs of DEAM initialized by IDD start from 35%, 3.5%, and 7.5%, and gradually decrease to 10%, 1.5%, 1.4% at 30, 60, and 150 keV, respectively, while DEAM initialized by MB-DECTNet-4 reaches 1% MAE at the 8th iteration for 60 and 150 keV. The errors at 30 keV are significantly larger than the errors at 60 and 150 keV. This disparity arises primarily because very few photons with energies near 30 keV are detected. Additionally, the low- and high-energy spectra greatly down-weight the contribution of the modeled LAC at extremely low energies.

Zoom In Zoom Out Reset image size