Image Deconvolution using Variational Method¶
This numerical tour explores the use of Sobolev and TV regularization to perform image deconvolution.
Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
options(warn=-1) # turns off warnings, to turn on: "options(warn=0)"
library(pracma)
library(imager)
# Importing the libraries
for (f in list.files(path="nt_toolbox/toolbox_general/", pattern="*.R")) {
source(paste("nt_toolbox/toolbox_general/", f, sep=""))
}
for (f in list.files(path="nt_toolbox/toolbox_signal/", pattern="*.R")) {
source(paste("nt_toolbox/toolbox_signal/", f, sep=""))
}
This tour is concerned with the deconvolution problem. The measurement are assumed to be blurry and noisy: $$y=\Phi f_0 + w = h \star f_0 + w$$
Where here |h| is the filter (low pass) and |w| some noise (here assumed to be white Gaussian).
We consider variational deconvolution methods, that finds a regularizer through a convex optimization: $$f^\star \in \text{argmin}_f \frac{1}{2}\|y-\Phi f\|^2 + \lambda J(f)$$
Where $J(f)$ is a prior energy. In this tour we consider a simple L2 prior (the image is assumed to have a bounded energy), a Sobolev prior (the image is uniformly smooth) and an approximate total variation (the image has edges of bounded perimeter).
Note that the parameter $\lambda$ should be carefully chosen to fit the noise level.
Image Blurring¶
Deconvolution corresponds to removing a blur from an image. We use here a Gaussian blur.
First we load the image to be inpainted.
n = 256
name = 'nt_toolbox/data/hibiscus.png'
f0 = load_image(name, n)
Initial image $f_0$.
options(repr.plot.width=4, repr.plot.height=4)
imageplot(f0, 'Image f0')
We build a convolution kernel. Since we are going to use Fourier to compute the convolution, we set the center of the kernel in the (1,1) pixel location.
Width $s$ of the kernel, in pixel.
s = 3
Define the convolution kernel $h$.
x = c(c(0:((n/2) - 1)), c(-(n/2):-1))
mesh = meshgrid(x)
X = mesh$X
Y = mesh$Y
h = exp( -(X**2 + Y**2)/(2.0 * s**2) )
h = h/sum(h)
Useful for later : the Fourier transform (should be real because of symmetry).
hF = Re(fft(h))
Display the kernel and its transform. We use |fftshift| to center the filter for display.
imageplot(fftshift(h), 'Filter', c(1,2,1))
imageplot(fftshift(hF), 'Fourier transform', c(1,2,2))
We use this short hand for the filtering. Note that this is a symmetric operator.
Phi = function(x,h){Re(fft(fft(x) * fft(h), inverse=TRUE) / length(h))}
Apply the filter.
y0 = Phi(f0[,], h)
Display the filtered observation.
imageplot(y0, 'Observation without noise')
Variance $\sigma^2$ of the noise $w$.
sigma = .02
Add some noise to obtain the measurements $y = \Phi f_0 + w$.
set.seed(123)
y = y0 + randn(n,n) * sigma
Display.
imageplot(y, 'Observation with noise')
# Normalize to 0-1
imageplot(y, paste('Observation, SNR = ',round(snr(f0, y),2),'dB'))
Deconvolution with L2 Regularization¶
Deconvolution is obtained by dividing the Fourier transform of $y$ by $\hat h$. $$f^\star(\omega) = \frac{\hat y(\omega)}{\hat h(\omega)} = \hat f_0(\omega) + \hat w(\omega)/{\hat h(\omega)}$$
Unfortunately this creates an explosion of the Noise.
To avoid this explosion, we consider a simple regularization. $$f^{\star} = \text{argmin}_f \: \|y-\Phi f\|^2 + \lambda \|f\|^2$$
Since the filtering is diagonalized over Fourier, the solution is simply computed over the Fourier domain as: $$\hat f^\star(\omega) = \frac{\hat y(\omega) \hat h(\omega)}{ \|\hat h(\omega)\|^2 + \lambda }$$
Useful for later: Fourier transform of the observations.
yF = fft(y)
Select a value for the regularization parameter.
Lambda = 0.02
Perform the inversion.
fL2 = Re(fft(yF * hF / (abs(hF)**2 + Lambda), inverse=TRUE)) / length(yF * hF / (abs(hF)**2 + Lambda))
Display.
imageplot(fL2, paste('L2 deconvolution, SNR = ', round(snr(f0, fL2),2), 'dB'))