Texture Synthesis and Inpainting using Patch Projections¶

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 texture synthesis and inpainting using patch averaging.

In [2]:

addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/graphics_6_patches')

Creating a Dictionary of Patches¶

Given an exemplar image, we extract many patch that are our library. We even perform PCA dimensionality reduction to speed up nearest neighbors search.

The main parameter is the width of the patches.

In [3]:

w = 5;

The other parameter is the number of patch in the exemplar dictionary.

In [4]:

q = 1000;

We load an exemplar image.

In [5]:

name = 'corral';
n = 128;
M0 = load_image(name,n);
M0 = rescale( crop(M0,n) );

We set up larges |(w,w,n-w+1,n-w+1)| matrices representing the X and Y position of the pixel to extract.

In [6]:

p = n-w+1;

location of pixels

In [7]:

[Y,X] = meshgrid(1:p,1:p);

offsets

In [8]:

[dY,dX] = meshgrid(0:w-1,0:w-1);

location of pixels to extract

In [9]:

X = reshape(X, [1 1 p p]);
Y = reshape(Y, [1 1 p p]);
X = repmat(X, [w w 1 1]) + repmat(dX, [1 1 p p]);
Y = repmat(Y, [w w 1 1]) + repmat(dY, [1 1 p p]);

We extract all the patches and reshape the matrix of patch so that |P(:,:,i)| is a patch

In [10]:

P0 = M0(X + (Y-1)*n);
P0 = reshape(P0,w,w,p*p);

Sub-sample the patches

In [11]:

sel = randperm(size(P0,3)); sel = sel(1:q);
P0 = P0(:,:,sel);

Image Patch-wise Projection¶

The basic step for synthesis or inpainting is to project each patch of an image to its nearest neighbor in the dictionary.

The initial image is just noise for instance.

In [12]:

n = 128;
M = rand(n);

We define an offset vector to shift the projected patch. This needs to be changed during the iteration of the algorithm (either synthesis or inpainting).

In [13]:

ofx = 2;
ofy = 1;

Patch locations.

In [14]:

[Y,X] = meshgrid(1:w:n, 1:w:n);
p = size(X,1);
[dY,dX] = meshgrid(0:w-1,0:w-1);
X = reshape(X, [1 1 p p]);
Y = reshape(Y, [1 1 p p]);
X = repmat(X, [w w 1 1]) + repmat(dX, [1 1 p p]);
Y = repmat(Y, [w w 1 1]) + repmat(dY, [1 1 p p]);

Shift location, with proper boundary condition (cyclic).

In [15]:

Xs = mod(X+ofx-1, n)+1;
Ys = mod(Y+ofy-1, n)+1;

Extract the patches.

In [16]:

P = M(Xs + (Ys-1)*n);

Replace each patch by its closest match.

In [17]:

for i=1:p*p
    % distance to current patch
    d = sum(sum( (P0 - repmat(P(:,:,i), [1 1 q])).^2 ) );
    % best match
    [tmp,s] = min(d);
    % replace the patch
    P(:,:,i) = P0(:,:,s);
end

Reconstruct the image.

In [18]:

Mp = M;
Mp(Xs + (Ys-1)*n) = P;

Display.

In [19]:

clf;
imageplot(M,'Input', 1,2,1);
imageplot(Mp,'Projected', 1,2,2);

Texture Synthesis¶

Texture synthesis is obtained by performing the projection for several offset and averaging the results.

To speed up performance, we consider only a subset of all possible |w*w| offsets.

In [20]:

noffs = 10;

Compute a randomized list of offsets

In [21]:

sel = randperm(w*w); sel = sel(1:noffs);
OffX = dX(sel); OffY = dY(sel);

We perform one step of synthesis by cycling through the offset.

In [22]:

M1 = zeros(n);
for j=1:noffs
    ofx = OffX(j);
    ofy = OffY(j);
    % shift locations
    Xs = mod(X+ofx-1, n)+1;
    Ys = mod(Y+ofy-1, n)+1;
    % extract patch
    P = M(Xs + (Ys-1)*n);
    % replace by closest patch
    for i=1:p*p
        d = sum(sum( (P0 - repmat(P(:,:,i), [1 1 q])).^2 ) );
        [tmp,s] = min(d);
        P(:,:,i) = P0(:,:,s);
    end
    % reconstruct the image.
    M1(Xs + (Ys-1)*n) = M1(Xs + (Ys-1)*n) + P;
end
M1 = M1 / noffs;

To enhance the synthesis result, we perform histogram equalization.

In [23]:

M1 = perform_hist_eq(M1,M0);

Display the result.

In [24]:

clf;
imageplot(M,'Input', 1,2,1);
imageplot(M1,'Projected', 1,2,2);

Exercise 1

Perform several step of synthesis.

In [25]:

exo1()