Variational Inpainting
This numerical tour explores the use of variational energies (Sobolev, total variation and sparsity) to regularize the image inpaiting problem.
Contents
Installing toolboxes and setting up the path.
You need to download the general purpose toolbox and the signal toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_general/ and toolbox_signal/ in your directory.
For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.
Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab commands you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.
Execute this line only if you are using Matlab.
getd = @(p)path(path,p); % scilab users must *not* execute this
Then you can add these toolboxes to the path.
% Add some directories to the path getd('toolbox_signal/'); getd('toolbox_general/');
Missing Pixels and Inpaiting
Inpainting corresponds to filling holes in images.
First we load the image to be inpainted.
name = 'lena'; name = 'cameraman'; n = 256; M = load_image(name, n); M = rescale(crop(M,n));
Then we construct a mask by removing random pixels.
% amount of removed pixels. rho = .7; % random mask, mask==1 for removed pixels mask = zeros(n,n); sel = randperm(n^2); sel = sel(1:round(rho*n^2)); mask(sel) = 1;
Compute the damaged observation (no noise).
% remove pixels y = M; y(mask==1) = 0; % display clf; imageplot(M, 'Image M', 1,2,1); imageplot(y, 'Observation y', 1,2,2);
Inpainting by Blurring
A simple way to recover the missing values of the image is to perform a blurring. To progressively fill the hole, the blurring must be repeated.
Initially, we put an average gray for the missing value. In fact this has no effect on the final inpainted result, but a clever initialization speeds up the convergence of the algorithms.
Mb = y;
% put average value inside the hole
Mb(mask==1) = .5;
We perform a blurring to fill missing values. A large kernel size s product smoother result but might destroy edges.
% blur the whole image s = 3; % number of pixels Mb = perform_blurring(Mb,3); % impose the known values Mb(mask==0) = y(mask==0);
You can iterate this several time.
niter = 50; err = []; for i=1:niter % blur the whole image Mb = perform_blurring(Mb,s); % impose the known values Mb(mask==0) = y(mask==0); % compute the error err(i) = snr(M, Mb); end % display error decay clf; plot(1:niter, err); axis('tight'); set_label('#iteration', 'SNR');
Display.
clf; imageplot(y, 'Obervations y', 1,2,1); imageplot(Mb, strcat(['Inpainted, SNR=' num2str(snr(M,Mb),3) 'dB']), 1,2,2);
Sobolev Impainting
The simplest regularization is simply the Sobolev energy \sum grad(M)^2. Its gradient is the Laplacian, and the gradient descent corresponds to an iterated blurring with a small 3x3 filter.
Parameters for the projected gradient descent of the Sobolev norm. The gradient step size tau must be less than 1/4, according to the CLF condition.
% gradient descent step size tau = .2; % use centered differences with reflecting boundary conditions options.bound = 'sym'; options.order = 1; % initialization Msob = y;
The gradient of the Sobolev energy is (minus) the Laplacian.
% compute Laplacian L = div( grad(Msob, options), options ); % compare with the Laplacian of the original image L0 = div(grad(M,options),options); % display clf; imageplot(L0, 'Laplacian of M', 1,2,1); imageplot(L, 'Laplacian of y', 1,2,2);
The first step of the inpainting iteration is a gradient descent of the Sobolev norm. A step of gradient descent performs a small blurring with a 3x3 filter.
Msob = Msob + tau*L;
The second step correspond to a projection on the measurement constraint, by imposing the value of the known pixels.
Msob(mask==0) = y(mask==0);
Exercice 1: (the solution is exo1.m) Perform the inpainting by a projected gradient descent of the Sobolev energy. Keep track of the Sobolev energy during the iterations. Also keep track of the decay of the error err(i)=norm(Msob-MsobC)/norm(MsobC) where MsobC is the solution at convergence.
exo1;
Display the decay of the convergence error, in log plot. This shows the linear convergence speed. This linear rate (the slope of the line) can be improved by carefully chosing the value of the step size tau.
clf; plot(log10(err)); axis('tight'); set_label('Iteration #', 'log_{10}(error)');
Display the result.
clf; imageplot(y, 'Obervations y', 1,2,1); imageplot(Msob, strcat(['Inpainted, SNR=' num2str(snr(M,Msob),3) 'dB']), 1,2,2);
Optimal Step Size
At each iteration of the inpainting process, one can look for the optimal step size tau so that Msob + tau*Lap(Msob) minimizes the Sobolev norm. This corresponds to a 1D linesearch along the gradient direction.
The optimal step size is given as
tau = -<Lap(Msob),Lap(Msob)> / Lap(Lap(Msob))
Initialize the projected gradient descent.
Msob = y;
Compute optimal step size.
L = div( grad(Msob, options), options ); LL = div( grad(L, options), options ); tau = -norm(L(:))^2 / sum(L(:).*LL(:));
Perfrom descent.
Msob = Msob + tau*L;
Exercice 2: (the solution is exo2.m) Compute the projected gradient descent with optimal step size. Compare the decay of the convergence error for the fixed step size strategy and for the optimal line search.
exo2;
Inpainting with TV Regularization
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.
Parameters for the projected gradient descent of the (smoothed) TV norm. The value of tau, the step size, should be smaller than epsilon/4 for convergence.
% regularization parameter for the TV norm epsilon = 1e-2; % gradient descent step size tau = epsilon/4; % initialization Mtv = y; % number of iteration (quite a large number is required) niter = 500;
Compute the gradient of the smoothed TV functional.
Gr = grad(Mtv, options); d = sqrt( epsilon^2 + sum3(Gr.^2,3) ); G = -div( Gr./repmat(d, [1 1 2]), options );
Compute the TV norm, usefull to keep track of its decay through iterations.
tv = sum(d(:));
The first step of the inpainting iteration is a gradient descent of the TV norm.
Mtv = Mtv - tau*G;
The second step correspond to a projection on the measurement constraint, by imposing the value of the known pixels.
Mtv(mask==0) = y(mask==0);
Exercice 3: (the solution is exo3.m) Perform the inpainting by a projected gradient descent of the TV norm. Keep track of the TV energy. Also keep track of the decay of the error err(i)=norm(Mtv-MtvC)/norm(MtvC) where MtvC is the solution at convergence.
exo3;
Display the decay of the convergence error, in log plot. This shows the approximate linear convergence speed. This linear rate (the slope of the line) can be improved by carefully chosing the value of the step size tau, but this is more difficult than with Sobolev inpainting because the energy is non quadratic.
clf; plot(log10(err)); axis('tight'); set_label('Iteration #', 'log_{10}(error)');
Display the result.
clf; imageplot(clamp(Msob), strcat(['Sobolev, SNR=' num2str(snr(M,Msob),3) 'dB']), 1,2,1); imageplot(clamp(Mtv), strcat(['TV, SNR=' num2str(snr(M,Mtv),3) 'dB']), 1,2,2);
Inpainting with non-random mask
Inpainting can be used to remove objects in pictures.
We load another image with a non-random mask.
n = 256; M = load_image('parrot', n); if nb_dims(M)==3 % flatten color images M = sum3(M,3); end M = rescale( M ); mask1 = load_image('parrot-mask', n); mask1 = double(rescale(mask1)<.5);
Exercice 4: (the solution is exo4.m) Perform Sobolev inpainting on the image 'parrot' using as mask the binary image 'perrot-mask'.
exo4;
Exercice 5: (the solution is exo5.m) Try other methods to solve this inpainting problem. You can for instance have a look on the numerical on sparsity for deconvolution and inpainting.
exo5;