<!-- DO NOT EDIT. -->
<!-- THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. -->
<!-- TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: -->
<!-- "auto_examples/blind-inverse-problems/demo_blind_richardsonlucy.py" -->
<!-- LINE NUMBERS ARE GIVEN BELOW. -->

#### NOTE
New to DeepInverse? Get started with the basics with the
[5 minute quickstart tutorial](https://deepinv.org/auto_examples/basics/demo_quickstart.html.md#sphx-glr-auto-examples-basics-demo-quickstart-py)..

<a id="sphx-glr-auto-examples-blind-inverse-problems-demo-blind-richardsonlucy-py"></a>

# Blind Richardson-Lucy deconvolution

This example introduces the Richardson-Lucy algorithm for deconvolution in the blind
setting, where both the underlying clean image and the blur kernel are unknown.

We first consider the non-blind problem

$$
y = h * x,
$$

where $h$ is the known convolution kernel.
Richardson-Lucy <sup>[1](#footcite-richardsonbayesianbasediterativemethod1972)</sup><sup>[2](#footcite-lucyiterativetechniquerectification1974)</sup>
is a deconvolution algorithm used when the data is corrupted by Poisson noise.
Starting from a nonnegative image $x^{(0)}$, it iterates

$$
x^{(k+1)} = \frac{x^{(k)}}{h^\dagger * \mathbf{1}} \odot h^\dagger * \left(\frac{y}{h * x^{(k)}}\right),
$$

where $h^\dagger$ is the spatially flipped kernel, such that $h^\dagger *$ is the adjoint of convolution by $h$. All products and divisions outside convolutions are pointwise.

When $h$ is unknown, blind Richardson-Lucy alternates these updates for the
image $x$ and the kernel $h$. This is a classical, simple baseline
for blind deconvolution. It can use regularization on the image and/or the kernel to
improve performances.

<!-- GENERATED FROM PYTHON SOURCE LINES 34-38 -->

## Setup the physics

Here we consider a simple Gaussian blur kernel and a medium amount of Poisson noise.

<!-- GENERATED FROM PYTHON SOURCE LINES 38-68 -->
```Python
import torch
import deepinv as dinv

torch.manual_seed(0)

device = dinv.utils.get_device()
img_size = 128 if torch.cuda.is_available() else 64

psnr = dinv.metric.PSNR()
mae = dinv.metric.MAE()

gaussian_psf = dinv.physics.functional.gaussian_blur(
    psf_size=(17, 17),
    sigma=2.0,
    device=device,
)
kernel_init = torch.ones_like(gaussian_psf)
kernel_init = kernel_init / kernel_init.sum()
gain = 1 / 100

physics = dinv.physics.BlurFFT(
    img_size=(1, img_size, img_size),
    filter=gaussian_psf,
    device=device,
    noise_model=dinv.physics.PoissonNoise(
        gain=gain, normalize=True, clip_positive=True
    ),
)
```

```none
Selected GPU 0 with 24381.25 MiB free memory
```

<!-- GENERATED FROM PYTHON SOURCE LINES 69-77 -->

## Non-blind Richardson-Lucy on a grayscale image

In the non-blind setting, the kernel is known and the only unknown is the image.
In deepinv, the non-blind Richardson-Lucy algorithm corresponds to the
[`deepinv.optim.MLEM`](https://deepinv.org/api/stubs/deepinv.optim.MLEM.html.md#deepinv.optim.MLEM) used on measurements involving a blur physics.
It works both on grayscale and RGB images, but for simplicity we first consider
grayscale images.

<!-- GENERATED FROM PYTHON SOURCE LINES 77-125 -->
```Python
x_gray = dinv.utils.load_example(
    "SheppLogan.png",
    img_size=img_size,
    grayscale=True,
    resize_mode="resize",
    device=device,
)

y_gray = physics(x_gray)

data_fidelity = dinv.optim.PoissonLikelihood(gain=gain)
mlem = dinv.optim.MLEM(
    data_fidelity=data_fidelity,
    prior=None,
    max_iter=70 if torch.cuda.is_available() else 20,
    early_stop=True,
    thres_conv=1e-6,
    crit_conv="residual",
    verbose=True,
)

x_rl, metrics = mlem(
    y_gray,
    physics,
    x_gt=x_gray,
    compute_metrics=True,
)

dinv.utils.plot(
    {
        "Ground truth": x_gray,
        "Blurred": y_gray,
        "Richardson-Lucy": x_rl.clamp(0, 1),
    },
    subtitles=[
        "Reference",
        f"PSNR: {psnr(x_gray, y_gray).item():.2f} dB",
        f"PSNR: {psnr(x_gray, x_rl).item():.2f} dB",
    ],
    figsize=(9, 3),
    rescale_mode="clip",
    vmin=0,
    vmax=1,
)

dinv.utils.plot_curves(metrics)
```

<!-- GENERATED FROM PYTHON SOURCE LINES 126-158 -->

## Blind Richardson-Lucy on clean data

We now assume that both the clean image $x$ and the blur kernel
$h$ are unknown. The convolution is written $h * x$ when updating
the image and $x * h$ when updating the kernel.
The blind Richardson-Lucy algorithm simply alternates the Richardson-Lucy updates
for these two convolutions.
The updates are given by:

$$
h^{(k+1)} = \Pi_{\Delta}\left[\frac{h^{(k)}}{(x^{(k)})^\dagger * \mathbf{1}} \odot (x^{(k)})^\dagger * \left(\frac{y}{x^{(k)} * h^{(k)}}\right)\right],
$$

and

$$
x^{(k+1)} = \frac{x^{(k)}}{(h^{(k+1)})^\dagger * \mathbf{1}} \odot (h^{(k+1)})^\dagger * \left(\frac{y}{h^{(k+1)} * x^{(k)}}\right).
$$

Here, $z^\dagger$ denotes the spatially flipped $z$, such that $z^\dagger *$ is the adjoint of convolution by $z$.
The operation $\Pi_{\Delta}$ keeps the kernel nonnegative and
normalized to unit sum. The [`deepinv.optim.BlindRL`](https://deepinv.org/api/stubs/deepinv.optim.BlindRL.html.md#deepinv.optim.BlindRL) class implements
these alternating updates.

This algorithm is implemented under the [`deepinv.optim.BlindRL`](https://deepinv.org/api/stubs/deepinv.optim.BlindRL.html.md#deepinv.optim.BlindRL) class.
The number of iterations of each of the two updates can be controlled by the
parameters `x_steps` and `k_steps`.
Here, we set `use_fft=True` to compute the filter adjoint with FFTs.

For simplicity we first test this algorithm on a clean measurement, i.e. without
noise. In this case, the algorithm performs fairly well, with only some oscillation artefacts near the edges known as [Gibbs phenomenon](https://en.wikipedia.org/wiki/Gibbs_phenomenon).

<!-- GENERATED FROM PYTHON SOURCE LINES 158-204 -->
```Python
physics_clean = dinv.physics.Blur(
    filter=gaussian_psf,
    padding="circular",
    device=device,
)
y_clean = physics_clean(x_gray)

blindrl = dinv.optim.BlindRL(
    max_iter=100 if torch.cuda.is_available() else 40,
    init=(y_clean.clamp_min(1e-12), kernel_init),
    x_steps=1,
    k_steps=1,
    normalize_kernel=True,
    use_fft=True,
    verbose=True,
)
(x_blindrl, k_blindrl), metrics_clean_blind = blindrl(
    y_clean,
    x_gt=x_gray,
    compute_metrics=True,
)
x_blindrl = x_blindrl.clamp(0, 1)

dinv.utils.plot(
    {
        "Ground truth": x_gray,
        "Measurement": y_clean,
        "Blind RL": x_blindrl,
        "True kernel": gaussian_psf,
        "Estimated kernel": k_blindrl,
    },
    subtitles=[
        "Reference",
        f"PSNR: {psnr(x_gray, y_clean).item():.2f} dB",
        f"PSNR: {psnr(x_gray, x_blindrl).item():.2f} dB",
        "Reference",
        f"MAE: {mae(gaussian_psf, k_blindrl).item():.4f}",
    ],
    figsize=(13, 3),
    rescale_mode="clip",
    vmin=0,
    vmax=1,
)
dinv.utils.plot_curves(metrics_clean_blind)
```

<!-- GENERATED FROM PYTHON SOURCE LINES 205-216 -->

## Blind Richardson-Lucy on noisy data

Blind Richardson-Lucy is much less stable when the data is noisy.
Without regularization, the multiplicative updates tend to amplify noise and
the blur kernel tends to converge very fast to a Dirac.
The Dirac is unfortunately never a valid solution, because it corresponds to the
situation where our input image $y$ was already sharp.

Here we show an example with a medium amount of Poisson noise.
The algorithm further degrades the image and is unable to recover the correct kernel.

<!-- GENERATED FROM PYTHON SOURCE LINES 216-258 -->
```Python
y_noisy = physics(x_gray)

physics_clean.update_parameters(filter=gaussian_psf)
blindrl = dinv.optim.BlindRL(
    max_iter=80 if torch.cuda.is_available() else 40,
    init=(y_noisy.clamp_min(1e-12), kernel_init),
    x_steps=1,
    k_steps=1,
    normalize_kernel=True,
    use_fft=True,
    verbose=True,
)
(x_blindrl, k_blindrl), metrics = blindrl(
    y_noisy,
    x_gt=x_gray,
    compute_metrics=True,
)
x_blindrl = x_blindrl.clamp(0, 1)

dinv.utils.plot(
    {
        "Ground truth": x_gray,
        "Measurement": y_noisy,
        "Blind RL, noisy": x_blindrl,
        "True kernel": gaussian_psf,
        "Estimated kernel": k_blindrl,
    },
    subtitles=[
        "Reference",
        f"PSNR: {psnr(x_gray, y_noisy).item():.2f} dB",
        f"PSNR: {psnr(x_gray, x_blindrl).item():.2f} dB",
        "Reference",
        f"MAE: {mae(gaussian_psf, k_blindrl).item():.4f}",
    ],
    figsize=(13, 3),
    rescale_mode="clip",
    vmin=0,
    vmax=1,
)
dinv.utils.plot_curves(metrics)
```

<!-- GENERATED FROM PYTHON SOURCE LINES 259-287 -->

## One-Step-Late regularization

A standard heuristic to extend EM methods to the regularized setting is called
One-Step-Late (OSL) regularization <sup>[3](#footcite-greenuseemalgorithm1990)</sup>.
Since Richardson-Lucy is an instance of the EM algorithm, it can be extended to the
regularized setting using OSL.
For an image prior or regularizer $R_x$ and a kernel regularizer $R_h$, the denominators are modified using the current gradients:

$$
h^{(k+1)} = \Pi_{\Delta}\left[\frac{h^{(k)}}{(x^{(k)})^\dagger * \mathbf{1} + \lambda_h \nabla R_h(h^{(k)})} \odot (x^{(k)})^\dagger * \left(\frac{y}{x^{(k)} * h^{(k)}}\right)\right].
$$

$$
x^{(k+1)} = \frac{x^{(k)}}{(h^{(k+1)})^\dagger * \mathbf{1} + \lambda_x \nabla R_x(x^{(k)})} \odot (h^{(k+1)})^\dagger * \left(\frac{y}{h^{(k+1)} * x^{(k)}}\right),
$$

In the non-smooth case, the gradients are replaced by subgradients.
This enables the use of any prior implementing the [`deepinv.optim.prior.Prior`](https://deepinv.org/api/stubs/deepinv.optim.Prior.html.md#deepinv.optim.Prior)
interface inside the [`deepinv.optim.BlindRL`](https://deepinv.org/api/stubs/deepinv.optim.BlindRL.html.md#deepinv.optim.BlindRL) class.
Remember that the OSL method is heuristic and in particular is not guaranteed
to converge.
Still it is fairly robust in most cases and computationally efficient compared
to more sophisticated regularized algorithms.
The image and kernel regularizers can be specified independently using the `x_prior` and `k_prior` arguments of the [`deepinv.optim.BlindRL`](https://deepinv.org/api/stubs/deepinv.optim.BlindRL.html.md#deepinv.optim.BlindRL) class.
Below, we use TV regularization on the image through the [`deepinv.optim.TVPrior`](https://deepinv.org/api/stubs/deepinv.optim.TVPrior.html.md#deepinv.optim.TVPrior) class and no regularization on the kernel.

The image is a bit sharper and we get a better estimate of the kernel.

<!-- GENERATED FROM PYTHON SOURCE LINES 287-332 -->
```Python
tv_prior = dinv.optim.TVPrior(n_it_max=10)

blindrl = dinv.optim.BlindRL(
    x_prior=tv_prior,
    k_prior=None,
    lambda_reg_x=0.015,
    lambda_reg_k=0.0,
    max_iter=100 if torch.cuda.is_available() else 40,
    init=(y_noisy.clamp_min(1e-12), kernel_init),
    x_steps=1,
    k_steps=1,
    normalize_kernel=True,
    use_fft=True,
    verbose=True,
)
(x_blindrl, k_blindrl), metrics = blindrl(
    y_noisy,
    x_gt=x_gray,
    compute_metrics=True,
)
x_blindrl = x_blindrl.clamp(0, 1)

dinv.utils.plot(
    {
        "Ground truth": x_gray,
        "Measurement": y_noisy,
        "Blind RL + TV": x_blindrl,
        "True kernel": gaussian_psf,
        "Estimated kernel": k_blindrl,
    },
    subtitles=[
        "Reference",
        f"PSNR: {psnr(x_gray, y_noisy).item():.2f} dB",
        f"PSNR: {psnr(x_gray, x_blindrl).item():.2f} dB",
        "Reference",
        f"MAE: {mae(gaussian_psf, k_blindrl).item():.4f}",
    ],
    figsize=(13, 3),
    rescale_mode="clip",
    vmin=0,
    vmax=1,
)
dinv.utils.plot_curves(metrics)
```

<!-- GENERATED FROM PYTHON SOURCE LINES 333-340 -->

## Blind Richardson-Lucy on an RGB image

The same algorithm also works on RGB images. The blur kernel is shared across
channels, and the kernel update aggregates information from all channels.
Here we show a noiseless example, but in the noisy setting, regularization
should again be used to avoid noise amplification and kernel collapse.

<!-- GENERATED FROM PYTHON SOURCE LINES 340-390 -->
```Python
x_rgb = dinv.utils.load_example(
    "butterfly.png",
    img_size=img_size,
    grayscale=False,
    resize_mode="resize",
    device=device,
)

physics_clean.update_parameters(filter=gaussian_psf)
y_rgb = physics_clean(x_rgb)

blindrl = dinv.optim.BlindRL(
    max_iter=400 if torch.cuda.is_available() else 40,
    init=(y_rgb.clamp_min(1e-12), kernel_init),
    x_steps=1,
    k_steps=2,
    normalize_kernel=True,
    use_fft=True,
    verbose=True,
)
(x_blindrl, k_blindrl), metrics = blindrl(
    y_rgb,
    x_gt=x_rgb,
    compute_metrics=True,
)
x_blindrl = x_blindrl.clamp(0, 1)

dinv.utils.plot(
    {
        "Ground truth": x_rgb,
        "Measurement": y_rgb,
        "Blind RL": x_blindrl,
        "True kernel": gaussian_psf,
        "Estimated kernel": k_blindrl,
    },
    subtitles=[
        "Reference",
        f"PSNR: {psnr(x_rgb, y_rgb).item():.2f} dB",
        f"PSNR: {psnr(x_rgb, x_blindrl).item():.2f} dB",
        "Reference",
        f"MAE: {mae(gaussian_psf, k_blindrl).item():.4f}",
    ],
    figsize=(13, 3),
    rescale_mode="clip",
    vmin=0,
    vmax=1,
)
dinv.utils.plot_curves(metrics)
```

<!-- GENERATED FROM PYTHON SOURCE LINES 391-394 -->
* **References:**

* <a id='footcite-richardsonbayesianbasediterativemethod1972'>**[1]**</a> William Hadley Richardson. Bayesian-based iterative method of image restoration. *Journal of the Optical Society of America*, 62(1):55, 1972. [doi:10.1364/JOSA.62.000055](https://doi.org/10.1364/JOSA.62.000055).
* <a id='footcite-lucyiterativetechniquerectification1974'>**[2]**</a> L. B. Lucy. An iterative technique for the rectification of observed distributions. *The Astronomical Journal*, 79:745, 1974. [doi:10.1086/111605](https://doi.org/10.1086/111605).
* <a id='footcite-greenuseemalgorithm1990'>**[3]**</a> Peter J. Green. On Use of the Em Algorithm for Penalized Likelihood Estimation. *Journal of the Royal Statistical Society: Series B (Methodological)*, 52(3):443–452, 1990. [doi:10.1111/j.2517-6161.1990.tb01798.x](https://doi.org/10.1111/j.2517-6161.1990.tb01798.x).

**Total running time of the script:** (0 minutes 2.511 seconds)

<a id="sphx-glr-download-auto-examples-blind-inverse-problems-demo-blind-richardsonlucy-py"></a>

[`Download Jupyter notebook: demo_blind_richardsonlucy.ipynb`](https://deepinv.org/auto_examples/blind-inverse-problems/_downloads/a57a05ae19bc0fef4b23d3f4f78ccb9f/demo_blind_richardsonlucy.ipynb)

[`Download Python source code: demo_blind_richardsonlucy.py`](https://deepinv.org/auto_examples/blind-inverse-problems/_downloads/bbd2028ff436b13d54b47f7801f23bfc/demo_blind_richardsonlucy.py)

[`Download zipped: demo_blind_richardsonlucy.zip`](https://deepinv.org/auto_examples/blind-inverse-problems/_downloads/40c9c3cab20619b85d7d120dc25d1dc1/demo_blind_richardsonlucy.zip)

[Gallery generated by Sphinx-Gallery](https://sphinx-gallery.github.io)