Expected Patch Log Likelihood (EPLL) for Denoising and Inpainting

In this example we use the expected patch log likelihood (EPLL) prior EPLL proposed in “From learning models of natural image patches to whole image restoration”. for denoising and inpainting of natural images. To this end, we consider the inverse problem \(y = Ax+\epsilon\), where \(A\) is either the identity (for denoising) or a masking operator (for inpainting) and \(\epsilon\sim\mathcal{N}(0,\sigma^2 I)\) is white Gaussian noise with standard deviation \(\sigma\).

import torch
from deepinv.optim import EPLL
from deepinv.physics import GaussianNoise, Denoising, Inpainting
from deepinv.loss.metric import PSNR
from deepinv.utils import plot
from deepinv.utils.demo import load_url_image, get_image_url

device = "cuda" if torch.cuda.is_available() else "cpu"

Load test image and model

As default EPLL loads pretrained weights for the Gaussian mixture model which where estimted based on 50 mio patches extracted from the BSDS500 dataset. An example how to estimate the parameters of GMM is included in the demo for limited-angle CT with patch priors.

url = get_image_url("CBSD_0010.png")
test_img = load_url_image(url, grayscale=False).to(device)
patch_size = 6
model = EPLL(channels=test_img.shape[1], patch_size=patch_size, device=device)

Downloading: "https://huggingface.co/deepinv/EPLL/resolve/main/GMM_BSDS_color2.pt?download=true" to /home/runner/.cache/torch/hub/checkpoints/GMM_BSDS_color2.pt

  0%|          | 0.00/35.7M [00:00<?, ?B/s]
  3%|▎         | 1.12M/35.7M [00:00<00:03, 11.1MB/s]
  6%|▋         | 2.25M/35.7M [00:00<00:03, 11.4MB/s]
  9%|▉         | 3.38M/35.7M [00:00<00:03, 10.4MB/s]
 13%|█▎        | 4.50M/35.7M [00:00<00:03, 10.9MB/s]
 16%|█▌        | 5.62M/35.7M [00:00<00:03, 10.3MB/s]
 19%|█▉        | 6.75M/35.7M [00:00<00:02, 10.7MB/s]
 22%|██▏       | 7.88M/35.7M [00:00<00:02, 10.2MB/s]
 25%|██▌       | 9.00M/35.7M [00:00<00:02, 10.6MB/s]
 28%|██▊       | 10.1M/35.7M [00:01<00:02, 10.3MB/s]
 31%|███       | 11.1M/35.7M [00:01<00:02, 10.3MB/s]
 34%|███▍      | 12.1M/35.7M [00:01<00:02, 10.4MB/s]
 37%|███▋      | 13.1M/35.7M [00:01<00:02, 10.4MB/s]
 40%|███▉      | 14.1M/35.7M [00:01<00:02, 10.4MB/s]
 42%|████▏     | 15.1M/35.7M [00:01<00:02, 10.4MB/s]
 46%|████▌     | 16.2M/35.7M [00:01<00:01, 10.8MB/s]
 49%|████▊     | 17.4M/35.7M [00:01<00:01, 10.3MB/s]
 52%|█████▏    | 18.5M/35.7M [00:01<00:01, 10.7MB/s]
 55%|█████▍    | 19.6M/35.7M [00:01<00:01, 10.3MB/s]
 58%|█████▊    | 20.6M/35.7M [00:02<00:01, 10.3MB/s]
 61%|██████    | 21.6M/35.7M [00:02<00:01, 10.4MB/s]
 63%|██████▎   | 22.6M/35.7M [00:02<00:01, 10.4MB/s]
 67%|██████▋   | 23.8M/35.7M [00:02<00:01, 10.8MB/s]
 70%|██████▉   | 24.9M/35.7M [00:02<00:01, 10.3MB/s]
 73%|███████▎  | 25.9M/35.7M [00:02<00:00, 10.4MB/s]
 75%|███████▌  | 26.9M/35.7M [00:02<00:00, 10.4MB/s]
 78%|███████▊  | 27.9M/35.7M [00:02<00:00, 10.4MB/s]
 81%|████████▏ | 29.0M/35.7M [00:02<00:00, 10.8MB/s]
 84%|████████▍ | 30.1M/35.7M [00:03<00:00, 10.3MB/s]
 87%|████████▋ | 31.1M/35.7M [00:03<00:00, 10.3MB/s]
 90%|█████████ | 32.1M/35.7M [00:03<00:00, 10.3MB/s]
 93%|█████████▎| 33.1M/35.7M [00:03<00:00, 10.3MB/s]
 96%|█████████▌| 34.1M/35.7M [00:03<00:00, 10.4MB/s]
 98%|█████████▊| 35.1M/35.7M [00:03<00:00, 10.4MB/s]
100%|██████████| 35.7M/35.7M [00:03<00:00, 10.5MB/s]

Denoising

In this setting, the operator \(A\) is the identity; we set the noise level to \(\sigma=25/255\). Define noise model, operator and generate observation

sigma = 0.1
noise_model = GaussianNoise(sigma)
physics = Denoising(device=device, noise_model=noise_model)
observation = physics(test_img)

We use the default choice of the betas in the half quadratic splitting given by \(\beta \in \sigma^{-2} \{1,4,8,16,32\}\). Generally, the betas are hyperparameters, which have to be chosen for each inverse problem separately.

# Reconstruction
with torch.no_grad():
    x_out = model(observation, physics, batch_size=5000)

# PSNR computation and plots.
psnr_obs = PSNR()(observation, test_img).item()
psnr_recon = PSNR()(x_out, test_img).item()

print("PSNRs for Denoising:")
print("Observation: {0:.2f}".format(psnr_obs))
print("EPLL: {0:.2f}".format(psnr_recon))

plot(
    [test_img, observation.clip(0, 1), x_out.clip(0, 1)],
    ["Ground truth", "Observation", "EPLL"],
)

PSNRs for Denoising:
Observation: 20.00
EPLL: 28.96

Inpainting

We now turn to a noisy inpainting problem, where the operator \(A\) is a masking operator. Define noise model, operator and generate observation

sigma = 0.01
physics = Inpainting(
    tensor_size=test_img[0].shape,
    mask=0.7,
    device=device,
    noise_model=GaussianNoise(sigma),
)
observation = physics(test_img)

Here, we need a different choice of beta. To this end, we extended the default choice by two values and optimized a constant factor via grid search.

betas = [1.0, 5.0, 10.0, 40.0, 80.0, 160.0, 320.0]

# Reconstruction
with torch.no_grad():
    x_out = model(observation, physics, betas=betas, batch_size=5000)

# PSNR computation and plots
psnr_obs = PSNR()(observation, test_img).item()
psnr_recon = PSNR()(x_out, test_img).item()

print("PSNRs for Inpainting:")
print("Observation: {0:.2f}".format(psnr_obs))
print("EPLL: {0:.2f}".format(psnr_recon))

plot(
    [test_img, observation.clip(0, 1), x_out.clip(0, 1)],
    ["Ground truth", "Observation", "EPLL"],
)

PSNRs for Inpainting:
Observation: 10.94
EPLL: 29.80