Inpainting using Sparse Regularization
This numerical tour explores the use of sparse energies 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 Scilabe) or numericaltour.m to write all the Scilab/Matlab command 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. 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;
M = load_image('lena');
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;
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 with Iterative Thresholding
We perform iterative soft thresholding to regularize the result of the inpainting. For an orthogonal wavelet transform, this corresponds to the minimization of the L1 norm of the wavelet coefficients.
Set up parameters for the wavele transform.
% we use orthogonal transform options.ti = 0; % we use only a few scales for the transform Jmax = log2(n)-1; Jmin = Jmax-3;
Since there is no noise, lambda should be chosen as small as possible (but the smaller lambda is the longer it will take to do the computations).
We will thus decay the value of lambda through the iterations
% number of iterations niter = 500; % list of lambda lambda_list = linspace(.02,0, niter);
Initialize the solution.
Mspars = y; Mspars(mask==1) = 1/2;
Number of the iteration and corresponding threshold.
i = 1; lambda = lambda_list(i);
Denoise the solution by thresholding.
MW = perform_wavelet_transf(Mspars, Jmin, +1,options);
MW = perform_thresholding( MW, lambda, 'soft' );
Mspars = perform_wavelet_transf(MW, Jmin, -1,options);
Project on the inpainting constraint.
Mspars(mask==0) = y(mask==0);
Exercice 1: (the solution is exo1.m) Compute the solution of inpainting using iterative projected thresholding (MCA).
exo1;
Display the result.
clf; imageplot(clamp(M), strcat('Original'), 1,2,1); imageplot(clamp(Mspars), strcat(['Sparsity, SNR=' num2str(snr(M,Mspars),3) 'dB']), 1,2,2);
Exercice 2: (the solution is exo2.m) Compare with the result Mb of Sobolev inpainting.
exo2;
Translation Invariant Wavelet Inpainting
The orthogonal wavelet transform is not really efficient for inpainting because it lacks translation invariant.
One can use instead a translation invariant wavelet tight frame.
options.ti = 1;
To speed up the compuation, we initialize the process using the result of the Sobolev inpainting.
MsparsTI = Mb;
Better results and faster algorithm is obtained by using only very high frequency wavelets.
Jmin = Jmax;
We decay the value of lambda through the iterations
% number of iterations niter = 500; % list of lambda lambda_list = linspace(.02,0, niter);
An iteration of the MCA algorithm compute a soft thresholding and then project on the value of the known pixels. To enhance the result, we do not threshold the low pass residual of the wavelet transform.
i = 1; % forward transform MW = perform_wavelet_transf(MsparsTI, Jmin, +1,options); % threshold MWT = perform_thresholding( MW, lambda_list(i), 'soft' ); % keep low pass MWT(:,:,1) = MW(:,:,1); % inverse transform MsparsTI = perform_wavelet_transf(MWT, Jmin, -1,options);
Note that in this case, the iterative thresholding does not corresponds to the minimization of L1 the wavelet coefficients. This algorihtm corresponds to the Morphological Component Analysis of Starck, Elad and Donoho.
Exercice 3: (the solution is exo3.m) Compute the solution of inpainting using iterative projected thresholding (MCA), with translation invariant tight frames.
exo3;
