Deconvolution of Lena

In this example, we deconvolve a noisy version of Lena using Wiener and unsupervised Wiener algorithms. This algorithms are based on linear models that can’t restore sharp edge as much as non-linear methods (like TV restoration) but are much faster.

Wiener filter

The inverse filter based on the PSF (Point Spread Function), the prior regularisation (penalisation of high frequency) and the tradeoff between the data and prior adequacy. The regularization parameter must be hand tuned.

Unsupervised Wiener

This algorithm has a self-tuned regularisation parameters based on data learning. This is not common and based on the following publication. The algorithm is based on a iterative Gibbs sampler that draw alternatively samples of posterior conditionnal law of the image, the noise power and the image frequency power.

[1]François Orieux, Jean-François Giovannelli, and Thomas Rodet, “Bayesian estimation of regularization and point spread function parameters for Wiener-Hunt deconvolution”, J. Opt. Soc. Am. A 27, 1593-1607 (2010)
../_images/plot_restoration_1.png

import numpy as np
import matplotlib.pyplot as plt

from skimage import color, data, restoration

lena = color.rgb2gray(data.lena())
from scipy.signal import convolve2d as conv2
psf = np.ones((5, 5)) / 25
lena = conv2(lena, psf, 'same')
lena += 0.1 * lena.std() * np.random.standard_normal(lena.shape)

deconvolved, _ = restoration.unsupervised_wiener(lena, psf)

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 5))

plt.gray()

ax[0].imshow(lena, vmin=deconvolved.min(), vmax=deconvolved.max())
ax[0].axis('off')
ax[0].set_title('Data')

ax[1].imshow(deconvolved)
ax[1].axis('off')
ax[1].set_title('Self tuned restoration')

fig.subplots_adjust(wspace=0.02, hspace=0.2,
                    top=0.9, bottom=0.05, left=0, right=1)

plt.show()

STDOUT


        

STDERR


        

Python source code: download (generated using skimage 0.11dev)

IPython Notebook: download (generated using skimage 0.11dev)

IyAtKi0gY29kaW5nOiB1dGYtOCAtKi0=
aW1wb3J0IG51bXB5IGFzIG5wCmltcG9ydCBtYXRwbG90bGliLnB5cGxvdCBhcyBwbHQKCmZyb20gc2tpbWFnZSBpbXBvcnQgY29sb3IsIGRhdGEsIHJlc3RvcmF0aW9uCgpsZW5hID0gY29sb3IucmdiMmdyYXkoZGF0YS5sZW5hKCkpCmZyb20gc2NpcHkuc2lnbmFsIGltcG9ydCBjb252b2x2ZTJkIGFzIGNvbnYyCnBzZiA9IG5wLm9uZXMoKDUsIDUpKSAvIDI1CmxlbmEgPSBjb252MihsZW5hLCBwc2YsICdzYW1lJykKbGVuYSArPSAwLjEgKiBsZW5hLnN0ZCgpICogbnAucmFuZG9tLnN0YW5kYXJkX25vcm1hbChsZW5hLnNoYXBlKQoKZGVjb252b2x2ZWQsIF8gPSByZXN0b3JhdGlvbi51bnN1cGVydmlzZWRfd2llbmVyKGxlbmEsIHBzZikKCmZpZywgYXggPSBwbHQuc3VicGxvdHMobnJvd3M9MSwgbmNvbHM9MiwgZmlnc2l6ZT0oOCwgNSkpCgpwbHQuZ3JheSgpCgpheFswXS5pbXNob3cobGVuYSwgdm1pbj1kZWNvbnZvbHZlZC5taW4oKSwgdm1heD1kZWNvbnZvbHZlZC5tYXgoKSkKYXhbMF0uYXhpcygnb2ZmJykKYXhbMF0uc2V0X3RpdGxlKCdEYXRhJykKCmF4WzFdLmltc2hvdyhkZWNvbnZvbHZlZCkKYXhbMV0uYXhpcygnb2ZmJykKYXhbMV0uc2V0X3RpdGxlKCdTZWxmIHR1bmVkIHJlc3RvcmF0aW9uJykKCmZpZy5zdWJwbG90c19hZGp1c3Qod3NwYWNlPTAuMDIsIGhzcGFjZT0wLjIsCiAgICAgICAgICAgICAgICAgICAgdG9wPTAuOSwgYm90dG9tPTAuMDUsIGxlZnQ9MCwgcmlnaHQ9MSkKCnBsdC5zaG93KCk=