Numerical Tours
← All Python tours
Download notebook Open in Colab

Sparse Image Inpainting

Reconstruct missing pixels by promoting sparsity in a wavelet representation. Alternate consistency with the observed pixels and coefficient shrinkage, then inspect how the sampling pattern affects recovery.

Run this tour

Run the cells in order with a Python 3 kernel. The first cell locates the companion data and toolbox and installs missing dependencies when needed. All worked examples include their implementation directly in this notebook. Random seeds make comparisons reproducible; you can change them to explore other samples.

# Locate the companion toolbox locally, or fetch it for a standalone/Colab copy.
from pathlib import Path
import importlib.util
import os
import subprocess
import sys

working = Path.cwd()
candidates = [working, working / "python", working.parent / "python"]
python_dir = next((p for p in candidates if (p / "nt_toolbox").is_dir()), None)
if python_dir is None:
    checkout = working / "numerical-tours-support"
    if not checkout.exists():
        subprocess.run(
            [
                "git",
                "clone",
                "--depth",
                "1",
                "--branch",
                "master",
                "https://github.com/gpeyre/numerical-tours.git",
                str(checkout),
            ],
            check=True,
        )
    python_dir = checkout / "python"
os.chdir(python_dir)
if str(python_dir) not in sys.path:
    sys.path.insert(0, str(python_dir))
requirements = python_dir / "requirements.txt"
if any(
    importlib.util.find_spec(name) is None
    for name in [
        "numpy",
        "scipy",
        "matplotlib",
        "skimage",
        "sklearn",
        "pywt",
        "ipywidgets",
        "cvxpy",
        "skfmm",
        "autograd",
        "progressbar",
        "celer",
    ]
):
    subprocess.run(
        [sys.executable, "-m", "pip", "install", "-r", str(requirements)], check=True
    )

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
plt.rcParams.update(
    {
        "figure.figsize": (8, 4),
        "figure.dpi": 100,
        "axes.spines.top": False,
        "axes.spines.right": False,
        "font.size": 11,
        "image.cmap": "gray",
    }
)
%matplotlib inline

$\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{\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{\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 sparse energies to regularize the image inpaiting problem.

import numpy as np
import scipy as scp
import pylab as pyl
import matplotlib.pyplot as plt

from nt_toolbox.general import clamp, np, plt, pylab, rescale
from nt_toolbox.signal import imageplot, load_image, np, plt, pylab, snr

import warnings

%matplotlib inline

Here we consider inpainting of damaged observation without noise.

Sparse Regularization

This tour consider measurements $y=\Phi f_0 + w$ where $\Phi$ is a masking operator and $w$ is an additive noise.

This tour is focused on using sparsity to recover an image from the measurements $y$. It considers a synthesis-based regularization, that compute a sparse set of coefficients $ (a_m^{\star})_m $ in a frame $\Psi = (\psi_m)_m$ that solves $$a^{\star} \in \text{argmin}_a \: \frac{1}{2}\|y-\Phi \Psi a\|^2 + \lambda J(a)$$

where $\lambda$ should be adapted to the noise level $\|w\|$. Since in this tour we consider damaged observation without noise, i.e. $w=0$, we use either a very small value of $\lambda$, or we decay its value through the iterations of the recovery process.

Here we use the notation $$\Psi a = \sum_m a_m \psi_m$$ to indicate the reconstruction operator, and $J(a)$ is the $\ell^1$ sparsity prior $$J(a)=\sum_m \|a_m\|.$$

Missing Pixels and Inpainting

Inpainting corresponds to filling holes in images. This corresponds to a linear ill posed inverse problem.

You might want to do first the numerical tour Variational image inpaiting that use Sobolev and TV priors to performs the inpainting.

First we load the image to be inpainted.

n = 128
f0 = load_image("nt_toolbox/data/lena.bmp")
f0 = rescale(f0[256 - n // 2 : 256 + n // 2, 256 - n // 2 : 256 + n // 2])

Display it.

plt.figure(figsize=(6, 6))
imageplot(f0, "Image f_0")
No description has been provided for this image

Amount of removed pixels.

rho = 0.7

Then we construct a mask $\Omega$ made of random pixel locations.

from numpy import random

Omega = np.zeros([n, n])
sel = random.permutation(n**2)
np.ravel(Omega)[sel[np.arange(int(rho * n**2))]] = 1

The damaging operator put to zeros the pixel locations $x$ for which $\Omega(x)=1$

Phi = lambda f, Omega: f * (1 - Omega)

The damaged observations reads $y = \Phi f_0$.

y = Phi(f0, Omega)

Display the observations.

plt.figure(figsize=(6, 6))
imageplot(y, "Observations y")
No description has been provided for this image

Soft Thresholding in a Basis

The soft thresholding operator is at the heart of $\ell^1$ minimization schemes. It can be applied to coefficients $a$, or to an image $f$ in an ortho-basis.

The soft thresholding is a 1-D functional that shrinks the value of coefficients. $$ s_T(u)=\max(0,1-T/|u|)u $$

Define a shortcut for this soft thresholding 1-D functional.

SoftThresh = lambda x, T: (
    x
    * np.maximum(
        1 - T / np.maximum(abs(x), 1e-10 * np.ones(np.shape(x))), np.zeros(np.shape(x))
    )
)

Display a curve of the 1D soft thresholding.

x = np.linspace(-1, 1, 1000)

plt.figure(figsize=(7, 5))
plt.plot(x, SoftThresh(x, 0.5))
plt.show()
No description has been provided for this image

Note that the function SoftThresh can also be applied to vector which defines an operator on coefficients: $$ S_T(a) = ( s_T(a_m) )_m. $$

In the next section, we use an orthogonal wavelet basis $\Psi$.

We set the parameters of the wavelet transform.

Jmax = np.log2(n) - 1
Jmin = Jmax - 3

Shortcut for $\Psi$ and $\Psi^*$ in the orthogonal case.

from nt_toolbox.perform_wavelet_transf import np, perform_wavelet_transf

Psi = lambda a: perform_wavelet_transf(a, Jmin, -1, ti=0)
PsiS = lambda f: perform_wavelet_transf(f, Jmin, +1, ti=0)

The soft thresholding opterator in the basis $\Psi$ is defined as $$S_T^\Psi(f) = \sum_m s_T( \langle f,\psi_m \rangle ) \psi_m $$

It thus corresponds to applying the transform $\Psi^*$, thresholding the coefficients using $S_T$ and then undoing the transform using $\Psi$. $$ S_T^\Psi(f) = \Psi \circ S_T \circ \Psi^*$$

SoftThreshPsi = lambda f, T: Psi(SoftThresh(PsiS(f), T))

This soft thresholding corresponds to a denoising operator.

plt.figure(figsize=(6, 6))
imageplot(clamp(SoftThreshPsi(f0, 0.1)))
No description has been provided for this image

Inpainting using Orthogonal Wavelet Sparsity

If $\Psi$ is an orthogonal basis, a change of variable shows that the synthesis prior is also an analysis prior, that reads $$f^{\star} \in \text{argmin}_f \: E(f) = \frac{1}{2}\|y-\Phi f\|^2 + \lambda \sum_m \|\langle f,\psi_m \rangle\|. $$

To solve this non-smooth optimization problem, one can use forward-backward splitting, also known as iterative soft thresholding.

It computes a series of images $f^{(\ell)}$ defined as $$ f^{(\ell+1)} = S_{\tau\lambda}^{\Psi}( f^{(\ell)} - \tau \Phi^{*} (\Phi f^{(\ell)} - y) ) $$

Set up the value of the threshold.

lambd = 0.03

In our setting, we have $ \Phi^* = \Phi $ which is an operator of norm 1.

For $f^{(\ell)}$ to converge to a solution of the problem, the gradient step size should be chosen as $$\tau < \frac{2}{\|\Phi^* \Phi\|} = 2$$

In the following we use: $$\tau = 1$$

Since we use $ \tau=1 $ and $ \Phi = \Phi^* = \text{diag}(1-\Omega) $, the gradient descent step is a projection on the inpainting constraint $$ C = \{ f \backslash \forall \Omega(x)=0, f(x)=y(x) \} $$ One thus has $$ f - \tau \Phi^{*} (\Phi f - y) = \text{Proj}_C(f) $$

For the sake of simplicity, we define a shortcut for this projection operator.

ProjC = lambda f, Omega: Omega * f + (1 - Omega) * y

Each iteration of the forward-backward (iterative thresholding) algorithm thus reads: $$ f^{(\ell+1)} = S_{\lambda}^\Psi( \text{Proj}_C(f^{(\ell)}) ). $$

Initialize the iterations.

fSpars = y

First step: gradient descent.

fSpars = ProjC(fSpars, Omega)

Second step: denoise the solution by thresholding.

fSpars = SoftThreshPsi(fSpars, lambd)

Worked example 1

Perform the iterative soft thresholding. Monitor the decay of the energy $E$ you are minimizing.

from numpy import linalg

fSpars = y
energy = []
niter = 1000
for i in range(niter):
    fSpars = SoftThreshPsi(ProjC(fSpars, Omega), lambd)
    # record the energy
    fW = PsiS(fSpars)
    energy = energy + [
        1 / 2 * linalg.norm(y - Phi(fSpars, Omega), "fro") ** 2
        + lambd * np.sum(abs(fW))
    ]

plt.figure(figsize=(7, 5))
plt.plot(energy, linewidth=2)
plt.xlabel("Iteration")
plt.ylabel("E")
plt.show()
No description has been provided for this image

Display the result.

plt.figure(figsize=(6, 6))
imageplot(clamp(fSpars))
No description has been provided for this image

Worked example 2

Since there is no noise, one should in theory take $\lambda \rightarrow 0$. To do this, decay the value of $\lambda$ through the iterations.

niter = 1000
lambda_list = np.linspace(0.03, 0, niter)
err = []

for i in range(niter):
    fSpars = SoftThreshPsi(ProjC(fSpars, Omega), lambda_list[i])

plt.figure(figsize=(6, 6))
imageplot(clamp(fSpars), "Sparsity inpainting, SNR = %.1f dB" % snr(f0, fSpars))
No description has been provided for this image

Inpainting using Translation Invariant Wavelet Sparsity

Orthogonal sparsity performs a poor regularization because of the lack of translation invariance. This regularization is enhanced by considering $\Psi$ as a redundant tight frame of translation invariant wavelets.

One thus looks for optimal coefficients $a^\star$ that solves $$a^{\star} \in \text{argmin}_a \: E(a) = \frac{1}{2}\|y-\Phi \Psi a\|^2 + \lambda J(a)$$

Important: The operator $\Psi^*$ is the forward translation invariant wavelet transform. It computes the inner product with the unit norm wavelet atoms: $$ (\Psi^* f)_m = \langle f,\psi_m \rangle \quad \text{with} \quad \|\psi_m\|=1. $$

The reconstruction operator $\Xi$ satisfies $ \Xi \Psi^* f = f $, and is the pseudo inverse of the analysis operator $ \Xi = (\Psi^*)^+ $.

For our algorithm, we will need to use $\Psi$ and not $\Xi$. Lukily, for the wavelet transform, one has $$ \Xi = \Psi \text{diag(U)} f $$ where $U_m$ account for the redundancy of the scale of the atom $\psi_m$.

Compute the scaling factor (inverse of the redundancy).

J = Jmax - Jmin + 1
u = np.hstack(([4 ** (-J)], 4 ** (-np.floor(np.arange(J + 2.0 / 3, 1, -1.0 / 3)))))
U = np.transpose(np.tile(u, (n, n, 1)), (2, 0, 1))

Choose a value of the regularization parameter.

lambd = 0.01

Shortcut for the wavelet transform and the reconstruction.

Important: Scilab users have to create files |Xi.m|, |PsiS.m| and |Psi.m| to implement this function.

Xi = lambda a: perform_wavelet_transf(a, Jmin, -1, ti=1)
PsiS = lambda f: perform_wavelet_transf(f, Jmin, +1, ti=1)
Psi = lambda a: Xi(a / U)

The forward-backward algorithm now compute a series of wavelet coefficients $a^{(\ell)}$ computed as $$a^{(\ell+1)} = S_{\tau\lambda}( a^{(\ell)} + \Psi^*\Phi( y - \Phi\Psi a^{(\ell)} ) ). $$

The soft thresholding is defined as: $$\forall m, \quad S_T(a)_m = \max(0, 1-T/\|a_m\|)a_m. $$

The step size should satisfy: $$\tau < \frac{2}{\|\Psi\Phi \|} \leq 2 \min( u ). $$

tau = 1.9 * np.min(u)

Initialize the wavelet coefficients with those of the previous reconstruction.

a = U * PsiS(fSpars)

Gradient descent.

fTI = Psi(a)
a = a + tau * PsiS(Phi(y - Phi(fTI, Omega), Omega))

Soft threshold.

a = SoftThresh(a, lambd * tau)

Worked example 3

Perform the iterative soft thresholding. Monitor the decay of the energy $E$.

niter = 1000
a = U * PsiS(fSpars)
E = []
for i in range(niter):
    fTI = Psi(a)
    d = y - Phi(fTI, Omega)
    E = E + [1 / 2 * linalg.norm(d, "fro") ** 2 + lambd * np.sum(abs(a))]
    # step
    a = SoftThresh(a + tau * PsiS(Phi(d, Omega)), lambd * tau)

plt.figure(figsize=(7, 5))
plt.plot(E)
plt.show()
No description has been provided for this image

Perform the reconstruction.

fTI = Psi(a)

Display the result.

plt.figure(figsize=(6, 6))
imageplot(clamp(fTI))
No description has been provided for this image

Worked example 4

Perform the iteration with a decaying value of $\lambda$

niter = 3000
lambda_list = np.linspace(0.03, 0, niter)
err = []

for i in range(niter):
    fTI = Psi(a)
    d = y - Phi(fTI, Omega)
    # step
    a = SoftThresh(a + tau * PsiS(Phi(d, Omega)), lambda_list[i] * tau)

plt.figure(figsize=(6, 6))
imageplot(clamp(fSpars), "Sparsity inpainting TI, SNR = %.1f dB" % snr(f0, fTI))
No description has been provided for this image

Inpainting using Iterative Hard Thresholding

To improve the sparsity of the solution, it is possible to replace the soft thresholding by a hard threshdoling. In this case, the resulting algorihtm does not perform anymore a variational minimization of an energy.

The hard thresholding is defined as $h_T(x)=0$ if $-T < x < T$ and $h_T(x)=x$ otherwise. It thus defines a thresholding operator of wavelet coefficients as $H_T(a)_m = h_T(a_m)$.

Define a shortcut for this vectorialized hard thresholding

Important: Scilab users have to create a file |HardThresh.m| to implement this function.

HardThresh = lambda x, t: x * (abs(x) > t)

Display a curve of the 1-D Hard thresholding.

x = np.linspace(-1, 1, 1000)

plt.figure(figsize=(7, 5))
plt.plot(x, HardThresh(x, 0.5))
plt.show()
No description has been provided for this image

The hard thresholding in the translation invariant wavelet basis $\Psi$ reads $$ H_T^\Psi(f) = \Xi \circ H_T \circ \Psi^* (f) $$ where $\Xi = (\Phi^*)^+$ is the reconstruction operator.

We follow the MCA paradigm of Jean-Luc Starck, that alternates between a gradient descent step and a hard thresholding denoising, using a decaying threshold. $$f^{(\ell+1)} = H_{\tau\lambda_\ell}^\Psi( f^{(\ell)} - \tau \Phi^*(\Phi f^{(\ell)} - y) ). $$

Number of iterations.

niter = 500

List of thresholds. One must start by a large enough initial threshold.

lambda_list = np.linspace(1, 0, niter)

Initialization.

fHard = y

Gradient descent.

fHard = ProjC(fHard, Omega)

Hard threshold (here $\lambda=\lambda_0$) is used).

fHard = Xi(HardThresh(PsiS(fHard), tau * lambda_list[1]))

Worked example 5

Perform the iteration with a decaying value of $\lambda$

lambda_list = np.linspace(1, 0, niter)
fHard = y

for i in range(niter):
    fHard = Xi(HardThresh(PsiS(ProjC(fHard, Omega)), lambda_list[i]))


plt.figure(figsize=(6, 6))
imageplot(clamp(fSpars), "Inpainting hard thresh., SNR = %.1f dB" % snr(f0, fHard))
No description has been provided for this image

References and further reading