#### 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)..
# Image deblurring with Total-Variation (TV) prior
This example shows how to use a standard TV prior for image deblurring. The problem writes as $y = Ax + \epsilon$
where $A$ is a convolutional operator and $\epsilon$ is the realization of some Gaussian noise. The goal is
to recover the original image $x$ from the blurred and noisy image $y$. The TV prior is used to regularize
the problem.
```Python
import deepinv as dinv
from pathlib import Path
import torch
from torchvision import transforms
from deepinv.optim.data_fidelity import L2
from deepinv.optim import PGD
from deepinv.utils.demo import load_dataset, load_degradation
from deepinv.utils.plotting import plot, plot_curves
```
## Setup paths for data loading and results.
```Python
BASE_DIR = Path(".")
DATA_DIR = BASE_DIR / "measurements"
RESULTS_DIR = BASE_DIR / "results"
DEG_DIR = BASE_DIR / "degradations"
```
## Load base image datasets and degradation operators.
In this example, we use the Set3C dataset and a motion blur kernel from Levin *et al.*[1](#footcite-levin2009understanding).
```Python
# Set the global random seed from pytorch to ensure reproducibility of the example.
torch.manual_seed(0)
device = dinv.utils.get_device()
# Set up the variable to fetch dataset and operators.
dataset_name = "set3c"
img_size = 256 if torch.cuda.is_available() else 64
val_transform = transforms.Compose(
[transforms.CenterCrop(img_size), transforms.ToTensor()]
)
# Generate a motion blur operator.
kernel_index = 1 # which kernel to chose among the 8 motion kernels from 'Levin09.mat'
kernel_torch = load_degradation("Levin09.npy", DEG_DIR / "kernels", index=kernel_index)
kernel_torch = kernel_torch.unsqueeze(0).unsqueeze(
0
) # add batch and channel dimensions
dataset = load_dataset(dataset_name, transform=val_transform)
```
```none
Selected GPU 0 with 28321.25 MiB free memory
Levin09.npy degradation downloaded in degradations/kernels
/local/jtachell/deepinv/deepinv/deepinv/datasets/base.py:495: DeprecationWarning: The tuple format for dataset outputs is deprecated and will be removed in a future version.It is recommended to set `use_dict_output=True` for better readability and flexibility in returned outputs.The default is currently `False` for backward compatibility, but will be switched to `True` in a future version.
super().__init__(use_dict_output=use_dict_output)
```
## Generate a dataset of blurred images and load it.
We use the BlurFFT class from the physics module to generate a dataset of blurred images.
```Python
noise_level_img = 0.05 # Gaussian Noise standard deviation for the degradation
n_channels = 3 # 3 for color images, 1 for gray-scale images
physics = dinv.physics.BlurFFT(
img_size=(n_channels, img_size, img_size),
filter=kernel_torch,
device=device,
noise_model=dinv.physics.GaussianNoise(sigma=noise_level_img),
)
# Select the first image from the dataset
x = dataset[0].unsqueeze(0).to(device)
# Apply the degradation to the image
y = physics(x)
```
## Exploring the total variation prior.
In this example, we will use the total variation prior, which can be done with the [`deepinv.optim.prior.Prior`](https://deepinv.org/api/stubs/deepinv.optim.Prior.html.md#deepinv.optim.Prior)
class. The prior object represents the cost function of the prior (TV in this case), as well as convenient methods,
such as its proximal operator $\text{prox}_{\tau g}$.
```Python
# Set up the total variation prior
prior = dinv.optim.prior.TVPrior(n_it_max=2000)
# Compute the total variation prior cost
cost_tv = prior(y).item()
print(f"Cost TV: g(y) = {cost_tv:.2f}")
# Apply the proximal operator of the TV prior
x_tv = prior.prox(y, gamma=0.1)
cost_tv_prox = prior(x_tv).item()
```
```none
Cost TV: g(y) = 19505.50
```
#### NOTE
The output of the proximity operator of TV is **not** the solution to our deblurring problem. It is only a
step towards the solution and is used in the proximal gradient descent algorithm to solve the inverse
problem.
```Python
# Plot the input and the output of the TV proximal operator
imgs = [y, x_tv]
plot(
{"Input": y, "Output": x_tv},
subtitles=[
f"TV cost: {int(cost_tv)}",
f"TV cost: {int(cost_tv_prox)}",
],
tight=False,
)
```
## Set up the optimization algorithm to solve the inverse problem.
The problem we want to minimize is the following:
$$
\begin{equation*}
\underset{x}{\operatorname{min}} \,\, \frac{1}{2} \|Ax-y\|_2^2 + \lambda \|Dx\|_{1,2}(x),
\end{equation*}
$$
where $1/2 \|A(x)-y\|_2^2$ is the a data-fidelity term, $\lambda \|Dx\|_{2,1}(x)$ is the total variation (TV)
norm of the image $x$, and $\lambda>0$ is a regularisation parameters.
We use a Proximal Gradient Descent (PGD) algorithm to solve the inverse problem.
```Python
# Select the data fidelity term
data_fidelity = L2()
# Specify the prior (we redefine it with a smaller number of iteration for faster computation)
prior = dinv.optim.prior.TVPrior(n_it_max=20)
# Logging parameters
verbose = True
plot_convergence_metrics = (
True # compute performance and convergence metrics along the algorithm.
)
# Algorithm parameters
stepsize = 1.0
lambda_reg = 1e-2 # TV regularisation parameter
max_iter = 300
early_stop = True
# Instantiate the algorithm class to solve the problem.
model = PGD(
prior=prior,
data_fidelity=data_fidelity,
stepsize=stepsize,
lambda_reg=lambda_reg,
early_stop=early_stop,
max_iter=max_iter,
verbose=verbose,
)
```
## Evaluate the model on the problem and plot the results.
The model returns the output and the metrics computed along the iterations.
For computing PSNR, the ground truth image `x_gt` must be provided.
```Python
x_lin = physics.A_adjoint(y) # linear reconstruction with the adjoint operator
# run the model on the problem.
x_model, metrics = model(
y, physics, x_gt=x, compute_metrics=True
) # reconstruction with PGD algorithm
# compute PSNR
psnr_input = dinv.metric.PSNR()(x, y)
psnr_lin = dinv.metric.PSNR()(x, x_lin)
psnr_model = dinv.metric.PSNR()(x, x_model)
# plot images. Images are saved in RESULTS_DIR.
imgs = [x, y, x_lin, x_model]
plot(
{"Ground Truth": x, "Input": y, "Linear Recon": x_lin, "Recons": x_model},
subtitles=[
"PSNR:",
f"{psnr_input.item():.2f} dB",
f"{psnr_lin.item():.2f} dB",
f"{psnr_model.item():.2f} dB",
],
)
# plot convergence curves
if plot_convergence_metrics:
plot_curves(metrics)
```
* **References:**
* **[1]** Anat Levin, Yair Weiss, Fredo Durand, and William T Freeman. Understanding and evaluating blind deconvolution algorithms. In *2009 IEEE conference on computer vision and pattern recognition*, 1964–1971. IEEE, 2009.
**Total running time of the script:** (0 minutes 3.719 seconds)
🏷 Tags: Deblurring
[`Download Jupyter notebook: demo_TV_minimisation.ipynb`](https://deepinv.org/auto_examples/optimization/_downloads/62d8fccb1dc20c16ea9832aeed2cdab6/demo_TV_minimisation.ipynb)
[`Download Python source code: demo_TV_minimisation.py`](https://deepinv.org/auto_examples/optimization/_downloads/d6d42a478eddf96e602235dd9d0410a9/demo_TV_minimisation.py)
[`Download zipped: demo_TV_minimisation.zip`](https://deepinv.org/auto_examples/optimization/_downloads/fe749e3d14ab646fb824bdec45530b59/demo_TV_minimisation.zip)
[Gallery generated by Sphinx-Gallery](https://sphinx-gallery.github.io)