PC-POCS Sampler to Reconstruct a Sparse-View Computer Tomography Image Consistent with Both the Prior and the Measurements

1. Introduction

With the advancement of machine learning, numerous methods have been proposed for reconstructing medical images from limited measurements [1,2,3,4,5,6]. Most of these approaches rely on supervised learning, where models are trained to directly map partial measurement data to reconstructed images using large datasets consisting of paired CT images and corresponding measurements. These measurements are typically generated by simulating a fixed physical model of the imaging process. However, when the measurement process changes — for example, when the number of CT projections is altered — it becomes necessary to regenerate the training dataset using the new acquisition settings and retrain the model from scratch. This process is both time-consuming and resource-intensive. Moreover, it limits the model’s ability to generalize across different acquisition conditions, and in some cases, may even lead to counterintuitive results, such as degraded performance with more measurements [7]. In this work, we address the aforementioned challenge by introducing unsupervised methods that eliminate the need for paired datasets during training, thus removing the dependency on a fixed measurement process. Our core idea is to model the prior distribution of medical images using a generative approach, which allows us to recover missing information resulting from incomplete measurements. Concretely, we train a score-based generative model[8,9] on a collection of high-quality medical images to serve as a data prior. Once trained, we employ a classifier-free [10] controllable sampling algorithm that generates images consistent with both the observed measurements and the learned prior, while incorporating the physical characteristics of the measurement process. Importantly, after training, the model can be flexibly applied to solve inverse problems within the same image domain, provided the measurement operator remains linear.

We evaluate the proposed method across multiple CT reconstruction tasks. Empirically, our approach achieves performance that is on par with or superior to existing unsupervised learning methods [11,12], even when all methods are trained and tested under the same measurement settings. Furthermore, when the measurement conditions are altered—for instance, by varying the number of projections in sparse-view CT—our method consistently outperforms all competing approaches. These empirical results collectively demonstrate that our method not only serves as a strong alternative to current unsupervised techniques in medical image reconstruction and artifact suppression but also shows great potential as a general-purpose framework for solving a wide range of inverse problems within the same image domain.

2. Background

2.1. Variational Diffusion Models

For many modalities, we can think of the data we observe as represented or generated by an associated unseen latent variable. The best intuition for expressing this idea is through Plato’s Allegory of the Cave. In this allegory, a group of people are chained inside a cave their entire life and can only see the two-dimensional shadows projected onto a wall in front of them, which are generated by unseen three-dimensional objects passed before a fire, as Figure 1 shows.

In generative modeling, observed data can be interpreted as a low-dimensional projection or manifestation of latent variables residing in a more complex, hidden space. These latent variables encapsulate high-level abstract features such as color, size, shape, and other semantic attributes. The objects we perceive in the real world can therefore be understood as realizations or projections of these latent concepts into observable form, analogous to how the cave dwellers perceive only shadows of unseen objects. While the individuals in the cave cannot observe the true sources directly, they may still reason about their existence. Likewise, in generative modeling, we aim to approximate these latent representations to understand and reconstruct the observed data.

2.1.1. Variational Autoencoders

In the default formulation of the Variational Autoencoder (VAE)[13,14], we directly maximize the ELBO. This approach is variational, because we optimize for the best $q_{\varphi }\left(z\right|x)$ amongst a family of potential posterior distributions parameterized by $\varphi$. It is called an autoencoder because it is reminiscent of a traditional autoencoder model, where input data is trained to predict itself after undergoing an intermediate bottlenecking representation step. To make this connection explicit, we rewrite the ELBO term:

$$\begin{array}{cc}\hfill & {\mathbb{E}}_{q_{\varphi }\left(z\right|x)}\left[ln{\displaystyle \frac{p(x,z)}{q_{\varphi }\left(z\right|x)}}\right]\hfill \\ \hfill & ={\mathbb{E}}_{q_{\varphi }\left(z\right|x)}\left[ln{p}_{\theta }\left(x\right|z)\right]-D_{KL}(q_{\varphi }\left(z\right|x)\left|\right|p\left(z\right))\hfill \end{array}$$

(1)

The first term is the reconstruction term, and the second term is the prior matching term. In this case, we learn an intermediate bottlenecking distribution $q_{\varphi }\left(z\right|x)$ that can be treated as an encoder; it transforms inputs into a distribution over possible latents. Simultaneously, we learn a deterministic function $p_{\theta }\left(x\right|z)$ to convert a given latent code z into an observation x, which can be interpreted as a decoder.

The KL divergence term (prior matching term) could be computed analytically, and the reconstruction term can be approximated using a Monte Carlo estimate. With this, we rewrite the objective function as:

$$\begin{array}{cc}\hfill & {\mathbb{E}}_{q_{\varphi }\left(z\right|x)}\left[ln{p}_{\theta }\left(x\right|z)\right]-D_{KL}(q_{\varphi }\left(z\right|x)\left|\right|p\left(z\right))\hfill \\ \hfill & ={\displaystyle \frac{1}{L}}{\Sigma }_{l=1}^{L}ln{p}_{\theta }\left(x\right|z^{\left(l\right)})-D_{KL}(q_{\varphi }\left(z\right|x)\left|\right|p\left(z\right))\hfill \end{array}$$

(2)

where latents ${\left\{z^{\left(l\right)}\right\}}_{l=1}^L$ are sampled from $q_{\varphi }\left(z\right|x)$, for every observation x in the dataset. After that, we use the reparameterization trick to conquer the non-differentiable problem, eventually getting:

$$z={\mu }_{\varphi }\left(x\right)+{\sigma }_{\varphi }\left(x\right)\ast ϵ,ϵ\sim \mathcal{N}(0,I)$$

(3)

In a VAE, each z is thus computed as a deterministic function of input x and auxiliary noise variable $ϵ$. Under this reparameterized version of z, gradients can then be computed for $\varphi$ as desired. See Figure 2.

2.1.2. Hierarchical Variational Autoencoders

A Hierarchical Variational Autoencoder (HVAE)[15,16] (See Figure 3) is a generalization of a VAE that extends to multiple hierarchies over latent variables. Under this formulation, latent variables themselves are interpreted as generated from other higher-level, more abstract latents. Whereas in the general HVAE with T hierarchical levels, each latent is conditioned on all previous latents, in this work, we focus on a special case, which we call a Markovian HVAE (MHVAE). In a MHVAE, the generative process is a Markov chain, i.e., each transition down the hierarchy is Markovian, where decoding each latent $z_t$ only conditions on the previous latent $z_{t+1}$. Mathematically, we represent the joint distribution and the posterior of a Markovian HVAE as:

$$\begin{array}{cc}\hfill p(x,z_{1:T})& =p\left(z_T\right)p_{\theta }\left(x\right|z_1){\Pi }_{t=2}^T{p}_{\theta }\left(z_{t-1}\right|z_t)\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill q_{\varphi }\left(z_{1:T}\right|x)& =q_{\varphi }\left(z_1\right|x){\Pi }_{t=2}^T{q}_{\varphi }\left(z_t\right|z_{t-1})\hfill \end{array}$$

(5)

Substituting the original ELBO expression, we get the new objective:

$$\begin{array}{cc}\hfill & {\mathbb{E}}_{q_{\varphi }\left(z_{1:T|x}\right)}\left[ln{\displaystyle \frac{p(x,z_{1:T})}{q_{\varphi }\left(z_{1:T|x}\right)}}\right]\hfill \\ \hfill & ={\mathbb{E}}_{q_{\varphi }\left(z_{1:T}\right|x)}\left[ln{\displaystyle \frac{p\left(z_T\right)p_{\theta }\left(x\right|z_1){\Pi }_{t=2}^T{p}_{\theta }\left(z_{t-1}\right|z_t)}{q_{\Phi }\left(z_1\right|x){\Pi }_{t=2}^T{q}_{\Phi }\left(z_t\right|z_{t-1})}}\right]\hfill \end{array}$$

(6)

2.1.3. Unconditional Score-based Generative Model

The Variational diffusion model[17,18,19] formulation has an equivalent Score-based Generative Modeling formulation. Arbitrary probability distributions can be written in the form:

$$p_{\theta }\left(x\right)={\displaystyle \frac{1}{Z_{\theta }}}{e}^{-f_{\theta }\left(x\right)}$$

(7)

where $f_{\theta }\left(x\right)$ is the energy function, often modeled by a neural network, and $Z_{\theta }$ is a normalizing constant to ensure that $\int p_{\theta }\left(x\right)dx=1$. One way to avoid modeling the normalization constant is to use a neural network $s_{\theta }\left(x\right)$ to learn the score function $\nabla lnp\left(x\right)$ of distribution $p\left(x\right)$ instead. This comes from taking the derivative of the ln of both sides of the equation above:

$$\begin{array}{cc}\hfill \nabla ln{p}_{\theta }\left(x\right)& =\nabla ln\left({\displaystyle \frac{1}{Z_{\theta }}}{e}^{-f_{\theta }\left(x\right)}\right)\hfill \\ \hfill & =\nabla ln{\displaystyle \frac{1}{Z_{\theta }}}+\nabla ln{e}^{-f_{\theta }\left(x\right)}\hfill \\ \hfill & =-\nabla f_{\theta }\left(x\right)\hfill \\ \hfill & =s_{\theta }\left(x\right)\hfill \end{array}$$

(8)

The equation above could be represented as a neural network without involving any normalization constants. The score model could be optimized by minimizing the Fisher Divergence with the ground truth score function:

$${\mathbb{E}}_{p\left(x\right)}\left[\right||s_{\theta }\left(x\right)-{\nabla lnp\left(x\right)\left|\right|}_2^2]$$

(9)

The score function defines a vector field over the entire space where the data x lies, pointing towards the modes. Then, by learning the score function of the true data distribution, we can generate samples by starting at an arbitrary point in the same space and iteratively following the score until a mode is reached. This sampling procedure is known as Langevin dynamics, and is mathematically described as:

$$x_{i+1}\leftarrow x_i+c\nabla lnp\left(x_i\right)+\sqrt{2c}ϵ,ϵ\sim \mathcal{N}(0,I)$$

(10)

The initial value $x_0$ is randomly sampled from a prior distribution. Eventually, we can choose a positive sequence of noise scale ${\left\{{\sigma }_t\right\}}_{t=1}^T$ and then define a sequence of progressively perturbed data distributions:

$$p_{{\sigma }_t}\left(x_t\right)=\int p\left(x\right)\mathcal{N}\left(x_t\right|x,{\sigma }_t^{2}I)dx$$

(11)

Then, a neural network $s_{\theta }(x,t)$ is learned using denoising score matching to learn the score function for all noise levels:

$${\Sigma }_{t=1}^T{\mathbb{E}}_{p_{{\sigma }_t}\left(x_t\right)}\left[\right||s_{\theta }(x,t)-\nabla ln{p}_{{\sigma }_t}\left(x_t\right){\left|\right|}_2^2]$$

(12)

where $\lambda \left(t\right)$ is a positive weighting function that conditions on noise level t. Furthermore, Song[9] proposes annealed Langevin Dynamics sampling as a generative procedure. The initialization is chosen from some hand-made prior, and each subsequent sampling step starts from the final samples of the previous simulation. This is analogous to the sampling procedure performed in the Markovian HVAE interpretation of a Variational Diffusion Model, where a randomly initialized data vector is iteratively refined over decreasing noise levels. Under an infinite number of noise scales, the perturbation of an image over continuous time can be represented as a stochastic process, and therefore described by a stochastic differential equation (SDE). Sampling is then performed by reversing the SDE, which naturally requires estimating the score function at each continuous-valued noise level. Different parameterizations of the SDE essentially describe different perturbation schemes over time. In this work, we use Variance Exploring (VE) SDE (See Figure 4 as a reference).

2.1.4. Classifier-Free Guidance

Because learning two separate diffusion models is expensive, in this work, we just consider the classifier-free guidance[10,20,21]. To derive the score function under Classifier-Free Guidance, first, we have:

$$\nabla lnp\left(y\right|x_t)=\nabla lnp\left(x_t\right|y)-\nabla lnp\left(x_t\right)$$

(13)

Then, substituting this into the classifier-guidance equation $\nabla lnp\left(x_t\right|y)=\nabla lnp\left(x_t\right)+\gamma \nabla lnp\left(y\right|x_t)$, get:

$$\nabla lnp\left(x_t\right|y)=\gamma \nabla lnp\left(x_t\right|y)+(1-\gamma )\nabla lnp\left(x_t\right)$$

(14)

The first term on the right-hand is the conditional score, and the second term is the unconditional score. $\gamma$ is a term that controls how much the learned conditional model cares about the conditioning information. When $\gamma =0$, the learned conditional model completely ignores the conditioner and makes an unconditional inference; when $\gamma =1$, the model explicitly learns the vanilla conditional distribution without guidance; when $\gamma >1$, the diffusion model not only prioritizes the conditional score function, but also moves in the direction away from the unconditional score function.

2.2. Sparse-View Computer Tomography Inverse Problems

An inverse problem aims to recover an unknown signal from a set of observed measurements. Formally, let $\mathbf{x}\in {\mathbb{R}}^n$ denote the unknown signal, and let $\mathbf{y}\in {\mathbb{R}}^m$ represent the observed measurements modeled by the linear relation $\mathbf{y}=\mathbf{A}\mathbf{x}+ϵ$, where $\mathbf{A}\in {\mathbb{R}}^{m\times n}$ is the measurement operator and $ϵ\in {\mathbb{R}}^m$ is the additive noise. The task of solving the inverse problem involves estimating $\mathbf{x}$ from the noisy observation $\mathbf{y}$. When $m<n$, the system is underdetermined, making the problem ill-posed unless additional prior knowledge about $\mathbf{x}$ is introduced. To address this, we assume that $\mathbf{x}$ follows a prior distribution $p\left(\mathbf{x}\right)$. In the probabilistic formulation, the relationship between measurements and the signal is captured by the likelihood function $p\left(\mathbf{y}\right|\mathbf{x})=q_{ϵ}(\mathbf{y}-\mathbf{A}\mathbf{x})$, where $q_{ϵ}$ denotes the noise distribution. Given both the likelihood $p\left(\mathbf{y}\right|\mathbf{x})$ and the prior $p\left(\mathbf{x}\right)$, the inverse problem can be solved by inferring the posterior distribution $p\left(\mathbf{x}\right|\mathbf{y})$ and sampling from it.

Figure 5. Some unconditional samples. They look very realistic, but in fact, they have no reasonable anatomical structure. This is precisely the reason why guiding signals are needed.

Figure 5. Some unconditional samples. They look very realistic, but in fact, they have no reasonable anatomical structure. This is precisely the reason why guiding signals are needed.

3. Methodology

3.1. Projection Onto Convex Sets (POCS)

In mathematics, projection onto convex sets (POCS), sometimes known as the alternating projection method, is a method to find a point in the intersection of two closed convex sets. The case when the sets are affine spaces is special, since the iterates not only converge to a point in the intersection, assuming the intersection is non-empty, but to the orthogonal projection of the point onto the intersection. POCS always has the following iteration form:

$$x_{k+1}=P_{prior}\left(P_{data}\left(x_k\right)\right)$$

(15)

where $P_{data}$ denotes projecting data x onto the data fidelity set, and $P_{prior}$ denotes projecting data x onto the prior set.

3.2. Imputation in the Radon Domain

Imputation is a special case of conditional sampling. The known and unknown dimensions of $\mathbf{y}$ are denoted by $\Omega \left(\mathbf{y}\right)$ and $\overline{\Omega }\left(\mathbf{y}\right)$ respectively, and let ${\mathbf{f}}_{\overline{\Omega }}(·,t)$ and $g_{\overline{\Omega }}(·,t)$ denote the drift coefficient and the diffusion coefficient restricted to the unknown dimensions respectively. For VE SDEs, the drift coefficient $\mathbf{f}(·,t)$ is element-wise, and the diffusion coefficient $g(·,t)$ is diagonal. So, when $\mathbf{f}(·,t)$ is element-wise, ${\mathbf{f}}_{\overline{\Omega }}(·,t)$ denotes the same element-wise function applied only to the unknown dimensions; when $g(·,t)$ is diagonal, $g_{\overline{\Omega }}(·,t)$ denotes the sub-matrix restricted to unknown dimensions.

In this setting, $P_{prior}$ specifically denotes the samples transformed by the forward model from the diffusion model at each timestep, and $P_{prior}$ denotes the constraint coming from the forward model. Because the Radon transform operator $\mathbf{A}$ is not orthogonal, we need an orthogonal term ${\left(\mathbf{A}{\mathbf{A}}^T\right)}^{†}$ analogous to the filtering step to implement the orthogonal projection, to make POCS work well.

3.3. PC-POCS Sampling

Combining POCS with the PC sampling method, we derive the PC-POCS algorithm. Specifically, we use the discretized timestep, from $t=1$ to $t=ϵ$, each timestep conducts two steps of update. The first step is used to predict the diffusion trajectory:

$$x_t^{\prime }=Predictor(x_t,t)$$

(16)

which always has the form of

$$x_{t-\Delta t}=x_t+\Delta t\times \mathbf{f}(x_t,t)+\sqrt{2g^2\left(t\right)\Delta t}\times ϵ$$

(17)

The second step is used to correct the trajectory. The first part of the second step is Langevin correct:

$$x_t^{″}=x_t^{\prime }+\alpha \times \nabla ln{p}_t\left(x_t^{\prime }\right)+z,z\sim \mathcal{N}(0,I)$$

(18)

After the corrector, come into the Kaczmarz iteration to satisfy the projection consistency:

$$x\leftarrow x+\lambda \times {\displaystyle \frac{{\left(\mathbf{A}{\mathbf{A}}^T\right)}^{†}(\mathbf{y}-\mathbf{A}\mathbf{x})}{c}}$$

(19)

where $\mathbf{A}$ denotes the forward radon projection operator, ${\mathbf{A}}^T$ denotes the back projection operator, ${\left(\mathbf{A}{\mathbf{A}}^T\right)}^{†}$ denotes the filter back projection operator. $\mathbf{y}$ denotes the measurement, i.e., sparse-view sinogram, and $\lambda$ denotes the stepsize scheduler coefficient, $c={\mathbf{A}}^T\left(\mathbf{A}\left(I\right)\right)$ denotes the normalization constant. The details of the whole algorithm are in the pseudocode.

Algorithm 1 PC-POCS

Require:  $\mathbf{y},\mathbf{A},{\left\{{\sigma }_i\right\}}_{i=1}^N,s_{\theta },{\left\{\lambda \right\}}_{i=1}^N$ 1: ${\mathbf{x}}_N\sim \mathcal{N}(0,I)$ 2: $c={\mathbf{A}}^T\mathbf{A}\left(I\right)$ 3: for $i=N$ to 1 do 4:       x← Predictor($s_{\theta }$, $x_t$, t) 5:       x← Corrector($s_{\theta }$, $x_t$, t) 6:       ${\lambda }_t\leftarrow {\left\{\lambda \right\}}_{i=1}^{\mathbf{N}}$ 7:       $x\leftarrow x+{\lambda }_t\times {\displaystyle \frac{{\left(\mathbf{A}{\mathbf{A}}^T\right)}^{†}(\mathbf{y}-\mathbf{A}\mathbf{x})}{c}}$ 8: end for 9: return x

4. Experiment

For all tasks, we aim to verify the superiority of our method against other diffusion model-based approaches.

4.1. Dataset

We take the dataset from the AAPM 2016 CT low-dose grand challenge, where the data are acquired in a fan-beam geometry with varying parameters. The data preparation steps follow those of [22]. From the helical cone beam projections, an approximation to fanbeam geometry is performed via a single-slice rebinning technique. Reconstruction is then performed via standard filtered backprojection (FBP), where the reconstructed axial images have a matrix size of $512\times 512$. We resize the axial slices to have the size $256\times 256$, and use these slices to train the score function. The whole dataset consists of 9 volumes (3142 slices) of training data, and 1 volume (823 slices) of testing data. To generate sparse-view measurements, we retrospectively employ the parallel-view geometry for simplicity.

4.2. Details of Network Training

For the CT score function, we train the ncsnpp network with setting $\lambda ={\sigma }^2\left(t\right)$, and $ϵ={10}^{-5}$. Training is conducted by using two RTX 3080 Ti for 100 epochs, and the average loss eventually converges to $2\times {10}^{-2}$.

4.3. Evaluation Metrics

We use three standard metrics: SSIM, PSNR, and FSIM for quantitative evaluation. FSIM may be a relatively novel metric, and makes an improvement on the SSIM, introducing phase consistency with position as a variable and weighting the SSIM based on this. The basic formula of FSIM is:

$$PC\left(x\right)={\mathbf{max}}_{\overline{\varphi }\left(x\right)\in [0,2\pi ]}{\displaystyle \frac{{\Sigma }_n{A}_{n}cos({\varphi }_n\left(x\right)-\overline{\varphi }\left(x\right))}{{\Sigma }_n{A}_n}}$$

(20)

The equation above expresses the image signal as a Fourier series expansion, where each component has a corresponding phase. $\overline{\varphi }\left(x\right)$ represents the average phase of all components, ${\varphi }_n\left(x\right)$ denotes the phase of each individual component, and $A_n$ is the amplitude.

4.4. Comparative Experiment with the State-of-the-Art Methods

To the best of our knowledge, [11,12] are the only two methods that solve CT reconstruction directly with diffusion models. So, we compare our method against theirs. We call the method from [11] as Song, and the approach from [12] as MCG. To control the variables, the experiments of the three methods were all conducted under the condition of sparsity of 6 (i.e., the number of projections was 30). The numerical results of the experiment are shown in Table 1, and the visualization results are shown in Figure 6.

From Table 1, we could see that the PC-POCS method outperforms the other two methods in terms of all three metrics, whether the max value or the second best value among all the slices of patient L096 shown. In addition, we also presented some reconstruction results to provide a more intuitive feeling. See supplement material. Visually, under the sparsity of 6 (30 projections), the PC-POCS method has the best visual effect.

4.5. Ablation Study

For different methods, determining the lower bound of sparsity is a very important task. Generally speaking, the higher the sparsity, the worse the quality of the reconstructed image. For instance, taking SSIM as the standard, the SSIM of reliably reconstructed images in medicine is generally greater than $0.7$. Reconstructed images below this lower bound, due to the loss of most of the true structure, will not have clinical value.

Due to this consideration, we conducted an ablation experiment to determine the sparsity corresponding to the minimum trusted image. The relationship between the sparsity and the number of projections is:

$$180=\mathrm{sparsity}\times \mathrm{number}\mathrm{of}\mathrm{projections}$$

(21)

From Figure 8, the relationship between the sparsity (number of projections) and metrics is obvious. For the integrity of the experiment, we also conducted experiments on FSIM and PSNR, and provided the results. According to Figure 8, when $\mathrm{sparsity}=26$, the SSIM reaches the lower bound of 0.7.

For the sake of demonstrating the performance of our method, we conduct an ablation study on the other two methods. The experiment result is shown in Figure 9. Besides that, the visualization results of this ablation study are also provided in Figure 7. Our proposed method needs less physical information than the other two methods and has the highest threshold sparsity for reconstruction.

Figure 7. The ablation study result of slice 290 of patient L096. From left to right, the images correspond to different sparsity levels. From upwards to downwards, the images correspond to ground truth, PC-POCS, MCG, and Song, respectively. For PC-POCS, the sparsity is set to $6,8,10,12,18,30$. Because the quality of the image is very low under the method of Song, we just set the sparsity of $2,3,5,6,10$ for reference. We find that under high sparsity, some subtle structure is lost, for example, the part in the white border is different in the ground truth and the reconstruction under high sparsity. You can zoom in to see more details of the image.

Figure 7. The ablation study result of slice 290 of patient L096. From left to right, the images correspond to different sparsity levels. From upwards to downwards, the images correspond to ground truth, PC-POCS, MCG, and Song, respectively. For PC-POCS, the sparsity is set to $6,8,10,12,18,30$. Because the quality of the image is very low under the method of Song, we just set the sparsity of $2,3,5,6,10$ for reference. We find that under high sparsity, some subtle structure is lost, for example, the part in the white border is different in the ground truth and the reconstruction under high sparsity. You can zoom in to see more details of the image.

Figure 8. The relationship between the sparsity and each metric. With the increase of sparsity, the number of projections decreases, so the true physical information decreases, and the quality of the reconstruction images deteriorates. According to the plot, when sparsity = 26, the SSIM reaches the threshold of 0.7.

Figure 8. The relationship between the sparsity and each metric. With the increase of sparsity, the number of projections decreases, so the true physical information decreases, and the quality of the reconstruction images deteriorates. According to the plot, when sparsity = 26, the SSIM reaches the threshold of 0.7.

Figure 9. MCG(left) and Song(right)’s ablation study results. The threshold sparsity is 18, 10, respectively. The bigger the numerical value is, the more physical information the method needs. So from this perspective, our proposed method also outperforms the other two methods.

Figure 9. MCG(left) and Song(right)’s ablation study results. The threshold sparsity is 18, 10, respectively. The bigger the numerical value is, the more physical information the method needs. So from this perspective, our proposed method also outperforms the other two methods.

5. Conclusion

In this work, we proposed a general framework that can greatly enhance the performance of the score-based generative model’s solver for solving the low-dose computer tomography inverse problem. We show that our method can outperform the current state-of-the-art method. In addition, we conduct an ablation study on the proposed method, the result shows that our method could work under high sparsity with a reliable reconstruction result, and also has a stronger ability than the other two methods.

5.1. Limitations

The proposed method is intrinsically stochastic, as it fundamentally relies on the diffusion model. When the measurement dimension m is significantly reduced, the method occasionally fails to yield high-quality reconstructions, although it still generally outperforms competing approaches.

Furthermore, we acknowledge that our method is relatively slow in generating samples, as it inherits the limitations of diffusion models. This limitation could potentially be mitigated by incorporating recent solvers designed to accelerate diffusion model inference.

5.2. Broader Impact

Aligned with the concerns associated with generative model-based approaches to inverse problems, our method’s strong reliance on the underlying diffusion model raises the potential for generating harmful content, such as deepfakes. In addition, it may amplify pre-existing social biases embedded in the training dataset.