# Score-Based Generative Models for PET Image Reconstruction

Imraj R D Singh  https://orcid.org/0000-0003-2186-0977
<EMAIL_ADDRESS>
Department of Computer Science, University College London, 66-72 Gower St,
WC1E 6EA, London, United Kingdom. Alexander Denker11footnotemark: 1
https://orcid.org/0000-0002-7265-261X<EMAIL_ADDRESS>
Center for Industrial Mathematics, University of Bremen, Bibliothekstr. 5,
28359 Bremen, Germany. Riccardo Barbano11footnotemark: 1
https://orcid.org/0000-0003-1863-2092<EMAIL_ADDRESS>
Department of Computer Science, University College London, 66-72 Gower St,
WC1E 6EA, London, United Kingdom. Željko Kereta
https://orcid.org/0000-0003-2805-0037<EMAIL_ADDRESS>
Department of Computer Science, University College London, 66-72 Gower St,
WC1E 6EA, London, United Kingdom. Bangti Jin
https://orcid.org/0000-0002-3775-9155<EMAIL_ADDRESS>
Department of Mathematics, The Chinese University of Hong Kong, Shatin, N.T.,
Hong Kong. Kris Thielemans  https://orcid.org/0000-0002-5514-199X
<EMAIL_ADDRESS>
Institute of Nuclear Medicine, University College London, London, United
Kingdom. Peter Maass  https://orcid.org/0000-0003-1448-8345 pmaass@uni-
bremen.de
Center for Industrial Mathematics, University of Bremen, Bibliothekstr. 5,
28359 Bremen, Germany. Simon Arridge  https://orcid.org/0000-0003-1292-0210
<EMAIL_ADDRESS>
Department of Computer Science, University College London, 66-72 Gower St,
WC1E 6EA, London, United Kingdom Equal contribution.

###### Abstract

Score-based generative models have demonstrated highly promising results for
medical image reconstruction tasks in magnetic resonance imaging or computed
tomography. However, their application to Positron Emission Tomography (PET)
is still largely unexplored. PET image reconstruction involves a variety of
challenges, including Poisson noise with high variance and a wide dynamic
range. To address these challenges, we propose several PET-specific
adaptations of score-based generative models. The proposed framework is
developed for both 2D and 3D PET. In addition, we provide an extension to
guided reconstruction using magnetic resonance images. We validate the
approach through extensive 2D and 3D in-silico experiments with a model
trained on patient-realistic data without lesions, and evaluate on data
without lesions as well as out-of-distribution data with lesions. This
demonstrates the proposed method’s robustness and significant potential for
improved PET reconstruction.

Keywords: Positron emission tomography, score-based generative models, image
reconstruction

##  1 Introduction

Positron Emission Tomography (PET) is a functional medical imaging technique
for quantifying and visualising the distribution of a radio-tracer within the
body, and is vital in clinical practice for accurate diagnosis, treatment
planning, and monitoring of diseases. In a PET scan, radio-tracers are
injected to interrogate a specific biological pathway of interest. Through the
decay of the radio-tracer a positron is released, which upon annihilating with
an electron produces a pair of coincident photons that travel in approximately
anti-parallel directions. These emitted photons are detected and are then used
to reconstruct the underlying radio-tracer distribution. The relationship
between the measured emissions and the radio-tracer can be approximated with
the Poisson noise model as

${\mathbf{y}}\sim{\text{Pois}}(\mathbf{\bar{y}}),\quad\mathbf{\bar{y}}={\mathbf{A}}{\mathbf{x}}+\mathbf{\bar{b}},$
(1)

where $\mathbf{\bar{y}}\in\mathbb{R}^{m}$ is the expected value of the
measurements ($m$ is the number of detector bins) and
${\mathbf{x}}\in\mathbb{R}^{n}$ is the discrete (voxel) basis representation
of the tracer distribution ($n$ is the number of voxels). The system matrix
${\mathbf{A}}\in\mathbb{R}^{m\times n}$ includes approximate line integrals
between detectors as well as physical phenomena such as photon attenuation,
positron range, and detector sensitivity. It should be noted that 3D
measurements detect pairs of photons between detector rings, i.e. they are not
a stack of 2D measurements. The expected background
$\mathbf{\bar{b}}\in\mathbb{R}^{m}$ are estimates of scatter and randoms
events Qi and Leahy (2006). The unique challenges that distinguish PET from
other imaging modalities, e.g. Magnetic Resonance Imaging (MRI) and Computed
Tomography (CT), include Poisson noise with low mean number of counts, and
widely varying dynamic range of images due to functional differences between
patients.

Most inverse problems in imaging are ill-posed, in the sense that the solution
may not exist, not be unique, or not depend continuously on the measurement
noise (Engl et al., 1996; Ito and Jin, 2015). To stabilise the reconstruction
process, prior knowledge is often leveraged through a penalising functional
that promotes solutions from a desirable image subset. The priors are
typically hand-crafted to promote desired features in the reconstructed image,
such as sparsity of edges (Rudin et al., 1992) or smoothness. Furthermore, if
an additional image is available, e.g. with high resolution structural
information, a suitable prior can promote common features between the two
images (Ehrhardt, 2021). This is often referred to as guided reconstruction.
In recent years, deep learning approaches have shown state-of-the-art
performance in PET image reconstruction, see surveys (Reader et al., 2021;
Pain et al., 2022). Existing approaches include post-processing (Kaplan and
Zhu, 2019), to synthesise high-dose images from low-dose ones (which is akin
to denoising), and deep unrolled optimisation (Mehranian and Reader, 2021;
Guazzo and Colarieti-Tosti, 2021). However, these supervised approaches
require large volumes of paired data that is often hard to acquire.

In contrast, generative models are unsupervised, requiring a dataset only of
images of the target domain. These can, for example, be high-quality
reconstructions acquired from prior scans. The aim of generative modelling is
to approximate the image manifold of a given dataset (Bengio et al., 2013).
There are a variety of methods for this task, e.g. generative adversarial
networks (Goodfellow et al., 2014), variational autoencoders (Kingma and
Welling, 2014) and recently Score-based Generative Models (SGMs), which aim to
generate high-quality samples, sample quickly, and have adequate mode coverage
(Xiao et al., 2022). Over recent years, SGMs have become the de facto method
for image generation due to the quality and diversity of generated images
(Dhariwal and Nichol, 2021). Generative models can be integrated into the
reconstruction process as data-driven priors that are independent of the
forward model, cf. (Dimakis, 2022). This modularity separates the classical
forward modelling problem from the generative image modelling problem. Where
the generative model is trained to generate images specific to the task at
hand, i.e. PET images of brains. The model can then be used across scanners
and noise levels given the learnt image manifold is still relevant.

SGMs have been applied to CT and MR image reconstruction (Song et al., 2022).
These reconstructions condition the SGM image generation on measurements, and
balance the consistency with measurements versus consistency with the SGM
learnt image manifold. From this perspective the SGM acts as a prior (Kobler
and Pock, 2023). There are different methods to enforce measurement
consistency of the reconstructions, which can be broadly classified into
gradient based methods (Jalal et al., 2021; Chung et al., 2023a) and
projection based methods (Song et al., 2022; Chung and Ye, 2022; Chung et al.,
2023b). Recently, denoising diffusion models (discrete variants of SGMs) were
used for PET image denoising (Gong et al., 2022). Instead, our work focuses on
PET image reconstruction, and we present the following contributions:

  * •

We develop a novel algorithmic framework building upon SGMs that carefully
addresses the challenges inherent to PET. To do so, we modify the conditional
sampling method (Chung et al., 2023b; Zhu et al., 2023), recently proposed for
inverse problems with Gaussian noise, for PET image reconstruction. This is
achieved with a penalised Maximum A Posteriori (MAP) estimator computed with
an accelerated algorithm that evaluates subsets of the measurements.

  * •

We leverage additional MR images to enhance the proposed framework, leading to
improved image quality that better agrees with the measured data.

  * •

We scale the approach to 3D PET reconstruction.

The proposed method is tested on multiple noise levels, radio-tracers, and in
both 2D and 3D settings with an SGM trained on patient-realistic BrainWeb data
without lesions (Collins et al., 1998). In addition to data without lesions,
we test on out-of-distribution (OOD) data with lesions to validate method
robustness. The rest of the paper is structured as follows. In Section 2 we
provide the background on PET reconstruction and SGMs. In particular, we
present different methods for using SGMs in image reconstruction. In Section 3
we propose modifications needed to apply SGMs for PET reconstruction. We
describe the experimental setting in Section 4, and present and discuss the
results in Section 5. The code is publicly available at
https://github.com/Imraj-Singh/Score-Based-Generative-Models-for-PET-Image-
Reconstruction

##  2 Background

###  2.1 Fundamentals of Positron Emission Tomography Reconstruction

PET measurements are the result of a low-count photon-counting process. The
true forward process, from tracer-distribution to photon detection, is
approximated by the forward model defined in Eq. (1). The likelihood of the
measured photon counts, for an unknown tracer distribution, can be modelled by
an independent Poisson distribution. One of the first methods developed for
estimating the tracer distribution through a Poisson model was maximum
likelihood. This selects an image ${\mathbf{x}}\in\mathbb{R}_{\geq 0}^{n}$ by
maximising the Poisson Log-Likelihood (PLL) function, given by

$L({\mathbf{y}}|{\mathbf{x}})=\sum_{i=1}^{m}y_{i}\log([{\mathbf{A}}{\mathbf{x}}+\mathbf{\bar{b}}]_{i})-[{\mathbf{A}}{\mathbf{x}}+\mathbf{\bar{b}}]_{i}-\log(y_{i}!).$
(2)

By maximising the PLL, the Maximum Likelihood Estimate (MLE) is obtained. A
particularly important algorithm for computing the MLE is Expectation
Maximisation (EM) (Shepp and Vardi, 1982). However, due to its slow
convergence, acceleration is sought through splitting the PLL into a sum of
$n_{\mathrm{sub}}\geq 1$ sub-objectives. This gives rise to the greatly sped-
up Ordered Subset Expectation Maximisation (OSEM) (Hudson and Larkin, 1994)
algorithm. Because of the ill-conditioning of PET reconstruction, the MLE
tends to overfit to measurement noise. To address the ill-conditioning and
improve reconstruction quality, it is common practice to regularise the
reconstruction problem via the use of an image-based prior. This gives rise to
the MAP objective

$\Phi(\mathbf{x})=L(\mathbf{y}|\mathbf{x})+\lambda R(\mathbf{x}),$ (3)

where $R({\mathbf{x}})$ is the log of a chosen image-based prior with penalty
strength $\lambda$. Block-Sequential Regularised Expectation Maximisation
(BSREM) (Pierro and Yamagishi, 2001; Ahn and Fessler, 2003) is an iterative
algorithm, globally convergent under mild assumptions, that applies the subset
idea of OSEM to the MAP objective. For
$\Phi({\mathbf{x}})=\sum_{j=1}^{n_{\mathrm{sub}}}\Phi_{j}({\mathbf{x}})$,
where $\Phi_{j}$ is the sub-objective
$\Phi_{j}({\mathbf{x}})=L_{j}({\mathbf{y}}|{\mathbf{x}})+\lambda
R({\mathbf{x}})/n_{\mathrm{sub}}$, and $L_{j}$ is the likelihood for a subset
of the measurements. The BSREM update iterations are given by

$\mathbf{x}^{i+1}=P_{\mathbf{x}\geq
0}\left[\mathbf{x}^{i}+\alpha_{i}\mathbf{D}(\mathbf{x}^{i})\nabla\Phi_{j}(\mathbf{x}^{i})\right]\quad
i\geq 0,$ (4)

where $P_{\mathbf{x}\geq 0}[\cdot]$ denotes the non-negativity projection, $i$
is the iteration number, and index $j=(i\mod n_{\mathrm{sub}})+1$ cyclically
accesses sub-objectives. The preconditioner is
$\mathbf{D}({\mathbf{x}}^{i})=\mathrm{diag}\left\\{{\max({\mathbf{x}}^{i},\delta)}/{\mathbf{A}^{\top}\mathbf{1}}\right\\}$,
where $\delta$ is a small positive constant to ensure positive definiteness,
$\mathbf{A}^{\top}\mathbf{1}$ is referred to as the sensitivity image, and
$\mathbf{A}^{\top}$ the matrix transpose. The step-sizes are
$\alpha_{i}=\alpha_{0}/(\zeta\lfloor i/n_{\mathrm{sub}}\rfloor+1)$, where
$\alpha_{0}=1$ and $\zeta$ is a relaxation coefficient. A common regulariser
for PET reconstruction is the Relative Difference Prior (RDP) (Nuyts et al.,
2002), see Appendix A.3 for details. The gradient of the RDP is scale-
invariant as it is computed using the ratio of voxel values. This partially
overcomes the issue with the wide dynamic range observed in emission
tomography images, helping to simplify the choice of the penalty strength
across noise levels.

###  2.2 Score-based Generative Models

SGMs have emerged as a state-of-the-art method for modelling, and sampling
from, high-dimensional image distributions (Song et al., 2021c). They
reinterpret denoising diffusion probabilistic modelling (Sohl-Dickstein et
al., 2015; Ho et al., 2020) and score-matching Langevin dynamics (Song and
Ermon, 2019) through the lens of Stochastic Differential Equations (SDE). SGMs
are often formulated by prescribing a forward diffusion process defined by an
Itô SDE

$d{\mathbf{x}_{t}}={\mathbf{f}}({\mathbf{x}_{t}},t)dt+g(t)d\mathbf{w}_{t},\quad{\mathbf{x}}_{0}\sim
p_{0}:={\pi},$ (5)

where $\\{{\mathbf{x}}_{t}\\}_{t\in[0,T]}$ is a stochastic process indexed by
time $t$ and ${\pi}$ is the image distribution. Each random vector
${\mathbf{x}}_{t}$ has an associated time-dependent density
$p({\mathbf{x}}_{t})$. To emphasise that the density is a function of $t$ we
write $p_{t}({\mathbf{x}}_{t}):=p({\mathbf{x}}_{t})$. The multivariate Wiener
process $\\{{\mathbf{w}}_{t}\\}_{t\geq 0}$ is the standard Brownian motion.
Starting at the image distribution ${\pi}$, the drift function
${\mathbf{f}}(\cdot,t):{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ and the diffusion
function $g:{\mathbb{R}}^{n}\to{\mathbb{R}}$ are chosen such that the terminal
distribution at $t=T$ approximates the standard Gaussian,
$p_{T}\approx\mathcal{N}(0,I)$. Thus, the forward diffusion process maps the
image distribution ${\pi}$ to a simple, tractable distribution. The aim of
SGMs is to invert this process, i.e. start at the Gaussian distribution and go
back to the image distribution ${\pi}$. Under certain conditions on
${\mathbf{f}}$ and $g$, a reverse diffusion process can be defined (Anderson,
1982)

$d{\mathbf{x}_{t}}=[{\mathbf{f}}({\mathbf{x}_{t}},t)-g(t)^{2}\nabla_{\mathbf{x}}\log{p_{t}}({\mathbf{x}_{t}})]dt+g(t)d\bar{\mathbf{w}}_{t},$
(6)

that runs backwards in time. The Wiener process
$\\{\bar{\mathbf{w}}_{t}\\}_{t\geq 0}$ is time-reversed Brownian motion, and
the term $\nabla_{\mathbf{x}}\log{p_{t}}({\mathbf{x}_{t}})$ is the score
function. Denoising Score Matching (DSM) (Vincent, 2011) provides a
methodology for estimating
$\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{x}_{t}})$ by matching the
transition densities $p_{t}({\mathbf{x}_{t}}|{\mathbf{x}_{0}})$ with a time-
conditional neural network $s_{\theta}({\mathbf{x}_{t}},t)$, called the score
model, parametrised by $\theta$. The resulting optimisation problem is given
by

$\min_{\theta}\left\\{L_{\text{DSM}}(\theta)=\mathbb{E}_{t\sim
U[0,T]}\mathbb{E}_{{\mathbf{x}_{0}}\sim{\pi}}\mathbb{E}_{{\mathbf{x}_{t}}\sim
p_{t}({\mathbf{x}_{t}}|{\mathbf{x}_{0}})}\left[\omega_{t}\|s_{\theta}({\mathbf{x}_{t}},t)-\nabla_{{\mathbf{x}}}\log
p_{t}({\mathbf{x}_{t}}|{\mathbf{x}_{0}})\|_{2}^{2}\right]\right\\},$ (7)

where $\omega_{t}>0$ are weighting factors, balancing the scores at different
time steps. For general SDEs, the loss $L_{\text{DSM}}(\theta)$ may still be
intractable, since it requires access to the transition density
$p_{t}({\mathbf{x}_{t}}|{\mathbf{x}_{0}})$. However, for SDEs with an affine
linear drift function, $p_{t}({\mathbf{x}_{t}}|{\mathbf{x}_{0}})$ is a
Gaussian and thus can be given in closed form (Särkkä and Solin, 2019).
Throughout the paper, we use $T=1$ and the variance preserving SDE (Ho et al.,
2020) given by

$d{\mathbf{x}_{t}}=-\frac{\beta(t)}{2}{\mathbf{x}_{t}}dt+\sqrt{\beta(t)}d\mathbf{w}_{t},$
(8)

where $\beta(t):[0,1]\to{\mathbb{R}}_{>0}$ is an increasing function defining
the noise schedule. We use
$\beta(t)=\beta_{\text{min}}+t(\beta_{\text{max}}-\beta_{\text{min}})$ giving
the transition kernel
$p_{t}({\mathbf{x}_{t}}|{\mathbf{x}_{0}})=\mathcal{N}({\mathbf{x}_{t}};\gamma_{t}{\mathbf{x}_{0}},\nu_{t}^{2}I)$
with coefficients $\gamma_{t},\nu_{t}\in{\mathbb{R}}$ computed from drift and
diffusion coefficients, see Appendix A.1 for details. Generating samples with
the score model $s_{\theta}({\mathbf{x}_{t}},t)$ as a surrogate requires
solving the reverse SDE (6), with the score model
$s_{\theta}({\mathbf{x}_{t}},t)$ in place of
$\nabla_{\mathbf{x}}\log{p_{t}}({\mathbf{x}_{t}})$:

$d{\mathbf{x}_{t}}=[{\mathbf{f}}({\mathbf{x}_{t}},t)-g(t)^{2}s_{\theta}({\mathbf{x}_{t}},t)]dt+g(t)d\bar{\mathbf{w}}_{t}.$
(9)

Drawing samples from the resulting generative model thus involves two steps.
First, drawing a sample from the terminal distribution
${\mathbf{x}}_{1}\sim\mathcal{N}(0,I)\approx p_{1}$, and second, initialising
the reverse SDE (9) with ${\mathbf{x}}_{1}$ and simulating backwards in time
until $t=0$. The latter can be achieved by Euler-Maruyama schemes or
predictor-corrector methods (Song et al., 2021c).

####  2.2.1 Denoising diffusion implicit models

Simulating the reverse SDE can be computationally expensive as a fine time
grid is often necessary to produce realistic samples. Denoising Diffusion
Implicit Models (DDIMs) (Song et al., 2021a) were introduced to allow faster
sampling, and build upon a result by Tweedie (Efron, 2011) to approximate the
expectation ${\mathbb{E}}[{\mathbf{x}_{0}}|{\mathbf{x}_{t}}]$ via the score
model $s_{\theta}({\mathbf{x}_{t}},t)$ as

${\mathbb{E}}[{\mathbf{x}_{0}}|{\mathbf{x}_{t}}]=\frac{{\mathbf{x}_{t}}+\nu_{t}^{2}\nabla_{\mathbf{x}}\log{p_{t}}({\mathbf{x}_{t}})}{\gamma_{t}}\approx\frac{{\mathbf{x}_{t}}+\nu_{t}^{2}s_{\theta}({\mathbf{x}_{t}},t)}{\gamma_{t}}:={\hat{\mathbf{x}}_{0}}({\mathbf{x}_{t}}).$
(10)

DDIM defines a non-Markovian sampling rule, which uses both the current sample
${\mathbf{x}_{t}}$ and Tweedie’s estimate
${\hat{\mathbf{x}}_{0}}({\mathbf{x}_{t}})$ to create an accelerated sampler.
Let $0=t_{k_{1}}\leq\dots\leq t_{k_{N}}=1$ be the time discretisation. The
DDIM sampling update can be written as

$\begin{split}{\mathbf{x}}_{t_{k-1}}&=\gamma_{t_{k-1}}{\hat{\mathbf{x}}_{0}}({\mathbf{x}}_{t_{k}})+\text{Noise}({\mathbf{x}}_{t_{k}},s_{\theta})+\eta_{t_{k}}{\mathbf{z}},\quad{\mathbf{z}}\sim\mathcal{N}(0,I)\\\
&\text{ with
}\text{Noise}({\mathbf{x}}_{t_{k}},s_{\theta}):=-\nu_{t_{k}}\sqrt{\nu_{t_{k-1}}^{2}-\eta_{t_{k}}^{2}}s_{\theta}({\mathbf{x}}_{t_{k}},t_{k}).\end{split}$
(11)

The sampling rule can be split into a denoising step (predicting
${\hat{\mathbf{x}}_{0}}({\mathbf{x}}_{t_{k}})$ using the score model), and
adding an appropriate amount of noise back. Thus, the sampling mimics an
iterative refinement process, as the prediction of the denoised estimate
${\hat{\mathbf{x}}_{0}}({\mathbf{x}}_{t_{k}})$ will be more accurate for
smaller $t_{k}$. Different choices of $\eta_{t}$ result in different sampling
schemes. We choose $\eta_{t_{k}}=\eta\beta_{t_{k}}$ with a hyperparameter
$\eta\in[0,1]$, controlling the amount of stochasticity in the sampling, and
$\beta_{t_{k}}=\nu_{t_{k-1}}/\nu_{t_{k}}\sqrt{1-\gamma_{t_{k}}/\gamma_{t_{k-1}}}$
(Song et al., 2021a).

###  2.3 Using Score-based Generative Models for Inverse Problems

The goal of the Bayesian framework of inverse problems is to estimate the
posterior ${p^{\text{post}}}({\mathbf{x}}|{\mathbf{y}})$, i.e. the conditional
distribution of images ${\mathbf{x}}$ given noisy measurements ${\mathbf{y}}$.
Using Bayes’ theorem the posterior can be factored into

${p^{\text{post}}}({\mathbf{x}}|{\mathbf{y}})\propto{p^{\text{lkhd}}}({\mathbf{y}}|{\mathbf{x}}){\pi}({\mathbf{x}}),\quad\nabla_{{\mathbf{x}}}\log{p^{\text{post}}}({\mathbf{x}}|{\mathbf{y}})=\nabla_{{\mathbf{x}}}\log{p^{\text{lkhd}}}({\mathbf{y}}|{\mathbf{x}})+\nabla_{{\mathbf{x}}}\log{\pi}({\mathbf{x}}),$
(12)

where ${p^{\text{lkhd}}}$ denotes the likelihood and ${\pi}$ the prior given
by the image distribution. We can set up a generative model for the posterior
in the same way as for the prior ${\pi}$ in Section 2.2 by defining a forward
SDE which maps the posterior to random noise. To generate a sample from the
posterior ${p^{\text{post}}}({\mathbf{x}}|{\mathbf{y}})$, we can simulate the
corresponding reverse SDE

$d{\mathbf{x}_{t}}=[{\mathbf{f}}({\mathbf{x}_{t}},t)-g(t)^{2}\nabla_{{\mathbf{x}}}\log
p_{t}({\mathbf{x}_{t}}|{\mathbf{y}})]dt+g(t)d\bar{\mathbf{w}}_{t},$ (13)

where we need access to the time-dependent posterior
$\nabla_{{\mathbf{x}}}\log p_{t}({\mathbf{x}_{t}}|{\mathbf{y}})$. Similar to
Eq. (12), we use Bayes’ theorem for the score of the posterior and decompose
$\nabla_{{\mathbf{x}}}\log p_{t}({\mathbf{x}_{t}}|{\mathbf{y}})$ into a prior
and a likelihood term, where the former is approximated with the trained score
model

$\begin{split}\nabla_{{\mathbf{x}}}\log
p_{t}({\mathbf{x}_{t}}|{\mathbf{y}})&=\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{x}_{t}})+\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}_{t}})\\\
&\approx s_{\theta}({\mathbf{x}_{t}},t)+\nabla_{{\mathbf{x}}}\log
p_{t}({\mathbf{y}}|{\mathbf{x}_{t}}).\end{split}$ (14)

Substituting the above approximation into (13), we obtain

$d{\mathbf{x}_{t}}=[{\mathbf{f}}({\mathbf{x}_{t}},t)-g(t)^{2}(s_{\theta}({\mathbf{x}_{t}},t)+\nabla_{{\mathbf{x}_{t}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}_{t}}))]dt+g(t)d\bar{\mathbf{w}}_{t}.$
(15)

We can recover approximate samples from the posterior
${p^{\text{post}}}({\mathbf{x}}|{\mathbf{y}})$, by simulating the reverse SDE
(15). Through iterative simulation of the reverse SDE with varying noise
initialisations, we can estimate moments of the posterior distribution. As is
common practice in the field (Song et al., 2022; Chung and Ye, 2022; Jalal et
al., 2021) we use one sample for the reconstruction, due to computational
costs of repeatedly solving the reverse SDE. In addition to the score model
$s_{\theta}({\mathbf{x}_{t}},t)$, we need the score of the time-dependent
likelihood $\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}_{t}})$.
At the start of the forward SDE (for $t=0$), it is equal to the true
likelihood ${p^{\text{lkhd}}}$. However, for $t>0$ the score
$\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}_{t}})$ is
intractable to compute exactly and different approximations have been
proposed. In (Jalal et al., 2021; Ramzi et al., 2020), this term was
approximated with the likelihood ${p^{\text{lkhd}}}$ evaluated at the noisy
sample ${\mathbf{x}_{t}}$ with time-dependent penalty strength $\lambda_{t}$

$\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}_{t}})\approx\lambda_{t}\nabla_{{\mathbf{x}}}\log{p^{\text{lkhd}}}({\mathbf{y}}|{\mathbf{x}_{t}}),$
(16)

We refer to Eq. (16) as the Naive approximation. The Diffusion Posterior
Sampling (DPS) (Chung et al., 2023a) uses Tweedie’s formula to obtain
$\hat{{\mathbf{x}}}_{0}({\mathbf{x}_{t}})\approx{\mathbb{E}}[{\mathbf{x}_{0}}|{\mathbf{x}_{t}}]$
and approximates
$\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}_{t}})$ by

$\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}_{t}})\approx\nabla_{{\mathbf{x}}}\log{p^{\text{lkhd}}}({\mathbf{y}}|{\hat{\mathbf{x}}_{0}}({\mathbf{x}_{t}})),$
(17)

where $\nabla_{{\mathbf{x}}}$ denotes taking derivative in ${\mathbf{x}}_{t}$
(instead of $\hat{{\mathbf{x}}}_{0}$). It was shown that this approximation
leads to improved performance for several image reconstruction tasks (Chung et
al., 2023a). However, DPS comes with a higher computational cost, due to the
need to back-propagate the gradient through the score model.

Recently, several works proposed modifying the DDIM sampling rule in Eq. (11)
for conditional generation (Zhu et al., 2023; Chung et al., 2023b). These
methods generally consist of three steps: (1) estimating the denoised image
${\mathbf{x}_{0}}$ using Tweedie’s estimate
${\hat{\mathbf{x}}_{0}}({\mathbf{x}}_{t_{k}})$; (2) updating
${\hat{\mathbf{x}}_{0}}({\mathbf{x}}_{t_{k}})$ with a data consistency step
using the measurements ${\mathbf{y}}$; and (3) adding the noise back,
according to the DDIM update rule, in order to get a sample for the next time
step $t_{k-1}$. Importantly, with this approach there is no need to estimate
the gradient of the time-dependent likelihood
$\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}_{t}})$ as data
consistency is only enforced on Tweedie’s estimate at $t=0$. These conditional
DDIM samplers differ most greatly in the implementation of the data
consistency update. Decomposed Diffusion Sampling (DDS) (Chung et al., 2023b)
proposes to align Tweedie’s estimate with the measurements by running $p$
steps of a Conjugate Gradient (CG) scheme for minimising the negative log-
likelihood at each sampling step. Let
$\text{CG}^{(p)}({\hat{\mathbf{x}}_{0}})$ denote the $p$-th CG update
initialised with ${\hat{\mathbf{x}}_{0}}({\mathbf{x}}_{t_{k}})$. This can be
seen as an approximation to the conditional expectation, i.e.
${\mathbb{E}}[{\mathbf{x}_{0}}|{\mathbf{x}_{t}},{\mathbf{y}}]\approx\text{CG}^{(p)}({\hat{\mathbf{x}}_{0}})$
(Ravula et al., 2023). Using this approximation, the update step for DDS can
be written as

${\mathbf{x}}_{t_{k-1}}=\gamma_{t_{k-1}}\text{CG}^{(p)}({\hat{\mathbf{x}}_{0}})+\text{Noise}({\mathbf{x}}_{t_{k}},s_{\theta})+\eta_{t_{k}}{\mathbf{z}},\text{
with }{\mathbf{z}}\sim\mathcal{N}(0,I),$ (18)

where the introduction of the conditional expectation offers us the
possibility to explore different approximations specific for PET image
reconstruction.

##  3 PET-specific Adaptations for SGMs

To apply SGMs to PET reconstruction, several key components of the pipeline in
Section 2.2 have to be modified in order to incorporate PET-specific
constraints. Namely, we introduce measurement-based normalisation of the input
to the score model, and explain how to apply a score model trained on 2D
slices for 3D reconstruction. Additionally, we adapt the sampling methods from
Section 2.3 to incorporate the Poisson noise model. Finally, we demonstrate
that the SGM framework allows for the incorporation of additional information,
e.g. MR images, by using classifier-free guidance (Ho and Salimans, 2022).

###  3.1 Measurement-based Normalisation

The intensity of the unknown tracer distribution in emission tomography can
significantly vary across different scans, resulting in a high dynamic range
that poses challenges for deep learning approaches. Neural networks may
exhibit bias toward intensity levels that appear more frequently in the
training set. Consequently, the network might struggle to handle new images
with unseen intensity levels, leading to instability in the learning and
evaluation process (Tran et al., 2020). For SGMs the intensity range of images
must be predefined to ensure that the forward diffusion process converges to a
standard Gaussian distribution and to stabilise the sampling process (Lou and
Ermon, 2023). Input normalisation is a standard deep learning methodology to
deal with intensity shifts and normalise the inputs to the network. In a
similar vein, we propose a PET-specific normalisation method to ensure that
the score model $s_{\theta}({\mathbf{x}_{t}},t)$ is able to estimate the score
function of images with arbitrary intensity values. Namely, we normalise each
training image ${\mathbf{x}_{0}}$ to ensure voxel intensities are within a
certain range. To do this we introduce a training normalisation factor
$c_{\text{train}}$ that when applied ensures the average emission per emission
voxels (a voxel with non-zero intensity value) is 1. This is computed as

$c_{\text{train}}=c({\mathbf{x}_{0}}):=\frac{\sum_{j=1}^{m}[{\mathbf{x}_{0}}]_{j}}{\\#\\{j:[{\mathbf{x}_{0}}]_{j}>0\\}},$
(19)

where the numerator captures the total emission in the image, and the
denominator is the number of emission voxels. The normalisation factor is
incorporated into the DSM training objective function by rescaling the initial
image, yielding the objective

$\mathbb{E}_{t\sim
U[0,T]}\mathbb{E}_{{\mathbf{x}_{0}}\sim{\pi}}\mathbb{E}_{{\mathbf{z}}\sim\mathcal{N}(0,I)}\mathbb{E}_{c\sim
U[\frac{c_{\rm train}}{2},\frac{3c_{\rm
train}}{2}]}\left[\omega_{t}\left\|s_{\theta}\left(\tilde{\mathbf{x}}_{t},t\right)-\nabla_{{\mathbf{x}}}\log
p_{t}(\tilde{\mathbf{x}}_{t}|{\mathbf{x}_{0}}/c)\right\|_{2}^{2}\right],$ (20)

with
$\tilde{\mathbf{x}}_{t}=\gamma_{t}{\mathbf{x}_{0}}/c+\nu_{t}{\mathbf{z}}$.
Compared with Eq. (7), the scale-factor in range $c\sim U[\frac{c_{\rm
train}}{2},\frac{3c_{\rm train}}{2}]$ is used to encourage the score model to
be more robust with respect to misestimations of the normalisation constant
during sampling.

An analogue of Eq. (19) is unavailable during the sampling, and thus a
surrogate is required. This is obtained through an approximate reconstruction,
computed using a single epoch of OSEM from a constant non-negative
initialisation. The resulting sampling normalisation factor is given by

$c_{\text{OSEM}}=\frac{\sum_{j=1}^{m}[{\mathbf{x}}_{\text{OSEM}}]_{j}}{\\#\\{j:[{\mathbf{x}}_{\rm
OSEM}]_{j}>Q_{0.01}\\}},$ (21)

where $Q_{0.01}$ defines the $1\%$ percentile of ${\mathbf{x}}_{\text{OSEM}}$
values. This threshold is heuristically chosen to ensure that noise and
reconstruction artefacts do not cause an over-estimation of the number of
emission voxels. In the reconstruction process the normalisation constant
$c_{\rm OSEM}$ is applied as a factor scaling the time-dependent likelihood,
giving

$d{\mathbf{x}_{t}}=[{\mathbf{f}}({\mathbf{x}_{t}},t)-g(t)^{2}(s_{\theta}({\mathbf{x}_{t}},t)+\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|c_{\text{OSEM}}{\mathbf{x}_{t}}))]dt+g(t)d\bar{\mathbf{w}}_{t}.$
(22)

At final time step $t=0$, the output ${\mathbf{x}}$ is rescaled by $c_{\rm
OSEM}$ to recover the correct intensity level.

###  3.2 Scaling to 3D Reconstruction

While some SGM studies deal with fully 3D image generation (Pinaya et al.,
2022), the majority of work is restricted to 2D images. This is largely due to
the fact that the learning of full 3D volume distributions is computationally
expensive and requires access to many training volumes. Therefore, we propose
to train the score model on 2D axial slices and use a specific decomposition
of the conditional sampling rules to apply the model for 3D reconstruction.
Upon simulating the conditional reverse SDE in Eq. (15) using the Euler-
Maruyama approach, we arrive at the iteration rule:

$\displaystyle\tilde{\mathbf{x}}_{t_{k-1}}={\mathbf{x}}_{t_{k}}+\big{[}{\mathbf{f}}({\mathbf{x}}_{t_{k}},t_{k})-g(t_{k})^{2}s_{\theta}({\mathbf{x}}_{t_{k}},t_{k})\big{]}\Delta
t+g(t_{k})\sqrt{|\Delta
t|}{\mathbf{z}},\quad{\mathbf{z}}\sim\mathcal{N}(0,I),$ (23)
$\displaystyle{\mathbf{x}}_{t_{k-1}}=\tilde{\mathbf{x}}_{t_{k-1}}-g(t_{k})^{2}\nabla_{{\mathbf{x}}}\log
p_{t_{k}}({\mathbf{y}}|{\mathbf{x}_{t}}_{k})\Delta t,$ (24)

using an equidistant time discretisation $0=t_{k_{1}}\leq\dots\leq
t_{k_{N}}=1$ for $N\in{\mathbb{N}}$, with a time step $\Delta t=-1/N$. We
split the Euler-Maruyama update into two equations to highlight the influences
of the score model and the measurements ${\mathbf{y}}$. First, Eq. (23) is the
Euler-Maruyama discretisation for the unconditional reverse SDE, see Eq. (9).
This update is independent of the measurements ${\mathbf{y}}$ and can be
interpreted as a prior update, increasing the likelihood of
${\mathbf{x}}_{t_{k}}$ under the SGM. The second step in Eq. (24) is a data
consistency update, pushing the current iterate to better fit the
measurements. Notably, this step is fully independent of the score model. This
strategy was developed for 3D reconstruction, focusing on sparse view CT and
MRI (Chung et al., 2023c). It was proposed to apply the prior update in Eq.
(23) to all slices in the volume independently and use the 3D forward model in
the data consistency step. Further, a regulariser in the direction orthogonal
to the slice was introduced, to improve consistency of neighbouring slices.
However, applying this approach to the Euler-Maruyama discretisation results
in slow sampling times as a small time step $|\Delta t|$ is necessary. To
accelerate the sampling of high quality samples, we propose to use the DDS
update in Eq. (18) that uses a similar decomposition of independent score
model updates to axial slices, and 3D data consistency updates. Additionally,
we accelerate data consistency updates by splitting the measurement data into
ordered subsets and applying the forward model block-sequentially. The details
are explained below.

###  3.3 Modifications of Sampling Methods

The sampling schemes and approximations in Section 2.3 were originally
proposed for inverse problems with Gaussian noise. The work on DPS (Chung et
al., 2023a) also considers inverse problems with Poisson noise, but utilises a
Gaussian approximation to the Poisson noise, which is not appropriate for PET
reconstruction due to the low photon count. To apply the Naive approximation
or the DPS approach to Poisson inverse problems, one could simply replace the
Gaussian log-likelihood with the PLL in Eq. (16) and Eq. (17). However, PLL
and its gradient are only defined for non-negative values. Therefore, we have
to introduce a non-negativity projection into the sampling to ensure that the
gradient of the PLL can be evaluated. In the context of guided diffusion, it
was proposed to project the iterates ${\mathbf{x}}_{t_{k}}$ to a specified
domain after each sampling step (Li et al., 2022; Saharia et al., 2022). In
our case this would require thresholding all negative values. However, this
creates a mismatch between the forward and reverse SDEs. It was observed that
this mismatch results in artefacts in the reconstructions and may even lead to
divergence of the sampling (Lou and Ermon, 2023). Therefore, we propose to
only threshold the input to the PLL, i.e. with $L$ being the PLL, see Eq. (2),
for the PET-Naive approximation we use

$\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}_{t}})\approx\lambda_{t}^{\text{{Naive}{}}}\nabla_{{\mathbf{x}}}L({\mathbf{y}}|c_{\text{OSEM}}P_{{\mathbf{x}}\geq
0}[{\mathbf{x}}_{t}]),$ (25)

and likewise for PET-DPS

$\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}_{t}})\approx\lambda_{t}^{\text{DPS}}\nabla_{{\mathbf{x}}}L({\mathbf{y}}|c_{\text{OSEM}}P_{{\mathbf{x}}\geq
0}[{\hat{\mathbf{x}}_{0}}({\mathbf{x}}_{t})]).$ (26)

Note that this leads to a perturbed likelihood gradient that is not computed
on the true iterate ${\mathbf{x}}_{t}$, but only on the projection. In order
to reconstruct the PET image we have to solve the reverse SDE using the
specific approximation (PET-Naive or PET-DPS) as the likelihood term. This
usually requires around $1000$ sampling steps to produce an appropriate
reconstruction and results in impractically long reconstruction times for 3D
volumes.

To reduce the reconstruction times we propose to modify the conditional DDIM
sampling rule, which we call PET-DDS, similar to the DDS framework (Zhu et
al., 2023; Chung et al., 2023b), cf. Section 3.3. This circumvents the usage
of $\nabla_{\mathbf{x}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}_{t}})$, instead
enforcing data consistency for Tweedie’s estimate
${\hat{\mathbf{x}}_{0}}({\mathbf{x}}_{t_{k}})$. For PET reconstruction we
propose to implement this data consistency with a MAP objective, leading to
the PET-DDS update

$\displaystyle{\mathbf{x}}^{0}_{t_{k}}$
$\displaystyle={\hat{\mathbf{x}}_{0}}({\mathbf{x}}_{t_{k}})$ (27)
$\displaystyle{\mathbf{x}}^{i+1}_{t_{k}}$ $\displaystyle=P_{\mathbf{x}\geq
0}\left[{\mathbf{x}}^{i}_{t_{k}}+\mathbf{D}({\mathbf{x}}^{i}_{t_{k}})\nabla_{\mathbf{x}}\Phi_{j}({\mathbf{x}}^{i}_{t_{k}})\right]$
(28) $\displaystyle\quad\quad i=0,\dots,p-1$
$\displaystyle{\mathbf{x}}_{t_{k-1}}$
$\displaystyle=\gamma_{t_{k-1}}{\mathbf{x}}^{p}_{t_{k}}+\text{Noise}({\mathbf{x}}_{t_{k}},s_{\theta})+\eta_{t_{k}}{\mathbf{z}},\quad{\mathbf{z}}\sim\mathcal{N}(0,I).$
(29)

where the sub-objective is

$\Phi_{j}({\mathbf{x}}^{i})=L_{j}({\mathbf{y}}|c_{\text{OSEM}}{\mathbf{x}}^{i})+(\lambda^{\text{RDP}}R_{z}({\mathbf{x}}^{i})-\lambda^{\text{DDS}}\|{\mathbf{x}}^{i}-{\hat{\mathbf{x}}_{0}}\|_{2}^{2})/n_{\mathrm{sub}}.$
(30)

The sub-objective index $j=j(i)$ is given by $j=(p(N-k)+i\mod
n_{\rm{sub}})+1$, which cyclically accesses sub-objectives. The RDP used for
3D data $R_{z}$ is applied in the $z$-direction, perpendicular to the axial
slice, see Appendix A.3 for more details. The prior in Eq. (30) consists of
two components: one anchoring to the Tweedie’s estimate
$\|{\mathbf{x}}-{\hat{\mathbf{x}}_{0}}\|_{2}^{2}$, and the other RDP in the
$z$-direction $R_{z}({\mathbf{x}})$. The components have associated penalty
strengths $\lambda^{\text{RDP}}$ and $\lambda^{\text{DDS}}$, respectively.

In a PET-DDS update we first independently compute Tweedie’s estimate based on
${\mathbf{x}}_{t_{k}}$ for each axial slice (Eq. 27). Tweedie’s estimate
${\hat{\mathbf{x}}_{0}}$ impacts the reconstruction in two ways: first through
the Tikhonov regulariser scaled with $\lambda^{\text{DDS}}$, and second as the
initial value for the projected gradient descent in Eq. (28). Through running
$p$ steps of projected gradient descent consistency is balanced between a PLL
on measurements, RDP in the $z$-direction, and Tweedie’s estimate (Eq. 30). To
speed up computation of the objective gradient, the objective is split into
sub-objectives and the gradient of the log-likelihood is evaluated using only
subsets of the measurements ${\mathbf{y}}$, similar to the BSREM update in Eq.
(4). The subsets are partitioned in a staggered configuration and are ordered
with a Herman-Meyer order (Herman and Meyer, 1993). Eq. (29) is the DDIM
update applied to the conditioned Tweedie estimate ${\mathbf{x}}^{p}$, where
the score update is again applied independently for each axial slice, here the
notation of $\text{Noise}({\mathbf{x}}_{t_{k}},s_{\theta})$ is overloaded. The
DDIM update gives ${\mathbf{x}}_{t_{k-1}}$, and these PET-DDS updates repeat
until $t_{0}=0$ giving reconstruction $\hat{{\mathbf{x}}}$.

###  3.4 MR Image Guided Reconstruction

In recent years, several regularisation methods have been proposed which
leverage the availability of additional MR images to improve PET image
reconstruction (Ehrhardt et al., 2016; Bai et al., 2013; Somayajula et al.,
2011). These studies often encode anatomical features of the MR image as edges
or level sets and build hand-crafted regularisers based on these encoded
features. This is commonly referred to as guided reconstruction (Ehrhardt,
2021), where the MR image is first reconstructed and is then used in the PET
reconstruction pipeline. The SGM approach can be modified for guided
reconstruction. In this setting we can use Bayes’ theorem to express the
posterior

$\nabla_{\mathbf{x}}\log{p^{\text{post}}}({\mathbf{x}}|{\mathbf{y}},{\mathbf{x}_{\text{MR}}})=\nabla_{\mathbf{x}}\log{p^{\text{lkhd}}}({\mathbf{y}}|{\mathbf{x}})+\nabla_{\mathbf{x}}\log{\pi}({\mathbf{x}}|{\mathbf{x}_{\text{MR}}}),$
(31)

assuming that ${\mathbf{y}}$ and ${\mathbf{x}_{\text{MR}}}$ are conditionally
independent given ${\mathbf{x}}$. Here, the likelihood
${p^{\text{lkhd}}}({\mathbf{y}}|{\mathbf{x}})$ is given by the Poisson noise
model and ${\pi}({\mathbf{x}}|{\mathbf{x}_{\text{MR}}})$ is a prior
conditioned on the MR image ${\mathbf{x}_{\text{MR}}}$, which will be learned
via a score model. Using this decomposition, the reverse SDE, given the MR
image, can be written as

$d{\mathbf{x}_{t}}=[{\mathbf{f}}({\mathbf{x}_{t}},t)-g(t)^{2}\left(\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}})+\nabla_{{\mathbf{x}}}\log{\pi}({\mathbf{x}}|{\mathbf{x}_{\text{MR}}})\right)]dt+g(t)d\bar{\mathbf{w}}_{t}.$
(32)

We can use PET-Naive or PET-DPS to approximate the score of the time dependent
likelihood $\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{y}}|{\mathbf{x}})$.
However, we have to train a score model, conditioned on the MR image, to
estimate the conditional score function

$s_{\theta}({\mathbf{x}_{t}};t,{\mathbf{x}_{\text{MR}}})\approx\nabla_{\mathbf{x}}\log{p_{t}}({\mathbf{x}}|{\mathbf{x}_{\text{MR}}}),$
(33)

where ${\mathbf{x}_{\text{MR}}}$ is an additional input to the score model. In
order to train such a score model, we need a dataset of PET images with
corresponding MR images. Learning such a conditional score model,
$s_{\theta}({\mathbf{x}};t,{\mathbf{x}_{\text{MR}}})\approx\nabla_{\mathbf{x}}\log{p_{t}}({\mathbf{x}}|{\mathbf{x}_{\text{MR}}})$,
was recently proposed and applied for PET image denoising (Gong et al., 2022).
However, learning a conditional score model requires a paired dataset
$\\{({\mathbf{x}}^{i},\mathbf{x}_{\text{MR}}^{i})\\}_{i=1}^{m}$ of PET images
and corresponding MR images. In contrast, using the Classifier Free Guidance
(CFG) framework (Ho and Salimans, 2022), we only need a partly paired dataset,
i.e. besides paired data
$\\{({\mathbf{x}}^{i},\mathbf{x}_{\text{MR}}^{i})\\}_{i=1}^{m_{1}}$ we can
also make use of unpaired data $\\{{\mathbf{x}}^{i}\\}_{i=1}^{m_{2}}$. In
particular, CFG trains both a conditional and unconditional score model
simultaneously and utilises their combination during the sampling process. CFG
uses a zero image $\mathbf{0}$ to distinguish between the conditional and
unconditional score model

$s_{\theta}({\mathbf{x}_{t}};t,{\mathbf{x}_{\text{MR}}})\approx\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{x}_{t}}|{\mathbf{x}_{\text{MR}}})\quad\mbox{and}\quad
s_{\theta}({\mathbf{x}_{t}};t,{\mathbf{x}_{\text{MR}}}=\mathbf{0})\approx\nabla_{{\mathbf{x}}}\log{p_{t}}({\mathbf{x}_{t}}),$
(34)

yielding a conditional DSM objective

$\begin{split}\\!{\mathbb{E}}_{t\sim
U[0,T]}{\mathbb{E}}_{{\mathbf{x}_{0}},{\mathbf{x}_{\text{MR}}}\sim{\pi}}{\mathbb{E}}_{{\mathbf{x}_{t}}\sim
p_{t}({\mathbf{x}_{t}}|{\mathbf{x}_{0}})}{\mathbb{E}}_{\rho\sim
B(q)}\\!\left[\omega_{t}\|s_{\theta}({\mathbf{x}_{t}},t;\rho~{}{\mathbf{x}_{\text{MR}}})\\!-\\!\nabla_{{\mathbf{x}}}\log
p_{t}({\mathbf{x}_{t}}|{\mathbf{x}_{0}})\|_{2}^{2}\right]\\!\\},\end{split}$
(35)

where $B(q)$ is a Bernoulli distribution with parameter $q$. Thus, if the
additional MR input is set to zero, the conditional DSM loss matches the
unconditional DSM loss defined in Eq. (7). After training, CFG defines a
combined score model

$\tilde{s}_{\theta}({\mathbf{x}_{t}};t,{\mathbf{x}_{\text{MR}}})=(1+w)s_{\theta}({\mathbf{x}_{t}};t,{\mathbf{x}_{\text{MR}}})-ws_{\theta}({\mathbf{x}_{t}};t,\mathbf{0}),$
(36)

as a linear combination with $w$ as the guidance strength. This combined score
model $s_{\theta}({\mathbf{x}_{t}};t,{\mathbf{x}_{\text{MR}}})$ can then be
used for any of the presented sampling methods.

##  4 Experimental Setup

###  4.1 Dataset and Evaluation Metrics

We use the BrainWeb dataset consisting of $20$ patient-realistic volumes
(Aubert-Broche et al., 2006). The tracer simulated was 18F-Fluorodeoxyglucose
(FDG) and the volumes were further perturbed by three realisations of random
distortions (Schramm, 2021). $19$ out of the $20$ volumes were used for
training. Axial slices with non-zero intensity were extracted, resulting in a
training dataset of $4569$ slices. For 2D evaluation, we used $20$ equidistant
axial slices from remaining volume (subject 04). An additional OOD dataset was
created by simulating ellipsoid hot lesions of random size and location within
soft-tissue. The noise level of simulated measurements was set by re-scaling
forward projected ground truth images, where the scale ensured that the total
counts divided by emission volume was 2.5 or 10. These rescaled measurements
are the clean measurements, which were then corrupted with Poisson noise and
constant background contamination. In addition, $10$ noise realisations were
obtained. Herein, we refer to the noise levels as $2.5$ and $10$, where the
total true counts averaged over evaluation dataset were $122\,808$ and
$491\,232$, respectively. Resolution modelling, attenuation, sensitivity, and
background contamination were modelled and subsequently included in the
forward model utilising ParallelProj (Schramm, 2022).

For the 3D evaluation, measurements of subject 04 were simulated with an
Siemens Biograph mMR scanner geometry (Karlberg et al., 2016). Measurements
with detector sensitivities and attenuation were simulated and included in the
forward model using SIRF and STIR (Ovtchinnikov et al., 2020; Thielemans et
al., 2012). The noise level was equivalent to 40 million counts without
background, and 5 noisy realisations were obtained. Both FDG and Amyloid
tracers were simulated with hot lesions. The projector and measurements were
split into 28 ordered subsets.

We evaluate the performance between reconstructions and ground truths using
two global metrics: Peak-Signal-to-Noise Ratio (PSNR) and Structural
Similarity Index Measure (SSIM) (Wang et al., 2004). Moreover, we compute two
local quality scores over a Region of Interest (ROI). First, to quantify the
detectability of lesions, we compute the Contrast Recovery Coefficient (CRC).
Second, the noise in reconstructions is estimated over background ROIs using
Standard Deviation (STD), which computes standard deviation across
realisations and then averages over ROI voxels (Tong et al., 2010). In 2D, we
evaluate the reconstruction consistency by computing the Kullback-Leibler
Divergence (KLDIV) between measurements ${\mathbf{y}}$ and estimated
measurements
$\mathbf{\bar{y}}={\mathbf{A}}\hat{{\mathbf{x}}}+\mathbf{\bar{b}}$, where
$\hat{{\mathbf{x}}}$ denotes the reconstruction. Furthermore, we include the
“mean KL” between the noisy and clean measurements across the 2D evaluation
dataset. More information about quality metrics can be found in the Appendix
A.4.

We present tables of best performing methods with optimal penalty strengths,
as well as qualitative figures of reconstructed images. Furthermore, to allow
direct comparison between methods, we give sensitivity plots of PSNR, SSIM,
KLDIV or CRC vs. STDs. Since STD gives an estimate of the noise in the image,
these plots can show the effect of varying penalty strength on reconstruction
quality or data-consistency. A lower STD typically corresponds to lower data-
fidelity (a higher prior strength), and the converse is true for higher STD.
In practice, with a generative model as a prior, higher penalty strengths do
not necessarily lead to a lower STD as there may be multiple reconstruction
with high likelihood under the model. Variations in STD are further
exacerbated by approximate nature and stochasticity of SGM sampling.

###  4.2 Comparison Methods

In the 2D setting, we compare against two established supervised learning
methods used in medical image reconstruction: the UNet post-processing
FBPConvNet (Jin et al., 2017) and unrolled iterative learned primal dual
(Adler and Öktem, 2018). We modify both models for PET reconstruction; the
post-processing method is referred to as PET-UNet, and the unrolled method as
PET-LPD. Additionally we compare against a state-of-the-art SGM approach for
PET image denoising (Gong et al., 2022), referred to as Naive (OSEM). This
denoising approach replaces the likelihood on the measurements with a
likelihood modelled as a Gaussian centred at the noisy reconstruction.
Therefore, Naive (OSEM) is able to use the same pre-trained score model as our
proposed PET-Naive, PET-DPS, and PET-DDS methods.

For 3D evaluation, Deep Image Prior (DIP) reconstruction was included as an
unsupervised comparison method with a 3D network architecture well-established
in literature (Gong et al., 2019; Ote et al., 2023; Singh et al., 2023). For
comparison, converged MAP solutions with an RDP regulariser were computed.
BSREM algorithm was used with a range of penalty strengths, cf. PET background
Sect. 2.1.

Further details on all comparison methods can be found in Appendix A.3.

##  5 Numerical Experiments

The first set of experiments investigates the performance of the SGM methods
(Naive (OSEM), PET-Naive, PET-DPS, PET-DDS) against one-another and against
established supervised methods (PET-UNet, PET-LPD). This is done in 2D and at
two noise levels, with and without lesions. In the second set of experiments
we present results with MR image guidance. The last set of experiments
investigates the best performing SGM method (PET-DDS) on 3D reconstruction,
and provides a comparison against classical MAP and state-of-the-art DIP
reconstructions with lesions and two simulated tracers. For all SGM results we
make use of a single score model trained on the dataset of axial BrainWeb
slices as discussed in Section 4.1. The details about the training process and
network architecture can be found in Appendix A.1. Further results can be
found in the Appendix B.

###  5.1 2D Reconstruction

The aim of 2D experiments is to benchmark the SGM and supervised methods, and
analyse the stability of SGM methods with respect to the choice of different
penalty strengths $\lambda_{t}^{\text{{Naive}}},\lambda_{t}^{\text{DPS}}$ and
$\lambda^{\text{DDS}}$. The penalty strengths for PET-Naive and PET-DPS
depends on the time step $t$, and the details about their specific choice can
be found in Appendix A.2.

####  5.1.1 Reconstruction without Lesion

The results in Fig. 1 show that the performance of the four SGM methods vary
greatly for data of noise level 2.5 with no lesions. PET-DPS is the best
performing method, consistently giving high PSNR, SSIM and low KLDIV values.
However, it is also computationally the most expensive, requiring $1000$ steps
with back-propagation through the score model. PET-DDS preforms competitively
with a much lower computational overhead of $100$ steps without score model
back-propagation. Naive (OSEM) performs well with regards to PSNR, but
performs poorly in terms of data-consistency (KLDIV) and SSIM. As Naive (OSEM)
computes the likelihood on an early-stopped OSEM image, increasing data-
consistency ensures the reconstruction approaches the OSEM image. The maximum
achievable likelihood of Naive (OSEM) does not give a KLDIV lower than the
“mean KL”. Hence it is not deemed a strong surrogate to the true likelihood
computed on measurements. The PET-Naive reconstructions have substantially
higher STD values. This is attributed to instability when computing the PLL
gradient due to non-negativity projection directly applied on
${\mathbf{x}_{t}}$.

Figure 1: Results for BrainWeb without lesions with noise level 2.5 for
different penalty parameters. Standard deviation is across reconstructions
from different realisations of measurements.

In Table 1 we show quantitative results of the optimal penalty strength choice
for each metric, and comparisons against PET-UNet and PET-LPD. These
supervised methods are trained on data with noise levels of 5, 10 and 50
without lesions. Using noise levels 2.5 and 10 in evaluation allows
investigating the effect of OOD noise levels on supervised methods. PET-LPD is
the best performing method, giving the best SSIM at noise level 10, and best
PSNR at both noise levels. Between noise levels 10 and 2.5 PET-LPD observes a
drop of 6.7% and 6.6% for PSNR and SSIM, whereas PET-DPS exhibits a drop of
3.4% and 3.8%, respectively. PET-DPS performs competitively across both noise
levels and metrics, and gives the best SSIM value at noise level 2.5. Given
this competitive performance and lower reduction in drop of quality metrics
between noise level, PET-DPS is deemed more robust to different noise level.
This may be attributed to the unsupervised nature of SGM methods. Namely, as
they are not trained on data of given noise levels they are less affected by
distributional differences in noise levels at evaluation and training stages.

Table 1: Results using the best hyperparameters for each method for BrainWeb without lesions for noise level 2.5 and 10. The best SGM is highlighted in grey, and overall best metric is underlined. Supervised methods are trained with data of noise level 10, but not 2.5, and are in-distribution when evaluated with noise level 10.  | Noise Level | 2.5 | 10
---|---|---|---
|  | PSNR, $\lambda$ | SSIM, $\lambda$ | PSNR, $\lambda$ | SSIM, $\lambda$
Score-based Model | Naive (OSEM) |  $22.38$, $0.527$ |  $0.770$, $3.08$ |  $23.40$, $0.2$ |  $0.792$, $0.9$
PET-Naive |  $21.52$, $12.0$ |  $0.781$, $12.0$ |  $22.81$, $10.0$ |  $0.815$, $10.0$
PET-DPS |  $22.80$, $650.$ |  $0.818$, $750.$ |  $23.70$, $400.$ |  $0.850$, $400.$
PET-DDS |  $22.46$, $0.25$ |  $0.789$, $0.2$ |  $23.55$, $0.025$ |  $0.849$, $0.025$
Supervised Models | PET-LPD |  $23.07$, N/A |  $0.813$, N/A |  $24.72$, N/A |  $0.87$, N/A
PET-UNet |  $22.80$, N/A |  $0.80$, N/A |  $24.52$, N/A |  $0.868$, N/A

####  5.1.2 2D Reconstruction with Lesion

As the score model was trained on data without lesions, testing on data with
simulated hot lesions gives an insight into generalisability to OOD data. The
quantitative results in Fig. 2 and Table 2 show results that are consistent
with those for data with no lesions in Fig. 1. CRC was computed to quantify
the detectability of hot lesions. The CRC results indicate that PET-DDS is
better at resolving lesions than other SGM methods. Further, Fig. 2 shows a
clear trade-off between reconstruction quality in terms of PSNR and SSIM and
visibility of lesions. Here, a lower regularisation results in a better
performance in terms of CRC. Results for noise level $10$ are shown in
Appendix B.1.

Figure 2: Results for BrainWeb with lesions with noise level 2.5 for different
penalty parameters. Standard deviation is across reconstructions from
different realisations of measurements.

Comparing the results between noise levels $2.5$ and $10$ in Table 2, we
observe that SGMs increase CRC values as compared to supervised methods. SGMs
also compare favourably with regards to PSNR and SSIM. CRC is local metric
that is more relevant than PSNR or SSIM in a clinical setting, as it
quantifies the detectability of lesions. Therefore, it is of greater interest
to improve this local metric rather than global metrics. With this
perspective, SGMs outperform supervised methods, and the best-preforming SGM
methods are PET-DPS and PET-DDS . Due to the performance observed and
computational overhead, PET-DDS is considered the most appropriate method to
test in guided reconstruction and in the 3D setting.

Table 2: Results using the best hyperparameters for each method for BrainWeb with lesions for noise level 2.5 and 10. The penalty strength for the SGMs methods is denoted by $\lambda$. The best score-based method is highlighted in grey. The overall best score per noise level is underlined. Noise Level |  |  | PSNR, $\lambda$ | SSIM, $\lambda$ | CRC, $\lambda$
---|---|---|---|---|---
2.5 |  Score-based Model | Naive (OSEM) |  $27.60$, $0.527$ |  $0.821$, $1.71$ |  $0.865$, $50.$
PET-Naive |  $26.82$, $12.0$ |  $0.817$, $12.0$ |  $0.761$, $50.$
PET-DPS |  $27.99$, $625.$ |  $0.855$, $650.$ |  $0.822$, $1500.$
PET-DDS |  $27.46$, $0.15$ |  $0.841$, $0.15$ |  $0.910$, $0.01$
Supervised Models | PET-LPD |  $28.30$, N/A |  $0.853$, N/A |  $0.865$, N/A
PET-UNet |  $27.74$, N/A |  $0.836$, N/A |  $0.787$, N/A
10 |  Score-based Model | Naive (OSEM) |  $28.87$, $0.25$ |  $0.847$, $0.9$ |  $0.898$, $4.$
PET-Naive |  $28.07$, $10.0$ |  $0.845$, $7.5$ |  $0.829$, $20.$
PET-DPS |  $29.01$, $400.$ |  $0.878$, $400.$ |  $0.907$, $550.$
PET-DDS |  $28.99$, $0.025$ |  $0.879$, $0.025$ |  $0.962$, $0.$111Regularised due to denoised score estimate initialisation.
Supervised Models | PET-LPD |  $30.07$, N/A |  $0.894$, N/A |  $0.904$, N/A
PET-UNet |  $29.41$, N/A |  $0.889$, N/A |  $0.856$, N/A

####  5.1.3 MR Guided Reconstruction

Experiments with and without additional MR image guidance were conducted to
illustrate the flexibility of the proposed approach, and tested at three
guidance strengths $w=0.25$, $0.5$, $1.0$, where the guidance strength $w$
closer to zero constitutes more guidance. The results with best hyper-
parameters are given in Table 3. It is observed that there are significant
improvements to PSNR ($>18\%$) and SSIM ($>13\%$) with guidance. On PET data
with lesions, the lesions were only simulated for PET and not MR images.
Therefore, the data was of a worse-case scenario where clinically important
features are only present in the PET image. The results with lesions show
increasing the guidance strength decreased of CRC values and the lesions were
more difficult to detect - see Fig. 12. Conversely, the PSNR and SSIM values
on with lesions data increased with $w$ closer to zero (more guidance). This
highlights the potential dangers of guidance, as well as the importance of
evaluating local and global quality metrics.

Table 3: Results using the best hyperparameters for SGM methods for noise level 2.5 with MR image guidance. The best method for each setting (with/out lesion, and performance metric) is highlighted in gray, where the penalty strength is tuned for each method individually. | without lesions | with lesions
---|---|---
| PSNR, $\lambda$ | SSIM, $\lambda$ | PSNR, $\lambda$ | SSIM, $\lambda$ | CRC, $\lambda$
DDS (w/o MR) |  $22.46$, $0.25$ |  $0.789$, $0.2$ |  $27.46$, $0.15$ |  $0.841$, $0.15$ |  $0.910$, $0.01$
DDS $w=0.25$ |  $30.22$, $0.35$ |  $0.950$, $0.35$ |  $31.21$, $0.15$ |  $0.954$, $0.25$ |  $0.726$, $0.0$
DDS $w=0.5$ |  $29.32$, $0.25$ |  $0.940$, $0.25$ |  $31.12$, $0.15$ |  $0.946$, $0.25$ |  $0.778$, $0.0$
DDS $w=1.0$ |  $26.66$, $0.15$ |  $0.899$, $0.15$ |  $29.31$, $0.1$ |  $0.906$, $0.15$ |  $0.939$, $0.0$

From Fig. 3 a reconstruction without guidance and with guidance of various
strengths is presented for data without lesions at noise level 2.5. The
reconstructions indicate that MR guidance helps to reconstruct the specific
anatomical boundaries and structure, i.e. white matter tracts. In Appendix B.2
we give additional qualitative slices with and without lesions, cf. Figs. 11
and 12, and the associated sensitivity plots in Figs. 9 and 10.

Figure 3: Comparisons of single slice reconstructions with the PET-DDS MR
guided vs. unguided at noise level 2.5 without lesions.

###  5.2 3D Reconstruction

Full 3D reconstructions were analysed for two tracers with simulated lesions.
We evaluate the performance of PET-DDS with additional RDP regularisation in
the $z$-direction perpendicular to axial slices (termed RDPz), and introduce
subset-based data consistency updates as in Eq. (28). Acceleration of PET-DDS
was obtained through the use of subset-based data consistency updates, see
Table 4. For further experiments 28 subsets were used. We compare against a
BSREM computed MAP solution with RDP, and DIP with RDP, similar to Singh et
al. (2023). In Fig. 4 we show sensitivity plots for the FDG tracer and in Fig.
5 we plot the axial, coronal and sagittal slices centred on the lesion
location. Additionally, sensitivity curves for the amyloid tracer are given in
Fig. 6, and the associated reconstructions are available in Appendix B, see
Fig. 15.

Table 4: Computational time for 3D PET-DDS with different numbers of subsets. Number of subsets | 1 | 4 | 7 | 14 | 28 | 42
---|---|---|---|---|---|---
Time for reconstruction (mins) | 47.8 | 13.6 | 8.6 | 5.1 | 3.4 | 2.8

The FDG tracer sensitivity plot in Fig. 5 shows that adding RDPz into PET-DDS
improves SSIM and CRC metrics, while classical RDP provides highest PSNR
values. Since PSNR is computed using a mean squared reconstruction error, the
resulting metric is biased toward blurrier reconstructions. This can be
observed in the qualitative images given in Fig. 5, where RDP gives high PSNR
values while the image insets show excessive blurring on the lesion. PET-DDS
without RDPz performs worse than with RDPz, since the score model only acts on
axial slices and, without RDPz, consistency in $z$-direction is only ensured
through data consistency. Qualitatively this can be observed in Fig. 5, where
coronal and sagittal slices display discontinuities in the $z$-direction
whereas the axial slice is smoother. DIP reconstructions give improvements in
SSIM and CRC as compared to classical RDP results, but fail to improve PSNR.
Results with OOD Amyloid tracer show milder improvements with PET-DDS, with
trends similar to those seen with the FDG tracer.

Figure 4: Results for 3D reconstruction using the FDG tracer for different penalty values. PET-DDS-RDPz $\beta=21.9$, and DIP+RDP $\beta=0.1$. Standard deviation is across reconstructions from different realisations of measurements. Table 5: Results using the best hyperparameters for each method for 3D BrainWeb data with FDG and Amyloid tracers.  |  | PSNR, $\lambda$ | SSIM, $\lambda$ | CRC, $\lambda$
---|---|---|---|---
FDG Tracer | RDP |  $25.74$, $1.81$ |  $0.911$, $2.77$ |  $0.994$, $0.5$
DIP+RDP |  $25.26$, $9,800$ |  $0.917$, $10,800$ |  $0.966$, $9,500$
PET-DDS |  $24.83$, $398$ |  $0.910$, $398$ |  $1.01$, $158$
PET-DDS+RDPz |  $25.70$, $158$ |  $0.922$, $63.1$ |  $0.996$, $158$
Amyloid | RDP |  $24.15$, $2.77$ |  $0.898$, $1.81$ |  $0.996$, $0.5$
DIP+RDP |  $24.10$, $10,200$ |  $0.894$, $10,800$ |  $0.964$, $9,500$
PET-DDS |  $23.08$, $1000$ |  $0.890$, $398$ |  $1.009$, $10$
PET-DDS+RDPz |  $24.15$, $398$ |  $0.906$, $158$ |  $0.999$, $10$
Figure 5: 3D reconstruction for the different methods with FDG tracer, and
metrics computed on the inset lesion. Figure 6: Results for 3D reconstruction
using the Amyloid tracer for different penalty values. PET-DDS-RDPz
$\beta=21.9$, and DIP+RDP $\beta=0.1$. Standard deviation is across
reconstructions from different realisations of measurements.

##  6 Conclusion

In this work we adapt SGMs for PET image reconstruction by incorporating PET
specific constraints, e.g. Poisson noise and non-negativity, into several
popular sampling techniques. We further introduce a measurement-based
normalisation technique, to improve the generalisability to different
intensity values by stabilising the dynamic range encountered by the score
model. In future work, reflected SGMs, recently proposed by Lou and Ermon
(2023), could be leveraged to introduce non-negativity into the sampling
procedure in a more principled manner. This work provides a first
investigation of the generalisation capabilities by training the score model
on patient-realistic slices without lesions and testing on slices with
lesions. However, further work is needed to comprehensively evaluate the
generalisation performance on in-vivo data, and investigate the biases of
SGMs, which is vitally important for clinical adoption. The proposed SGM
sampling methods can produce multiple samples from the posterior
$p({\mathbf{x}}|{\mathbf{y}})$, and in this vein one can draw multiple samples
from the posterior for empirical uncertainty estimation; this is left for
future work. This work proposes guided SGM reconstruction with an additional
MR guidance image using CFG. The preliminary results are promising and further
validation is required. A clinically pertinent investigation into robustness
to misregistration of the MR image could be investigated. Furthermore,
guidance could be extended to a joint PET-MRI reconstruction. Recently, Levac
et al. (2023) used similar ideas for a joint reconstruction of multi-contrast
MR images.

Acknowledgments

I.R.D. Singh and R. Barbano are supported by the EPSRC-funded UCL Centre for
Doctoral Training in Intelligent, Integrated Imaging in Healthcare (i4Health)
(EP/S021930/1) and the Department of Health’s NIHR-funded Biomedical Research
Centre at University College London Hospitals. A. Denker acknowledges the
support by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) - Project number 281474342/GRK2224/2. Ž. Kereta was supported by
the UK EPSRC grant EP/X010740/1. B. Jin and S. Arridge were supported by the
UK EPSRC EP/V026259/1. Software used in this project is partially maintained
by CCP SyneRBI (EPSRC EP/T026693/1). P. Maass acknowledges support by DFG-NSFC
project M-0187 of the Sino-German Center mobility programme. The authors thank
Georg Schramm for his help with ParallelProj.

Ethical Standards

The work follows appropriate ethical standards in conducting research and
writing the manuscript, following all applicable laws and regulations
regarding treatment of animals or human subjects.

Conflicts of Interest

We declare we don’t have conflicts of interest.

## References

  * Adler and Öktem (2018) Jonas Adler and Ozan Öktem.  Learned primal-dual reconstruction.  _IEEE Transactions on Medical Imaging_ , 37(6):1322–1332, 2018.
  * Ahn and Fessler (2003) Sangtae Ahn and Jeffrey A Fessler.  A globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms.  _IEEE Transactions on Medical Imaging_ , 22(5):623–626, 2003.
  * Anderson (1982) Brian D O Anderson.  Reverse-time diffusion equation models.  _Stochastic Processes and their Applications_ , 12(3):313–326, 1982.
  * Aubert-Broche et al. (2006) Bérengère Aubert-Broche, M Griffin, G Bruce Pike, Alan C Evans, and D Louis Collins.  Twenty new digital brain phantoms for creation of validation image data bases.  _IEEE Transactions on Medical Imaging_ , 25(11):1410–1416, 2006.
  * Bai et al. (2013) Bing Bai, Quanzheng Li, and Richard M Leahy.  Magnetic resonance-guided positron emission tomography image reconstruction.  In _Seminars in Nuclear Medicine_ , volume 43, pages 30–44. Elsevier, 2013.
  * Bengio et al. (2013) Yoshua Bengio, Aaron Courville, and Pascal Vincent.  Representation learning: A review and new perspectives.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 35(8):1798–1828, 2013.
  * Chung and Ye (2022) Hyungjin Chung and Jong Chul Ye.  Score-based diffusion models for accelerated MRI.  _Medical Image Analysis_ , 80:102479, 2022.
  * Chung et al. (2023a) Hyungjin Chung, Jeongsol Kim, Michael T McCann, Marc L Klasky, and Jong Chul Ye.  Diffusion posterior sampling for general noisy inverse problems.  In _Proceedings of the Eleventh International Conference on Learning Representations_. ICLR, 2023a.
  * Chung et al. (2023b) Hyungjin Chung, Suhyeon Lee, and Jong Chul Ye.  Fast diffusion sampler for inverse problems by geometric decomposition.  _CoRR_ , abs/2303.05754, 2023b.
  * Chung et al. (2023c) Hyungjin Chung, Dohoon Ryu, Michael T McCann, Marc L Klasky, and Jong Chul Ye.  Solving 3D inverse problems using pre-trained 2D diffusion models.  In _Proceedings of the Thirty Seventh IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 22542–22551. CVPR, 2023c.
  * Collins et al. (1998) D Louis Collins, Alex P Zijdenbos, Vasken Kollokian, John G Sled, Noor Jehan Kabani, Colin J Holmes, and Alan C Evans.  Design and construction of a realistic digital brain phantom.  _IEEE Transactions on Medical Imaging_ , 17(3):463–468, 1998.
  * Dhariwal and Nichol (2021) Prafulla Dhariwal and Alexander Q Nichol.  Diffusion models beat GANs on image synthesis.  In _Proceedings of the Thirty Fourth Conference on Advances in Neural Information Processing Systems_ , pages 8780–8794. NeurIPS, 2021.
  * Dimakis (2022) Alexandros G Dimakis.  _Mathematical Aspects of Deep Learning_ , chapter Deep Generative Models and Inverse Problems, page 400–421.  Cambridge University Press, Cambridge, UK, 2022.
  * Efron (2011) Bradley Efron.  Tweedie’s formula and selection bias.  _Journal of the American Statistical Association_ , 106(496):1602–1614, 2011.
  * Ehrhardt (2021) Matthias J Ehrhardt.  _Handbook of Mathematical Models and Algorithms in Computer Vision and Imaging: Mathematical Imaging and Vision_ , chapter Multi-modality Imaging with Structure-Promoting Regularizers, pages 1–38.  Springer International Publishing, Cham, Switzerland, 2021.  ISBN 978-3-030-03009-4.
  * Ehrhardt et al. (2016) Matthias J Ehrhardt, Pawel Markiewicz, Maria Liljeroth, Anna Barnes, Ville Kolehmainen, John S Duncan, Luis Pizarro, David Atkinson, Brian F Hutton, Sébastien Ourselin, Kris Thielemans, and Simon R Arridge.  PET reconstruction with an anatomical MRI prior using parallel level sets.  _IEEE Transactions on Medical Imaging_ , 35(9):2189–2199, 2016.
  * Engl et al. (1996) Heinz W Engl, Martin Hanke, and Andreas Neubauer.  _Regularization of inverse problems_.  Mathematics and its Applications. Kluwer Academic Publishers Group, Dordrecht, 1996.
  * Gong et al. (2019) Kuang Gong, Ciprian Catana, Jinyi Qi, and Quanzheng Li.  PET image reconstruction using deep image prior.  _IEEE Transactions on Medical Imaging_ , 38(7):1655–1665, 2019.
  * Gong et al. (2022) Kuang Gong, Keith A. Johnson, Georges El Fakhri, Quanzheng Li, and Tinsu Pan.  PET image denoising based on denoising diffusion probabilistic models.  _CoRR_ , abs/2209.06167, 2022.
  * Goodfellow et al. (2014) Ian J Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron C Courville, and Yoshua Bengio.  Generative adversarial nets.  In _Proceedings of the Second Conference on Advances in Neural Information Processing Systems_ , pages 2672–2680. NeurIPS, 2014.
  * Guazzo and Colarieti-Tosti (2021) Alessandro Guazzo and Massimiliano Colarieti-Tosti.  Learned primal dual reconstruction for PET.  _Journal of Imaging_ , 7(12), 2021.
  * Herman and Meyer (1993) Gabor T Herman and Lorraine B Meyer.  Algebraic reconstruction techniques can be made computationally efficient (positron emission tomography application).  _IEEE Transactions on Medical Imaging_ , 12(3):600–609, 1993.
  * Ho and Salimans (2022) Jonathan Ho and Tim Salimans.  Classifier-free diffusion guidance.  _CoRR_ , abs/2207.12598, 2022.
  * Ho et al. (2020) Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  In _Proceedings of the Thirty Third Conference on Advances in Neural Information Processing Systems_. NeurIPS, 2020.
  * Hudson and Larkin (1994) H Malcolm Hudson and Richard S Larkin.  Accelerated image reconstruction using ordered subsets of projection data.  _IEEE Transactions on Medical Imaging_ , 13(4):601–609, 1994.
  * Ito and Jin (2015) Kazufumi Ito and Bangti Jin.  _Inverse problems: Tikhonov theory and algorithms_.  World Scientific, Hackensack, NJ, 2015.  ISBN 978-981-4596-19-0.
  * Jalal et al. (2021) Ajil Jalal, Marius Arvinte, Giannis Daras, Eric Price, Alexandros G. Dimakis, and Jonathan I. Tamir.  Robust compressed sensing MRI with deep generative priors.  In _Proceedings of the Thirty Fourth Conference on Advances in Neural Information Processing Systems_ , pages 14938–14954. NeurIPS, 2021.
  * Jin et al. (2017) Kyong H Jin, Michael T McCann, Emmanuel Froustey, and Michael Unser.  Deep convolutional neural network for inverse problems in imaging.  _IEEE Transactions on Image Processing_ , 26(9):4509–4522, 2017.
  * Kaplan and Zhu (2019) Sydney Kaplan and Yang-Ming Zhu.  Full-dose PET image estimation from low-dose PET image using deep learning: a pilot study.  _Journal of Digital Imaging_ , 32(5):773–778, 2019.
  * Karlberg et al. (2016) Anna M Karlberg, Oddbjørn Sæther, Live Eikenes, and Pål Erik Goa.  Quantitative comparison of PET performance—siemens biograph mCT and mMR.  _European Journal of Nuclear Medicine and Molecular Imaging: Physics_ , 3(1), 2016.
  * Kingma and Welling (2014) Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  In _Proceedings of the Second International Conference on Learning Representations_. ICLR, 2014.
  * Kobler and Pock (2023) Erich Kobler and Thomas Pock.  Learning gradually non-convex image priors using score matching.  _CoRR_ , abs/2302.10502, 2023.
  * Leuschner et al. (2021) Johannes Leuschner, Maximilian Schmidt, Daniel Otero Baguer, David Erzmann, and Mateus Baltazar.  DIV$\alpha$l library, 2021.  URL https://doi.org/10.5281/zenodo.4428220.
  * Levac et al. (2023) Brett Levac, Ajil Jalal, Kannan Ramchandran, and Jonathan I Tamir.  MRI reconstruction with side information using diffusion models.  _arXiv_ , 2023.  URL https://arxiv.org/abs/2303.14795.
  * Li et al. (2022) Xiang Li, John Thickstun, Ishaan Gulrajani, Percy Liang, and Tatsunori B Hashimoto.  Diffusion-LM improves controllable text generation.  In _Proceedings to the Thirty Fifth Conference on Advances in Neural Information Processing Systems_ , pages 4328–4343. NeurIPS, 2022.
  * Lou and Ermon (2023) Aaron Lou and Stefano Ermon.  Reflected diffusion models.  In _Proceedings of the Fortieth International Conference on Machine Learning_ , volume 202, pages 22675–22701. PMLR, 2023.
  * Mehranian and Reader (2021) Abolfazl Mehranian and Andrew J Reader.  Model-based deep learning PET image reconstruction using forward–backward splitting expectation–maximization.  _IEEE Transactions on Radiation and Plasma Medical Sciences_ , 5(1):54–64, 2021.
  * Nuyts et al. (2002) Johan Nuyts, Dirk Bequé, Patrick Dupont, and Luc Mortelmans.  A concave prior penalizing relative differences for maximum-a-posteriori reconstruction in emission tomography.  _IEEE Transactions on Nuclear Science_ , 49(1):56–60, 2002.
  * Ote et al. (2023) Kibo Ote, Fumio Hashimoto, Yuya Onishi, Takashi Isobe, and Yasuomi Ouchi.  List-mode PET image reconstruction using deep image prior.  _IEEE Transactions on Medical Imaging_ , 42(6):1822–1834, 2023.
  * Ovtchinnikov et al. (2020) Evgueni Ovtchinnikov, Richard Brown, Christoph Kolbitsch, Edoardo Pasca, Casper O da Costa-Luis, Ashley Gillman, Benjamin A Thomas, Nikos Efthimiou, Johannes Mayer, Palak Wadhwa, Matthias J Ehrhardt, Sam Ellis, Jakob S Jørgensen, Julian C Matthews, Claudia Prieto, Andrew J Reader, Charalampos Tsoumpas, Martin J Turner, and Kris Thielemans.  SIRF: synergistic image reconstruction framework.  _Computer Physics Communications_ , 249:107087, 2020.
  * Pain et al. (2022) Cameron D Pain, Gary F Egan, and Zhaolin Chen.  Deep learning-based image reconstruction and post-processing methods in positron emission tomography for low-dose imaging and resolution enhancement.  _European Journal of Nuclear Medicine and Molecular Imaging_ , 49(9):3098–3118, 2022.
  * Pierro and Yamagishi (2001) Alvaro R De Pierro and Michel E B Yamagishi.  Fast EM-like methods for maximum ‘a posteriori’ estimates in emission tomography.  _IEEE Transactions on Medical Imaging_ , 20(4):280–288, 2001.
  * Pinaya et al. (2022) Walter H L Pinaya, Petru-Daniel Tudosiu, Jessica Dafflon, Pedro F Da Costa, Virginia Fernandez, Parashkev Nachev, Sébastien Ourselin, and M Jorge Cardoso.  Brain imaging generation with latent diffusion models.  In _Proceedings of the Twenty Fifth Conference on Medical Image Computing and Computer Assisted Intervention_ , pages 117–126. MICCAI, 2022.
  * Qi and Leahy (2006) Jinyi Qi and Richard M Leahy.  Iterative reconstruction techniques in emission computed tomography.  _Physics in Medicine & Biology_, 51(15):R541, 2006.
  * Ramzi et al. (2020) Zaccharie Ramzi, Benjamin Remy, François Lanusse, Jean-Luc Starck, and Philippe Ciuciu.  Denoising score-matching for uncertainty quantification in inverse problems.  _CoRR_ , abs/2011.08698, 2020.
  * Ravula et al. (2023) Sriram Ravula, Brett Levac, Ajil Jalal, Jonathan I Tamir, and Alexandros G Dimakis.  Optimizing sampling patterns for compressed sensing MRI with diffusion generative models.  _CoRR_ , abs/2306.03284, 2023.
  * Reader et al. (2021) Andrew J Reader, Guillaume Corda, Abolfazl Mehranian, Casper da Costa-Luis, Sam Ellis, and Julia A Schnabel.  Deep learning for PET image reconstruction.  _IEEE Transactions on Radiation and Plasma Medical Sciences_ , 5(1):1–25, 2021.
  * Ronneberger et al. (2015) Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _Proceedings of the Eighteenth Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 234–241. MICCAI, 2015.
  * Rudin et al. (1992) Leonid I Rudin, Stanley Osher, and Emad Fatemi.  Nonlinear total variation based noise removal algorithms.  _Physica D: Nonlinear Phenomena_ , 60(1-4):259–268, 1992.
  * Saharia et al. (2022) Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily L Denton, Seyed Kamyar Seyed Ghasemipour, Raphael Gontijo Lopes, Burcu Karagol Ayan, Tim Salimans, Jonathan Ho, David J Fleet, and Mohammad Norouzi.  Photorealistic text-to-image diffusion models with deep language understanding.  In _Proceedings of the Thirty Fifth Conference on Advances in Neural Information Processing Systems_ , pages 36479–36494. NeurIPS, 2022.
  * Schramm (2021) Georg Schramm.  Simulated brainweb PET/MR data sets for denoising and deblurring, 2021\.  URL https://zenodo.org/record/4897350.
  * Schramm (2022) Georg Schramm.  PARALLELPROJ – an open-source framework for fast calculation of projections in tomography, 2022.  URL https://arxiv.org/abs/2212.12519.
  * Shepp and Vardi (1982) Lawrence A Shepp and Yehuda Vardi.  Maximum likelihood reconstruction for emission tomography.  _IEEE Transactions on Medical Imaging_ , 1(2):113–122, 1982.
  * Singh et al. (2023) Imraj R D Singh, Riccardo Barbano, Željko Kereta, Bangti Jin, Kris Thielemans, and Simon Arridge.  3D PET-DIP reconstruction with relative difference prior using a SIRF-based objective.  In _Proceedings to the Seventeenth International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine_. Fully3D, 2023.
  * Sohl-Dickstein et al. (2015) Jascha Sohl-Dickstein, Eric A Weiss, Niru Maheswaranathan, and Surya Ganguli.  Deep unsupervised learning using nonequilibrium thermodynamics.  In _Proceedings of the Thirty Second International Conference on Machine Learning_ , pages 2256–2265, 2015.
  * Somayajula et al. (2011) Sangeetha Somayajula, Christos Panagiotou, Anand Rangarajan, Quanzheng Li, Simon R Arridge, and Richard M Leahy.  PET image reconstruction using information theoretic anatomical priors.  _IEEE Transactions on Medical Imaging_ , 30(3):537–549, 2011.
  * Song et al. (2021a) Jiaming Song, Chenlin Meng, and Stefano Ermon.  Denoising diffusion implicit models.  In _Proceedings of the Ninth International Conference on Learning Representations_. ICLR, 2021a.
  * Song and Ermon (2019) Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  In _Proceedings of the Thirty Second Conference on Advances in Neural Information Processing Systems_ , pages 11895–11907. NeurIPS, 2019.
  * Song et al. (2021b) Yang Song, Conor Durkan, Iain Murray, and Stefano Ermon.  Maximum likelihood training of score-based diffusion models.  In _Proceedings of the Thirty Fourth Conference on Advances in Neural Information Processing Systems_ , pages 1415–1428. NeurIPS, 2021b.
  * Song et al. (2021c) Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  In _Proceedings of the Ninth International Conference on Learning Representations_. ICLR, 2021c.
  * Song et al. (2022) Yang Song, Liyue Shen, Lei Xing, and Stefano Ermon.  Solving inverse problems in medical imaging with score-based generative models.  In _Proceedings of the Tenth International Conference on Learning Representations_. ICLR, 2022.
  * Särkkä and Solin (2019) Simo Särkkä and Arno Solin.  _Applied Stochastic Differential Equations_.  Institute of Mathematical Statistics Textbooks. Cambridge University Press, Cambridge, 2019.  ISBN 978-131-6649-46-6.
  * Thielemans et al. (2012) Kris Thielemans, Charalampos Tsoumpas, Sanida Mustafovic, Tobias Beisel, Pablo Aguiar, Nikolaos Dikaios, and Matthew W Jacobson.  STIR: software for tomographic image reconstruction release 2.  _Physics in Medicine & Biology_, 57(4):867, 2012\.
  * Tong et al. (2010) Shan Tong, Adam M Alessio, and Paul E Kinahan.  Noise and signal properties in PSF-based fully 3D PET image reconstruction: an experimental evaluation.  _Physics in Medicine & Biology_, 55(5):1453, 2010.
  * Tran et al. (2020) Dustin Tran, Jasper Snoek, and Balaji Lakshminarayanan.  Practical uncertainty estimation and out-of-distribution robustness in deep learning.  NeurIPS Tutorial, Google Brain, 2020.  URL https://neurips.cc/media/neurips-2020/Slides/16649.pdf.
  * Ulyanov et al. (2018) Dmitry Ulyanov, Andrea Vedaldi, and Victor S Lempitsky.  Deep image prior.  In _Proceedings of the Thirty First IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 9446–9454. CVPR, 2018.
  * Vincent (2011) Pascal Vincent.  A connection between score matching and denoising autoencoders.  _Neural Computation_ , 23(7):1661–1674, 2011\.
  * Wang et al. (2004) Zhou Wang, Alan C Bovik, Hamid R Sheikh, and Eero P Simoncelli.  Image quality assessment: from error visibility to structural similarity.  _IEEE Transactions on Image Processing_ , 13(4):600–612, 2004.
  * Xiao et al. (2022) Zhisheng Xiao, Karsten Kreis, and Arash Vahdat.  Tackling the generative learning trilemma with denoising diffusion GANs.  In _Proceedings of the Tenth International Conference on Learning Representations_. ICLR, 2022.
  * Zhu et al. (2023) Yuanzhi Zhu, Kai Zhang, Jingyun Liang, Jiezhang Cao, Bihan Wen, Radu Timofte, and Luc Van Gool.  Denoising diffusion models for plug-and-play image restoration.  _CoRR_ , abs/2305.08995, 2023.

##  A Appendix

###  A.1 Training of Score-based Models

##### Forward SDE

In our experiments, we make use of the variance preserving SDE (Ho et al.,
2020)

$\displaystyle
d{\mathbf{x}_{t}}=-\frac{\beta(t)}{2}{\mathbf{x}_{t}}dt+\sqrt{\beta(t)}d\mathbf{w},$
(37)

were we employ
$\beta(t)=\beta_{\text{min}}+t(\beta_{\text{max}}-\beta_{\text{min}})$ as a
linear schedule with $\beta_{\text{min}}=0.1$ and $\beta_{\text{max}}=10$ with
a terminal time $T=1$. The coefficients were chosen such that the terminal
distribution approximates a Gaussian, i.e.
$p_{1}({\mathbf{x}})\approx\mathcal{N}(0,I)$. We also tested the Variance
Exploding (VE) SDE (Song et al., 2021c); it was found that VE-SDE was more
unstable than VP-SDE for PET image reconstruction. The transition kernel for
the variance preserving SDE is a Gaussian, i.e.
$p_{t}({\mathbf{x}_{t}}|{\mathbf{x}_{0}})=\mathcal{N}({\mathbf{x}_{t}};\gamma_{t}{\mathbf{x}_{0}},\nu_{t}^{2}I)$,
with coefficients

$\displaystyle\gamma_{t}=\exp{\left(-\frac{1}{2}\int_{0}^{t}\beta(s)ds\right)},\quad\nu_{t}^{2}=1-\exp{\left(-\int_{0}^{t}\beta(s)ds\right)}.$

Using this closed form expression for the transition kernel, the denoising
score matching loss can be rewritten as

$\displaystyle L_{\text{DSM}}(\theta)=\mathbb{E}_{t\sim
U[0,1]}\mathbb{E}_{{\mathbf{x}_{0}}\sim{\pi}}\mathbb{E}_{{\mathbf{z}}\sim\mathcal{N}(0,I)}\left[\omega_{t}\left\|s_{\theta}({\mathbf{x}_{t}},t)+\frac{{\mathbf{z}}}{\nu_{t}}\right\|_{2}^{2}\right],$
(38)

with ${\mathbf{x}}_{t}=\gamma_{t}{\mathbf{x}}_{0}+\nu_{t}{\mathbf{z}}$. The
weighting $\omega_{t}$ is chosen as $\omega_{t}=\nu_{t}^{2}$ to approximate
maximum likelihood training (Song et al., 2021b).

##### Model Architecture

We use the architecture proposed by Dhariwal and Nichol (2021)222available at
https://github.com/openai/guided-diffusion. The architecture is based on the
popular U-Net architecture (Ronneberger et al., 2015) consisting of a decoder
implemented as a stack of residual blocks and downsampling operations and an
encoder of residual blocks and upsampling operations. At the lowest resolution
($8\times 8$), additional global attention layers are used. To incorporate the
timestep into each residual block, the authors use adaptive group
normalisation (AdaGN) layers defined as
$\text{AdaGN}(h,e)=e_{s}\text{GroupNorm}(h)+e_{b}$, where $h$ are intermediate
features and $e=[e_{s},e_{b}]$ is the encoded time step. The specific
implementation and the choice of our hyperparameters can be seen in our
github. For the MRI guided model we apply the clean MRI image as an additional
channel to the input of the network.

###  A.2 Experimental Details

The sampling methods presented in Section 3 use different penalty strengths in
order to scale the likelihood term for PET-Naive and PET-DPS or to set the
strength of the additional Tikhonov regularisation for PET-DDS. For Naive it
is recommended to choose $\lambda_{t}^{\text{naive}}$ s.t. the penalty is zero
at the start of sampling and increased as $t\to 0$ (Jalal et al., 2021). We
use $\lambda_{t}=\lambda(1-t)$ in all our experiments. For the PET-DPS
approach (Chung et al., 2023a) define the sampling iteration as

$\begin{split}&\tilde{{\mathbf{x}}}_{t_{k-1}}={\mathbf{x}}_{t_{k}}+[{\mathbf{f}}({\mathbf{x}}_{t_{k}},t_{k})-g(t_{k})^{2}s_{\theta}({\mathbf{x}}_{t_{k}},t_{k})]\Delta
t+g(t_{k})\sqrt{|\Delta
t|}{\mathbf{z}}\quad{\mathbf{z}}\sim\mathcal{N}(0,I),\\\
&{\mathbf{x}}_{t_{k-1}}=\tilde{{\mathbf{x}}}_{t_{k-1}}-\lambda_{t_{k}}^{\text{DPS}}\nabla_{\mathbf{x}}L({\mathbf{y}}|{\hat{\mathbf{x}}_{0}}({\mathbf{x}}_{t_{k}})).\end{split}$
(39)

which is equivalent to the classical Euler-Maruyama scheme, when
$\lambda_{t_{k}}^{\text{DPS}}$ is chosen in such a way that it incorporates
the step size $\Delta t$ and the diffusion function $g(t_{k})^{2}$. We follow
Chung et al. (2023a) and define
$\lambda_{t}^{\text{DPS}}=\frac{\lambda}{D_{KL}(A{\hat{\mathbf{x}}_{0}}||{\mathbf{y}})}$.
For PET-DDS a constant penalty $\lambda^{\text{DDS}}$, without time
dependency, is used. Heuristically, it was found that the number of iterations
used for data-consistency projection were adjust such that the results, with
$\lambda^{\text{DDS}}=0$, overfit to noise. The penalty strength
$\lambda^{\text{DDS}}$ was then increased to regularise the reconstruction
more. In 2D the number of projection steps for PET-DDS were set to $4$ for
noise level $2.5$ and $15$ for noise level $10$.

###  A.3 Baseline Methods

##### Classical Methods

Relative Difference Prior (RDP) is a common penalty for PET reconstruction
(Nuyts et al., 2002), defined by

$\displaystyle R({\mathbf{x}})=-\sum_{j=1}^{n}\sum_{k\in
N_{j}}\frac{(x_{j}-x_{k})^{2}}{x_{j}+x_{k}+\xi\lvert x_{j}-x_{k}\rvert},$

where $N_{j}$ is a pre-defined neighbourhood around $x_{j}$, typically
$3\times 3$ in 2D or $3\times 3\times 3$ in 3D. $R_{z}({\mathbf{x}})$ is a
variant of RDP whereby the neighbourhood is defined in the axial dimensions in
3D, i.e. $3\times 1\times 1$. A Neumann boundary was chosen where
neighbourhoods that were outside of the domain. The tunable parameter $\xi>0$
controls the degree of edge-preservation ($\xi=1$, in-line with clinical
practice), and with its gradient given by

$\displaystyle\frac{\partial R({\mathbf{x}})}{\partial x_{j}}=\sum_{k\in
N_{j}}-\frac{(r_{jk}-1)(\xi|r_{jk}-1|+r_{jk}+3)}{(r_{jk}+1+\xi|r_{jk}-1|)^{2}},\quad\text{with
}r_{jk}:=\frac{x_{j}}{x_{k}}.$ (40)

The penalisation is scale-invariant since the gradient is computed using the
ratio of voxel values $r_{jk}$. This partially overcomes the issue with the
wide dynamic range observed in emission tomography images. For BSREM algorithm
the convergence criteria was set based on the change of voxel values within
the reconstruction between iterates. Specifically, the change in mean voxel
values across non-zero voxel values was less than $0.01\%$, we set the
relaxation coefficient to $\zeta=0.1$.

##### PET Image Denoising with SGM

In PET image denoising, the goal is to sample from the posterior
$p({\mathbf{x}}|\mathbf{x}_{\text{noisy}})$ of the true image ${\mathbf{x}}$
given an initial (low-count) reconstruction $\mathbf{x}_{\text{noisy}}$. This
is differs from PET reconstruction, where the goal is to sample from the
posterior ${p^{\text{post}}}({\mathbf{x}}|{\mathbf{y}})$ conditioned on the
measurements ${\mathbf{y}}$. In this framework the denoising likelihood is
given by Gaussian noise, i.e.

$\displaystyle
p(\mathbf{x}_{\text{OSEM}}|{\mathbf{x}})=\mathcal{N}(\mathbf{x}_{\text{noisy}};{\mathbf{x}},\sigma_{d}^{2}I),$
(41)

with the noise level $\sigma_{d}$ to be specified. Using the Naive
approximation, we get the following reverse SDE for the PET denoising
likelihood

$\displaystyle
d{\mathbf{x}_{t}}=[{\mathbf{f}}({\mathbf{x}_{t}},t)-g(t)^{2}(s_{\theta}({\mathbf{x}_{t}},t)-1/\sigma_{d}^{2}(\mathbf{x}_{\text{noisy}}-{\mathbf{x}_{t}})]dt+g(t)d\bar{\mathbf{w}}_{t}.$
(42)

In our implementation we estimate the initial reconstruction using OSEM with
34 subsets and iterations (i.e. 1 epoch). The same score model
$s_{\theta}({\mathbf{x}_{t}},t)$ is used for both PET denoising and
reconstruction. The noise level $\sigma_{d}$ is chosen based on a held-out
evaluation dataset.

##### Supervised Learning

We are using two popular supervised learning techniques: post-processing and
learned iterative methods. For the post-processing method we used a variant of
the FBPConvNet (Jin et al., 2017), modified to PET reconstruction. The input
to the FBPConvNet was changed to an OSEM with 34 subsets and iterations, this
variation is denoted by PET-UNet. For the learned iterative method, we adopt
Learned Primal Dual (LPD) (Adler and Öktem, 2018), referred to as PET-LPD. For
PET-LPD we use the same OSEM reconstruction as initialisation for the primal
channels and include the affine forward model with sample specific attenuation
maps. Note that these sample specific factors were not included in previous
implementation of learned iterative methods for PET image reconstruction
(Guazzo and Colarieti-Tosti, 2021). Three primal and dual unrolled iterations
were used. Both of these networks were implemented using Div$\alpha$L
(Leuschner et al., 2021) with only minimal changes to the architecture; PET-
UNet was a UNet with $1\,783\,249$ parameters, and PET-LPD used a block of
convolutional filters for each primal and dual network with a total of
$132\,300$ parameters. Both networks were trained using the dataset in Section
4.1 without lesions and noise levels of $5$, $10$, and $50$. The dataset was
split into training and evaluation, and training was terminated when over-
fitting was observed. Additionally, data-corrected mean normalisation was
included to promote generalisability between noise levels. The code for these
supervised learning models is publicly available at https://github.com/Imraj-
Singh/pet_supervised_normalisation.

##### Deep Image Prior

The Deep Image Prior (DIP) (Ulyanov et al., 2018) is a popular framework for
unsupervised image reconstruction, relying only on a single measurement. A
common problem of the DIP is its tendency to overfit to noise. Therefore some
regularisation has to be used. We included RDP into the objective function to
elevate the need for early-stopping and prevent over-fitting to noise. The
architecture used was a three-scale U-Net (Ronneberger et al., 2015) with
$1\,606\,899$ parameters, with a rectified linear unit on the output to ensure
non-negativity. This architecture is minimally changed from previous
applications of DIP to PET (Gong et al., 2019; Ote et al., 2023; Singh et al.,
2023). DIP results are computed on reconstructions along the optimisation
trajectory, every 100 iterations from 6,600 iterations to 11,600.

###  A.4 Evaluation metrics

In addition to peak-signal-to-noise ratio (PSNR) and structural similarity
index measure (SSIM) (Wang et al., 2004), we compute two local quality scores
over a Region of Interest (ROI). For reconstructions with lesions a Contrast
Recovery Coefficient (CRC) was computed to quantify detectability of these
local features. This was computed between lesion $L$ and background $B$ ROIs,
these have $N_{L}$ and $N_{B}$ number of elements respectively333We use $L$ to
denote the lesion ROI in this section only; in the main manuscript $L$ is the
likelihood.. Additionally, there are $R$ realisations of the measured data.
Given an ROI $Z$, we include subscript indices for element and realisation
$Z_{r,k}$, where $r$ is the realisation index, and $k$ is the element index.
An average over the elements of the ROI is denoted as
$\bar{Z}_{r}=\frac{1}{N_{Z}}\sum_{k=1}^{N_{Z}}Z_{r,k}$. The CRC is defined by

$\mathrm{CRC}:=\sum_{r=1}^{R}\left(\frac{\bar{L}_{r}}{\bar{B}_{r}}-1\right)/\left(\frac{L_{\mathrm{t}}}{B_{\mathrm{t}}}-1\right),$
(43)

where the subscript $\rm t$ denotes the ground truth ROIs. We study the noise
over realisations of the measured data using normalised STD (also referred to
as ensemble noise, see Tong et al. (2010), and is reported to give a true
estimate of noise in the image). We define an average over realisations of the
ROI as $\bar{Z}_{k}=\frac{1}{R}\sum_{r=1}^{R}Z_{r,k}$ where STD is computed on
background ROIs it is given by:

$\mathrm{STD}:=\frac{1}{N_{B}}\sum_{k=1}^{N_{B}}\sqrt{\frac{1}{R-1}\sum^{R}_{r=1}\frac{(B_{r,k}-\bar{B}_{k})^{2}}{\bar{B}_{k}}}.$
(44)

For reconstructions with lesions the background ROI was used, and without
lesions a background of the whole emission volume (defined on reference
images) was used. In 2D $R=10$ noise realisations of acquisition data were
used, and $R=5$ in 3D. To evaluate the consistency of our reconstructions to
the true measurements, we compute the Kullback-Leibler divergence (KLDIV)

$\displaystyle\mathrm{KLDIV}:=\sum_{j=1}^{m}\bar{y}_{j}\log\left(\frac{\bar{y}_{j}}{y_{j}}\right)-\bar{y}_{j}+y_{j},$
(45)

between measurements ${\mathbf{y}}$ and estimated measurements
$\mathbf{\bar{y}}={\mathbf{A}}\hat{{\mathbf{x}}}+\mathbf{\bar{b}}$ where
$\hat{{\mathbf{x}}}$ denotes the reconstruction.

##  B Additional Results

###  B.1 2D Reconstruction

We show additional sensitivity plots for 2D reconstruction. For noise level
$10$ these results are presented in Fig. 7 and Fig. 8 without and with
lesions, respectively. The results are similar to the settings for noise level
$2.5$, as we see a clear trade-off between reconstruction quality in terms of
PSNR/SSIM and visibility of lesions in terms of CRC in Fig. 8. Here, a higher
regularisation leads to better PSNR/SSIM scores and a lower regularisation to
a better recovery of lesions. A high regularisation, i.e. a high influence of
the score model, may lead to a worse reconstruction of the lesions, as the
score model was trained on images without lesions.

Figure 7: Results for BrainWeb without lesions with noise level 10 for
different penalty parameters. The Standard Deviation is computed over
reconstructions of different noise realisations ${\mathbf{y}}$.  Figure 8:
Results for BrainWeb with lesions with noise level 10 for different penalty
parameters. The Standard Deviation is computed over reconstructions of
different noise realisations ${\mathbf{y}}$.

###  B.2 MR guidance

We show additional results for the MR guided model. Sensitivity plots without
and with lesions are presented in Fig. 9 and Fig. 10. These results support
the findings of the paper, as the MR guided models achieve better
reconstruction quality w.r.t. PSNR and SSIM. However, the CRC is similar to
the unguided model. As the lesions were not visible in the MR image, no
additional information about the lesions are introduced through guidance. We
show two more reconstruction examples without lesions in Fig. 11 and examples
with lesions in Fig. 12.

Figure 9: Results for 2D reconstruction guided vs unguided without lesions for
noise level 2.5. Figure 10: Results for 2D reconstruction guided vs unguided
with lesion for noise level 2.5. Figure 11: Comparisons of single central
slice reconstructions with the PET-DDS MR guided vs. unguided at noise level
2.5 without lesions.

Figure 12: Comparison of the PET-DDS MR guided vs. unguided with a noise level
2.5 with lesions.

###  B.3 3D results RDPz sweeps

We show the sensitivity plots for different penalty values of the additional
RDP regularizer in $z$-direction for PET-DDS in Fig. 13 and 14 for the two
different tracers. In addition, we show axial, coronal and saggital slices of
the reconstruction with the Amyloid tracer in 15.

Figure 13: Results for 3D reconstruction using the FDG tracer for different
penalty values. Figure 14: Results for 3D reconstruction using the Amyloid
tracer for different penalty values. Figure 15: 3D reconstruction for the
different method with Amyloid tracer, and metrics computed on inset lesion.