# Photon-Limited Blind Deconvolution using Unsupervised Iterative Kernel
Estimation

Yash Sanghvi, , Abhiram Gnanasambandam, , Zhiyuan Mao, , and Stanley H. Chan
Y. Sanghvi, Z. Mao and S. Chan are with the School of Electrical and Computer
Engineering, Purdue University, West Lafayette, IN 47907, USA. The work of A.
Gnanasambandam was completed when he was a graduate student at Purdue
University. Email: {ysanghvi, mao114<EMAIL_ADDRESS><EMAIL_ADDRESS>The work is supported, in part, by the US National
Science Foundation under the grants IIS-2133032, and ECCS-2030570.

###### Abstract

Blind deconvolution is a challenging problem, but in low-light it is even more
difficult. Existing algorithms, both classical and deep-learning based, are
not designed for this condition. When the photon shot noise is strong,
conventional deconvolution methods fail because (1) the image does not have
enough signal-to-noise ratio to perform the blur estimation; (2) While deep
neural networks are powerful, many of them do not consider the forward
process. When the noise is strong, these networks fail to simultaneously
deblur and denoise; (3) While iterative schemes are known to be robust in the
classical frameworks, they are seldom considered in deep neural networks
because it requires a differentiable non-blind solver.

This paper addresses the above challenges by presenting an _unsupervised_
blind deconvolution method. At the core of this method is a reformulation of
the general blind deconvolution framework from the conventional image-kernel
alternating minimization to a purely kernel-based minimization. This kernel-
based minimization leads to a new iterative scheme that backpropagates an
unsupervised loss through a pre-trained non-blind solver to update the blur
kernel. Experimental results show that the proposed framework achieves
superior results than state-of-the-art blind deconvolution algorithms in low-
light conditions.

###### Index Terms:

photon-limited, low-light, deconvolution, inverse problems, deblurring, shot
noise

##  I Introduction

In low-light imaging applications such as microscopy [1, 2, 3], and astronomy
[4, 5], a consistent question being asked is how to estimate the blur kernel
and deblur the image in the presence of photon shot noise. A natural
formulation of the problem is the Poisson _blind_ deconvolution where the goal
is to simultaneously recover the blur kernel $\mathbf{h}$ and the latent image
$\mathbf{x}$ from the Poisson forward model

$\mathbf{y}=\text{Poisson}(\alpha\mathbf{h}\circledast\mathbf{x}),$ (1)

where $\circledast$ denotes the convolution. The constant $\alpha$ here is a
parameter that determines the mean photon level of the image
$\mathbf{h}\circledast\mathbf{x}$. For low-light photography problems,
$\alpha$ can be as low as a few photons per pixel, assuming that the latent
image $\mathbf{x}$ is normalized to the range of $[0,1]$.

Figure 1: Proposed blind deconvolution framework. The proposed method is
unsupervised and iterative. At the $k$th iteration, the algorithm takes the
current estimate of the blur kernel $\mathbf{h}^{k}$ and feeds forward through
a nonblind solver to obtain an image estimate $F(\mathbf{y},\mathbf{h}^{k})$.
An unsupervised loss is then used to backpropagate and update the kernel.
Since the nonblind solver is differentiable, the gradient of the loss can be
computed.

Because of the low-light condition (i.e., a small $\alpha$), a brute-force
implementation of existing blind deconvolution algorithms, including fine-
tuning state-of-the-art deep-learning models using the Poisson data, would
face three challenges:

  1. (i)

Heavy noise makes kernel estimation hard. Even in the noiseless case, the
joint optimization of the blur kernel and the latent image is known to cause
degenerate solutions according to Levin et al. [6]. This finding is consistent
with many previous algorithms that propose to find the blur kernel. For
example, many propose to obtain the pilot estimates of the images, such as
applying shock filters and extracting edges, and using these pilot estimates
to update the kernel [7, 8, 9]. However, when the signal-to-noise ratio is
low, the kernel estimation would be very challenging because no pilot
estimates of sufficient quality can be constructed.

  2. (ii)

Blur agnostic deep neural networks seldom utilize the forward model to guide
the restoration. Generic deep-learning methods that are agnostic to the actual
forward equation are known to perform worse than those explicitly taking into
account the blurring process, as evident in learning-based linear inverse
solvers [10]. This gap is further widened in heavy-noise conditions. For
Poisson distributions, there are currently very few deep-learning methods that
explicitly model them in the inversion process [11, 12].

  3. (iii)

Iterative methods are important but not used. Many classical methods estimate
the kernel iteratively because the alternating minimization strategy is known
to be effective. Even for non-blind deconvolution, people found iterative
methods perform well [13, 14]. Delbracio et al. [15] argued that iteratively
deblurring the image leads to better and more stable solutions. However, this
concept is missing in today’s deep-learning based blind deconvolution method
because there is no fast, effective, and _differentiable_ non-blind solvers.

|  |
---|---|---
Blurred noisy image | MPR-Net [16] | Proposed
Figure 2:  Demonstration of the proposed method. A Blurred image corrupted
with Poisson noise (photon level $\alpha=20$) is provided on the left. The
blur kernel, unknown to the reconstruction methods, is provided in the inset.
We take a state-of-the-art deblurring supervised method, i.e., MPR-Net [16],
retrain it using Poisson blurred images, and show the corresponding
reconstructions using MPR-Net (middle) and technique proposed in this paper
(right).

Overcoming the challenges mentioned above requires a new design of the blind
deconvolution algorithm. The main idea presented in this paper is a
differentiable framework that back-propagates the loss to update the _blur
kernel_ (instead of the model parameters), as illustrated in Figure 1. The
framework is unsupervised because the non-blind solver shown in the middle of
the figure is pre-trained and fixed. Once the non-blind solver is plugged in,
there is no training. During testing, the blur kernel is updated by back-
propagating an unsupervised loss which is computed by propagating the gradient
of the non-blind solver. The blur kernel is updated iteratively instead of
estimated directly by a separate module.

The proposed method has two advantages compared to the existing state-of-the-
art deep neural network blind deconvolution algorithms: (1) It is
unsupervised. This helps the algorithm to be adaptive to different blur
kernels. (2) It is flexible as any differentiable non-blind solver can be
plugged into this framework.

A snapshot of the main results of the proposed method is shown in Figure 2. In
this example, the photon level is set to be as low as $\alpha=20$ (equivalent
to Gaussian noise of standard deviation $\sigma\approx 60/255$ in terms of the
signal-to-noise ratio). This is a substantially noisier problem compared to
the ones reported in the deconvolution literature. Compared to state-of-the-
art methods such as the MPR-Net [16] that has been fine-tuned using the
Poisson data, the performance of the proposed method is better both in terms
of the peak signal-to-noise ratio (PSNR) and the visual quality.

##  II Related Work

Non-blind deconvolution. Non-blind deconvolution algorithms assume a known
blur kernel [17]. For the Poisson case, methods have been studied in as early
as the 1980s. For example, the Richardson-Lucy algorithm [18, 19] is iterative
method that converges to the maximum likelihood estimate (MLE) of the Poisson
deconvolution problem. A variety of Maximum A Posteriori (MAP) based solutions
with different regularization methods have been used to solve the problem [20,
21, 22]. PURE-LET [23] estimates the clean image by minimizing the Poisson
unbiased risk estimate. Iterative solutions using the ADMM scheme have also
been proposed to incorporate Total Variation (TV), and complex image priors
[24, 20, 25].

Traditional blind deconvolution. Single image blind deconvolution has a long
list of work since the 1990s [26]. A common theme is to jointly optimize for
both the kernel and the image iteratively [9, 8, 7, 27]. Under this joint
optimization framework, many image priors and kernel priors are utilized, such
as total variation (TV) [26], local gradient maxima [28], internal patch
recurrence [29], and sparse wavelet approximation [30]. Among those, a mixture
of $\ell_{1}$-norm and Total Variation (TV) have been proposed in [31] to
deblur images under a high noise regime.

As far as kernel estimation is concerned, people observe that it is more
robust to use image gradients instead of the image itself [9, 32, 8, 33].
Methods such as [7, 34] improve the joint optimization scheme by identifying
the salient gradients of the image and using them in the kernel estimation
process. These methods are often employed in a multi-scale manner, starting
from a coarse kernel and image estimate, and using these estimates at a finer
scale [35, 8].

Despite the empirical success, there is an issue in the framework. Levin et
al. [36, 6] showed that the joint minimization can theoretically converge to a
degenerate solution pair, i.e., it will recover an identity kernel instead of
the true kernel. Levin et al.’s argument was that the joint minimization needs
to simultaneously estimate the kernel and the image. The search space is too
large for the available amount of observations. Their suggestion was to
exploit the asymmetry of the problem by estimating the kernel first because
the kernel has fewer variables than the image. So, it is easier to recover the
kernel than jointly recovering the kernel and the image.

Learning based blind deconvolution. Since the last decade, deep learning
methods have been used to solve the blind deconvolution problem. Different
CNN-architectures such as Scale Recurrent Network (SRN) [37], Deep-Deblur
[38], MPR-Net [16] have achieved impressive performance on single image
deblurring without considering the forward model. Some learning based methods
such as [39, 40, 41, 42] do the opposite by taking the forward imaging model
into consideration. For example, [39] learns the Fourier coefficients of the
deconvolution filter from patches of image and deconvolves the blurred image
by applying the patchwise average of the predicted filter. [43] predicts the
non-uniform motion blur field using a CNN and Markov Random Field (MRF) model.
[44] designs its network structures by unfolding the classic Richardson-Lucy
algorithm. Another line of work is to use deep neural network as implicit
prior for images or kernels [45], [46], which also achieves competitive
performance in blind or non-blind image deblurring tasks.

##  III Iterative Kernel Estimation

###  III-A Main Idea

To explain the proposed approach, it would be useful to start with the
classical blind deblurring algorithm for Gaussian noise. In the classical
setting, the joint optimization of the blur kernel $\mathbf{h}$ and the image
$\mathbf{x}$ is

$(\widehat{\mathbf{x}},\widehat{\mathbf{h}})=\mathop{\underset{\mathbf{x},\mathbf{h}}{\mbox{argmin}}}\;\bigg{\\{}\|\mathbf{y}-\mathbf{h}\circledast\mathbf{x}\|_{2}^{2}+\lambda
R(\mathbf{x})+\gamma S(\mathbf{h})\bigg{\\}},$ (2)

for some regularization functions $R(\mathbf{x})$ and $S(\mathbf{h})$, and
$\circledast$ denotes the convolution operator. The standard strategy to solve
the joint optimization is to alternate between $\mathbf{x}$ and $\mathbf{h}$
by fixing one and updating the other.

For Poisson noise, the likelihood function changes from the $\ell_{2}$-norm
squares to the Poisson-likelihood:

$\mathcal{P}(\mathbf{y},\mathbf{h}\circledast\mathbf{x})\overset{\text{def}}{=}\boldsymbol{1}^{T}(\mathbf{h}\circledast\mathbf{x})-\mathbf{y}^{T}\log(\mathbf{h}\circledast\mathbf{x})+\boldsymbol{1}^{T}\log\mathbf{y}!$
(3)

where the last term $\boldsymbol{1}^{T}\log\mathbf{y}!$ can be dropped because
it does not depend on $\mathbf{x}$ and $\mathbf{h}$. Following the same
alternating minimization principle, the algorithm for the Poisson blind
deblurring consists of two steps:

$\displaystyle\mathbf{x}^{k+1}$ $\displaystyle=\underset{\text{non-blind
Poisson
deblurring}}{\underbrace{\mathop{\underset{\mathbf{x}}{\mbox{argmin}}}\;\;\mathcal{P}(\mathbf{y},\mathbf{h}^{k}\circledast\mathbf{x})+\lambda
R(\mathbf{x})}}\overset{\text{def}}{=}F(\mathbf{y},\mathbf{h}^{k}),$ (4)
$\displaystyle\mathbf{h}^{k+1}$
$\displaystyle=\mathop{\underset{\mathbf{h}}{\mbox{argmin}}}\;\;\mathcal{P}(\mathbf{y},\mathbf{h}\circledast\mathbf{x}^{k+1})+\gamma
S(\mathbf{h}).$ (5)

A key observation here is that the $\mathbf{x}$-subproblem in (4) is a Poisson
proximal map using some regularization $R(\mathbf{x})$. We define this mapping
as a function $F(\mathbf{y},\mathbf{h}^{k})$ which takes the corrupted image
$\mathbf{y}$ and the current estimate of the kernel $\mathbf{h}^{k}$ to
produce a deblurred image $\mathbf{x}^{k+1}$. Since $\mathbf{h}$ is known and
fixed at $\mathbf{h}^{k}$, the mapping $F(\mathbf{y},\mathbf{h}^{k})$ is a
_non-blind_ Poisson deblurring method. As will be elaborated later, this non-
blind Poisson deblurring step is implemented via a deep neural network.

The first step of the proposed method is to substitute
$F(\mathbf{y},\mathbf{h}^{k})$ into the $\mathbf{h}$-subproblem in (5) and
merge the pair of alternating equations (4) and (5) into one. This will give
us

$\mathbf{h}^{k+1}=\mathop{\underset{\mathbf{h}}{\mbox{argmin}}}\;\;\bigg{\\{}\mathcal{P}(\mathbf{y},\mathbf{h}\circledast
F(\mathbf{y},\mathbf{h}^{k}))+\gamma S(\mathbf{h})\bigg{\\}}.$ (6)

Notice here (6) has no difference with the equation pair (4) and (5) as far as
implementation is concerned: we fix $\mathbf{h}^{k}$ and compute
$F(\mathbf{y},\mathbf{h}^{k})$, and then we estimate the kernel by assuming a
fixed $F(\mathbf{y},\mathbf{h}^{k})$.

The second and a very important step is to recognize that (6) is still a two-
variable optimization because $F(\mathbf{y},\mathbf{h}^{k})$ is the non-blind
solver: when we run $F(\mathbf{y},\mathbf{h}^{k})$, we are solving the
optimization (4). The innovation here is to turn this two-variable alternating
optimization in $(\mathbf{x},\mathbf{h})$ into a one-variable optimization in
$\mathbf{h}$. To do so, we recognize (6) takes the form of a fixed-point
equation — given $\mathbf{h}^{k}$, plug it into an equation, obtain
$\mathbf{h}^{k+1}$, and repeat. We propose to consider the _equilibrium_ of
this fixed-point iteration. The equilibrium is obtained by dropping the
indices $k$ and $k+1$ from the equation. This will lead to an optimization
that does not involve the iteration index $k$:

$\mathbf{h}=\mathop{\underset{\mathbf{h}}{\mbox{argmin}}}\;\;\bigg{\\{}\mathcal{P}(\mathbf{y},\mathbf{h}\circledast
F(\mathbf{y},\mathbf{h}))+\gamma S(\mathbf{h})\bigg{\\}}.$ (7)

Although in terms of notations, the difference between (7) and (6) is subtle,
the physical significance is huge. Figure 3 shows a pictorial illustration. In
(7), the optimization is completely in $\mathbf{h}$. The variable $\mathbf{x}$
is never involved. The way to visualize is that while $F$ theoretically takes
the form of an optimization which is (4), it is nonetheless just a neural
network. As long as we know the input-output relationship of this neural
network, it can be absorbed as a part of the forward image formation model.
Because it is now part of the image formation model, the optimization in
$\mathbf{h}$ just needs to know how to take gradients of this nonlinear
forward model.

Figure 3: Conceptual illustration of the proposed method. In conventional
alternating minimization, we alternatingly update the image and the blur
kernel through two separate modules. In the proposed approach, we leverage the
differentiability of the non-blind Poisson solver and directly estimate the
blur kernel.

At this point, one may anticipate that a reasonable algorithm would be to run
a gradient descent step to minimize the objective function. However, there is
a caveat of adopting a gradient-based algorithm — the Poisson-likelihood is
not differentiable at zero and hence it will cause difficulties. Besides, (7)
asks us to estimate the kernel $\mathbf{h}$ directly from the noisy image
$\mathbf{y}$. This is not a good approach because the Poisson noise is often
very strong relative to the signal. Even in the classical literature, people
have proposed to estimate $\mathbf{h}$ from some variants of $\mathbf{y}$
(such as pixels around the strong edges) instead of the original noisy-blurry
image $\mathbf{y}$.

Following the above argument, the proposed method introduces the third step by
applying a Poisson _denoiser_ $G$ to the input noisy image. That is, we define

$G(\mathbf{y})=\mathop{\underset{\mathbf{x}}{\mbox{argmin}}}\;\;\bigg{\\{}\boldsymbol{1}^{T}\mathbf{x}-\mathbf{y}^{T}\log(\mathbf{x})+\lambda
R(\mathbf{x})\bigg{\\}},$ (8)

for some regularization function $R(\mathbf{x})$. Just like the non-blind
deblurring step, $G(\mathbf{y})$ can be implemented via a neural network.

The next big observation here is that while $\mathbf{y}$ is the Poisson count
that ranges from zero to infinity, the estimate $\mathbf{x}$ is normalized to
the range of $[0,1]$. Therefore, while $\mathbf{y}$ exhibits the Poisson
characteristics, the estimate $G(\mathbf{y})$ does not have any Poissonian
structure and only contains the blur. The “noise” remaining in $G(\mathbf{y})$
is algorithm-dependent because a different denoising algorithm will generate a
different residue pattern. From a statistical point of view, we can
empirically find out the distribution of the residue but practically this
would require access to the ground-truth. However, since $G(\mathbf{y})$ is
distributed in the range of $[0,1]$, it would match with the range of the
prediction $\mathbf{h}\circledast F(\mathbf{y},\mathbf{h})$. This gives us an
approximated objective function:

$\mathbf{h}=\mathop{\underset{\mathbf{h}}{\mbox{argmin}}}\;\;\bigg{\\{}\|G(\mathbf{y})-\mathbf{h}\circledast
F(\mathbf{y},\mathbf{h}))\|_{2}^{2}+\gamma S(\mathbf{h})\bigg{\\}}.$ (9)

The choice of the $\ell_{2}$-norm squares is more of convenience than a
rigorous statistical reasoning because we do not know the distribution of the
residue $G(\mathbf{y})-\mathbf{h}\circledast F(\mathbf{y},\mathbf{h})$.
However, this slackness does not seem to cause significant drawbacks in our
experimental results.

The optimization in (9) is the core of the proposed idea. Comparing this with
the alternating minimization in (4) and (5), the new formulation has two
advantages:

  1. 1.

Optimization in Kernel Space. The new formulation pushes all the efforts to
estimating the kernel. This is consistent with the literature such as [8, 7]
that also spend most of the time trying to obtain a good kernel. Also, as
suggested in Levin et al. [6], we perform estimation in the smaller kernel
space $\mathbf{h}$ instead of joint estimation which has shown to converge to
the no-blur degenerate solution.

  2. 2.

Unsupervised Loss (9) enables _unsupervised_ estimation of $\mathbf{h}$. Note
that the neural network $F$ and $G$ and pre-trained and fixed under the
proposed framework. No ground truth is needed as far as the estimation of the
kernel is concerned.

###  III-B Choice of $F$ and $G$

Since the premise of the problem is low-light, the Poissonian structure has to
be handled. In the proposed approach, the Poisson part is realized through the
choice of $F$ and $G$.

For the choice of $F$, its goal is to take a blur kernel $\mathbf{h}$ and a
noisy-blurry image $\mathbf{y}$ to return an deblurred image
$F(\mathbf{y},\mathbf{h})$, which is essentially (4). There are many non-blind
Poisson deconvolution solvers. This paper uses a recent method by Sanghvi et
al. [11]. In this approach, the minimization is solved via an algorithm
unrolling that unrolls the plug-and-play alternating direction method of
multiplier (PnP-ADMM) into a chain of repeated blocks of steps. The advantage
of the unfolded network is that $F(\cdot)$ constructed in this way will be
_differentiable_ , a key property required to solve (9).

The denoiser $G$ follows equation (8). $G$ is a special case of $F$ where
there is no blur kernel. Unlike a generic denoiser, $G$ is specifically
trained to remove the Poisson noise from the noisy-blurred input and keep the
blur intact so it can serve as a target for the estimated blurred image i.e.
$\mathbf{h}\circledast F(\mathbf{y},\mathbf{h})$. The impact of $G$ is
illustrated in Figure 4. We consider a noisy-blurry image
$\mathbf{y}=\text{Poisson}(\alpha\mathbf{x})$ and take a cross-section of the
pixels. Because of the noise, the cross-section plot is extremely noisy and
blurry. The Poisson denoiser $G$ will remove the noise but the blur along the
edges is preserved. As a result, estimating the blur kernel from
$G(\mathbf{y})$ would become easier.

Figure 4: Why use $G(\mathbf{y})$? (Left) Using the denoiser $G(\mathbf{y})$,
we aim to remove the Poisson noise from the blurred image. This image serves
as a target for the data fidelity term in the cost function we optimize later.
(right) 1D scan of the middle row of the 3 images provided on the left.
$G(\mathbf{y})$ removes the noise from the image $\mathbf{y}$ but preserves
the blurred edges.

The regularization function $\lambda R(\mathbf{x})$ used to define $F$ and $G$
is a symbolic place-holder. During the actual implementation, since both $F$
and $G$ are the unfolded PnP-ADMM, the regularization function $\lambda
R(\mathbf{x})$ is realized via a denoising network, as documented in several
prior work [47, 48, 49]. The exact value of the parameter $\lambda$ and the
explicit form of the regularization $R(\mathbf{x})$ are unimportant.

###  III-C Iterative Scheme

We now discuss how to solve the minimization in (9). Define the following as
the data fidelity term:

$\mathcal{L}(\mathbf{h})\overset{\text{def}}{=}\|G(\mathbf{y})-\mathbf{h}\circledast
F(\mathbf{y},\mathbf{h}))\|_{2}^{2}$ (10)

so that the overall optimization becomes

$\mathbf{h}=\mathop{\underset{\mathbf{h}}{\mbox{argmin}}}\;\underset{\text{Overall
Loss}}{\underbrace{\mathcal{L}(\mathbf{h})+\gamma S(\mathbf{h})}},$

where we define the overall loss as the sum of the data fidelity term
$\mathcal{L}(\mathbf{h})$ and the regularization $\gamma S(\mathbf{h})$.

To decouple the minimization of the first term $\mathcal{L}(\mathbf{h})$ from
the minimization of $\mathcal{S}(\mathbf{h})$ (which can potentially be non-
differentiable such as $\mathcal{S}(\mathbf{h})=\|\mathbf{h}\|_{1}$), we
perform variable splitting as done in a standard manner in half-quadratic
splitting. More specifically, we introduce an additional variable $\mathbf{v}$
and enforce the constraint $\mathbf{h}=\mathbf{v}$ using a quadratic penalty.
This leads us to the following iterative steps to minimize the cost function
in (10)

$\displaystyle\mathbf{h}^{k+1}$
$\displaystyle=\mathop{\underset{\mathbf{h}}{\mbox{argmin}}}\;\mathcal{L}(\mathbf{h})+\frac{\mu}{2}\|\mathbf{h}-\mathbf{v}^{k}\|_{2}^{2},$
(11) $\displaystyle\mathbf{v}^{k+1}$
$\displaystyle=\mathop{\underset{\mathbf{v}}{\mbox{argmin}}}\;\;\gamma
S(\mathbf{v})+\frac{\mu}{2}\|\mathbf{v}-\mathbf{h}^{k+1}\|_{2}^{2},$ (12)

for some hyperparameter $\mu$.

The $\mathbf{h}$-subproblem (11) is solved using gradient descent. To save
computation, $\mathbf{h}$ can be updated _inexactly_ via one gradient descent
step:

$\mathbf{h}^{k+1}=\mathbf{h}^{k}-\delta\cdot\Bigg{\\{}\underset{\text{backpropagate
}F}{\underbrace{\nabla_{\mathbf{h}}\big{\\{}L(\mathbf{h})\big{\\}}\Big{|}_{\mathbf{h}=\mathbf{h}^{k}}}}+\mu(\mathbf{h}^{k}-\mathbf{v}^{k})\Bigg{\\}},$
(13)

where $\delta$ is the gradient descent step size. The gradient can be found
using automatic differentiation, i.e., backpropagation through the non-blind
deblurring algorithm. In popular deep learning packages such as `PyTorch` and
`TensorFlow`, the implementation is done by the `autograd` function in these
packages. In one line, namely `loss.backward()`, one can compute the gradients
of all the variables involved in the calculation of the variable `loss`.

For the $\mathbf{v}$-subproblem (12), a natural choice of the regularization
function is $S(\mathbf{v})=\|\mathbf{v}\|_{1}$. Other regularization functions
can also be used; however, if the problem of interest mainly concerns about
the support of kernel such as motion blur, then $\ell_{1}$-norm is a
reasonable choice. Computationally, if $S(\mathbf{v})=\|\mathbf{v}\|_{1}$,
there exists a closed-form solution using the shrinkage formula:

$\displaystyle\mathbf{v}^{k+1}$
$\displaystyle=\max\left(\left|\mathbf{h}^{k+1}\right|-\gamma/\mu,0\right)\cdot\text{sign}(\mathbf{h}^{k+1})$
$\displaystyle\overset{\text{def}}{=}\mathcal{S}_{\gamma/\mu}(\mathbf{h}^{k+1}).$
(14)

###  III-D Overall Algorithm and Initialization

The overall algorithm is shown in Algorithm 1. The algorithm contains a few
additional modifications to make the algorithm more robust. For example, the
parameter $\mu$ and $\gamma$ are updated according to the following
heuristics. $\mu$ is increased by a constant factor to ensure that the two
variables representing the kernel estimate $\mathbf{h},\mathbf{v}$ converge to
the same value. $\gamma$, which represents the strength of the $\ell_{1}$ norm
prior, is gradually decreased over the iterations to allow for sparse
solutions initially and then for less sparse solutions in later iterations. In
our implementation, the operators $G$ and $F$ also requires a photon level
estimate (equivalent to the noise level in the Gaussian case.)

Due to the ill-posedness of the blind deconvolution problem, the
initialization of the kernel estimate plays an important role. The goal of our
initialization schemes is to obtain a reasonable starting point which is
computationally inexpensive. For this paper, we use the method provided in
[15]. Specifically, the method assumes the blur kernel as a tilted anisotopic
Gaussian kernel parametrized by 3 parameters - the major axis, minor axis, and
tilt of the kernel. We use the denoised image $G(\mathbf{y})$ as input to the
scheme and use the resulting anisotropic Gaussian kernel as initial kernel
estimate $\mathbf{h}^{0}$.

1:Input: Noisy-blurry $\mathbf{y}$, denoiser $G(\cdot)$, and non-blind solver
$F(\cdot)$.

2:Denoise $\mathbf{y}$ to obtain a noiseless blurred image $G(\mathbf{y})$

3:Initialize $\mathbf{h}^{0}$ and $\mathbf{v}^{0}$.

4:$\mu\leftarrow 2.0$, $\gamma\leftarrow 10^{-3}$

5:for $k=0,1,2,\cdot\cdot\cdot$ do

6: $\mathcal{L}(\mathbf{h})\leftarrow\|G(\mathbf{y})-\mathbf{h}\circledast
F(\mathbf{y},\mathbf{h})\|_{2}^{2}$

7: Calculate $\nabla_{\mathbf{h}}\mathcal{L}(\mathbf{h}^{k})$ using automatic
differentiation

8:
$\mathbf{h}^{k+1}\leftarrow\mathbf{h}^{k}-\delta\big{(}\nabla_{\mathbf{h}}\mathcal{L}(\mathbf{h}^{k})+\mu(\mathbf{h}^{k}-\mathbf{v}^{k})\big{)}$

9: $\mathbf{v}^{k+1}\leftarrow\mathcal{S}_{\gamma/\mu}(\mathbf{h}^{k+1})$

10: $\mu\leftarrow 1.01\mu$, $\gamma\leftarrow\gamma/1.01$

11:end for

12:return $\mathbf{h}^{(\infty)}$ and
$\mathbf{x}^{(\infty)}=F(\mathbf{y},\mathbf{h}^{(\infty)})$

Algorithm 1 Iterative Poisson Deconvolution Scheme

##  IV Experiments

Method $\rightarrow$ | Two-Phase | SRN | DMPHN | Deep-Deblur | MPRNet | Poisson PnP | PURE-LET | Ours | P4IP
---|---|---|---|---|---|---|---|---|---
Photon lvl $\downarrow$ | [7] | [37] | [50] | [38] | [16] | [25] | [23] |  |  [11]
$\alpha=10$ | 16.42 | 20.33 | 20.25 | 20.92 | 21.03 | 19.83 | 21.63 | 22.10 | 22.45
0.511 | 0.511 | 0.509 | 0.523 | 0.533 | 0.464 | 0.607 | 0.598 | 0.639
$\alpha=20$ | 17.40 | 20.46 | 20.43 | 21.11 | 21.34 | 19.20 | 21.82 | 22.52 | 22.80
0.558 | 0.523 | 0.524 | 0.537 | 0.552 | 0.441 | 0.622 | 0.622 | 0.665
$\alpha=40$ | 17.84 | 20.56 | 20.51 | 21.21 | 21.53 | 17.34 | 21.76 | 22.67 | 22.99
0.565 | 0.533 | 0.531 | 0.545 | 0.566 | 0.382 | 0.634 | 0.638 | 0.690
Blind? | ✓ | ✓ | ✓ | ✓ | ✓ | ✕ | ✕ | ✓ | ✕
Network? | ✕ | ✓ | ✓ | ✓ | ✓ | ✓ | ✕ | ✓ | ✓
Unsupervised? | ✓ | ✕ | ✕ | ✕ | ✕ | ✓ | ✓ | ✓ | ✕
TABLE I: Performance on Levin et al. Dataset [6]: (Top) Average PSNR in dB (Bottom) SSIM. The last column, P4IP, is a non-blind deconvolution method, used as $F(\cdot)$ in the iterative blind deconvolution method and serves as an upper bound of the iterative scheme. The gray-colored columns represent the unsupervised methods which cannot be trained end-to-end. Photon lvl $\rightarrow$ | $\alpha=10$ | $\alpha=20$ |  $\alpha=40$
---|---|---|---
Method $\downarrow$
DMPHN [50] | 23.88 | 24.26 | 24.20
SRN [37] | 23.89 | 24.05 | 24.00
Deep-Deblur [38] | 24.57 | 24.45 | 24.65
MPRNet [16] | 25.50 | 25.93 | 25.97
Ours | 26.49 | 27.43 | 26.67
TABLE II: Performance on Real-Blur Dataset [51]: Average PSNR in dB.

###  IV-A Training $F(\cdot)$ and $G(\cdot)$

In this subsection, we describe the training procedure for the non-blind
solver $F(\cdot)$ and the denoiser $G(\cdot)$. These training processes, along
with other experiments described in paper, are implemented in PyTorch 1.7.0,
and use an NVIDIA Titan Xp GP102 GPU.

For the non-blind solver $F(\cdot)$, we use a similar process as that of [11].
We use 2650 clean images from the Flickr2K dataset [52] and divide them into a
80:20 training and validation dataset. We generate 60 motion kernels described
using the code provided in [53]. To generate synthetically blurred-noisy
image, clean images from the Flickr2K training subset are blurred by the
motion kernels followed by the Poisson noise corruption at photon level
$\alpha$. The photon level $\alpha$ is uniformly sampled from $[1,60]$. Using
the pairs of noisy-blurred images, corresponding blur kernels and the clean
images, we train PhD-Net ($K=8$ iterations) from [11] using the $\ell_{1}$
loss.

For training the blur-noisy to blur-only denoiser $G(\cdot)$, we use a similar
method. We use synthetically generated noisy-blurred images from the training
process of non-blind solver $F(\cdot)$. However, there are a few key
distinctions here. We use the input as noisy-blurred images and don’t use the
relevant blur kernel. The target for training is _not the clean image but the
blur-only image_. For the network architecture, we use the same configuration
as that of PhD-Net, but change we fix the kernel input to be the identity
operator. We use the $\ell_{1}$ loss function to train the denoiser.

###  IV-B Quantitative Comparison

We compare our proposed iterative scheme with five other blind deconvolution
methods in this paper: a classical Two-Phase Estimation [7], Scale-Recurrent
Network [37], Deep-Deblur [38], Deep-Hierarchical Multi-Patch Network [50],
and MPR-Net [16]. When the ground-truth blur kernel is available, we also
include the following non-blind deconvolution algorithms in the comparison:
PURE-LET [23], Poisson PnP [25]. Since Two-Phase estimation [7] is an
unsupervised method designed for noiseless images, we use the denoiser output
$G(\mathbf{y})$ as the input to the algorithm instead of Poisson corrupted
image $\mathbf{y}$.

For the supervised end-to-end trainable methods i.e., [37, 38, 50, 16], we
retrain the networks for deconvolving Poissonian images in the same manner as
training $F(\cdot)$ and $G(\cdot)$. Specifically, we use the artificially
generated clean and noisy-blurred image pairs from Fickr2K dataset as
described in the previous subsection, are retrain them with Adam optimizer of
learning rate $10^{-4}$, $\ell_{1}$ loss functions until training converges.

We remark that the training procedure outlined above is a fair setting for all
the methods we consider in this paper. The reason is that supervised methods
such as [16, 38, 50, 37] have a an unfair advantage over our scheme because
they can see the ground truth. The ground truths allow them to handle blurs
beyond a simple spatially invariant blur we assume (which is also assumed by
all other unsupervised blind deconvolution methods). Since in this paper we
are focused on a spatially invariant blur, we argue that it is more useful to
compare all methods using the same forward image formation model. Future work
could address how the scheme can be modified to take into account more complex
forward models such as spatially varying blur and blur due to object motion.

Levin Dataset [6]: Once the supervised methods are retrained using Poisson
data, we compare the different schemes on the Levin et al. dataset [6] which
consists of 4 images, each blurred by 8 different motion kernels. The Poisson
noise is artificially added at different light levels $\alpha=10,20,40$. The
results of the comparison are summarized in Table I. Note that in Table I, we
also compare the blind deconvolution methods with the non-blind deconvolution
network P4IP [11]. While this network is used iteratively in our scheme as
$F(\cdot)$, we provide its performance using the ground truth blur kernel, to
serve as an upper bound to the performance of blind deconvolution algorithms.

Real-Blur Dataset [51]: In addition to the Levin dataset, we perform a
quantitative evaluation on the Real-Blur dataset [51] which contains pairs of
sharp and blurred images for 232 different scenes with 20 image pairs per
scene. These pair of images are captured using a specially designed image
acquisition system that contains a beam splitter and two cameras - one for
ground-truth operating at $1/80$s exposure time and the other at $1/2$s for
the blurred images.

For our experiment, we evaluate different blind-deconvolution methods on a
randomly chosen $256\times 256$ patch from a ground-truth, blurred pair from
50 different scenes in RealBlur-J subset. We degrade the blurred images with
Poisson noise at photon level $\alpha=10,20,\text{ and }40$ and the resulting
PSNR values are shown in Table II. Since the blur kernel is not available for
this dataset, unlike Levin, the non-blind deconvolution methods i.e. Poisson-
PnP [25], PURE-LET [23], P4IP [11] cannot be evaluated on this dataset.

Input 16.35 / 0.280 |  Two-Phase [7] 17.82 / 0.640 |  SRN [37] 18.89 / 0.459 |  DHMPN [50] 18.70 / 0.452 |  Deep-Deblur [38] 18.85 / 0.455
---|---|---|---|---
MPR-Net [16] 19.11 / 0.477 |  Ours 19.86 / 0.599 |  Non-Blind [11] 21.93 / 0.621 |  Ground-Truth

Input 15.72 / 0.278 |  Two-Phase [7] 12.99 / 0.599 |  SRN [37] 17.84 / 0.486 |  DMPHN [50] 17.63 / 0.483 |  Deep-Deblur [38] 17.69 / 0.475
---|---|---|---|---
MPR-Net [16] 18.48 / 0.551 |  Ours 19.76 / 0.639 |  Non-Blind [11] 20.74 / 0.652 |  Ground-Truth

Input 14.73 / 0.191 |  Two-Phase [7] 14.73 / 0.191 |  SRN [37] 18.91 / 0.358 |  DMPHN [50] 18.87 / 0.356 |  Deep-Deblur [38] 19.11 / 0.365
---|---|---|---|---
MPR-Net [16] 19.00 / 0.358 |  Ours 20.03 / 0.459 |  Non-Blind [11] 21.99 / 0.549 |  Ground-Truth

Figure 5: Visual comparisons. Qualitative results of the competing methods on
synthetically blurred and Poisson corrupted images from the BSD100 dataset.
The specific blur kernel used in each experiment is shown in the inset of the
input. The photon level for car, viking, speed-limit is $\alpha=40,40$ and
$20$ respectively.

Input |  Input |  SRN [37] |  Deep-Deblur [38]
---|---|---|---
MPR-Net [16] |  Ours |  Ground-Truth

Input |  Input |  SRN [37] |  Deep-Deblur [38]
---|---|---|---
MPR-Net [16] |  Ours |  Ground-Truth

Figure 6: Visual comparisons on Realistic Blur. Qualitative results on realistically blurred images from the RealBlur dataset [51]. Photon levels for the first and second images are $\alpha=20$ and $\alpha=10$ respectively.  |  |  |
---|---|---|---
Blurred and Noisy | MPR-Net [16] | Ours | Non-Blind [11]
|  |  |
Blurred and Noisy | MPR-Net [16] | Ours | Non-Blind [11]
Figure 7: Real world data from [11]. Note that the the blur kernel is unknown
to the reconstruction methods except the non-blind solver in the last column.

###  IV-C Qualitative Comparison

A qualitative comparison of grayscale and color reconstructions using
different schemes is provided in Figure 5. We also provide the corresponding
result using the non-blind solver from [11] as a reference since this method
assumes the ground truth kernel to be known. For reconstructing color images,
we first describe how to simulate synthetic blur and noise and then how to
modify the scheme from deconvolving grayscale to color images.

To simulate blur and Poisson noise in color images, we convert the RGB image
to the corresponding Bayer pattern image - which itself can be viewed as 4
channels, namely R, G1, G2, B, interleaved with each other. To simulate
photon-limited blur, we blur each channel (R, G1, G2, and B) individually and
add Poisson noise at the same photon level.

For our deconvolution scheme, we pick a single channel, say R, and apply the
kernel estimation process described in Section 3 to it. Using the kernel
estimate obtained, we deconvolve each of the four channels with the non-blind
solver $F(\mathbf{y},\mathbf{k})$. The deconvolved channels are combined into
an RGB image using an off-the-shelf demosaicing method. For other schemes such
as [7, 37, 38, 16], we perform deconvolution for each channel followed by
demosaicing.

###  IV-D Qualitative Evaluation on Real World Data

We also demonstrate that our iterative scheme can reconstruct real-world
blurred and noisy images. First we demonstrate our qualitative reconstructions
on examples from the RealBlur dataset in Figure 6. The blurred images in the
dataset are degraded by Poisson shot noise at photon level $\alpha=20$ and
$\alpha=10$.

Next, we demonstrate the iterative scheme proposed in this paper on patches
from the dataset provided in [11]. This dataset consists of 30 images
collected using a DSLR camera with blur in the images generated using handheld
motion. Due to the photon-limited setting of these scenes, the images are
naturally corrupted by the photon shot noise, which follows a Poisson
distribution. Hence, unlike the RealBlur dataset, this contains both realistic
blur and photon-shot noise. While the blur kernel is also estimated using a
point source, it is not used in the reconstruction scheme in this paper. For
further details about the dataset, we refer the reader to [11]. The
reconstructions using the dataset are reported in Figure 7 along with the
reconstruction using MPR-Net [16] and the non-blind solver
$F(\mathbf{y},\mathbf{k})$ which uses the ground truth kernel.

Note that for images with realistic blur, circular boundary conditions for
convolution lead to artifacts and symmetric boundary conditions are more
accurate to model real-world blur. As a result, the input to $F(\cdot)$ is
padded symmetrically on both sides before the iterative scheme starts. The
relevant center portion of the image is cropped out of the output
$F(\mathbf{y},\mathbf{h})$.

###  IV-E Kernel Estimation

We quantitatively compare the estimated kernel. In Table III, we show the mean
absolute error for each of the eight kernels in the Levin dataset [54]. The
mean absolute error (MAE) is defined as follows:

$\displaystyle\text{MAE}\overset{\text{def}}{=}\|\widehat{\mathbf{h}}-\mathbf{h}\|_{1}/M$
(15)

where $\widehat{\mathbf{h}}$, $\mathbf{h}$ are the estimated and the ground-
truth kernels respectively, and $M$ represents the number of pixels of the
kernel. We compare the estimated kernel with that of the Two-Phase estimation
[7] which takes $G(\mathbf{y})$ as the input and gives the estimated kernel.
We observe a consistently better kernel estimate using our proposed method.

Kernel $\rightarrow$ Photon Level $\downarrow$ |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---
| [7] | Ours | [7] | Ours | [7] | Ours | [7] | Ours | [7] | Ours | [7] | Ours | [7] | Ours | [7] | Ours
$\alpha=10$ | 1.40 | 1.08 | 1.36 | 1.11 | 1.03 | 0.84 | 1.53 | 1.47 | 1.13 | 0.92 | 1.38 | 1.16 | 1.47 | 1.27 | 1.46 | 1.29
$\alpha=20$ | 1.37 | 1.15 | 1.28 | 1.12 | 0.91 | 0.86 | 1.51 | 1.40 | 1.07 | 0.99 | 1.34 | 1.09 | 1.45 | 1.28 | 1.47 | 1.31
$\alpha=40$ | 1.25 | 1.27 | 1.33 | 1.26 | 0.94 | 0.88 | 1.52 | 1.42 | 1.08 | 1.01 | 1.26 | 1.12 | 1.42 | 1.26 | 1.42 | 1.33
TABLE III: Mean absolute error (MAE) on kernels estimated on the Levin
dataset.

##  V Ablation Study

###  V-A Non-Blind Solver

Since a photon-limited non-blind deconvolution solver is at the core of the
proposed scheme, we evaluate the performance of the scheme using different
non-blind solvers. For this scheme, the non-blind solver $F(\cdot)$ can be any
differentiable solver which takes as input a noisy-blurred image and a kernel
input. Therefore, we compare the performance of the method on 3 different non-
blind solvers in Table IV \- P4IP [11], DWDN [55], and RGDN [13]. We use the
pre-trained P4IP model for the experiment and the latter non-blind solvers are
retrained using Poisson data.

The results for this comparison are provided in Table IV. We can see that the
performance on the non-blind and blind deconvolution task are correlated to
each other. Since P4IP has the best performance on non-blind deconvolution
task, it outperforms other methods on the blind-deconvolution problem as well.
However, while DWDN outperforms RGDN in the non-blind task, it only
outperforms RGDN in blind deconvolution at photon level $\alpha=40$.

Photon Level $\rightarrow$ | Blind / | $\alpha=10$ | $\alpha=20$ |  $\alpha=40$
---|---|---|---|---
Solver $\downarrow$ | Non-Blind |
P4IP [11] | Non-Blind | 24.70 | 25.89 | 27.10
Blind | 22.26 | 22.92 | 22.99
DWDN [55] | Non-Blind | 23.93 | 24.87 | 25.75
Blind | 19.03 | 19.87 | 23.20
RGDN [13] | Non-Blind | 22.41 | 22.97 | 23.44
Blind | 21.95 | 22.46 | 22.68
TABLE IV: Ablation Study for Effect of Non-Blind Solver The output PSNR (in
dB) on Levin dataset for different non-blind solvers (retrained for Poisson
noise) when plugged into the iterative blind deconvolution scheme.

###  V-B Hyperparameters

In addition to studying the effect of non-blind solvers, we also evaluate the
effect of different hyperparameters on the performance of the scheme in Table
V. The experiments in this subsection can be divided into 3 different
categories. (i) We evaluate the effect of the denoiser $G(\cdot)$ and
$\ell_{1}$ kernel prior on the performance of the iterative method. (ii) We
then evaluate how the performance varies with different strength of the
$\ell_{1}$ kernel prior by varying the parameter $\gamma$. (iii) We also vary
the maximum number of iterations of the iterative scheme.

From Table V, we observe that denoising the image $G(\mathbf{y})$ is an
important feature of the scheme and improves the performance by approximately
$4$dB when used as a target for the $\ell_{2}$ loss. Compared to denoiser
$G(\cdot)$, the $\ell_{1}$ kernel prior doesn’t have as significant
improvement on the scheme in lower photon levels i.e. $\alpha=10,20$. However,
at $\alpha=40$, it adds another $1.5$dB to the blind-deconvolution
performance. From the third experiment, we can infer that performance of the
scheme slightly degrades after $I=20$ iterations.

Photon Level $\rightarrow$ | $\alpha=10$ | $\alpha=20$ | $\alpha=40$ |  $\alpha=60$
---|---|---|---|---
Ablation Study $\downarrow$ | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
w/o $\ell_{1}$ kernel prior, w/o denoiser $G(\cdot)$ | 17.32 | 0.375 | 18.23 | 0.406 | 18.83 | 0.447 | 19.24 | 0.470
w/o denoiser $G(\cdot)$ | 18.41 | 0.402 | 19.75 | 0.450 | 20.90 | 0.508 | 21.92 | 0.558
w/o $\ell_{1}$ kernel prior | 22.36 | 0.614 | 23.03 | 0.645 | 21.40 | 0.447 | 23.05 | 0.657
with $\ell_{1}$ kernel prior, $\gamma=10^{-3}$ | 22.38 | 0.590 | 22.78 | 0.613 | 22.93 | 0.630 | 23.03 | 0.643
with $\ell_{1}$ kernel prior, $\gamma=10^{-4}$ | 22.50 | 0.600 | 22.95 | 0.626 | 23.02 | 0.644 | 23.04 | 0.654
with $\ell_{1}$ kernel prior, $\gamma=10^{-5}$ | 22.44 | 0.606 | 22.98 | 0.631 | 22.64 | 0.660 | 23.04 | 0.656
with $\ell_{1}$ kernel prior, $\gamma=10^{-6}$ | 22.29 | 0.610 | 22.96 | 0.639 | 21.59 | 0.642 | 23.06 | 0.658
Initialization only, no iterations | 18.20 | 0.466 | 20.07 | 0.563 | 19.64 | 0.561 | 22.73 | 0.657
maximum iterations $I=10$ | 22.38 | 0.613 | 22.99 | 0.642 | 22.89 | 0.663 | 23.04 | 0.654
maximum iterations $I=20$ | 22.56 | 0.608 | 23.03 | 0.636 | 23.08 | 0.660 | 23.04 | 0.654
maximum iterations $I=100$ | 22.29 | 0.597 | 22.98 | 0.630 | 23.03 | 0.644 | 23.04 | 0.654
TABLE V: Ablation Study we perform ablation studies to better understand the
impact of different hyperparameters on the iterative scheme. First, we test
the performance of the scheme without using the denoiser $G(y)$ as a target
and without using the $\ell_{1}$ prior. We study the effect of varying the
strength of the $\ell_{1}$ prior i.e. hyperparameter $\gamma$ and changing the
number of iteration of the scheme on the performance.

##  VI Conclusion

This paper presented a new iterative scheme for photon-limited blind
deconvolution. The success of the method depends on the following
characteristics of the proposed scheme which were absent in previous methods:
(1) A differentiable and powerful non-blind solver as the backbone. While
classical methods only use hand-crafted image priors, we use the latest deep
learning based non-blind Poisson solver to estimate the latent image. This new
solver is differentiable and thus allows back-propagation. (2) Alternating
minimization is known to have difficulties when noise is strong. We re-
formulate the problem as an equilibrium statement by absorbing the non-blind
Poisson solver into the kernel estimation process. Experimental results
confirm the effectiveness of the method.

Two aspects of the algorithm can be improved: (1) The prior of the kernel can
be better chosen, as the $\ell_{1}$-norm prior is suitable only for a subset
of motion blurs. A more powerful learning-based prior can be considered. (2)
The computational cost needs to be reduced as backpropagation through a large
non-blind solver requires a lot of memory. Some approximation schemes could be
useful to speed up the computation.

## References

  * [1] P. Pankajakshan, B. Zhang, L. Blanc-Féraud, Z. Kam, J.-C. Olivo-Marin, and J. Zerubia, “Blind deconvolution for thin-layered confocal imaging,” _Applied Optics_ , vol. 48, no. 22, pp. 4437–4448, 2009.
  * [2] F. Soulez, L. Denis, Y. Tourneur, and É. Thiébaut, “Blind deconvolution of 3D data in wide field fluorescence microscopy,” in _Proceedings of the IEEE International Symposium on Biomedical Imaging_ , 2012, pp. 1735–1738.
  * [3] J. Chen, R. Lin, H. Wang, J. Meng, H. Zheng, and L. Song, “Blind-deconvolution optical-resolution photoacoustic microscopy in vivo,” _Optics Express_ , vol. 21, no. 6, pp. 7316–7327, 2013.
  * [4] S. M. Jefferies and J. C. Christou, “Restoration of astronomical images by iterative blind deconvolution,” _The Astrophysical Journal_ , vol. 415, p. 862, 1993.
  * [5] T. J. Schulz, “Multiframe blind deconvolution of astronomical images,” _Journal of the Optical Society of America A_ , vol. 10, no. 5, pp. 1064–1073, 1993.
  * [6] A. Levin, Y. Weiss, F. Durand, and W. T. Freeman, “Understanding blind deconvolution algorithms,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 33, no. 12, pp. 2354–2367, 2011.
  * [7] L. Xu and J. Jia, “Two-phase kernel estimation for robust motion deblurring,” in _European Conference on Computer Vision_. Springer, 2010, pp. 157–170.
  * [8] S. Cho and S. Lee, “Fast motion deblurring,” _ACM Transactions on Graphics_ , pp. 1–8, 2009.
  * [9] Q. Shan, J. Jia, and A. Agarwala, “High-quality motion deblurring from a single image,” _ACM Transactions on Graphics_ , vol. 27, no. 3, pp. 1–10, 2008.
  * [10] V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,” _IEEE Signal Processing Magazine_ , vol. 38, no. 2, pp. 18–44, 2021.
  * [11] Y. Sanghvi, A. Gnanasambandam, and S. H. Chan, “Photon limited non-blind deblurring using algorithm unrolling,” _IEEE Transactions on Computational Imaging_ , vol. 8, pp. 851–864, 2022.
  * [12] A. L. Mur, P. Bataille, F. Peyrin, and N. Ducros, “Deep expectation-maximization for image reconstruction from under-sampled poisson data,” in _Proceedings of the IEEE International Symposium on Biomedical Imaging_ , 2021, pp. 1535–1539.
  * [13] D. Gong, Z. Zhang, Q. Shi, A. van den Hengel, C. Shen, and Y. Zhang, “Learning deep gradient descent optimization for image deconvolution,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 31, no. 12, pp. 5468–5482, 2020.
  * [14] Y. Nan, Y. Quan, and H. Ji, “Variational-EM-based deep learning for noise-blind image deblurring,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 3626–3635.
  * [15] M. Delbracio, I. Garcia-Dorado, S. Choi, D. Kelly, and P. Milanfar, “Polyblur: Removing mild blur by polynomial reblurring,” _IEEE Transactions on Computational Imaging_ , vol. 7, pp. 837–848, 2021.
  * [16] S. W. Zamir, A. Arora, S. Khan, M. Hayat, F. S. Khan, M.-H. Yang, and L. Shao, “Multi-stage progressive image restoration,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 14 821–14 831.
  * [17] M. Banham and A. Katsaggelos, “Digital image restoration,” _IEEE Signal Processing Magazine_ , vol. 14, no. 2, pp. 24–41, 1997.
  * [18] W. H. Richardson, “Bayesian-based iterative method of image restoration,” _Journal of the Optical Society of America_ , vol. 62, no. 1, pp. 55–59, 1972\.
  * [19] L. B. Lucy, “An iterative technique for the rectification of observed distributions,” _The Astronomical Journal_ , vol. 79, p. 745, 1974.
  * [20] M. A. Figueiredo and J. M. Bioucas-Dias, “Restoration of Poissonian images using alternating direction optimization,” _IEEE Transactions on Image Processing_ , vol. 19, no. 12, pp. 3133–3145, 2010.
  * [21] Z. T. Harmany, R. F. Marcia, and R. M. Willett, “This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms— theory and practice,” _IEEE Transactions on Image Processing_ , vol. 21, no. 3, pp. 1084–1096, 2011.
  * [22] R. D. Nowak and E. D. Kolaczyk, “A statistical multiscale framework for Poisson inverse problems,” _IEEE Transactions on Information Theory_ , vol. 46, no. 5, pp. 1811–1825, 2000.
  * [23] J. Li, F. Luisier, and T. Blu, “Pure-let image deconvolution,” _IEEE Transactions on Image Processing_ , vol. 27, no. 1, pp. 92–105, 2017.
  * [24] M. A. Figueiredo and J. M. Bioucas-Dias, “Deconvolution of poissonian images using variable splitting and augmented lagrangian optimization,” in _Proceedings of the IEEE/SP Workshop on Statistical Signal Processing_ , 2009, pp. 733–736.
  * [25] A. Rond, R. Giryes, and M. Elad, “Poisson inverse problems by the plug-and-play scheme,” _Journal of Visual Communication and Image Representation_ , vol. 41, pp. 96–108, 2016.
  * [26] T. F. Chan and C.-K. Wong, “Total variation blind deconvolution,” _IEEE Transactions on Image Processing_ , vol. 7, no. 3, pp. 370–375, 1998.
  * [27] Z. Mao, N. Chimitt, and S. H. Chan, “Image reconstruction of static and dynamic scenes through anisoplanatic turbulence,” _IEEE Transactions on Computational Imaging_ , vol. 6, pp. 1415–1428, 2020.
  * [28] L. Chen, F. Fang, T. Wang, and G. Zhang, “Blind image deblurring with local maximum gradient prior,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 1742–1750.
  * [29] T. Michaeli and M. Irani, “Blind deblurring using internal patch recurrence,” in _European Conference on Computer Vision_. Springer, 2014, pp. 783–798.
  * [30] J.-F. Cai, H. Ji, C. Liu, and Z. Shen, “Blind motion deblurring from a single image using sparse approximation,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2009, pp. 104–111.
  * [31] J. Anger, M. Delbracio, and G. Facciolo, “Efficient blind deblurring under high noise levels,” in _2019 11th International Symposium on Image and Signal Processing and Analysis (ISPA)_ , 2019, pp. 123–128.
  * [32] N. Joshi, R. Szeliski, and D. J. Kriegman, “PSF estimation using sharp edge prediction,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2008, pp. 1–8.
  * [33] L. Xu, S. Zheng, and J. Jia, “Unnatural l0 sparse representation for natural image deblurring,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2013, pp. 1107–1114.
  * [34] D. Gong, M. Tan, Y. Zhang, A. Van den Hengel, and Q. Shi, “Blind image deconvolution by automatic gradient activation,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2016, pp. 1827–1836.
  * [35] R. Fergus, B. Singh, A. Hertzmann, S. T. Roweis, and W. T. Freeman, “Removing camera shake from a single photograph,” in _ACM SIGGRAPH 2006 Papers_ , 2006, pp. 787–794.
  * [36] A. Levin, Y. Weiss, F. Durand, and W. T. Freeman, “Understanding and evaluating blind deconvolution algorithms,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2009, pp. 1964–1971.
  * [37] X. Tao, H. Gao, X. Shen, J. Wang, and J. Jia, “Scale-recurrent network for deep image deblurring,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 8174–8182.
  * [38] S. Nah, T. H. Kim, and K. M. Lee, “Deep multi-scale convolutional neural network for dynamic scene deblurring,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , July 2017.
  * [39] A. Chakrabarti, “A neural approach to blind motion deblurring,” in _European Conference on Computer Vision_. Springer, 2016, pp. 221–235.
  * [40] D. Gong, J. Yang, L. Liu, Y. Zhang, I. Reid, C. Shen, A. Van Den Hengel, and Q. Shi, “From motion blur to motion flow: A deep learning solution for removing heterogeneous motion blur,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 2319–2328.
  * [41] C. J. Schuler, M. Hirsch, S. Harmeling, and B. Schölkopf, “Learning to deblur,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 38, no. 7, pp. 1439–1451, 2015.
  * [42] Z. Wang, Z. Wang, Q. Li, and H. Bilen, “Image deconvolution with deep image and kernel priors,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision Workshops_ , 2019, pp. 980–989.
  * [43] J. Sun, W. Cao, Z. Xu, and J. Ponce, “Learning a convolutional neural network for non-uniform motion blur removal,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2015, pp. 769–777.
  * [44] C. Agarwal, S. Khobahi, A. Bose, M. Soltanalian, and D. Schonfeld, “Deep-url: A model-aware approach to blind deconvolution based on deep unfolded richardson-lucy network,” in _2020 IEEE International Conference on Image Processing (ICIP)_ , 2020, pp. 3299–3303.
  * [45] D. Ren, K. Zhang, Q. Wang, Q. Hu, and W. Zuo, “Neural blind deconvolution using deep priors,” in _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2020, pp. 3338–3347.
  * [46] Y. Fang, H. Zhang, H. S. Wong, and T. Zeng, “A robust non-blind deblurring method using deep denoiser prior,” in _2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)_ , 2022, pp. 734–743.
  * [47] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg, “Plug-and-play priors for model based reconstruction,” in _Proceedings of the IEEE Global Conference on Signal and Information Processing_ , 2013, pp. 945–948.
  * [48] S. Sreehari, S. V. Venkatakrishnan, B. Wohlberg, G. T. Buzzard, L. F. Drummy, J. P. Simmons, and C. A. Bouman, “Plug-and-play priors for bright field electron tomography and sparse interpolation,” _IEEE Transactions on Computational Imaging_ , vol. 2, no. 4, pp. 408–423, 2016.
  * [49] R. Ahmad, C. A. Bouman, G. T. Buzzard, S. Chan, S. Liu, E. T. Reehorst, and P. Schniter, “Plug-and-play methods for magnetic resonance imaging: Using denoisers for image recovery,” _IEEE Signal Processing Magazine_ , vol. 37, no. 1, pp. 105–116, 2020.
  * [50] H. Zhang, Y. Dai, H. Li, and P. Koniusz, “Deep stacked hierarchical multi-patch network for image deblurring,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 5978–5986.
  * [51] J. Rim, H. Lee, J. Won, and S. Cho, “Real-world blur dataset for learning and benchmarking deblurring algorithms,” in _European Conference on Computer Vision_. Springer, 2020, pp. 184–201.
  * [52] B. Lim, S. Son, H. Kim, S. Nah, and K. M. Lee, “Enhanced deep residual networks for single image super-resolution,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , July 2017\.
  * [53] G. Boracchi and A. Foi, “Modeling the performance of image restoration from motion blur,” _IEEE Transactions on Image Processing_ , vol. 21, no. 8, pp. 3502–3517, 2012.
  * [54] A. Levin, “Blind motion deblurring using image statistics,” _Advances in Neural Information Processing Systems_ , vol. 19, 2006.
  * [55] J. Dong, S. Roth, and B. Schiele, “Deep wiener deconvolution: Wiener meets deep learning for image deblurring,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 1048–1059, 2020.