Accelerating direct material decomposition via diffusion probabilistic model for Sparse-view spectral computed tomography

Abstract

Accurate material decomposition constitutes the foundation of Spectral Computed Tomography (Spectral CT) applications across diverse domains. Nevertheless, conventional model-based material decomposition methods face significant limitations including sparse-view sampling artifacts, slow convergence rates, noise amplification, and inherent ill-posedness—challenges that are particularly pronounced in geometrically inconsistent imaging. To overcome these constraints, we propose an unsupervised deep learning framework that synergistically optimizes virtual monochromatic images (VMIs) through the probabilistic diffusion model for direct material decomposition in sparse-view spectral CT. The proposed methodology introduces VMIs as critical differentiation enhancers for polychromatic projections, effectively addressing convergence limitations in iterative reconstruction algorithms. By incorporating probabilistic diffusion priors into the optimization process, we achieve superior refinement of material-specific representations. Our framework systematically enforces dual constraint: 1) data fidelity term ensuring measurement consistency, and 2) probabilistic regularization suppressing unwanted structures, thereby guaranteeing anatomically plausible material image reconstruction. Comprehensive validation on preclinical data demonstrates that our method achieves a 10 dB improvement in the peak-signal-to-noise ratio (PSNR) and a 4.31% increase in structural similarity (SSIM) for soft-tissue reconstructions compared to the optimal comparison algorithm with 90 projections. Experimental results confirm the algorithm's robustness under challenging conditions, maintaining reconstruction fidelity even with geometric inconsistency and sparse sampling.

Keywords

spectral CT material decomposition diffusion probabilistic models virtual monochromatic image

Introduction

Spectral Computed Tomography (Spectral CT) leverages multi-energy X-ray spectra to enable material identification, a critical capability for applications such as bone density measurement,¹ liver fibrosis assessment,² abdominal imaging,³ and breast cancer risk assessment.⁴ Material decomposition algorithms in spectral CT are broadly classified into direct and indirect methods. Direct material decomposition methods (DMDM),⁵ contrast with indirect methods, which encompass image-domain and projection-domain techniques. Image-domain methods involve reconstructing multi-channel spectral CT images followed by material decomposition.⁶ These methods rely on simplified linear imaging model that inadequately align with actual polychromatic scanning condition. Projection-domain methods decompose projection data into material-specific projections prior to image reconstruction,⁷ leveraging spectral and polychromatic models for nonlinear decomposition.⁸ While continuously refined,^9,10 these methods require precise spectral calibration and strict geometrical consistency across energy channels, posing practical implementation challenges.

DMDM directly reconstruct material images from projection data, accurately modeling spectral CT physics to mitigate beam-hardening artifacts and minimize cascaded information loss. In 2015, Zhao et al. proposed the extended algebraic reconstruction technique (EART) for dual energy CT, achieving promising results at the cost of low computational efficiency.¹¹ Subsequent improved algorithms are developed, such as the extended simultaneous algebraic reconstruction technique (ESART)¹² and virtual monochromatic images (VMIs)-guided iterative methods.¹³ To improve reconstruction quality, DMDM incorporate regularization terms exploiting intrinsic image properties, such as the total variation (TV)¹⁴ and block matching strategies.¹⁵ Chang et al.¹⁶ unified basis material reconstruction and spectral estimation into a TV-regularized non-convex optimization. However, traditional regularization designs remain empirically driven, limited in characterizing noise correlations and insufficient in describing the mapping relationships between spectral CT images and basis materials.

Deep convolutional neural networks (CNNs) have revolutionized spectral CT material decomposition. Zhang et al.¹⁷ pioneered an end-to-end butterfly network architecture, later enhanced through loss functions optimization¹⁸ and the generative adversarial network (GAN) structure for multi-material decomposition (MMD) in.¹⁹ Moreover, the decomposition of three materials under dual-energy spectra have been realized by combining constraint condition,²⁰ while projection-domain decomposition and DMDM now incorporate sparse-view reconstruction and end-to-end mappings.²¹ Li et al.²² proposed interpretable decomposition iterative network for low-dose CT denoising. Despite overcoming material recognition via hierarchical features, supervised CNNs depend critically on high-quality labels—a significant limitation in medical image. Moreover, insufficient integration of physical imaging mechanisms often compromises data consistency and reconstruction fidelity.

Diffusion models (DMs) have gained prominence for their ability to learn image probability distributions and generate high-fidelity samples without label.^23,24 Comprising a forward diffusion process and a backward sampling process,²⁵ DM excel in inverse problem solving under noisy conditions, exemplified by diffusion posterior sampling (DPS)²⁶ and manifold constraint gradient (MCG).²⁷ Liu et al.²⁸ applied DM to sparse spectral CT material decomposition, achieving accurate contrast agent distribution recovery. To address the slow sampling of DM, various approaches have been proposed, including jump sampling strategies,²⁹ denoising diffusion implicit models (DDIM),³⁰ and medical segmentation optimizations³¹ — improve efficiency. DM-based spectral CT decomposition remains challenged by computational complexity, slow convergence and poor adaptability to material-specific nonlinearities, particularly under sparse sampling conditions.

This work proposes an unsupervised spectral CT direct material decomposition framework synergizing data-driven DM with physical imaging models. To address the slow convergence problem caused by geometric inconsistencies under fast kVp switching imaging, VMIs are employed to augment angular separation between polychromatic projection curves. Critically, DM iteratively refines VMIs via probability distribution learning and superior feature representation. To enhance the iteration efficiency, we introduce a truncated reverse diffusion method, allowing the DM to participate only in the latter half of the VMIs-guided iterative process. By bypassing early high-noise phase detrimental to CT image generation, this approach effectively mitigates the severe ill-posedness of spectral CT material decomposition, primarily caused by high correlation across different spectral channels, material attenuation curve similarities, and sparse-view protocols.

The structure of this paper is as follows: the second section introduces the background, primarily focusing on physics model for spectral CT and DM. The third section introduces the proposed method and the solution of the model. The fourth section is the experiment and analysis of the experimental results. The paper concludes with a brief discussion and conclusion.

Background

Physics model for spectral CT

This work takes into account the high and low energy projections to direct reconstruct two basis material maps. The dual-energy CT projection model can be expressed as follows:

\begin{aligned} p_{k, L} (b_{1}, b_{2}) = & - \ln \sum_{m = 1}^{M_{k}} S_{k, m} \exp (- [μ_{1} (e) A b_{1} + μ_{2} (e) A b_{2}]), \\ k = h, l \end{aligned}

(1)

where

k = h, l

denotes the spectrum index for high/low energy spectrum,

L \in L_{k}

represents the X-ray path,

L_{k}

indicates the index set of the X-ray paths corresponding to the k -th X-ray spectrum, and

p_{k, L} (b_{1}, b_{1})

is the polychromatic projection of the k -th X-ray spectrum along the L -th X-ray path based on the basis material images

b = (b_{1}, b_{2})

S_{k, m}

denotes the normalized X-ray spectrum at the k -th energy level satisfying

\sum_{m} S_{k, m} = 1

μ_{t} (e) (t = 1, 2)

represent the linear attenuation coefficient (LAC) of the t -th material at energy e.

A \in R^{M \times N}

represents the projection matrix, where M represents the number of detector elements, and N denotes the number of pixels in the reconstructed image. Assuming the current estimation is

b^{(n)} = (b_{1}^{(n)}; b_{2}^{(n)})

, and the first-order Taylor expansion of the right-hand side of (1) at this point can be expressed as follows,

p_{k, L} \approx p_{k, L}^{(n)} + \frac{Φ_{k, L}^{(n)}}{q_{k, L}^{(n)}} A (b_{1} - b_{1}^{(n)}) + \frac{Θ_{k, L}^{(n)}}{q_{k, L}^{(n)}} A (b_{2} - b_{2}^{(n)})

(2)

For the specific expression of the coefficients, please refer to appendix I. Eq. (2) can be briefly written as the following linear system,

p \approx H b, H = (\begin{matrix} \frac{Φ_{h, L}^{(n)}}{q_{h, L}^{(n)}} & \frac{Θ_{h, L}^{(n)}}{q_{h, L}^{(n)}} \\ \frac{Φ_{l, L}^{(n)}}{q_{l, L}^{(n)}} & \frac{Θ_{l, L}^{(n)}}{q_{l, L}^{(n)}} \end{matrix}) A

(3)

For the specific expression of $p$ , please refer to appendix I. Subsequently, the material image can be updated via error back-projection within the ESART algorithm.¹²

Diffusion probability model

Generative models are a class of models capable of learning data distributions and generating new samples, replicating, and even creating new samples similar to the training data.³² Diffusion models,^33,34 as one of the most popular probabilistic generative models, learn complex data distribution by estimating the gradient of the prior probability density $\nabla_{x} \log p_{t} (x)$ . Specifically, the forward diffusion process of diffusion model injects the noise into the data continuously until it becomes white noise, and the reverse diffusion process removes the noise progressively to recover clean data. The forward processes can be described by stochastic differential equation (SDE) as follows,

d x = f (x, t) d t + g (t) d w,

(4)

where

f (\cdot, t) : R^{n} \to R^{n}, g (t) : R \to R

are the drift coefficient and the diffusion coefficient, respectively, and vector

w

is the standard Brownian motion. The reverse diffusion process can be formulated as follows,

d x = [f (x, t) - g^{2} (t) \nabla_{x_{t}} \log p_{t} (x)] d t + g (t) d \bar{w},

(5)

where

d t

denotes an infinitesimal time step and the vector

\bar{w}

represents the standard Brownian motion. Following the approach in reference,²⁴ the drift coefficient and diffusion coefficient are set to

f (\cdot, t) = 0, g (t) = \sqrt{d [σ^{2} (t)] / d t}

, respectively, where

σ (t) > 0

is a monotonically increasing function with t. Accordingly, Eq.(5) can be reformulated as follows,

d x = \frac{d [σ^{2} (t)]}{d t} \nabla_{x_{t}} \log p_{t} (x) + \sqrt{\frac{d [σ^{2} (t)]}{d t}} d \bar{w} .

(6)

The score $\nabla_{x_{t}} \log p_{t} (x)$ is estimated by a neural network $s_{θ}$ parameterized by $θ$ , which will be trained with score matching,³⁵

θ * = \underset{θ}{\arg min} E_{x_{0}} E_{x_{t} | x_{0}} [‖ s_{θ} (x_{t}, t) - \nabla_{x_{t}} \log p (x_{t} | x_{0}) ‖_{2}^{2}],

(7)

where

x_{t} \sim p (x_{t} | x_{0})

. When the training iteration converges, the estimation of score function

s_{θ *} (x, t)

is obtained, where

θ *

is the optimal parameter.

Method

Direct material decomposition based on probability diffusion priors

The solution of the following inverse problem can be examined from a probabilistic point of view,

p = H b + ε .

(8)

The problem of reconstruction material images based on measurement projections $p$ can be approached from a statistical modeling perspective, which means estimating the material image $b$ given the measurements $p$ from the observation space. Specifically, by maximizing the posterior distribution $p (b | p)$ , we can obtain the maximum a posteriori (MAP) estimate of the material image, that is,

\hat{b} = \underset{b}{\arg max} \log p (b | p) .

(9)

According to Bayes’ rule, the MAP estimate can be rewritten as follows,

\hat{b} = \underset{b}{\arg max} \log p (b | p) \propto \underset{b}{\arg max} \log p (p | b) + \log p (b),

(10)

where

\log p (p | b)

denotes the “likelihood” term and

\log p (b)

delivers the prior. Under the assumption that

ε = H b - p

follows zero-mean Gaussian noise, we have

\log p (p | b)

\propto - ‖ H b - p ‖_{2}^{2}

.Then Eq. (10) can be rewritten in the following,

b * = \underset{b}{\arg min} \frac{1}{2} ‖ H b - p ‖_{2}^{2} + η (- \log p (b)),

(11)

where

η > 0

is the regularization parameter, and the term

- \log p (b)

denotes the regularization prior, which can leverage the inverse diffusion sampling process. During the minimization of Eq. (11), the gradient of the logarithm of the unknown data distribution

- \log p (b)

can be estimated with the trained neural network

s_{θ *} (b, t)

, satisfying

‖ s_{θ *} (b, t) - \nabla_{b} \log p (b) ‖ \leq τ^{2}

with

τ > 0

is the threshold value. The architectural design of

s_{θ} (b, t)

critically influences the performance of the DMDM framework. The reconstruction algorithm achieves improved convergence stability through the synergistic interaction between the data fidelity term and the regularizers.

Our methodology integrates pre-trained probabilistic diffusion models as implicit probability regularizer within the Plug-and-Play (PnP) architecture³⁶ framework, combining the theoretical rigor of model-based approaches with the representational power of deep learning techniques. This unified approach enables simultaneous optimization of generative models and direct material decomposition, establishing a novel paradigm that synergistically integrates data-driven and model-driven methodologies. The diffusion-based regularizer overcomes the representational limitations inherent in conventional regularization approaches, creating a new framework for addressing inverse problems in spectral CT material decomposition.

Direct material decomposition framework integrates VMIs guidance with diffusion probability priors

The polychromatic nature of broad-spectrum X-ray sources leads to spectral overlap in acquired projection. This inherent mixing of multi-energy photons introduces high cross-correlation among attenuation measurements at different energy bins, resulting in the system severely ill-posed system. Additionally, the LAC of different materials can be very similar. For example, in breast imaging, the LAC differences between adipose (fat) and glandular/fibrous tissues across different energy are small. Sparse sampling further exacerbates the insufficiency of the acquired projection data. Collectively, these factors increase the ill-posedness of the reconstruction problem. This causes the spectral attenuation cureves to lie extremely close to one another, leading to slow convergence of the ESART algorithm.¹²

To accelerate convergence, Zhang et al.¹³ proposed a VMIs-guided method. This method increases the angle between polychromatic projection curves by introducing monochromatic images, and thus accelerating the iterative process. Specifically, two distinct energies $E_{1}$ and $E_{2}$ are introduced, and the corresponding VMIs are denoted as $I_{E 1}, I_{E 2}$ . The material image reconstruction is transformed into a problem of reconstruction the VMIs $I_{E 1}, I_{E 2}$ as following

\nabla_{I_{E}} F (I_{E}, t) = \nabla_{I_{E}} ‖ H I_{E} - p ‖_{2}^{2} + \frac{η}{2} (s_{θ} (I_{E 1}, t) + s_{θ} (I_{E 2}, t)) .

(12)

These VMIs can be represented as linear combinations of the basis materials,

(\begin{matrix} I_{E 1} \\ I_{E 2} \end{matrix}) = (\begin{matrix} μ_{1} (E_{1}) & μ_{2} (E_{1}) \\ μ_{1} (E_{2}) & μ_{2} (E_{2}) \end{matrix}) (\begin{matrix} b_{1} \\ b_{1} \end{matrix}) .

(13)

Eq. (13) can be equivalently represented as the following,

(\begin{matrix} b_{1} \\ b_{1} \end{matrix}) = (\begin{matrix} c_{11} & c_{12} \\ c_{21} & c_{22} \end{matrix}) (\begin{matrix} I_{E 1} \\ I_{E 2} \end{matrix}),

(14)

where

c_{11}, c_{12}, c_{21}, c_{22}

denote the corresponding elements in the inverse of the coefficient matrix in (13).

The introduction of VMIs for material decomposition increases the angle between polychromatic attenuation curves, thereby accelerating the iterative convergence. More importantly, it facilitates the implementation of deep prior regularization using probability diffusion models. VMIs exhibit highly consistent distribution features with CT images, sharing similarities in anatomical structure edges and tissue contrast, and they also have the same grayscale scaling. In previous work, our research group successfully trained a probabilistic diffusion model using the publicly available LIDC-IDRI lung CT dataset.³⁷ Leveraging this, we can directly employ this pre-trained model without requiring retraining on the material distribution image dataset, providing substantial convenience for the smooth progress of our research. The flowchart of the proposed DMDM framework with virtual monochromatic guidance, incorporating a fused probability diffusion prior, is illustrated in Figure 1.

Figure 1.

Visual representation of the proposed method. Under inconsistent geometry, the material residual projections of the VMIs are first updated, followed by VMIs reconstruction using the ESART as defined in Eq.(17). Subsequently, denoising is performed on the updated VMIs via the probabilistic diffusion model. Finally, material decomposition is executed based on the optimized VMIs following Eq.(14) to obtain material images.

Solution of the direct material decomposition framework

The objection function based on VMIs guidance is rewritten as follows,

F (I_{E}, t) := \frac{1}{2} ‖ H^{'} I_{E} - p ‖_{2}^{2} - η \log p (I_{E}) \overset{Δ}{=} m (I_{E}) + h (I_{E}) .

(15)

The specific expression of $H^{'}$ is detailed in the appendix I. To maintain consistency with the index notation in the reverse diffusion process, the range of the iteration index i here is also from T to 0.

The iteration of the proximal gradient descent algorithm consist of two steps:

Step 1, Gradient descent step: Perform gradient updates on the fidelity term $m (I_{E})$ with step size $α$ .

I_{E}^{(i - 1)} = I_{E}^{(i)} - α \nabla ‖ H_{i}^{'} I_{E}^{(i)} - p ‖_{2}^{2} .

(16)

The solution process in Step 1 can utilize the approach of ESART,¹³ first updating the material residual projections, and then updating the VMIs. With the inconsistent geometry, the X-ray paths of the dual energy passing through the object to be scanned are no longer exactly the same, and the weights need to be recalculated to update the error. By the equation (25) in,¹³ the material projection residuals are obtained. To enable the iteration to converge rapidly, the entire projection dataset is divided into T parts, and the index set of the rays within each projection is $Γ_{t} (t = 1, 2, \dots, T)$ . Then, an ordered-subset accelerated version of the dual-energy ESART algorithm, namely the OS-ESART algorithm, is obtained. In this case, the iterative formula for the t-th operation in the n-th iteration of the OS-ESART algorithm denoted as follows,

\begin{aligned} I_{E 1}^{(n + 1)} = I_{E 1}^{(n)} + \frac{η}{\sum_{L \in Γ_{t}} A_{L}} \sum_{L \in Γ_{t}} A_{L} ((Δ F_{L, 1}^{^{'} (n + 1)}) / (\sum_{i = 0}^{I - 1} a_{L, i})), \\ I_{E 2}^{(n + 1)} = I_{E 2}^{(n)} + \frac{η}{\sum_{L \in Γ_{t}} A_{L}} \sum_{L \in Γ_{t}} A_{L} ((Δ G_{L, 2}^{^{'} (n + 1)}) / (\sum_{i = 0}^{I - 1} a_{L, i})) . \end{aligned}

(17)

where

A_{L} = (a_{L, 1}, \dots, a_{L, I})

denotes the LAC weights corresponding to the i-th pixel along the L-th X-ray,

Δ F_{L, 1}^{'}

and

Δ G_{L, 1}^{'}

are the material projection residuals in.¹³ For a more detailed derivation and analysis, please refer to Appendix I. To optimize the objective function (15) using the proximal operator step, we need to employ the reverse diffusion process of the DM as a first-order approximation of the proximal operator. We have the following proposition.

Proposition 1.

Assume $x_{i} = \sqrt{\bar{α}} x_{0} + \sqrt{1 - \bar{α}} ε,$ $ε \sim N (0, I)$ , be a noisy point generated from $x_{0} \sim p (x)$ through the forward diffusion $p (x)$ . Then the $x_{i - 1}$ obtained by reverse sampling from $x_{i}$ , i.e., $x_{i - 1} = x_{i} + ξ_{i} s_{θ *} (x_{i}, i) + \sqrt{2 ξ_{i}} ε_{i}, ε_{i} \sim N (0, I)$ , is a first-order approximation to the proximal operator of $h (x)$ at $x_{i}$ ,where ${ξ_{i} \in (0, 1)}_{i = 1}^{T}$ is a pre-defined variance schedule.

Remarks.

Reference³⁸ presented this proposition, and its explanation is similar to that in,³⁹ but neither provided an explicit proof. This paper offers a complete proof process for the proposition, and the details can be found in appendix II.

Step 2, Proximal operator step: Apply the proximal operator to the non-smooth term $h (I_{E})$ ,

I_{E}^{(i - 1)} = {prox}_{α η h} (I_{E}^{(i)}) .

(18)

The definition of the proximal operator is provided as follows, for a function f on $R^{n}$ and $λ > 0$ , define the proximal operator of f as

{prox}_{λ f} (x) := \underset{z \in R^{n}}{\arg min} {f (z) + \frac{1}{2 λ} ‖ z - x ‖^{2}} .

(19)

According to Proposition 1, the solution in Step 2 can be approximated by the sampling process of the DM, as represented below,

I_{E}^{(i - 1)} = I_{E}^{(i)} + ξ_{i} s_{θ *} (I_{E}^{(i)}, i) + \sqrt{2 ξ_{i}} ε_{i}, ε_{i} \sim N (0, I) .

(20)

After the above iterative process, the optimization $x *$ can be obtained, that is, the VMIs $I_{E 1} *, I_{E 2} *$ are derived, and we obtain the material image $(b_{1}, b_{2})$ with Eq.(14).

Summarizing the optimization process of the algorithm is listed in Algorithm 1.

Convergence analysis of the direct material decomposition framework

In this section, we present conditions for the convergence of DMDM to ensure the algorithm is convergent. First, under appropriate assumption, the monotonic decrease of the objective function is established.

Proposition 2.

Assume that: (1) the gradient of $m (x) = ‖ H^{'} x - y ‖_{2}^{2}$ is Lipschitz continuous, and $p (x)$ follows an exponential family distribution; (2) The step size $γ$ and the Lipschitz constant L satisfy $γ < 1 / L$ . (3) the sequence $- {x^{k}}$ is generated by proximal gradient descent, and the corresponding object function value sequence $- {F (x^{k})}$ . Then, the function sequence $- {F (x^{k})}$ is monotonically decrease and bounded below, and thus converges to a finite value.

Remarks.

The detailed proof is provided in Appendix II. This Proposition 2 ensures the convergence of the sequence of objective function values, thereby paving the way for the convergence of algorithm 1.

Now we are ready to establish the convergence analysis of Algorithm 1. By combining the result of Proposition 2 and some assumptions, we can derive the theoretical results as follows.

Theorem 1.

Assume that

$f (x)$ and $g (x)$ satisfy the assumption in the proposition 2 .

Existence of solution to subproblem and boundedness of the generated sequence $- {x^{k}}$ .

Then the Algorithm.1 based on proximal gradient descent converges.

Remarks.

The detailed proof can be seen in Appendix II. Theorem 1 indicates that Algorithm 1 designed in this paper is convergent, and the solution deviation is controllable under certain conditions.

Theorem 2.

Assume that the generated sequence $- {x^{k}}$ in Theorem 1 converge to $x *$ , that is

x * = \underset{x}{\arg min} ‖ H x - y ‖_{2}^{2} - η \log p (x) .

(21)

Replace the score function $\nabla_{x} \log p (x)$ with $s_{θ} (x, t)$ satisfying $‖ s_{θ} (x, t) - \nabla_{x} \log p (x) ‖^{2} \leq τ^{2}$ with $τ$ is positive, and the optimal point of the corresponding objective function is denoted as $\tilde{x}$ , Then

‖ x * - \tilde{x} ‖^{2} \leq \frac{η τ}{σ_{min}^{2} - η L},

(22)

where

σ_{min}

is the minimum singular value of the matrix

H

and L is the Lipschitz constant of

\nabla_{x} \log p (x)

, which satisfying

σ_{min}^{2} - η L > 0

Remarks.

The detailed proof can be seen in Appendix II. Theorem 2 establishes that when neural networks are employed to approximate the score function with approximation errors constrained with $τ$ , the solution deviation can be rigorously bounds. Specifically, the theoretical upper limit of such deviation is explicitly formulated by $η τ / (σ_{min}^{2} - η L)$ .

4 Experiment and results

Experiment design

In this section, we introduce the comparison algorithms and the datasets, as well as the key details of training and evaluating the model. We describe the comparison algorithms employed, including the two-step decomposition approach (FBP-decomp), the method based on TV (ESART-TV), IRM-MI,¹³ OPMT,⁴⁰ and GAN-MD. The FBP-decomp algorithm involves first reconstructing sparse images using FBP and then directly inverting them in the image domain. The GAN-MD algorithm, adapted from reference,¹⁹ replaces the multi-material decomposition network with a dual-material decomposition network and shifts the focus from low-dose denoising to sparse artifact suppression. The GAN-MD network was trained on 1095 data samples constructed from dataset of the local hospital, partitioned into training, validation, and test sets at an 8:1:1 ratio. The trained model was then directly applied to real mouse scan data for material decomposition. The modified GAN-MD model takes dual-energy CT images as input and outputs two basis material images (bone and soft tissue). However, two key limitations affect model performance: (1) insufficient training data impedes model convergence, and (2) the lack of ground truth material labels restricts the construction of high-quality training dataset. These challenges are further compounded by the sparse-view scanning protocol employed in this study, which significantly increase the problem complexity. Consequently, the GAN-MD model exhibits constrained performance, particularly impairing its artifact suppression capability in the output material images. Furthermore, since the GAN-MD model was trained solely on simulated data, its generalization capability to real data may be limited, potentially affecting artifact suppression efficacy in practical applications.

We have optimized the implementation of each comparison algorithm to ensure the fairness of the experiments. To quantitatively assess the performance of the proposed algorithm, we use evaluation metrics such as Peak Signal-to-Noise Ratio (PSNR), Structural Similarity (SSIM),⁴¹ and Root Mean Square Error (RMSE). The reconstructed image exhibits lower distortion when demonstrating either reduced RMSE or increased PSNR values. The RMSE and PSNR are calculated as follows:

RMSE = \sqrt{\frac{\sum_{i = 1}^{N} {| f (i) - f^{r e f} (i) |}^{2}}{N}}, PSNR = 10 \log_{10} \frac{‖ f^{r e f} (i) ‖_{\infty}^{2}}{RMSE},

where

f (i), f^{r e f} (i)

denote the reconstructed and reference images at i pixel, respectively, and N represents the total pixel of the image.

This model is constructed on the public PyTorch platform with the NVIDIA RTX A6000 GPU and Inter(R) Xeon(R) Silver 4216 CPU @ 2.10 GHz and 64GB RAM. We employed a pre-trained model provided by our team. Our diffusion model network is trained by Adam optimizer⁴² and the learning rate is set to 1 × 10⁻⁴. The number of epochs is 300 and batch size is 1. The total training time for the diffusion model is about 300 h. The distance between the scanner rotation center and the X-ray source is set to 300 mm. The projections in each view are collected using a linear detector, which consists of 1024 bins with a size of 0.124 mm. The spectra of the simulation data under tube voltages of 80 kVp and 140 kVp were obtained by using the ImaSim software. The corresponding spectrum distributions are shown in Figure 2. Through simulating the fast kVp switching scan, the ray paths corresponding to the two spectrum were inconsistent. In the experiment, 90 views or 120 views evenly distributed among 360 views are selected to validate the performance of the suggested algorithm. In this work, parameters $λ, β$ were empirically selected.

Figure 2.

The experimental data include images from the LIDC dataset slices in (a) and (b), with (c) is synthesized CT image from segmented image by local hospital with display window [0.001 0.03], and (d) is the real mouse data reconstructed from full energy spectra with display window [0 0.08]. Figure (e) denotes normalized X-ray spectra and (f) the attenuations coefficients of bone and tissue.

Data preparation

To evaluate the performance of the algorithm proposed in this paper, extensive experimental verifications were conducted using preclinical data as shown in Figure 2. This preclinical data includes two datasets: one is a thoracic dataset pre-segmented and labelled by a local hospital, covering 1591 thoracic images from 5 patients. The primary basis materials within the images include bones and soft tissues, with the segmentation tasks mainly assisted by radiologists from the local hospital. The radiologists manually corrected the labelled images for each pair of bone and tissue using the bone removal function in the Syngo. Via software, equipped with the Siemens SOMATOM Definition

Flash CT scanner. The LAC of bone and tissue materials can be obtained from the National Institute of Standards and Technology. The other is the publicly available LIDC-IDRI lung dataset, which contains lung images from 20 patients. Due to the lack of accurate basis material image segmentation within the LIDC dataset, it was unfeasible to use them to validate the algorithm performance. Instead, we selected 50 images from the pre-segmented thoracic dataset of the local hospital for verification.

The real experimental mouse data as shown in Figure 2(d) were acquired using a spectral CT imaging system at the Institute of High Energy Physics, Chinese Academy of Sciences. Post euthanasia, mouse-tail specimens injected with iodine contrast agent were scanned with fixed source-to-object and source-to-detector distances of $200.8 mm$ and $362.8 mm$ , respectively. The detector comprised $512 \times 15$ elements with a pixel size of $0.4 mm \times 0.4 mm$ . Scanning parameters included a tube voltage of $60 keV$ and current of $72 mAs$ , with energy thresholds set at $12, 26, 34, 42 keV$ . Prior to testing the real data, the initial spectrum was estimated through linear combination of model spectra, followed by optimization via the expectation-maximization (EM) algorithm⁴³ for spectral estimation.

Experiment results

In the numerical simulation experiments, we first validated the proposed algorithm using the pre-segmented and labelled dataset by the local hospital. The reconstructed image results and numerical evaluation metrics are presented in Figure 3 and Table 1, respectively. Figure 3 illustrates the reconstructed image of the second slice under 120 sparse views. It is evident that the basis material image reconstructed by the FBP-decomp exhibits significant noise and numerous streak artifacts. Although the OPMT method improves the quality of the material decomposition image, edge blurring issues remain. The IRM-MI method reduces artifacts compared to the two-step decomposition method and provides clearer structures compared to the OPMT method; however, soft tissue parts remain difficult to discern clearly. The GAN-MD model exhibits limited space artifact suppression capability, preserving most structural details in reconstructed images but introducing noticeable streak artifacts. Compared to the aforementioned methods, the ESART-TV method shows remarkable effectiveness in noise suppression and artifact removal, yet when examining the magnified regions of interest (ROI) in soft tissue basis materials, the details appear overly smooth, and the decomposed soft tissue material image contains noise. The proposed algorithm in this study achieves superior decomposition results compared to the ESART-TV method, producing synthesized VMIs closest to the reference images with minimal noise and more realistic textures.

Figure 3.

Reconstructed material images on the local hospital simulated clinical data from 120 projections. The 1^st-7^th columns stand for the reference, FBP-decomposition, ESART-TV, OPMT, IRM-MI, GAN-MD and our proposed methods. The 1^st row represent bone material and the corresponding ROIs, and the 2^nd row represent soft tissue material and the corresponding ROIs. The 3^rd row represent virtual monochromatic images at 75 keV.The display windows are [0.0001 0.8], [0.1 0.88], and [0.001 0.035], respectively.

Table 1.

Statistical reconstruction PSNR/SSIM (10⁻³) for different algorithms on the local hospital clinical data from 120 projections (50 test slices).

Material	Metrics	FBP-decomp	ESART-TV	OPMT	IRM-MI	GAN-MD	Ours
Bone	PSNR	27.6412 $\pm$ 1.2	42. 9141 $\pm$ 1.8	38.7519 $\pm$ 1.1	42.0109 $\pm$ 0.7	38.8213 $\pm$ 1.9	44.9020 $\pm$ 1.3
	SSIM	0.7349 $\pm$ 1.3	0.9715 $\pm$ 1.1	0.9501 $\pm$ 0.7	0.9607 $\pm$ 0.8	0.9452 $\pm$ 1.3	0.9931 $\pm$ 0.8
	RMSE	0.0298 $\pm$ 1.4	0.0134 $\pm$ 0.8	0.0251 $\pm$ 0.5	0.0205 $\pm$ 0.5	0.0232 $\pm$ 0.9	0.0042 $\pm$ 0.7
Soft	PSNR	14.4318 $\pm$ 1.6	31.8421 $\pm$ 1.6	26.8245 $\pm$ 1.1	31.1109 $\pm$ 0.5	29.3125 $\pm$ 1.7	34.3033 $\pm$ 1.5
	SSIM	0.5081 $\pm$ 1.5	0.9439 $\pm$ 0.9	0.8105 $\pm$ 0.7	0.8218 $\pm$ 0.7	0.8054 $\pm$ 0.8	0.9503 $\pm$ 0.9
	RMSE	0.2209 $\pm$ 1.3	0.0384 $\pm$ 0.7	0.1149 $\pm$ 0.5	0.0502 $\pm$ 0.9	0.0917 $\pm$ 0.9	0.0209 $\pm$ 1.0

The bolded data in this table represent the corresponding evaluation metrics of the algorithm proposed in this paper, highlighted solely for the purpose of emphasizing its performance indicators.

As indicated in Table 1, the numerical evaluation metrics for the bone material images significantly outperformed those of the soft tissue basis material images, primarily due to the simpler structure of bone materials leading to better decomposition results. In contrast, the complex structure of soft tissues, especially when using the FBP-decomp, results in ineffective artifact removal and poorer evaluation metrics. The proposed algorithm achieved the best evaluation metrics, consistent with the visual observation results.

The results of material image reconstruction on LIDC dataset are presented in Figure 4. The FBP-decomp results still exhibit noticeable artifacts. The OPMT produces soft tissue material images with blurred boundaries. The IRM-MI also retains some noise, which is particularly evident in the magnified ROI. The GAN-MD algorithm exhibits significant noise, resulting in degraded texture clarity. Visually, both the ESART-TV and the proposed method exhibit better performance in material image, effectively preserving tissue structures and image details. However, pixel line profiles drawn along the yellow dashed lines in Figure 4 reveal significant deviations between the pixel intensity values of ESART-TV and the ground truth, whereas the values from the proposed method closely match the true values. The proposed algorithm accurately maintains edge and texture information while removing sparse artifacts. The deviations of pixel values from the material ground truth for other comparison algorithms can also be visualized in Figure 5.

Figure 4.

Reconstructed material images on LIDC data from 120 projections. The 1^st-7^th columns stand for the reference, FBP-decomposition, ESART-TV, OPMT, IRM-MI, GAN-MD and our method. The 1^st row represent bone material and the corresponding ROIs, and the 2^nd row represent soft tissue material and the corresponding ROIs, with the gray-scale value ranges [0.0001 0.9], and [0.1 0.99], respectively.

Figure 5.

Line profile of different reconstruction methods corresponding to the yellow dotted line in figure 4.

As shown in Table 2, the proposed algorithm achieves the optimal values for all three evaluation metrics. Compared to the ESART-TV method, the PSNR of the reconstructed bone material images improved by 7.01%, and the RMSE decreased by 58.54%, demonstrating the superior performance and effectiveness of the proposed algorithm in material image reconstruction.

Table 2.

Reconstruction PSNR/SSIM of LIDC data by different methods from 120 views.

Material	Metrics	FBP-decomp	ESART-TV	OPMT	IRM-MI	GAN-MD	Ours
Bone	PSNR	27.3933	42.9058	36.1887	41.3784	36.4251	45.9120
	SSIM	0.7531	0.9695	0.9421	0.9512	0.9487	0.9926
	RMSE	0.0427	0.0123	0.0219	0.0135	0.0194	0.0051
Soft	PSNR	14.9169	31.8105	25.7093	31.1203	30.0122	34.8133
	SSIM	0.5860	0.9393	0.8055	0.8198	0.8093	0.9511
	RMSE	0.1900	0.0524	0.1640	0.0582	0.1120	0.0189

The bolded data in this table represent the corresponding evaluation metrics of the algorithm proposed in this paper, highlighted solely for the purpose of emphasizing its performance indicators.

To comprehensively validate the performance of the proposed algorithm, we conducted experiments using 90 sparse views. The basis material images under 90 sparse views on the LIDC dataset were presented in Figure 6. As shown in Figure 6, the FBP-decomp method produces poor-quality material images, making it difficult to clearly discern tissue structures and details. The OPMT and IRM-MI both reduce artifacts compared to the FBP-decomp but still exhibits significant noise. The GAN-MD algorithm yields favorable error images, yet exhibits blurred structural details that obscure texture discernibility, as visually confirmed in the corresponding ROI images. The ESART-TV algorithm provides relatively clear reconstructed images; however, it suffers from excessive smoothing, as evidenced in the magnified ROI. Additionally, noise residues are clearly visible in the error images of the second and fourth rows. In contrast, the proposed algorithm effectively eliminates artifacts and suppresses noise. The magnified ROI show that the reconstructed images well preserve structural details and remain sharp edges. Table 3 lists the quantitative evaluation metrics of different decomposition methods under 90 sparse angles, which are consistent with Figure 6.

Figure 6.

Reconstructed material images on the LIDC data from 90 projections. The 1^st-7^th columns stand for the reference, FBP-decomposition, ESART-TV, OPMT, IRM-MI, GAN-MD and our methods. The 1^st row and the 3^rd row represent bone material and soft tissue material and the corresponding ROIs. The 2^nd and 4^th row represent the error maps. The gray-scale value ranges for bone and soft maps are [0.0001 0.9], and [0.1 0.99], respectively.

Table 3.

Reconstruction PSNR/SSIM of preclinical LIDC data by different algorithms from 90 views.

Material	Metrics	FBP-decomp	ESART-TV	OPMT	IRM-MI	GAN-MD	Ours
Bone	PSNR	26.0559	35.8349	34.3153	34.9142	37.4549	42.5235
	SSIM	0.6721	0.9857	0.8957	0.9303	0.9665	0.9882
	RMSE	0.0498	0.0162	0.0193	0.0180	0.0119	0.0075
Soft	PSNR	14.4118	24.1138	19.9901	23.6503	20.1090	34.2754
	SSIM	0.5368	0.9072	0.6980	0.7519	0.8812	0.9463
	RMSE	0.1903	0.0623	0.1011	0.0657	0.0988	0.0193

The bolded data in this table represent the corresponding evaluation metrics of the algorithm proposed in this paper, highlighted solely for the purpose of emphasizing its performance indicators.

Real data experiments

This section evaluates algorithm performance using real mouse data. Figure 2(d) displays the full-spectrum CT reconstruction obtained through FBP. Comparative material decomposition results for bone and soft tissue are shown in Figure 7, where all algorithms demonstrate effective streak artifacts suppression. Our ROI analysis provides qualitative assessment of the artifact suppression and detail preservation capabilities of the different algorithms. The FBP-decomp reconstruction exhibits the most severe noise contamination with poorly resolved structures. Both OPMT and IRM-MI reconstructions exhibit significant edge blurring, while ESART-TV introduces noticeable blocky artifacts that compromise anatomical detail. Unlike the simulation training data, real scanned data present additional challenges including pronounced noise and beam hardening artifacts. These factors contribute to substantial residual artifacts in the deep learning-based GAN-MD reconstruction. As evidenced in Figure 7, despite being trained exclusively on simulation data, our proposed method achieves significantly better performance than GAN-MD on real mouse data. While GAN-MD reconstructions contain noticeable noise and artifacts, our method demonstrates superior artifact suppression while maintaining excellent preservation of texture details and edge structures. This performance gap highlight our model's enhanced generalization capability in real-world imaging scenarios.

Figure 7.

Reconstructed material images on the real data from 120 projections. The 1^st-6^th columns stand for FBP-decomposition, ESART-TV, OPMT, IRM-MI, GAN-MD and our methods. The first row and the second row represent bone material and soft tissue material and the corresponding ROIs. The gray-scale value ranges for bone and soft maps are [0.01 0.8], and [0.1 0.7], respectively.

The red boxes (A-D) in Figure 7 denote the ROIs for quantitative evaluation of reconstructed image quality. Regions A and B contain rich structural details, while regions C and D are smooth and homogeneous. Since the ground truth of real data reconstruction is unavailable, we employed two metrics for algorithm performance assessment: mean equivalent attenuation coefficient (for reconstruction accuracy) and noise standard deviation (SD, for noise suppression capability). Table 4 presents the comparative results of these metrics across different algorithms. Our method achieved mean equivalent attenuation coefficients comparable to other approaches, demonstrating its accuracy in spectral CT reconstruction of real data. Regarding noise standard deviation, while GAN-MD yielded the highest values, our algorithm maintained finer image details while achieving lower noise levels. The proposed method achieved optimal balance between reconstruction fidelity and noise reduction.

Table 4.

Mean values and SDs of ROI A, B, C, and D in figure 7.

Metric	ROI	FBP-decomp	ESART-TV	OPMT	IRM-MI	GAN-MD	Ours
Mean values $\pm$ STD	A	0.1059 $\pm$ 0.16	0.0982 $\pm$ 0.15	0.1015 $\pm$ 0.16	0.1076 $\pm$ 0.16	0.1053 $\pm$ 0.23	0.0966 $\pm$ 0.14
	B	0.1689 $\pm$ 0.22	0.1696 $\pm$ 0.22	0.1768 $\pm$ 0.22	0.1701 $\pm$ 0.22	0.1776 $\pm$ 0.30	0.1690 $\pm$ 0.20
	C	0.8083 $\pm$ 0.04	0.8034 $\pm$ 0.03	0.8017 $\pm$ 0.06	0.8103 $\pm$ 0.04	0.7998 $\pm$ 0.17	0.7981 $\pm$ 0.03
	D	0.7527 $\pm$ 0.16	0.7432 $\pm$ 0.13	0.7501 $\pm$ 0.11	0.7547 $\pm$ 0.12	0.7494 $\pm$ 0.23	0.7380 $\pm$ 0.11

The bolded data in this table represent the corresponding evaluation metrics of the algorithm proposed in this paper, highlighted solely for the purpose of emphasizing its performance indicators.

Ablation experiments

To optimize material image reconstruction fidelity, this study innovatively integrates deep priors from the probability diffusion model (DM) into the iterative reconstruction framework. The experimental results from comprehensive simulations demonstrate the algorithm's superior capability in material decomposition. Ablation studies contrasting implementations with and without DM integration reveal critical performance difference. As illustrated in Figure 8, the baseline method (without DM) exhibits pronounced artifacts in bone material decomposition and blurred soft-tissue boundaries, whereas the DM-enhanced approach effectively suppresses noise while preserving structural edges, particularly in arrow-indicated regions. Quantitative analysis in Table 5 confirms these improvements, showing an 11.76% reduction in bone material RMSE and 4.85% increase in soft-tissue PSNR compared to the algorithm without the probability diffusion model implementations, thereby robustly validating the DM's capability to enhance decomposition accuracy while maintaining projection data consistency.

Figure 8.

Reconstruction results from the ablation experiment from 90views. From left to right are: the reference, with the diffusion model, and no diffusion model. The gray-scale value ranges are [0.0001 0.9], and [0.1 0.99], respectively.

Table 5.

Reconstruction PSNR/SSIM of ablation experiment.

Methods	Material	PSNR	SSIM	RMSE
With-diffusion	Bone	42.5235	0.9882	0.0075
With-diffusion	Soft	34.2754	0.9463	0.0193
No-diffusion	Bone	41.4442	0.9847	0.0085
No-diffusion	Soft	32.6902	0.9341	0.0232

Discussion

To address the limitations of model-based material decomposition methods, we propose an unsupervised direct material decomposition framework that leverages deep priors. This approach synergistically optimizes VMIs through a diffusion probabilistic model in sparse-view spectral CT. This deep prior mechanism facilitates structural detail preservation by exploiting a comprehensive statistical distribution analysis of data. Through precise modeling of intrinsic image statistics, our framework simultaneously accomplishes two critical objectives: (1) robust image restoration through artifact suppression, and (2) fundamental recovery of essential physical attributes with preserved edge integrity. The proposed method employs VMIs as critical discriminative enhancers for polychromatic projections, effectively overcoming convergence limitations in iterative reconstruction algorithms. By integrating probabilistic diffusion priors into the optimization process, our framework achieves superior refinement of material-specific representations. It systematically enforces a data fidelity term to ensure measurement consistency, while applying probabilistic regularization to suppress unwanted structures–thereby guaranteeing anatomically plausible material image reconstruction. The results of the experiments on preclinical data and real data demonstrate the algorithm's robustness under challenging conditions, including geometric inconsistencies and sparse sampling, while maintaining enhanced reconstruction fidelity.

However, while diffusion models (DM) demonstrate remarkable reconstruction capabilities, they inherently suffer from low sampling efficiency due to their iterative refinement mechanism. Conventional implementations typically require hundreds to thousands of sequential steps to gradually denoise Gaussian noise into high-fidelity images, with each step demanding computationally intensive neural network predictions. As evidenced in Table 6, this cumulative computational burden results in slightly higher time costs for our proposed algorithm in this work compared with ESART-TV, posing practical challenges for real-time clinical applications. To address this efficiency bottleneck, we implement a truncation diffusion model activation strategy within the iterative reconstruction workflow. This approach strategically restricts diffusion model-based optimization program to the final stages of the reconstruction process. Empirical analysis in Table 7 reveals that material-specific RMSE values for both bone and soft tissue exhibit minimal fluctuations and insignificant changes during these final rounds. With a reduction in the number of truncated iterations, computational efficiency increases, but the evaluation metrics decline. Specifically, the PSNR and SSIM values peak when the iteration count is set to 50, indicating optimal reconstruction performance at this point. Therefore, this study determines to introduce the optimization strategy of the DM in the last 50 rounds of iterations.

Table 6.

Iteration time of different methods (s).

Views	ESART-TV	IRM-MI	Ours
90-views	3.3	1.1	3.6
120-views	3.8	1.4	4.4

Table 7.

Reconstruction PSNR/SSIM of different truncation iterations.

Material	Metrics	10	20	50	100	200
Bone	PSNR	54.1089	54.1093	54.1776	54.1056	54.1076
	SSIM	0.9985	0.9987	0.9989	0.9988	0.9986
	RMSE	0.0020	0.0020	0.0020	0.0020	0.0020
Soft	PSNR	43.2564	43.2736	43.3393	43.2846	43.2846
	SSIM	0.9929	0.9931	0.9943	0.9930	0.9930
	RMSE	0.0069	0.0068	0.0065	0.0066	0.0068

This paper proposes a novel framework for material image reconstruction that leverages direct material decomposition in the projection domain. This approach unifies the nonlinear observation model of polychromatic projections with the fundamental physics of material attenuation. By integrating spectral sensitivity modeling and decomposition priors into a single optimization paradigm, the proposed methodology eliminates cascading approximations inherent to traditional two-step approaches (e.g., projection-domain decomposition followed by image reconstruction). Consequently, it mitigates cumulative information loss while enhancing decomposition accuracy. Notably, the current DMDM implementation assumes precise a priori spectral calibration—a condition rarely satisfied in practical scenarios due to spectral drift caused by detector imperfections, source variability, or object scatter. Experimental validation on real-world datasets reveals that spectral estimation errors induce non-negligible biases in decomposition outcomes. This spectral dependency underscores a critical research frontier: advancing robust spectral self-calibration techniques or hybrid learning-physics architectures to relax stringent spectral knowledge requirements while preserving decomposition fidelity.

Conclusion

In this study, we propose a novel spectral CT direct material decomposition framework enhanced by a diffusion probability model to address the challenges of sparse-view spectral CT imaging, particularly in scenarios involving fast kVp switching protocols. To overcome the inherent limitations of slow convergence in iterative spectral CT reconstruction, we innovatively integrate virtual monochromatic images as spectral discriminators, which expand the angular coverage of the spectral surface equations and dynamically refine the reconstruction manifold. By synergistically optimizing virtual monochromatic images through the probability diffusion model-guided denoising, our framework achieves improvement in material decomposition accuracy compared to conventional two-step approaches. Comprehensive validation with preclinical data across diverse sparse-view configurations (90–120 projections) demonstrates consistent performance for both bone and soft-tissue materials. These results not only validate the framework's robustness against spectral inconsistency and noise amplification but also establish a paradigm for physics-informed diffusion models in spectral imaging. Future work will focus on relaxing spectral calibration dependencies through self-supervised spectral estimation, further bridging the gap between simulation-driven innovation and clinical applicability.

Footnotes

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 62271504, 62101596) and Natural Science Foundation of Henan (Grant No.252300420395). This work was also supported by Technology Innovation Leading Talent Project of Zhongyuan (Grant No. 244200510015).

ORCID iDs

Jie Guo

Ailong Cai

Zhizhong Zheng

Bin Yan

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Research and development of early defect insitu detection technology and intelligent warning equipment for key components of ultra-high/extra high voltage electrical equipment Key Research and Development Special Project of Henan Province, China under Grant number 251111220600.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Appendix I

In this section, we provide the specific expressions for the coefficients in the first-order Taylor expansion, and the linearized representation of the DECT imaging model is as follow,

(23)

p_{k, L} \approx p_{k, L}^{(n)} + \frac{Φ_{k, L}^{(n)}}{q_{k, L}^{(n)}} A (b_{1} - b_{1}^{(n)}) + \frac{Θ_{k, L}^{(n)}}{q_{k, L}^{(n)}} A (b_{2} - b_{2}^{(n)}) .

where (24)

\begin{aligned} p_{k, L}^{(n)} = - \ln (\sum_{m = 1}^{M_{k}} S_{k, m} \cdot δ_{E} \cdot \exp (- μ_{1} (e) A b_{1}^{(n)} - μ_{2} (e) A b_{2}^{(n)})) \\ q_{k, L}^{(n)} = \sum_{m = 1}^{M_{k}} S_{k, m} \cdot δ_{E} \cdot \exp (- μ_{1} (e) A b_{1}^{(n)} - μ_{2} (e) A b_{2}^{(n)}) \\ Φ_{k, L}^{(n)} = \sum_{m = 1}^{M_{k}} S_{k, m} \cdot δ_{E} \cdot μ_{1} (e) \cdot \exp (- μ_{1} (e) A b_{1}^{(n)} - μ_{2} (e) A b_{2}^{(n)}) \\ Θ_{k, L}^{(n)} = \sum_{m = 1}^{M_{k}} S_{k, m} \cdot δ_{E} \cdot μ_{2} (e) \cdot \exp (- μ_{1} (e) A b_{1}^{(n)} - μ_{2} (e) A b_{2}^{(n)}) . \end{aligned}

For the sake of simplicity, Eq. (24) is rewritten as a linear system in the following form, (25)

p \approx H_{(n)} b, H_{(n)} = (\begin{matrix} \frac{Φ_{h, L}^{(n)}}{q_{h, L}^{(n)}} & \frac{Θ_{h, L}^{(n)}}{q_{h, L}^{(n)}} \\ \frac{Φ_{l, L}^{(n)}}{q_{l, L}^{(n)}} & \frac{Θ_{l, L}^{(n)}}{q_{l, L}^{(n)}} \end{matrix}) A,

where the column vector

p

can be indicated as follows with the

2 N \times 1

dimension, (26)

p = (\begin{matrix} p_{h, L} \\ p_{l, L} \end{matrix}) - (\begin{matrix} p_{h, L}^{(n)} \\ p_{l, L}^{(n)} \end{matrix}) + H_{(n)} (\begin{matrix} b_{1}^{(n)} \\ b_{1}^{(n)} \end{matrix}),

The VMIs-guided method increase the angles between polychromatic projection curves by introducing monochromatic images, and thus accelerating the iterative process. Specifically, two distinct energies $E_{1}$ and $E_{2}$ are introduced, and the corresponding VMIs are denoted as $I_{E 1}, I_{E 2}$ . These VMIs can be represented as linear combinations of the basis materials, (27)

(\begin{matrix} I_{E 1} \\ I_{E 2} \end{matrix}) = (\begin{matrix} μ_{1} (E_{1}) & μ_{2} (E_{1}) \\ μ_{1} (E_{2}) & μ_{2} (E_{2}) \end{matrix}) (\begin{matrix} b_{1} \\ b_{1} \end{matrix}),

Eq. (13) can be equivalently represented as the following, (28)

(\begin{matrix} b_{1} \\ b_{1} \end{matrix}) = (\begin{matrix} c_{11} & c_{12} \\ c_{21} & c_{22} \end{matrix}) (\begin{matrix} I_{E 1} \\ I_{E 2} \end{matrix}),

where

c_{11}, c_{12}, c_{21}, c_{22}

denote the corresponding elements in the inverse of the coefficient matrix in (13).

Applying the projection matric A to the above Eq.(28), and denote $F = A b_{1}, G = A b_{2}, F^{'} = A (I_{E 1}), G^{'} = A (I_{E 2}),$ then we substitute the formula $F, G$ into(1), the polychromatic projections can be represented as follows,

(29)

p_{k, L} = - \ln (\sum_{m = 1}^{M_{k}} S_{k, m} \cdot δ_{E} \cdot \exp (- μ_{1}^{'} (e) F^{'} - μ_{2}^{'} (e) G^{'})),

where

μ_{1}^{'} (e) = c_{11} μ_{1} (e) + c_{21} μ_{2} (e), μ_{2}^{'} (e) = c_{12} μ_{1} (e) + c_{22} μ_{2} (e)

represent the equivalent LAC. At this point, the material image reconstruction is transformed into a problem of reconstructing the VMIs

I_{E 1}

and

I_{E 1}

from Eq. (29) with the known projections

p_{k, L}

and equivalent LAC

μ_{1}^{'} (e), μ_{2}^{'} (e)

, which can similarly be solved using an iterative strategy.

Similarly, if the n -th iteration result of the VMIs is known as $(I_{E 1}^{(n)}, I_{E 1}^{(n)})$ , then the first-order Taylor expansion of the right-hand side of (29) at this point can be expressed as follows: (30)

p_{k, L} \approx p_{k, L}^{^{'} (n)} + \frac{Φ_{k, L}^{^{'} (n)}}{q_{k, L}^{^{'} (n)}} A (I_{E 1} - I_{E 1}^{(n)}) + \frac{Θ_{k, L}^{^{'} (n)}}{q_{k, L}^{^{'} (n)}} A (I_{E 2} - I_{E 2}^{(n)}),

where (31)

\begin{aligned} p_{k, L}^{^{'} (n)} = - \ln (\sum_{m = 1}^{M_{k}} S_{k, m} \cdot δ_{E} \cdot \exp (- μ_{1}^{'} (e) F^{^{'} (n)} - μ_{2}^{'} (e) G^{^{'} (n)})) \\ q_{k, L}^{^{'} (n)} = \sum_{m = 1}^{M_{k}} S_{k, m} \cdot δ_{E} \cdot \exp (- μ_{1}^{'} (e) F^{^{'} (n)} - μ_{2}^{'} (e) G^{^{'} (n)}) \\ Φ_{k, L}^{^{'} (n)} = \sum_{m = 1}^{M_{k}} S_{k, m} \cdot δ_{E} \cdot μ_{1}^{'} (e) \cdot \exp (- μ_{1}^{'} (e) F^{^{'} (n)} - μ_{2}^{'} (e) G^{^{'} (n)}) \\ Θ_{k, L}^{^{'} (n)} = \sum_{m = 1}^{M_{k}} S_{k, m} \cdot δ_{E} \cdot μ_{2}^{'} (e) \cdot \exp (- μ_{1}^{'} (e) F^{^{'} (n)} - μ_{2}^{'} (e) G^{^{'} (n)}) . \end{aligned}

Eq. (30) can be equivalently represented as the following linear system, (32)

p \approx H_{(n)}^{'} I_{E}, H_{(n)}^{'} = (\begin{matrix} \frac{Φ_{h, L}^{^{'} (n)}}{q_{h, L}^{^{'} (n)}} & \frac{Θ_{h, L}^{^{'} (n)}}{q_{h, L}^{^{'} (n)}} \\ \frac{Φ_{l, L}^{^{'} (n)}}{q_{l, L}^{^{'} (n)}} & \frac{Θ_{l, L}^{^{'} (n)}}{q_{l, L}^{^{'} (n)}} \end{matrix}) A,

Appendix II

References

Faul

Couch

Cann

, et al. Composition-selective reconstruction using dual-energy CT for bone mineral quantification. J. Comput. Assist. Tomography 1982; 6: 202–204.

Lamb

Sahani

Fuentes-Orrego

. Stratification of patients with liver fibrosis using dual-energy CT. IEEE Trans Med Imaging 2014; 34: 807–815.

Kalender

Klotz

Kostaridou

. An algorithm for noise suppression in dual energy CT material density images. IEEE Trans Med Imaging 1988; 7: 218–224.

Crawley

Evans

Owen

. A theoretical analysis of the accuracy of single-energy CT bone-mineral measurements. Physics in Medicine & Biology 1988; 33: 1113.

Cai

Rodet

Legoupil

, et al. A full-spectral Bayesian reconstruction approach based on the material decomposition model applied in dual-energy computed tomography. Med Phys 2013; 40: 111916.

Zhao

Niu

Xing

, et al. Using edge-preserving algorithm with non-local mean for significantly improved image-domain material decomposition in dual-energy CT. Phys Med Biol 2016; 61: 1332.

Kalender

. X-ray computed tomography. Phys Med Biol 2006; 51: 29–43.

Ducros

Abascal

Sixou

, et al. Regularization of nonlinear decomposition of spectral x-ray projection images. Med Phys 2017; 44: 174–187.

Alvarez

Macovski

. Energy-selective reconstructions in X-ray computerised tomography. Phys Med Biol 1976; 21: 733–744.

10.

Zhao

, et al. Iterative dual energy material decomposition from spatial mismatched raw data sets. J Xray Sci Technol 2014; 22: 745–762.

11.

Zhao

Zhang

. An extended algebraic reconstruction technique (E-ART) for dual spectral CT. IEEE Trans Med Imaging 2015; 34: 761–768.

12.

Zhao

Wang

. An extended simultaneous algebraic reconstruction technique (E-SART) for X-ray dual spectral computed tomography. Scanning 2016; 38: 599–611.

13.

Zhang

Zhao

Pan

, et al. An iterative reconstruction method based on monochromatic images for dual energy CT. Med Phys 2021; 48: 6437–6452.

14.

Sidky

Pan

. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Phys Med Biol 2008; 53: 4777–4807.

15.

Wang

Liu

, et al. Block matching frame based material reconstruction for spectral CT. Phys Med Biol 2019; 64: 235011.

16.

Chang

, et al. Spectrum estimation-guided iterative reconstruction algorithm for dual energy CT. IEEE Trans Med Imaging 2020; 39: 246–258.

17.

Zhang

Wang

, et al. Image domain dual material decomposition for dual-energy CT using butterfly network. Med Phys 2019; 46: 2037–2051.

18.

Ren

Zhang

Wang

, et al. A dual-energy CT reconstruction method based on anchor network from dual quarter scans. J Xray Sci Technol 2024; 32: 229–252.

19.

Ren

Zheng

Wang

, et al. Prior image-based generative adversarial learning for multi-material decomposition in photon counting computed tomography. Comput Biol Med 2024; 180: 108854.

20.

Zhu

Zhang

, et al. Feasibility study of three-material decomposition in dual-energy cone-beam ct imaging with deep learning. Phys Med Biol 2022; 67: 145012.

21.

Wang

Xia

, Z. et al. An End-to-End deep network for reconstructing CT images directly from sparse sinograms, IEEE Transactions on Computational Imaging 6(2020), 1548–1560.

22.

Liu

Zhang

, et al. Decomposition iteration strategy for low-dose CT denoising. J Xray Sci Technol 2024; 32: 493–512.

23.

Jain

Abbeel

. Denoising diffusion probabilistic models. arXiv: 2020, 2006.11239.

24.

Song

Shen

Xing

, et al. Solving inverse problems in medical imaging with score-based generative models. In International Conference on Learning Representations 2021: 186–203. https://doi.org/10.48550/arXiv.2111.08005

25.

Wang

. Time-reversion fast-sampling score-based model for limited-angle CT reconstruction. IEEE Trans Med Imaging 2024; 43: 3449–3460.

26.

Chung

Kim

Mccann

, et al. Diffusion posterior sampling for general noisy inverse problems, 2023, arXiv: 2209.14687.

27.

Chung

Kim

Mccann

, et al. Improving diffusion models for inverse problems using manifold constraints, 2022, arXiv:2206.00941.

28.

Liu

Zhou

Wei

, et al. Sparse-view spectral CT reconstruction and material decomposition based on multi-channel SGM. IEEE Trans Med Imaging 2024; 43: 3425–3435.

29.

Jiang

Gang

Stayman

. Multi-Material Decomposition Using Spectral Diffusion Posterior Sampling. 2024, arXiv, 2408.01519.

30.

Song

Meng

Ermon

. Denoising diffusion implicit models, 2020, arXiv, 2010.02502.

31.

Guo

Xiang

Yang

, et al. Accelerating denoising diffusion probabilistic model via truncated inverse processes for medical image segmentation. Comput Biol Med 2024; 180: 108933.

32.

Ruthotto

Haber

. An introduction to deep generative modeling. GAMM-Mitteilungen 2021; 44: e202100008.

33.

Sohl-Dickstein

Weiss

Maheswaranathan

, et al. Deep unsupervised learning using non-equilibrium thermodynamics. In International Conference on Machine Learning 2015; 2015: 2256–2265.

34.

Jain

Abbel

. Denoising diffusion probabilistic models. Adv Neural Inf Process Syst 2020; 33: 6840–6851.

35.

Song

Ermon

. Generative modeling by estimating gradients of the data distribution. Adv Neural Inf Process Syst 2019; 32: 65–77. https://doi.org/10.48550/arXiv.1907.05600

36.

Zhang

Zuo

, et al. Plug-and-play image restoration with deep denoiser prior. IEEE Trans Pattern Anal Mach Intell 2021; 44: 6360–6376.

37.

https://opendatalab.org.cn/OpenDataLab/LIDC-IDRI .

38.

Zhu

Zhang

Liang

, et al. Denoising diffusion models for plug-and-play image restoration. Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, Vancouver, BC, Canada, 2023, pp. 1219–1229. doi: https://doi.org/10.1109/CVPRW59228.2023.00129

39.

Zhang

Liu

, et al. Decoupling Training-Free Guided Diffusion by ADMM. arXiv 2024, 2411.12773.

40.

Zhao

Pan

Zhang

, et al. An oblique projection modification technique (OPMT) for fast multispectral CT reconstruction. Phys Med Biol 2021; 66: 065003.

41.

Wang

Bovik

Sheikh

, et al. Image quality assessment: from error visibility to structural similarity. IEEE Trans. Image Process 2004; 13: 600–612.

42.

Diederik

. Adam: A method for stochastic optimization, (No Title). 2014.

43.

Sidky

Pan

, et al. A robust method of x-ray source spectrum estimation from transmission measurements: Demonstrated on computer simulated, scatter-free transmission data. J Appl Phys 2005; 97(12): 4701–4713. https://doi.org/10.1063/1.1928312