# Complex-valued Retrievals From Noisy Images Using Diffusion Models

Nadav Torem , Roi Ronen11footnotemark: 1 , Yoav Y. Schechner
Viterbi Faculty of Electrical & Computer Eng.
Technion, Haifa, Israel
<EMAIL_ADDRESS>Equal contribution. Michael Elad
Department of Computer Science
Technion, Haifa, Israel

###### Abstract

In diverse microscopy modalities, sensors measure only real-valued
intensities. Additionally, the sensor readouts are affected by Poissonian-
distributed photon noise. Traditional restoration algorithms typically aim to
minimize the mean squared error (MSE) between the original and recovered
images. This often leads to blurry outcomes with poor perceptual quality.
Recently, deep diffusion models (DDMs) have proven to be highly capable of
sampling images from the a-posteriori probability of the sought variables,
resulting in visually pleasing high-quality images. These models have mostly
been suggested for real-valued images suffering from Gaussian noise. In this
study, we generalize annealed Langevin Dynamics, a type of DDM, to tackle the
fundamental challenges in optical imaging of complex-valued objects (and real
images) affected by Poisson noise. We apply our algorithm to various optical
scenarios, such as Fourier Ptychography, Phase Retrieval, and Poisson
denoising. Our algorithm is evaluated on simulations and biological empirical
data.

##  1 Introduction

Light is a complex-valued field. During imaging, both the phase and intensity
of the field change by the captured objects and the propagation within the
optical system. Typically, the examined objects are assumed to have only real
values. Recovering the phase of complex-valued objects is often crucial [9,
1]. However, image sensors can only measure real-valued non-negative
intensities. This limitation gives rise to inverse problems such as phase-
retrieval [22, 51, 53, 15, 50] ptychography [37, 8], Fourier ptychography [41,
47, 40, 56, 16], coded diffraction imaging [2] and holography. A common
characteristic of all these problems is the nonlinear relationship between the
measurements and the complex-valued unknowns.

Furthermore, the particle nature of light leads to Poissonian-distributed
photon noise, in contrast to the commonly assumed Gaussian noise [21, 29, 55,
7, 35, 27, 34]. Poisson noise complicates the recovery of phase-related
problems as well as linear problems involving real-valued images. Therefore,
addressing Poisson noise during restoration is of great importance, especially
in scientific tasks, such as cell imaging (see Fig. 1).

Figure 1: Poisson denoising. Results of BM3D + I + VST [3] and our algorithm
on real data [55], obtained via two-photon microscopy. [Bottom] A zoom-in. Our
recovery has high-perceptual quality (sharp) and low distortion.  Figure 2:
Imaging of a complex-valued objects. The imaging setup captures only real
valued noisy intensities. Given the noisy measurements, our method
stochastically recovers the objects amplitude and phase.

Over the years, restoration algorithms have varied based on their objective.
These algorithms, whether classical (NLM [6], KSVD [17], EPLL [57], WNNM [20],
and BM3D [30]) or based on deep-learning (TNRD [10], MALA [19], DnCNN [54],
Noise2Void [28] and NLRN [30]), have been typically designed to minimize the
mean squared error (MSE) between true and estimated images. These are MMSE
estimators. Alternatively, other methods pursued maximum a-posteriori
probability (MAP)111Interestingly, also works on these methods tend to report
performance using MSE..

The distortion-perception trade-off [4, 48] states that in image recovery, the
resulting distortion (MSE score) is reciprocal to the resulting perceptual
quality. Hence, there is a growing interest in recovery algorithms that yield
high perceptual quality, rather than minimizing the MSE. Some prior works
along this line use generative adversarial networks (GANs) for image recovery.
GANs can be turned into samplers from the posterior distribution. However,
training difficulties and lack of solid mathematical guarantees hindered these
techniques.

Newcomer generative models are deep diffusion models (DDMs). DDMs are a
highly-performing class of iterative probabilistic generative models. DMMs
have shown success in solving various inverse problems [45, 11, 24, 13, 23,
25, 14]. However, a DMM requires a likelihood term, which is often
analytically intractable. This is exacerbated by the iterative nature of DDMs.
As a result, restoration algorithms that utilize DDMs have primarily
concentrated on linear inverse problems affected by stationary Gaussian noise
[24, 25, 23]. An exception is [12] which proposes DDMs for nonlinear inverse
problems. It has two limitations: (i) Being restricted to real-valued images,
and (ii) Requiring costly back-propagation through a DNN in each iteration.
This lengthens computation time during recovery.

One type of DDM is annealed Langevin Dynamics. In this study, we generalize
this DDM type, to tackle fundamental challenges in optical imaging: recovering
complex-valued objects (and real images) affected by Poisson noise (see Fig.
2). To handle the untraceable likelihood term, we analytically develop a novel
relaxed model, derived from the optical setting and the Poisson noise
distribution. We show that the gradient of the prior can be expressed via an
expectation, which is conveniently approximated by a deep neural network
(DNN). We apply our algorithm to various optical scenarios, such as Fourier
ptychography, phase retrieval, and Poisson denoising. Our method samples from
the posterior distribution, thus exposing uncertainties in the recovery, via
the spread of the results. The uncertainty is valuable when imaging true
objects. We show results on simulations and real empirical biological data.

##  2 Theoretical Background

###  2.1 Noise

Empirical light measurements are random, mainly due to the integer nature of
photons and electric charges. This randomness is often regarded as photon
noise. Dark noise, quantization noise, and photon noise are the dominant noise
sources in common cameras [36]. Denote incorporation of noise into the
expected noiseless signal I by the operator $\cal N$. The noisy measurement is
$y={\cal N}({\texttt{I}})$.

Photon noise is modeled by a Poisson process [36]. A pixel is characterized by
its full well capacity (FWC), that is, the maximal amount of charge that can
be accumulated in the pixel before saturating. Let $x^{\rm e}\in[0,{\rm FWC}]$
be the expected (noiseless) value of a measured signal (photo-electrons). The
corresponding noisy measurement is

$y^{\rm e}\sim{\rm Poisson}(x^{\rm e}).$ (1)

Photon noise is fundamental. Its crucial nature of signal-dependency (1) is
valid both for photon-hungry scenes, and even more so for well-lit scenes, as
the noise variance is higher in brighter pixels.

Modeling the nature of noise, being signal-dependent, is important when
solving scientific inverse problems in optical imaging. However, data analysis
under Poisson noise is difficult to deal with in a closed form. Therefore, in
many inverse problems, Poisson noise is approximated by a Gaussian
distribution, $n^{\rm e}\sim N(0,x^{\rm e})$. The variance still depends on
the noiseless signal. Define normalized variables

$x=\frac{x^{e}}{{\rm FWC}}\in[0,1]\;,~{}~{}y=\frac{y^{\rm e}}{\rm
FWC}\;,~{}~{}\sigma_{0}=\frac{1}{\sqrt{{\rm FWC}}}\;.$ (2)

The Gaussian approximation of a measurement is

$n^{\rm meas}\sim N(0,\sigma_{0}^{2}x)\;.$ (3)

The normalized approximation of Eq. (1) is

$y\approx x+n^{\rm meas}\;.$ (4)

In Eqs. (1,4), the noise variance is linear in $x$, and the signal-to-noise-
ratio (SNR) is the same. Hence, in this work, the parameter $\sigma_{0}$
controls the measurement noise level.

###  2.2 Imaging of Complex-valued Objects

Let $\mathbf{r}$ represent the 2D spatial coordinate of an object. Let
$a(\mathbf{r})$, ${\phi(\mathbf{r})}$ be its amplitude and phase,
respectively, and $j=\sqrt{-1}$. An optically-thin object has transmittance
${o}(\mathbf{r})={a(\mathbf{r})}\exp\left[{j{\phi(\mathbf{r})}}\right]$.
Cameras only measure intensity. Therefore, optical imaging systems can often
be formulated by an imaging model, which is nonlinear in $o$. Let $\rho$ be
the object irradiance. Let $H$ be a linear transformation, which depends on
the optical system. A noisy measurement $y$ is modeled by

$y={\cal N}\left(\rho|{H}\\{o(\mathbf{r})\\}|^{2}\right)\;.$ (5)

For example, in the phase retrieval problem [42], ${H}$ is equivalent to the
Fourier transform ${\mathcal{F}}$.

Fourier ptychography uses a set of $M$ linear transformations
${H}=\\{{H}_{m}\\}_{m=1}^{M}$. This method enables to expand the diffraction
limit set by the numerical aperture of the optical system, thus increasing the
optical resolution. Let ${\mathcal{F}}^{-1}$ be the inverse Fourier transform.
In the spatial frequency domain $\mathbf{k}$, the pupil function
${P}(\mathbf{k}-\mathbf{k}_{m})$ is a 2D bandpass filter, centered at spatial
frequency $\mathbf{k}_{m}$ with a spatial bandwidth set by the system
numerical aperture. Overall,

${H}_{m}\\{o(\mathbf{r})\\}={\mathcal{F}}^{-1}\\{{P}(\mathbf{k}-\mathbf{k}_{m}){\mathcal{F}}\\{o(\mathbf{r})\\}\\}.$
(6)

Discretize $o(\mathbf{r})$ to $q$ elements and the measurements
$y(\mathbf{r})$ to $q^{\rm meas}=Mq$ sampled pixels, creating vectors
$\mathbf{o}\in\mathbb{C}^{q}$ and $\mathbf{y}\in\mathbb{R}^{q^{\rm meas}}$.
The discrete version of the linear transformation $H$ is
$\mathbf{H}\in\mathbb{C}^{q^{\rm meas}\times q}$. The discrete versions of the
Fourier transform, its inverse, and the pupil function are denoted by
$\mathbf{F},\mathbf{F}^{-1}$ and $\mathbf{P}_{m}$, respectively. Then,
$\mathbf{H}_{m}\mathbf{o}=\sqrt{\rho}\mathbf{F}^{-1}\left[\mathbf{P}_{m}\odot\left[\mathbf{Fo}\right]\right]$,
where $\odot$ is an element-wise multiplication. Assume that the dominant
noise in Eq. (5) is photon noise. Following Sec. (2.1), during data analysis,
we approximate this noise as Gaussian $\mathbf{n}^{\rm meas}\sim
N(0,\sigma_{0}^{2}|\mathbf{H}\mathbf{o}|^{2})$. Then, Eq. (5) takes the form

$\mathbf{y}=|\mathbf{H}\mathbf{o}|^{2}+\mathbf{n^{\rm meas}}\;.$ (7)

###  2.3 Langevin Dynamics

Denoising is ill-posed: a noisy input may have multiple possible solutions. To
address this, _Langevin Dynamics_ can be used for denoising, converging
stochastically to a single sharp possible solution. Here we detail it shortly.

Let $t\in[1,\dots,T]$ be an iteration index and $\alpha_{t}>0$ be an
appropriately chosen small function of $t$, having $\alpha_{t\rightarrow
T}=0$. To maximize the probability distribution $p(\mathbf{x}|\mathbf{y})$,
the following iterations can be apply

$\mathbf{x}_{t+1}=\mathbf{x}_{t}+\alpha_{t}\nabla_{\mathbf{x}_{t}}\log
p(\mathbf{x}_{t}|\mathbf{y})\;.$ (8)

Eq. (8) may converge to an undesired local maximum. To address that, Langevin
Dynamics [49] adds to Eq. (8) _synthetic_ noise $\mathbf{n}_{t}^{\rm
Langevin}\sim N(0,\mathbf{I})$ at iteration $t$,

$\mathbf{x}_{t+1}=\mathbf{x}_{t}+\alpha_{t}\nabla_{\mathbf{x}_{t}}\log
p(\mathbf{x}_{t}|\mathbf{y})+\sqrt{2\alpha_{t}}\mathbf{n}^{\rm
Langevin}_{t}\;.$ (9)

Langevin Dynamics is a type of gradient descent. It relies on stochastic steps
to converge to a probable solution.

In general, when solving an inverse problem, neither
$p(\mathbf{x}_{t}|\mathbf{y})$ nor $\nabla_{\mathbf{x}_{t}}\log
p(\mathbf{x}_{t}|\mathbf{y})$ are known. Still, to use Eq. (9), there is a
need to assess $\nabla_{\mathbf{x}_{t}}\log p(\mathbf{x}_{t}|\mathbf{y})$. To
address that, Ref. [43] extends Eq. (9) into _annealed_ Langevin Dynamics.
Hence, iterations would _not_ perform Eq. (9), and would not attempt to update
a vector $\mathbf{x}_{t}$ whose probability cannot be computed. Instead,
iterations would update a vector, that we denote $\mathbf{\tilde{x}}_{t}$, for
which we can analytically assess the gradient
$\nabla_{\mathbf{\tilde{x}}_{t}}\log{p(\mathbf{\tilde{x}}_{t}|\mathbf{{y}})}$.
We thus need to properly define $\mathbf{\tilde{x}}_{t}$. Analogously to Eq.
(9), annealed Langevin Dynamics has this iterative rule

$\displaystyle\mathbf{\tilde{x}}_{t+1}=\mathbf{\tilde{x}}_{t}+\alpha_{t}\nabla_{\mathbf{\tilde{x}}_{t}}\log{p(\mathbf{\tilde{x}}_{t}|\mathbf{y})}+\sqrt{2\alpha_{t}}\mathbf{n}^{\rm
Langevin}_{t}\;.$ (10)

Any $\mathbf{\tilde{x}}_{t}$ differs from $\mathbf{x}$ by an error
$\mathbf{e}_{t}^{\rm annea}$,

$\mathbf{\tilde{x}}_{t}=\mathbf{{x}}+\mathbf{e}_{t}^{\rm annea}\;.$ (11)

We stress that this induced error $\mathbf{e}_{t}^{\rm annea}$ is only a
conceptual mathematical maneuver, and is not actually added during the
iterative process. It is important that the probability distribution of
$\mathbf{e}_{t}^{\rm annea}$ would be known, despite not knowing the actual
value of $\mathbf{e}_{t}^{\rm annea}$. Moreover, the error _is designed_ such
that $\mathbf{e}_{t\rightarrow T}^{\rm annea}\approx 0$.

Let us proceed with the description of annealed Langevin Dynamics Eq. 10.
Using Bayes rule,

$\nabla_{\mathbf{\tilde{x}}_{t}}\log{p(\mathbf{\tilde{x}}_{t}|\mathbf{y})}=\nabla_{\mathbf{\tilde{x}}_{t}}\left[\log{p(\mathbf{y}|\mathbf{\tilde{x}}_{t})}+\log{p(\mathbf{\tilde{x}}_{t})}\right]\;.$
(12)

There is a need to derive tractable expressions to the right-hand side of Eq.
(12). The term
$\nabla_{\mathbf{\tilde{x}}_{t}}\log{p(\mathbf{y}|\mathbf{\tilde{x}}_{t})}$ in
Eq. (12) can sometimes be derived analytically using the known distribution of
$\mathbf{e}_{t}^{\rm annea}$, as we show in Sections 3.1 and 3.2. To obtain
$\nabla_{\mathbf{\tilde{x}}_{t}}\log p(\mathbf{\tilde{x}}_{t})$, Refs. [25,
26] follow [44]: Before imaging takes place, they pre-train a DNN denoiser to
approximate $\nabla_{\mathbf{\tilde{x}}_{t}}\log p(\mathbf{\tilde{x}}_{t})$.

##  3 Langevin Dynamics for Optical Imaging

As described in Sections 2.1 and 2.2, real-world measurements have: (a)
Poisson distributed noise; (b) Nonlinear transforms over complex-valued
objects. There is a need to extend annealed Langevin Dynamics for both (a) and
(b).

As discussed in Section 2.3, in general inverse problems
$\nabla_{\mathbf{x}}\log{p(\mathbf{x}|\mathbf{y})}$ is unknown. Therefore, in
annealed Langevin Dynamics,
$\nabla_{\mathbf{\tilde{x}}_{t}}\log{p(\mathbf{\tilde{x}}_{t}|\mathbf{y})}$ is
exploited, while using the connection

$\lim\limits_{t\to
T}\nabla_{\mathbf{\tilde{x}}_{t}}\log{p(\mathbf{\tilde{x}}_{t}|\mathbf{y})}\rightarrow\nabla_{\mathbf{x}}\log{p(\mathbf{x}|\mathbf{y})}.$
(13)

In a similar fashion, to analytically handle signal dependency of noise, we
model the noise as if the signal dependency is proportional to
$\mathbf{\tilde{x}}_{t}$ rather then $\mathbf{x}$. This allows us to
analytically develop
$\nabla_{\mathbf{\tilde{x}}_{t}}\log{p(\mathbf{y}|\mathbf{\tilde{x}}_{t})}$ as
described later on in Sections 3.1 and 3.3. From Eq. 11, notice that our
approximation becomes accurate as
${\tilde{\mathbf{x}}}_{t}\rightarrow{\mathbf{x}}$. The iterative algorithm is
described in Algorithm 1.

Initialize $\mathbf{\tilde{o}}_{0}$, set a step size $\epsilon>0$

for  _$t\leftarrow 1$ to T_ do

$\alpha_{t}\leftarrow\epsilon({\sigma^{2}_{t}}/{\sigma^{2}_{T}})$

Draw $\mathbf{n}^{\rm Langevin}_{t}$ from a Gaussian distribution

if  _$\mathbf{\tilde{o}}$ is complex-valued_ then

From Eq. 39 compute gradient

$\mathbf{\Delta}_{t}\leftarrow\\!\nabla_{\mathbf{{\tilde{o}}}_{t-1}}\log
p(\mathbf{y}|\mathbf{{\tilde{o}}}_{t-1})+h_{\hat{\Theta}}^{\rm
complex}(\mathbf{\tilde{o}}_{t-1}|\sigma_{t})$

else

From Eq. 26 compute gradient

$\mathbf{\Delta}_{t}\leftarrow\\!\nabla_{\mathbf{{\tilde{o}}}_{t-1}}\log
p(\mathbf{y}|\mathbf{{\tilde{o}}}_{t-1})+h_{\hat{\Theta}}^{\rm
Poisson}(\mathbf{\tilde{o}}_{t-1}|\sigma_{t})$

end if

$\mathbf{\tilde{o}}_{t}\leftarrow\mathbf{\tilde{o}}_{t-1}+\alpha_{t}\mathbf{\Delta}_{t}+\sqrt{2\alpha_{t}}\mathbf{n}^{\rm
Langevin}_{t}$

end for

$\hat{\mathbf{o}}\leftarrow\mathbf{\tilde{o}}_{T}$

Algorithm 1 Langevin Dynamics for complex-valued objects. For real-valued
images, $\mathbf{\tilde{o}}$ is replaced by $\mathbf{\tilde{x}}$.

###  3.1 Denoising Poissonian Image Intensities

We now attempt to recover an image $\mathbf{x}$, given its noisy version
$\mathbf{y}$. As discussed in Section 2.1, $\mathbf{y}$ contains several types
of noise. Indeed data is _generated_ using a rich noise model: true
Poissonian, compounded by graylevel quantization. For _data analysis_ ,
however, we act as if $\mathbf{y}$ is related by Eq. 4. The parameter
$\sigma_{0}$ is known from the camera properties by Eq. 2, as described in
Section 2.1. Here, all vector operations are _element-wise_. Analogously to
[25, 26], we set a sequence of noise levels
$\sigma_{0}>\sigma_{1}\geq\sigma_{2}\geq\dots\geq\sigma_{T}>\sigma_{T+1}=0$.
Let $\mathbf{e}^{\rm part}_{t}$ be a statistically independent synthetic
_partial_ dummy error term, distributed as

$\mathbf{e}^{\rm part}_{t}\sim
N\left[0,(\sigma^{2}_{t}-\sigma^{2}_{t+1})\mathbf{I}\right]\;.$ (14)

A sum of zero-mean Gaussian variables is a Gaussian, whose variance is the sum
of the variances. So,

$\sum^{t_{2}}_{\tau=t_{1}}\mathbf{e}^{\rm part}_{\tau}\sim\
N\left[0,(\sigma^{2}_{t_{1}}-\sigma^{2}_{t_{2}})\mathbf{I}\right]\;.$ (15)

Using Eqs. 3, 15, LABEL: and 3 , we design $\mathbf{e}_{t}^{\rm part}$ such
that

$\mathbf{{{n}}}^{\rm meas}=\sqrt{\mathbf{x}}\sum^{T}_{\tau=0}\mathbf{e}^{\rm
part}_{\tau}\;.$ (16)

Moreover, recall that $\mathbf{e}_{t}^{\rm annea}$ is a dummy error term. We
set $\mathbf{e}_{t}^{\rm annea}$ in this task as

$\mathbf{e}_{t}^{\rm annea}=\sum^{T}_{\tau=t}\mathbf{e}^{\rm part}_{\tau}\;.$
(17)

We do _not_ set $\mathbf{e}_{t}^{\rm annea}$ as Refs. [25, 26] because the
noise in Refs. [25, 26] has stationary Gaussian distribution and here it is
Poissonian (signal-dependent).

We generalize Eq. 11 to Poisson noise, by defining an annealed variable

$\mathbf{\tilde{x}}_{t}=\mathbf{x}+\sqrt{\mathbf{x}}\mathbf{e}^{\rm
annea}_{t}=\mathbf{x}+\sqrt{\mathbf{x}}\sum^{T}_{\tau=t}\mathbf{e}^{\rm
part}_{\tau}\;.$ (18)

Eq. 18 allows us to derive $\nabla_{\mathbf{\tilde{x}}_{t}}\log
p(\mathbf{y}|\mathbf{\tilde{x}}_{t})$ analytically for Poisson denoising.
Inserting Eqs. 18, LABEL: and 16 in Eq. 4 yields

$\mathbf{y}=\mathbf{\tilde{x}}_{t}-\sqrt{\mathbf{x}}\sum^{T}_{\tau=t}\mathbf{e}^{\rm
part}_{\tau}+\sqrt{\mathbf{x}}\sum^{T}_{\tau=0}\mathbf{e}^{\rm
part}_{\tau}=\mathbf{\tilde{x}}_{t}+\sqrt{\mathbf{{x}}}\sum^{t-1}_{\tau=0}\mathbf{e}^{\rm
part}_{\tau}.$ (19)

Notice that $\mathbf{\tilde{x}}_{t}$ is statistically independent of the
additional dummy error $\sum^{t-1}_{\tau=0}\mathbf{e}^{\rm part}_{\tau}$. We
apply the following steps. From Eq. 19,

$\displaystyle{p(\mathbf{y}|\mathbf{\tilde{x}}_{t})}={p(\mathbf{y}-\mathbf{\tilde{x}}_{t}|\mathbf{\tilde{x}}_{t})}={p\left(\left[\sqrt{\mathbf{{x}}}\sum^{t-1}_{\tau=0}\mathbf{e}^{\rm
part}_{\tau}\right]|\mathbf{\tilde{x}}_{t}\right)}.$ (20)

Recall that ${\mathbf{x}}$ in Eq. 20 is unknown. Therefore, to set an explicit
analytical expression of Eq. 20, we approximate the remaining measurement
noise as if the signal dependency is proportional to $\mathbf{\tilde{x}}_{t}$
rather222From Eq. 18, this approximation becomes accurate as
${\tilde{\mathbf{x}}}_{t}\rightarrow{\mathbf{x}}$. then $\mathbf{x}$. Overall,
Eq. 20 takes the form of

$p(\mathbf{y}|\mathbf{{\tilde{x}}}_{t})\approx{p\left(\left[\sqrt{{\mathbf{{\tilde{x}}}}}\sum^{t-1}_{\tau=0}\mathbf{e}^{\rm
part}_{\tau}\right]|\mathbf{\tilde{x}}_{t}\right)}.$ (21)

Finally, the expression
$\sqrt{\mathbf{{\tilde{x}}}_{t}}\sum^{t-1}_{\tau=0}\mathbf{e}^{\rm
part}_{\tau}$ is both known and statistically tractable. This allows us to
find an explicit expression of the probability distribution
$p(\mathbf{y}|\mathbf{\tilde{x}}_{t})$

$\\!\\!\\!\\!{{{p}}(\mathbf{y}|\mathbf{\tilde{x}}_{t})}\\!\approx\\!\frac{1}{\sqrt{2\pi(\sigma{{}^{2}_{0}}-\sigma{{}^{2}_{t}})\mathbf{\tilde{x}}}_{t}}\exp\left[{\\!-\frac{({\mathbf{y}-\mathbf{\tilde{x}}_{t}})^{2}}{2(\sigma{{}^{2}_{0}}-\sigma{{}^{2}_{t}}){\mathbf{\tilde{x}}}_{t}}}\right].$
(22)

Computing the gradient of the logarithm of Eq. 22 yields

$\displaystyle\nabla_{\mathbf{\tilde{x}}_{t}}\log{p}(\mathbf{y}|\mathbf{\tilde{x}}_{t})=\frac{1}{2(\sigma{{}^{2}_{0}}-\sigma{{}^{2}_{t}})}\left(\frac{\mathbf{y}^{2}}{\mathbf{\tilde{x}}^{2}_{t}}-1\right)-\frac{1}{2\mathbf{\tilde{x}}_{t}}\;.$
(23)

In analogy to Section 2.3, we compute
$\nabla_{\mathbf{\tilde{x}}_{t}}\log{p(\mathbf{\tilde{x}}_{t})}$. The
derivation is detailed in the supplementary material. Here we bring the
result: The gradient $\nabla_{\mathbf{{\tilde{x}}}_{t}}\log
p(\mathbf{{\tilde{x}}}_{t})$ is given by,

$\nabla_{\mathbf{{\tilde{x}}}_{t}}\log
p(\mathbf{{\tilde{x}}}_{t})=\operatorname{\mathbb{E}}\left[{\frac{\mathbf{{x}}-\mathbf{{\tilde{x}}}_{t}}{{\sigma^{2}}\mathbf{{x}}}}|\mathbf{\tilde{x}}_{t}\right].$
(24)

The expectation in Eq. 24 cannot be calculated analytically. The reason is
that the probability distribution of $\mathbf{{x}}$ is not accessible. So,
before imaging takes place, we pre-train a DNN $h_{\Theta}^{\rm Poisson}$
having learnable parameters $\Theta$ to estimate
$\nabla_{\mathbf{{\tilde{x}}}_{t}}\log p(\mathbf{{\tilde{x}}}_{t})$, using
synthetic pairs of corresponding clean and noisy images
$(\mathbf{x},\mathbf{x}^{\prime})$,

$\hat{\Theta}=\arg\min_{\Theta}\operatorname{\mathbb{E}}\left[\left\lVert\frac{\mathbf{{x}}-\mathbf{{x}}^{\prime}}{{\sigma^{2}}\mathbf{{x}}}-h_{\Theta}^{\rm
Poisson}(\mathbf{{x}^{\prime}}|\sigma)\right\rVert^{2}_{2}\right]\;.$ (25)

The expectation $\operatorname{\mathbb{E}}$ is approximated by an empirical
average of all pairs in the training set of corresponding clean and noisy
images $(\mathbf{x},\mathbf{x}^{\prime})$. Overall, Eqs. (23,24,25) yeild

$\displaystyle\nabla_{\mathbf{\tilde{x}}_{t}}\log
p(\mathbf{\tilde{x}}_{t}|\mathbf{y})=\nabla_{\mathbf{\tilde{x}}_{t}}\log
p(\mathbf{\mathbf{y}|\tilde{x}}_{t})+h_{\hat{\Theta}}^{\rm
Poisson}(\mathbf{\tilde{x}}_{t}|\sigma_{t})=$
$\displaystyle\frac{1}{2(\sigma{{}^{2}_{0}}-\sigma{{}^{2}_{t}})}\left(\frac{\mathbf{y}^{2}}{\mathbf{\tilde{x}}^{2}_{t}}-1\right)-\frac{1}{2\mathbf{\tilde{x}}_{t}}+h_{\hat{\Theta}}^{\rm
Poisson}(\mathbf{\tilde{x}}_{t}|\sigma_{t})\;.$ (26)

###  3.2 Recovering Complex-Valued Objects

This section formulates annealed Langevin Dynamics for complex-valued,
nonlinear imaging models, as described in Section 2.2. Let ${\Re}(\cdot)$,
$\Im(\cdot)$ be real and imaginary arguments of a complex-valued number,
respectively. Analogously to Eq. 9, for Langevin Dynamics,
$\nabla_{\mathbf{{o}}_{t}}\log{p(\mathbf{{o}}_{t}|\mathbf{y})}$ is needed.
However, $p(\mathbf{{o}}_{t}|\mathbf{y})$ is not known, thus analogously to
Eq. 10, annealed Langevin Dynamics for a complex-valued object is

$\displaystyle\mathbf{\tilde{o}}_{t+1}=\mathbf{\tilde{o}}_{t}+\alpha_{t}\nabla_{\mathbf{\tilde{o}}_{t}}\log{p(\mathbf{\tilde{o}}_{t}|\mathbf{y})}+\sqrt{2\alpha_{t}}\mathbf{n}^{\rm
Langevin}_{t}\;.$ (27)

Here, $\mathbf{n}^{\rm Langevin}$ is a random complex-valued vector, where
$\Re(\mathbf{n}^{\rm Langevin}),\Im(\mathbf{n}^{\rm Langevin})\sim
N(0,\mathbf{I})$ and independent. For simplicity, we present here the special
case $\mathbf{H=I}$ in Eq. 7. Thus, most derivations can be given in scalar
formulation, that is, ${{y}}=|{{o}}|^{2}+{{{n}^{\rm meas}}}$. The
generalization of the model to a general $\mathbf{H}$ is detailed in the
supplementary material. Analogously to Eq. 11, we set

$\tilde{{o}}_{t}={{o}}+{e}_{t}^{\rm annea}\;.$ (28)

Using Eq. 28 in Eq. 7 yields

$\displaystyle{{y}}=|\tilde{o}_{t}|^{2}-{e}_{t}^{\rm couple}+|{e}^{\rm
annea}_{t}|^{2}+{{{n}^{\rm meas}}}\;,$ (29)

where

$\displaystyle|{o}_{t}|^{2}=\Re({o}_{t})^{2}+\Im({o}_{t})^{2}\;,$ (30)
${e}_{t}^{\rm couple}=2\left[\Re({\tilde{o}}_{t})\Re({e}_{t}^{\rm
annea})+\Im({\tilde{o}}_{t})\Im({e}_{t}^{\rm annea})\right]\;,$ (31)

and $|{e}_{t}^{\rm annea}|^{2}=\Re({e}_{t}^{\rm annea})^{2}+\Im({e}_{t}^{\rm
annea})^{2}$. The non-linearity of Eq. 7 yields two terms in Eq. 29 that pose
analytic challenges: a coupling term of ${e}_{t}^{\rm couple}$; and a squared
annealing error term $|{e}^{\rm annea}_{t}|^{2}$. We handle them by
introducing a _relaxed_ model, creating annealed Langevin Dynamics for
complex-valued objects, as

${{y}}=|\tilde{o}_{t}|^{2}-{e}_{t}^{\rm couple}+|{e}^{\rm
annea}_{t}|^{2}+{\tilde{n}}_{t}^{\rm meas}-\sigma^{2}_{t}{z}_{t}\;.$ (32)

We detail Eq. 32 and specifically ${{e}}_{t}^{\rm annea},{\tilde{n}}_{t}^{\rm
meas}$ and ${z}_{t}$ in Section 3.3.

| $\sigma_{0}=0.025$ | $\sigma_{0}=0.05$ | $\sigma_{0}=0.1$ | $\sigma_{0}=0.2$
---|---|---|---|---
CelebA 64x64 |  PSNR$\uparrow$ |  SSIM% $\uparrow$ |  FID$\downarrow$ |  PSNR$\uparrow$ |  SSIM% $\uparrow$ |  FID$\downarrow$ |  PSNR$\uparrow$ |  SSIM% $\uparrow$ |  FID$\downarrow$ |  PSNR$\uparrow$ |  SSIM% $\uparrow$ |  FID$\downarrow$
I + VST + BM3D [3] |  ${39.1\\!\pm\\!1.0}$ |  ${99\\!\pm\\!0.1}$ |  2.80 |  ${35.8\\!\pm\\!1.1}$ |  ${98\\!\pm\\!0.2}$ |  8.90 |  ${30.9\\!\pm\\!1.2}$ |  ${96\\!\pm\\!1.0}$ |  25.4 |  ${27.2\\!\pm\\!1.1}$ |  ${93\\!\pm\\!2.0}$ |  63.2
PnP RED [38] |  ${39.8\\!\pm\\!0.9}$ |  ${99\\!\pm\\!0.1}$ |  4.00 |  ${36.9\\!\pm\\!1.0}$ |  ${98\\!\pm\\!0.3}$ |  7.50 |  ${31.8\\!\pm\\!0.9}$ |  ${96\\!\pm\\!1.0}$ |  12.4 |  ${27.8\\!\pm\\!1.2}$ |  ${93\\!\pm\\!2.0}$ |  28.4
Bahjat et al. [26] |  $34.2\\!\pm\\!1.1$ |  $97\\!\pm\\!0.9$ |  2.10 |  $34.1\\!\pm\\!0.7$ |  $97\\!\pm\\!1.0$ |  6.70 |  $30.6\\!\pm\\!0.9$ |  $95\\!\pm\\!2.0$ |  11.1 |  $27.1\\!\pm\\!1.2$ |  $91\\!\pm\\!4.0$ |  15.2
MMSE |  ${42.5\\!\pm\\!1.0}$ |  ${99\\!\pm\\!0.1}$ |  1.32 |  ${38.4\\!\pm\\!1.1}$ |  ${99\\!\pm\\!0.2}$ |  3.60 |  ${34.5\\!\pm\\!1.1}$ |  ${98\\!\pm\\!0.5}$ |  7.50 |  ${30.9\\!\pm\\!1.2}$ |  ${96\\!\pm\\!1.0}$ |  11.9
Ours |  ${41.5\\!\pm\\!1.1}$ |  ${99\\!\pm\\!0.0}$ |  1.25 |  ${37.1\\!\pm\\!1.0}$ |  ${99\\!\pm\\!0.4}$ |  2.70 |  ${32.8\\!\pm\\!1.1}$ |  ${98\\!\pm\\!0.5}$ |  6.10 |  ${28.6\\!\pm\\!1.3}$ |  ${94\\!\pm\\!1.0}$ |  13.0
FFHQ 256x256 |  |  |  |  |  |  |  |  |  |  |  |
DPS [12] |  ${36.8\\!\pm\\!4.7}$ |  ${99\\!\pm\\!0.1}$ |  2.00 |  ${33.8\\!\pm\\!3.8}$ |  ${98\\!\pm\\!0.4}$ |  5.47 |  ${31.5\\!\pm\\!2.7}$ |  ${98\\!\pm\\!2.0}$ |  12.3 |  ${28.6\\!\pm\\!2.0}$ |  ${93\\!\pm\\!2.0}$ |  23.8
Ours |  ${41.9\\!\pm\\!1.0}$ |  ${99\\!\pm\\!0.0}$ |  0.72 |  ${37.9\\!\pm\\!1.3}$ |  ${99\\!\pm\\!0.1}$ |  2.50 |  ${33.5\\!\pm\\!1.1}$ |  ${96\\!\pm\\!1.0}$ |  10.0 |  ${29.0\\!\pm\\!1.0}$ |  ${93\\!\pm\\!1.0}$ |  33.6

Table 1: Results on 1000 images of CelebA and FFHQ datasets. For
$\sigma_{0}\leq 0.1$, our algorithm has better FID and similar SSIM with the
MMSE denoiser and DPS [12]. The MMSE denoiser inherently achieves superior
PSNR.

###  3.3 Relaxation of the Nonlinear Problem

Let $\\{{\sigma}_{t}\\}^{T}_{t=0}$ be a sequence of standard deviations, set
to
$\sigma_{0}>\sigma_{1}\geq\sigma_{2}\geq\dots\geq\sigma_{T}>\sigma_{T+1}=0$.
Let

${\Re}({e}^{\rm part}_{t}),\Im({e}^{\rm part}_{t})\sim
N[0,\frac{1}{4}(\sigma^{2}_{t}-\sigma^{2}_{t+1})]$ (33)

be statistically independent Gaussian variables. As described in Section 2.2
and in Eq. 3, for $\mathbf{H=I}$, ${n}^{\rm meas}\sim
N(0,\sigma_{0}^{2}|{o}|^{2})$. Using $|o|^{2}$ from Eq. 30, ${n}^{\rm
meas}\sim N[0,\sigma_{0}^{2}\Re({o})^{2}+\sigma_{0}^{2}\Im({o})^{2}]$. In
analogy to Eq. 16, we use Eq. 33 and design dummy variables ${\Re}({e}^{\rm
part}_{t}),\Im({e}^{\rm part}_{t})$ such that,

${{n}^{\rm meas}}=2\Re({{o}})\sum^{T}_{\tau=0}{\Re}({e}^{\rm
part}_{\tau})+2\Im({{o}})\sum^{T}_{\tau=0}\Im({e}^{\rm part}_{\tau})\;.$ (34)

Notice that Eq. 34 now has a form similar to the coupling term ${e}_{t}^{\rm
couple}$ in Eq. 31. To make computations tractable, we assume
${o}\approx{\tilde{o}}_{t}$ in Eq. 34. We do not use ${\tilde{o}}_{t}$ as a
denoised version yet. An approximate333From Eq. 28, the approximation becomes
accurate as ${e}_{t}^{\rm annea}\rightarrow 0$, where
${\tilde{o}}_{t}\rightarrow{o}$ and ${{\tilde{n}^{\rm
meas}}}_{t}\rightarrow{n^{\rm meas}}$. measurement noise is

${\tilde{n}^{\rm
meas}}_{t}=2\Re(\tilde{{{o}}}_{t})\sum^{T}_{\tau=0}{\Re}({e}^{\rm
part}_{\tau})+2\Im(\tilde{{{o}}}_{t})\sum^{T}_{i=0}{\Im}({e}^{\rm
part}_{\tau})\;.$ (35)

Using Eqs. 31 and 35, define

$\displaystyle{e}_{t}^{\rm remain}(\tilde{{{o}}}_{t})={\tilde{n}^{\rm
meas}}_{t}-{e}_{t}^{\rm
couple}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle=2\Re(\tilde{{{o}}}_{t})\sum^{t-1}_{\tau=0}\Re{({e}^{\rm
part}_{\tau})}+2\Im(\tilde{{{o}}}_{t})\sum^{t-1}_{\tau=0}\Im{({e}^{\rm
part}_{\tau})}\;.$ (36)

The term ${z}_{t}$ in Eq. 32 is designed to address the challenge of a squared
annealing error term $|{e}^{\rm annea}_{t}|^{2}$. We detail ${z}_{t}$ in the
supplementary material, here we give the final result,

$\displaystyle{z}_{t}^{\rm remain}=|{e}^{\rm
annea}_{t}|^{2}-\sigma_{t}^{2}{z}_{t}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle={\left[\sum^{t-1}_{\tau=0}\Re{({e}^{\rm
part}_{\tau})}\right]^{2}}\\!\\!\\!+{\left[\sum^{t-1}_{\tau=0}\Im{({e}^{\rm
part}_{\tau})}\right]^{2}}\;.$ (37)

Insert Eqs. 36 and 37 into the relaxed model in Eq. 32. This yields

${{y}}=|{\tilde{o}}_{t}|^{2}+{e}_{t}^{\rm
remain}(\tilde{{{o}}}_{t})+{z}_{t}^{\rm remain}\;.$ (38)

Using Eq. 38, $\nabla_{{\tilde{o}}_{t}}\log p({y}|{\tilde{o}}_{t})$ can be
derived analytically. Here we give the results: the derivation is detailed in
the supplementary material. Let $I_{0},I_{1}$ be the modified Bessel functions
of the first kind [5]. Then,

$\nabla_{{\tilde{o}}_{t}}\log
p({y}|{\tilde{o}}_{t})={\frac{{\tilde{o}}_{t}}{2(\sigma^{2}_{0}-\sigma^{2}_{t})}\left[{\cal
I}({\tilde{o}}_{t},{y})\frac{\sqrt{{y}}}{|{\tilde{o}}_{t}|}-1\right]}\;,$ (39)

where ${\cal
I}({\tilde{o}}_{t},{y})={I_{1}\left({\frac{|{\tilde{o}}_{t}|\sqrt{{y}}}{\sigma^{2}_{0}-\sigma^{2}_{t}}}\right)}/{I_{0}{\left({\frac{|{\tilde{o}}_{t}|\sqrt{{y}}}{\sigma^{2}_{0}-\sigma^{2}_{t}}}\right)}}$.
The gradient $\nabla_{\mathbf{{\tilde{o}}}_{t}}\log
p(\mathbf{{\tilde{o}}}_{t})$ is given by

$\nabla_{\mathbf{{\tilde{o}}}_{t}}\log
p(\mathbf{{\tilde{o}}}_{t})=\operatorname{\mathbb{E}}\left[\frac{\mathbf{o}-\mathbf{{\tilde{o}}}_{t}}{\sigma^{2}}|\mathbf{\tilde{o}}_{t}\right]\;.$
(40)

See the supplemental material for the proof of Eq. 40.

Hence again, we pre-train444The network input is a two-channel tensor of
$\Re(\mathbf{{o}^{\prime}}),\Im(\mathbf{{o}^{\prime}})$. a DNN
$h_{\Theta}^{\rm complex}$ having learnable parameters $\Theta$ to estimate
$\nabla_{\mathbf{{\tilde{o}}}_{t}}\log p(\mathbf{{\tilde{o}}}_{t})$

$\hat{\Theta}=\arg\min_{\Theta}\operatorname{\mathbb{E}}\left[\left\lVert\frac{\mathbf{{o}}-\mathbf{{o}}^{\prime}}{{\sigma^{2}}}-h_{\Theta}^{\rm
complex}(\mathbf{{o}}^{\prime}|\sigma)\right\rVert^{2}_{2}\right]\;.$ (41)

Overall, the complex-valued derivative is

$\displaystyle\nabla_{\mathbf{\tilde{o}}_{t}}\log
p(\mathbf{\tilde{o}}_{t}|\mathbf{y})=\nabla_{\mathbf{\tilde{o}}_{t}}\log
p(\mathbf{\mathbf{y}|\tilde{o}}_{t})+h_{\hat{\Theta}}^{\rm
complex}(\mathbf{\tilde{o}}_{t}|\sigma_{t})=$
$\displaystyle{\frac{{\tilde{\mathbf{o}}}_{t}}{2(\sigma^{2}_{0}-\sigma^{2}_{t})}\odot\left[{\cal
I}({\tilde{\mathbf{o}}}_{t},{\mathbf{y}})\odot\frac{\sqrt{{\mathbf{y}}}}{|{\tilde{\mathbf{o}}}_{t}|}-1\right]}+h_{\hat{\Theta}}^{\rm
complex}(\mathbf{\tilde{o}}_{t}|\sigma_{t})\;.$ (42)

##  4 Experiments

To evaluate recoveries we use the FID, SSIM and PSNR555For complex-valued
objects, the amplitude is evaluated by standard PSNR metric. The recovered
object’s phase ${\angle}\mathbf{o}$ has an inherent ambiguity due to the
$|\cdot|^{2}$ operation and the periodicity of $2\pi$. Therefore, we evaluate
phase recovery by ${\rm PSNR}^{\rm
phase}({\angle}\hat{\mathbf{o}})=10\log_{10}\left[{(2\pi)^{2}}\left[{\rm
MSE}^{\rm phase}({\angle}\hat{\mathbf{o}})\right]^{-1}\right]$, where ${\rm
MSE}^{\rm
phase}({\angle}\hat{\mathbf{o}})=\arg\min_{\theta}\left\lVert{\angle}\mathbf{o}-({\angle}\hat{\mathbf{o}}+\theta)\right\rVert^{2}_{2}/q$.
criteria. In all experiments we used the NCSNv2 architecture [44] as the DNN.
The bottleneck of the proposed algorithm is iteratively running this DNN
$T\sim 1000$ times, where time complexity increases with image size and
linearly with $T$. We apply our algorithm on both real biological data [55,
47] and on synthetic datasets, CelebA ($64\times 64$) and FFHQ ($256\times
256$). We stress that during testing on synthetic scenes, we synthesized data
using actual Poisson noise (Eq. 1). Moreover, 8-bit quantization was applied
on the simulated electron count.

###  4.1 Denoising Poissonian Image Intensities

Figure 3: Poisson denoising results. From top to bottom: CelebA
$\sigma_{0}=0.1$, $\sigma_{0}=0.2$. FFHQ $\sigma_{0}=0.2$.

We compare our algorithm with [3, 12, 26, 38] and an MMSE solution. Fig. 3 and
Table 1 present Poisson denoising results. Across all tests, our method
outperforms [3], a non-learning-based state-of-the-art method for Poisson
denoising. Furthermore, we outperform Ref. [26], which utilizes annealed
Langevin Dynamics for _Gaussian_ noise. This indicates the benefit of
addressing the natural property of imaging noise. PnP RED [38] is a method to
solve general inverse problems. It uses Gaussian denoiser as its prior666The
data term derives from a Poissonian model.. Our method surpasses it,
emphasizing the significance of employing a customized denoiser for Poisson
denoising.

Typically, the SNR in standard optical systems satisfies $\sigma_{0}\leq 0.1$.
In that regime, we have better FID results and similar SSIM results than the
MMSE denoiser and DPS [12]. MMSE denoisers are designed to minimize the MSE
criterion, that is, maximizing the PSNR. Thus, unsurprisingly, an MMSE
denoiser yields a superior PSNR than other algorithms, however, with blur
outputs (see Fig. 3).

####  4.1.1 Poisson Denoising of Real Data

Fig. 1 shows our method on real data [55]. It contains noisy colored images of
fixed BPAE cells, obtained using two-photon microscopy. The imaged sample was
split into 20 non-overlapping patches, 19 for training, and the remaining for
testing. In [55], each patch was repeatedly captured 50 times, as 50 noise
realizations. The clean ground-truth image is based on temporal averaging of
the 50 images. The DNN trains only on 19 images with synthetic Poisson noise,
while testing on raw data involves real noise. The authors in Ref. [55]
modeled the imaging system noise by two additive components with zero means: A
signal-dependent Poisson noise, and a negligible777Let $\mathbf{x}\in[0,1]$ be
the true normalized image. Let $n^{G}$, $n^{P}$ be Gaussian noise and Gaussian
approximation for the Poisson noise, with corresponding standard deviations
$\sigma_{G}=5\times 10^{-4}$ and $\sigma_{P}=\sqrt{x}(2.5\times 10^{-2})$.
Then, the measured image $\mathbf{y}$ is approximated by
$\mathbf{y}=\mathbf{x}+n^{G}+n^{P}$.  signal-independent Gaussian noise.
Therefore, we apply our denoising algorithm as the measured image $\mathbf{y}$
can be described by Eq. 4. The PSNR and SSIM of the noisy raw data are 26 and
0.67, respectively, while ours are 33 and 0.92.

###  4.2 Phase Retrieval

As described in Section 2.2, phase retrieval can be modeled using
$\mathbf{H}=\mathbf{F}$, that is, $\mathbf{y}={\cal
N}\left(|\mathbf{F}\mathbf{x}|^{2}\right)$. Here $\cal N$ contains photon
noise and quantization, the latter using 8 bits per pixel. Note that even if
$\mathbf{x}$ is a real-valued image, the measurement $\mathbf{y}$ has no
phase. This task is ill-posed. The Hybrid Input-Output (HIO) [18] algorithm is
frequently used in phase retrieval. Let $\boldsymbol{\psi}$ be a random
Gaussian vector. HIO initializes hundreds of possible candidate solutions
$\\{\mathbf{x}_{l}\\}$ by
$\mathbf{F}\mathbf{x}_{l}=\sqrt{\mathbf{y}}\exp(j2\pi\boldsymbol{\psi}_{l})$.
Then, HIO iteratively uses negative feedback in Fourier space, in order to
progressively force the candidate solutions to conform to the problem
constraints. The candidate that holds best these constraints is the evaluated
HIO solution.

Figure 4: Phase retrieval results on CelebA images. A variance image is
obtained by using 50 stochastic solutions.  Figure 5: Phase retrieval
qualitative results for the dataset suggested in [32] $(\sigma_{0}=4/255)$.

$\sigma_{0}=$ | $2/255$ | $3/255$ | $4/255$
---|---|---|---
HIO [18] | 22.2 | 19.9 | 17.9
BM3D-prGAMP [33] | 24.5 | 23.2 | 21.0
BM3D-ADMM [32] | 27.5 | 24.5 | 21.8
DnCNN-ADMM [32] | 29.3 | 26.0 | 22.0
prDeep [32] | 28.4 | 28.5 | 26.4
Ours | 33.1 | 32.6 | 28.9

Table 2: Phase retrieval PSNR results for the “natural” images dataset in
[32], using different noise levels $\sigma_{0}$.

We compare our method (Algorithm 1) for phase retrieval, vs. prior art. Here,
we initialize our algorithm with an extremely noisy version of the HIO
solution. We found that the noisy initialization gives us better robustness
during iterations, but still stochastically converges to a variety of possible
solutions (see Fig. 4). Additionally, we follow the experiment suggested in
[32]. First, we trained $h_{\Theta}^{\rm
Complex}(\mathbf{{o}}^{\prime}|\sigma)$ with overlapping patches drawn from
400 images in the Berkeley Segmentation Dataset [31]. Then, we evaluate our
method on a test set of six “natural” images[32], three images are presented
in Fig. 5. Table 2 compares the results. In the supplementary material, we
present more results and further details.

###  4.3 Fourier Ptychography: Real Experiment

We further demonstrate our approach on the Fourier ptychography problem,
described in Section 2.2. We use data of real-world publicly available
complex-valued biological samples [46]. Using a Fourier ptychography imaging
system, the amplitude and phase of a sample of human bone cells are recovered
[47]. We randomly divided the recovered sample into non-spatially overlapping
patches, each of size $256\times 256$ pixels. Then, the patches were split
randomly into training and test sets. We train our denoiser over the training
set, according to Eq. 41. For each patch $\mathbf{o}$ in the test set, we
simulate $M$ noisy measured images, each corresponding to a single LED in the
optical system, according to Eqs. 5 and 6. Figs. 7 and 6 show the results
using Algorithm 1. Further details about this test and more results are in the
supplementary material.

Figure 6: Fourier ptychography. Average PSNR for 10 complex-valued objects:
[Top] Varying Poisson noise levels $\sigma_{0}$ and fixed $M=89$. [Bottom]
Varying number of LEDs $M$ and a fixed noise level $\sigma_{0}=0.1$.
[Baseline] Algorithm 1 excluding the prior term $h_{\hat{\Theta}}^{\rm
complex}$. Yeh et al. [52] uses quasi-Newton method. The bars stand for PSNR
standard deviations.  Figure 7: Fourier ptychography. Our results of a
complex-valued restoration, for $M=89$ LEDs, $\sigma_{0}=0.1$ and 8 bit
quantization. [Top] The amplitude and phase images’ dynamic range $\in[0,1]$
and $[-\pi,\pi]$ respectively. [Bottom] Scatter plots for the amplitude and
phase values, the dashed red lines represent a perfect recovery.

##  5 Discussion

We generalize annealed Langevin Dynamics to realistic optical imaging, which
involves complex-valued objects and Poisson noise. Our method is generic and
can be applied to a variety of imaging transformations $\mathbf{H}$.
Furthermore, the pre-trained DNN in Eq. 41 is independent of specific imaging
system parameters $\sigma_{0}$ and $\mathbf{H}$ of Eq. 7. Hence, the same DNN
can be used for various optical systems.

Dealing with non-Gaussian noises, such as quantization noise or dead-time
noise in SPAD-based imager, is essential in real-world imaging. Hence, noise
statistics create new and interesting challenges, specifically in low-light
conditions. There are more challenging forward models in imaging. These
include imaging by coherent illumination (with speckles), partial-coherence,
and models that are nonlinear with respect to the illumination (e.g. two-
photon microscopy). Other computer vision problems are highly nonlinear in the
unknowns, such as tomography scatter fields [39]. Therefore, motivated by this
paper, we hope that more research expands annealed Langevin Dynamics to
additional challenging models in real-world imaging.

## References

  * [1] Kazunori Akiyama, Antxon Alberdi, Walter Alef, Keiichi Asada, Rebecca Azulay, Anne-Kathrin Baczko, David Ball, Mislav Baloković, John Barrett, Dan Bintley, et al.  First M87 Event Horizon Telescope results. IV. imaging the central supermassive black hole.  The Astrophysical Journal Letters, 875(1):L4, 2019.
  * [2] Benjamin Attal and Matthew O’Toole.  Towards mixed-state coded diffraction imaging.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2022\.
  * [3] Lucio Azzari and Alessandro Foi.  Variance stabilization for noisy+ estimate combination in iterative poisson denoising.  IEEE signal processing letters, 23(8):1086–1090, 2016.
  * [4] Yochai Blau and Tomer Michaeli.  The perception-distortion tradeoff.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 6228–6237, 2018.
  * [5] Frank Bowman.  Introduction to Bessel Functions.  Courier Corporation, 2012.
  * [6] Antoni Buades, Bartomeu Coll, and J-M Morel.  A non-local algorithm for image denoising.  In 2005 IEEE computer society conference on computer vision and pattern recognition (CVPR’05), volume 2, pages 60–65. Ieee, 2005.
  * [7] Jaeseok Byun, Sungmin Cha, and Taesup Moon.  FBI-denoiser: Fast blind image denoiser for poisson-gaussian noise.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5768–5777, 2021.
  * [8] Eunju Cha, Chanseok Lee, Mooseok Jang, and Jong Chul Ye.  Deepphasecut: Deep relaxation in phase for unsupervised fourier phase retrieval.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(12):9931–9943, 2021.
  * [9] Dorian Chan, Srinivasa G Narasimhan, and Matthew O’Toole.  Holocurtains: Programming light curtains via binary holography.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 17886–17895, 2022.
  * [10] Yunjin Chen and Thomas Pock.  Trainable nonlinear reaction diffusion: A flexible framework for fast and effective image restoration.  IEEE transactions on pattern analysis and machine intelligence, 39(6):1256–1272, 2016.
  * [11] Jooyoung Choi, Sungwon Kim, Yonghyun Jeong, Youngjune Gwon, and Sungroh Yoon.  Ilvr: Conditioning method for denoising diffusion probabilistic models.  arXiv preprint arXiv:2108.02938, 2021.
  * [12] Hyungjin Chung, Jeongsol Kim, Michael T Mccann, Marc L Klasky, and Jong Chul Ye.  Diffusion posterior sampling for general noisy inverse problems.  arXiv preprint arXiv:2209.14687, 2022.
  * [13] Hyungjin Chung, Byeongsu Sim, Dohoon Ryu, and Jong Chul Ye.  Improving diffusion models for inverse problems using manifold constraints.  Advances in Neural Information Processing Systems, 35:25683–25696, 2022.
  * [14] Hyungjin Chung, Byeongsu Sim, and Jong Chul Ye.  Come-closer-diffuse-faster: Accelerating conditional diffusion models for inverse problems through stochastic contraction.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 12413–12422, 2022.
  * [15] Angélique Drémeau, Antoine Liutkus, David Martina, Ori Katz, Christophe Schülke, Florent Krzakala, Sylvain Gigan, and Laurent Daudet.  Reference-less measurement of the transmission matrix of a highly scattering material using a dmd and phase retrieval techniques.  Optics Express, 23(9):11898–11911, 2015.
  * [16] Regina Eckert, Zachary F Phillips, and Laura Waller.  Efficient illumination angle self-calibration in fourier ptychography.  Applied Optics, 57(19):5434–5442, 2018.
  * [17] Michael Elad and Michal Aharon.  Image denoising via sparse and redundant representations over learned dictionaries.  IEEE Transactions on Image processing, 15(12):3736–3745, 2006.
  * [18] James R Fienup.  Phase retrieval algorithms: A comparison.  Applied Optics, 21(15):2758–2769, 1982.
  * [19] Jan Funke, Fabian Tschopp, William Grisaitis, Arlo Sheridan, Chandan Singh, Stephan Saalfeld, and Srinivas C Turaga.  Large scale image segmentation with structured loss based deep learning for connectome reconstruction.  IEEE transactions on pattern analysis and machine intelligence, 41(7):1669–1680, 2018.
  * [20] Shuhang Gu, Lei Zhang, Wangmeng Zuo, and Xiangchu Feng.  Weighted nuclear norm minimization with application to image denoising.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2862–2869, 2014.
  * [21] Tao Huang, Songjiang Li, Xu Jia, Huchuan Lu, and Jianzhuang Liu.  Neighbor2neighbor: Self-supervised denoising from single noisy images.  In Proceedings of the IEEE/CVF conference on Computer Vision and Pattern Recognition (CVPR), pages 14781–14790, 2021.
  * [22] Rakib Hyder, Zikui Cai, and M Salman Asif.  Solving phase retrieval with a learned reference.  In European Conference on Computer Vision (ECCV), pages 425–441. Springer, 2020.
  * [23] Zahra Kadkhodaie and Eero Simoncelli.  Stochastic solutions for linear inverse problems using the prior implicit in a denoiser.  Advances in Neural Information Processing Systems (NIPS), 34:13242–13254, 2021.
  * [24] Bahjat Kawar, Michael Elad, Stefano Ermon, and Jiaming Song.  Denoising diffusion restoration models.  Advances in Neural Information Processing Systems, 35:23593–23606, 2022.
  * [25] Bahjat Kawar, Gregory Vaksman, and Michael Elad.  SNIPS: Solving noisy inverse problems stochastically.  Advances in Neural Information Processing Systems (NIPS), 34:21757–21769, 2021.
  * [26] Bahjat Kawar, Gregory Vaksman, and Michael Elad.  Stochastic image denoising by sampling from the posterior distribution.  In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), pages 1866–1875, 2021.
  * [27] Wesley Khademi, Sonia Rao, Clare Minnerath, Guy Hagen, and Jonathan Ventura.  Self-supervised poisson-gaussian denoising.  In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pages 2131–2139, 2021.
  * [28] Alexander Krull, Tim-Oliver Buchholz, and Florian Jug.  Noise2void-learning denoising from single noisy images.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 2129–2137, 2019.
  * [29] Prashanth G Kumar and Rajiv Ranjan Sahay.  Low rank poisson denoising (LRPD): A low rank approach using split bregman algorithm for poisson noise removal from images.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (WCVPR), pages 0–0, 2019.
  * [30] Ding Liu, Bihan Wen, Yuchen Fan, Chen Change Loy, and Thomas S Huang.  Non-local recurrent network for image restoration.  Advances in neural information processing systems, 31, 2018.
  * [31] David Martin, Charless Fowlkes, Doron Tal, and Jitendra Malik.  A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics.  In Proceedings Eighth IEEE International Conference on Computer Vision (ICCV), volume 2, pages 416–423. IEEE, 2001.
  * [32] Christopher Metzler, Phillip Schniter, Ashok Veeraraghavan, and Richard Baraniuk.  Prdeep: Robust phase retrieval with a flexible deep network.  In International Conference on Machine Learning, pages 3501–3510. PMLR, 2018.
  * [33] Christopher A Metzler, Arian Maleki, and Richard G Baraniuk.  BM3D-PRGAMP: Compressive phase retrieval based on BM3D denoising.  In 2016 IEEE International Conference on Image Processing (ICIP), pages 2504–2508. IEEE, 2016.
  * [34] Ben Mildenhall, Jonathan T Barron, Jiawen Chen, Dillon Sharlet, Ren Ng, and Robert Carroll.  Burst denoising with kernel prediction networks.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, pages 2502–2510, 2018.
  * [35] Ben Moseley, Valentin Bickel, Ignacio G López-Francos, and Loveneesh Rana.  Extreme low-light environment-driven image denoising over permanently shadowed lunar regions with a physical noise model.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6317–6327, 2021.
  * [36] Netanel Ratner and Yoav Y Schechner.  Illumination multiplexing within fundamental limits.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 1–8. IEEE, 2007.
  * [37] John Rodenburg and Andrew Maiden.  Ptychography.  Springer Handbook of Microscopy, pages 819–904, 2019.
  * [38] Yaniv Romano, Michael Elad, and Peyman Milanfar.  The little engine that could: Regularization by denoising (RED).  SIAM Journal on Imaging Sciences, 10(4):1804–1844, 2017.
  * [39] Roi Ronen, Vadim Holodovsky, and Yoav Y Schechner.  Variable imaging projection cloud scattering tomography.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2022\.
  * [40] Atreyee Saha, Salman S Khan, Sagar Sehrawat, Sanjana S Prabhu, Shanti Bhattacharya, and Kaushik Mitra.  Lwgnet-learned wirtinger gradients for fourier ptychographic phase retrieval.  In European Conference on Computer Vision (ECCV), pages 522–537. Springer, 2022.
  * [41] Fahad Shamshad, Asif Hanif, Farwa Abbas, Muhammad Awais, and Ali Ahmed.  Adaptive ptych: Leveraging image adaptive generative priors for subsampled fourier ptychography.  In Proceedings of the IEEE/CVF International Conference on Computer Vision Workshops (WICCV), pages 0–0, 2019.
  * [42] Yoav Shechtman, Yonina C Eldar, Oren Cohen, Henry Nicholas Chapman, Jianwei Miao, and Mordechai Segev.  Phase retrieval with application to optical imaging: A contemporary overview.  IEEE Signal Processing Magazine, 32(3):87–109, 2015.
  * [43] Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  Advances in Neural Information Processing Systems (NIPS), 32, 2019\.
  * [44] Yang Song and Stefano and Ermon.  Improved techniques for training score-based generative models.  Advances in Neural Information Processing Systems (NIPS, 33:12438–12448, 2020.
  * [45] Yang Song, Liyue Shen, Lei Xing, and Stefano Ermon.  Solving inverse problems in medical imaging with score-based generative models.  arXiv preprint arXiv:2111.08005, 2021.
  * [46] Lei Tian, Xiao Li, Kannan Ramchandran, and Laura Waller.  (Dataset) LED array 2D Fourier Ptychography. Available online.  http://gigapan.com/gigapans/170956.
  * [47] Lei Tian, Xiao Li, Kannan Ramchandran, and Laura Waller.  Multiplexed coded illumination for Fourier ptychography with an LED array microscope.  Biomedical Optics Express, 5(7):2376–2389, 2014.
  * [48] Zhou Wang and Alan C Bovik.  Mean squared error: Love it or leave it? a new look at signal fidelity measures.  IEEE Signal Processing Magazine, 26(1):98–117, 2009.
  * [49] Max Welling and Yee W Teh.  Bayesian learning via stochastic gradient langevin dynamics.  In Proceedings of the International Conference on Machine Learning (ICML), pages 681–688, 2011.
  * [50] Yicheng Wu, Manoj Kumar Sharma, and Ashok Veeraraghavan.  Wish: wavefront imaging sensor with high resolution.  Light: Science & Applications, 8(1):44, 2019.
  * [51] Duoduo Xue, Ziyang Zheng, Wenrui Dai, Chenglin Li, Junni Zou, and Hongkai Xiong.  On the convergence of non-convex phase retrieval with denoising priors.  IEEE Transactions on Signal Processing, 70:4424–4439, 2022.
  * [52] Li-Hao Yeh, Jonathan Dong, Jingshan Zhong, Lei Tian, Michael Chen, Gongguo Tang, Mahdi Soltanolkotabi, and Laura Waller.  Experimental robustness of fourier ptychography phase retrieval algorithms.  Optics express, 23(26):33214–33240, 2015.
  * [53] Feilong Zhang, Xianming Liu, Cheng Guo, Shiyi Lin, Junjun Jiang, and Xiangyang Ji.  Physics-based iterative projection complex neural network for phase retrieval in lensless microscopy imaging.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 10523–10531, 2021.
  * [54] Kai Zhang, Wangmeng Zuo, Yunjin Chen, Deyu Meng, and Lei Zhang.  Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising.  IEEE transactions on image processing, 26(7):3142–3155, 2017.
  * [55] Yide Zhang, Yinhao Zhu, Evan Nichols, Qingfei Wang, Siyuan Zhang, Cody Smith, and Scott Howard.  A Poisson-Gaussian denoising dataset with real fluorescence microscopy images.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 11710–11718, 2019.
  * [56] Guoan Zheng, Roarke Horstmeyer, and Changhuei Yang.  Wide-field, high-resolution fourier ptychographic microscopy.  Nature Photonics, 7(9):739–745, 2013.
  * [57] Daniel Zoran and Yair Weiss.  From learning models of natural image patches to whole image restoration.  In 2011 international conference on computer vision, pages 479–486. IEEE, 2011.