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}$
This numerical tour explores the use of variational energies (Sobolev, total variation) to regularize the image inpaiting problem.
Here we consider inpainting of damaged observation without noise.
from __future__ import division
import nt_toolbox as nt
from nt_solutions import inverse_4_inpainting_variational as solutions
%matplotlib inline
%load_ext autoreload
%autoreload 2
Inpainting corresponds to filling holes in images.
First we load the image $f_0 \in \RR^N$ of $N=n\times n$ to be inpainted.
name = 'cameraman'
n = 256
f0 = rescale(load_image(name, n))
Display the original image.
imageplot(f0)
Ratio of removed pixels.
rho = .7
Construct a random mask $\Ga = \chi_{\Om}$ so that $\Ga_i=0$ for removed pixels $i \notin \Om$, and $\Ga_i=1$ for kept pixels $i \in \Om$.
Gamma = rand(n) >rho
We create the masking operator $\Phi$ which is a diagonal operator: $$ (\Phi f)_i = \Ga_i f_i $$
Phi = lambda f: f.*Gamma
Compute the damaged observation $y=\Phi(f_0)$ (no noise is added).
y = Phi(f0)
Display the observations.
imageplot(y)
We solve the inpainting problem by minimzing the Sobolev norm of the image under the constraint of matching the observation $$ f^\star = \uargmin{ \Phi(f) = y } E(f) = \norm{\nabla f}^2 $$ where $\nabla$ is a finite difference approximation of the gradient.
It can be shown that the solution to this problem is an harmonic function with prescribed boundary condition $$ \forall i \notin \Om, \quad (\Delta f^\star)_i=0 \qandq \forall i \in \Om, \quad f^\star_i = y_i. $$
This problem requires the constrained minimization of a smooth function, it can thus be solved using a projected gradient descent $$ f^{(\ell+1)} = \Pi \pa{ f^{(\ell)} + \tau \Delta(f^{(\ell)}) }$$ where $ \Pi $ is the orthogonal projector on the constraint $y=\Phi f$ $$ (\Pi f)_i = \choice{ y_i \qifq i \in \Om, \\ f_i \qifq i \notin \Om, \\ } $$
Pi = lambda f: f.*(1-Gamma) + y.*Gamma
Here $ \Delta = -\nabla^* \circ \nabla = \text{div} \circ \nabla $ is the gradient of the Sobolev energy $ E $.
Delta = lambda f: div(grad(f))
For convergence, the gradient descent step size should satisfy: $$ \tau<\frac{2}{\norm{\Delta}}=\frac{1}{4} $$
tau = .8/ 4
Exercise 1
Perform the projected gradient descent. Record in a variable |E| the evolution of the Sobolev energy $E$.
solutions.exo1()
## Insert your code here.
Display the decay of the energy $E(f^{(\ell)})$ with the iterations.
plot(E); axis('tight')
set_label('Iteration #', 'E')
Display the result.
imageplot(f, strcat(['Inpainted, SNR = ' num2str(snr(f0, f), 3) 'dB']))
A non-linear prior replaces the Sobolev energy by the TV norm, that tends to better reconstruct edges. Here we use a smoothed TV norm to avoid convergence issue with gradient descent algorithms.
The smoothed TV norm reads: $$ J_\epsilon(f) = \sum_x \sqrt{\norm{ \nabla f(x) }^2+\epsilon^2} $$
We use a projected gradient descent to solve this problem $$ f^{(\ell+1)} = \Pi \pa{ f^{(\ell)} + \tau G_\epsilon(f^{(\ell)}) }$$ where $ G_\epsilon $ is the gradient of $J_\epsilon$, that is defined as $$ G_\epsilon(f) = -\text{div} N_\epsilon( \nabla f ) $$ where $ N_\epsilon $ is the following normalization operator $$ N_\epsilon(u)_i = \frac{u_i}{ \sqrt{\norm{u_i}^2 + \epsilon^2} } $$ that is applied to any vector field $u=(u_i)_i \in \RR^{N \times 2} $ for $u_i \in \RR^2$.
Regularization parameter $\epsilon$ for the TV norm
epsilon = 1e-2
Define the normalization operator.
Amplitude = lambda u: sqrt(sum(u.^2, 3) + epsilon^2)
Neps = lambda u: u./ repmat(Amplitude(u), [1 1 2])
The step size $\tau$, should satisfy $$ \tau<\frac{\epsilon}{4}. $$
tau = .9*epsilon/ 4
Define the gradient of $J$
G = lambda f: -div(Neps(grad(f)))
Exercise 2
Perform the projected gradient descent. Record in a variable |J| the evolution of the TV energy $J_\epsilon$.
solutions.exo2()
## Insert your code here.
Display the result.
imageplot(clamp(f), strcat(['SNR = ' num2str(snr(f0, f), 3) 'dB']))
Display the evolution of the TV norm $J_\epsilon$.
plot(J)
axis('tight')
set_label('Iteration #', 'J_\epsilon')
Inpainting can be used to remove objects in pictures.
Load an image.
n = 256
f0 = load_image('parrot', n)
f0 = rescale(sum(f0, 3))
Display it.
imageplot(f0)
Load the mask.
Gamma = load_image('parrot-mask', n)
Gamma = double(rescale(Gamma) >.5)
Masking operator $\Phi$.
Phi = lambda f: f.*Gamma
Observation $y=\Phi(f_0)$.
y = Phi(f0)
Display it.
imageplot(y)
Exercise 3
Perform Sobolev inpainting on this image.
solutions.exo3()
## Insert your code here.
Exercise 4
Try other methods to solve this inpainting problem. You can for instance have a look on the numerical on sparsity for deconvolution and inpainting.
solutions.exo4()
## Insert your code here.