Homec4science
Diffusion iCAPs public (master) Browse

iCAPs/packages/gspbox/demos/gsp_demo_tv_inpainting.m
2cb18aae5832master

Last Change
Raw File

gsp_demo_tv_inpainting.m
View Options

%GSP_DEMO_GRAPH_TV Reconstruction of missing sample on a graph using TV
%
% In this demo, we try to reconstruct missing sample of a piece-wise
% smooth signal on a graph. To do so, we will minimize the well-known TV
% norm defined on the graph.
%
% For this example, you need the unlocbox. You can download it here:
% http://unlocbox.sourceforge.net/download
%
% We express the recovery problem as a convex optimization problem of the
% following form:
%
% argmin ||grad(x)||_1 s. t. ||Mx-b||_2 < epsilon
%
% Where b represents the known measurements, M is an operator
% representing the mask and epsilon is the radius of the l2 ball.
%
% We set
%
% f_1(x)=||nabla x _1
% We define the prox of f_1 as:
%
% prox_{f1,gamma} (z) = argmin_{x} 1/2 ||x-z||_2^2 + gamma ||grad(z)||_1
%
% f_2 is the indicator function of the set S define by Mx-b||_2 < epsilon
% We define the prox of f_2 as
%
% prox_{f2,gamma} (z) = argmin_{x} 1/2 ||x-z||_2^2 + gamma i_S( x ),
%
% with i_S(x) is zero if x is in the set S and infinity otherwise.
% This previous problem has an identical solution as:
%
% argmin_{z} ||x - z||_2^2 s.t. ||b - M z||_2 < epsilon
%
% It is simply a projection on the B2-ball.
%
% Results
% -------
%
% Figure 1: Original signal on graph
%
% This figure shows the original signal on graph.
%
% Figure 2: Depleted signal on graph
%
% This figure shows the signal on graph after the application of the
% mask and addition of noise. Half of the vertices are set to 0.
%
% Figure 3: Reconstructed signal on graph usign TV
%
% This figure shows the reconstructed signal thanks to the algorithm.
%
% Comparison with Tikhonov regularization
% ---------------------------------------
%
% We can also use the Tikhonov regularizer that will promote smoothness.
% In this case, we solve:
%
% argmin ||grad(x)||_2^2 s. t. ||Mx-b||_2 < epsilon
%
% The result is presented in the following figure:
%
% Figure 4: Reconstructed signal on graph using Tikhonov
%
% This figure shows the reconstructed signal thanks to the algorithm.
%
%
% Url: http://lts2research.epfl.ch/gsp/doc/demos/gsp_demo_tv_inpainting.php
% Copyright (C) 2013-2016 Nathanael Perraudin, Johan Paratte, David I Shuman.
% This file is part of GSPbox version 0.7.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% If you use this toolbox please kindly cite
% N. Perraudin, J. Paratte, D. Shuman, V. Kalofolias, P. Vandergheynst,
% and D. K. Hammond. GSPBOX: A toolbox for signal processing on graphs.
% ArXiv e-prints, Aug. 2014.
% http://arxiv.org/abs/1408.5781
% Author: Nathanael Perraudin
% Date: 4th March 2014
%% Initialisation
clear;
close all;
% Load toolbox
init_unlocbox();
%% Important parameters
% size of the graph for the demo
N = 250;
% probability of having a label on a vertex.
p = 0.1;
verbose = 1; % verbosity level
sigma = 0.0;
%% Create a graph
% paramgraph.distribute = 1;
% G = gsp_sensor(N, paramgraph);
G = gsp_sensor(N);
G = gsp_community(N);
% for a path with irregular weights and distances
% G = gsp_path(N);
% G = gsp_update_coordinates(G, G.coords + [rand(G.N, 1), zeros(G.N, 1)]);
% G = gsp_update_weights(G, G.W .* exp(-gsp_distanz(G.coords').^2));
G = gsp_adj2vec(G);
G = gsp_estimate_lmax(G);
G = gsp_compute_fourier_basis(G);
% x_0 = (1 + sign(G.U(:,4))) / 2;
x_0 = round((1 + sign(G.U(:,2))) * 2+(1 + sign(G.U(:,3))) +(1 + sign(G.U(:,4))) / 2);
% x_0 = round(linspace(0, 5, N))';
%% Set the labels
% create the mask
% ind_obs = boolean(full(sparse(1:4:G.N, 1, 1, G.N, 1)));
ind_obs = rand(G.N, 1) < p;
ind_unobs = not(ind_obs);
%applying the Mask to the data
x_obs = ind_obs.*(x_0+sigma*randn(G.N,1));
%% Method 1: Graph TV
% setting different parameter for the simulation
param_solver.verbose = verbose; % display parameter
param_solver.tol = 1e-5;
param_solver.maxit = 4000;
tic;
sol_tv = gsp_regression_tv(G, ind_obs, x_obs, 0, param_solver);
toc
%% Method 2: Tikhonov
tic;
G2 = G;
G2.L = G.L^3;
%sol_tik = gsp_regression_tik(G, ind_obs, x_obs, 0, param_solver);
sol_tik = gsp_regression_tik(G2, ind_obs, x_obs, 0, param_solver);
toc
%% Method 3: My idea: inpaint using graph TV on one variable that is l-2
%% away from another one that is smooth
%%
% x, y are primal variables: z is the dual
% min_x,y,z ||grad(x_all)||_1 + a/2||x-y||^2 + b/2 ||grad(y_all)||^2
%
% grad(x_all) = grad([x;x_obs]) = A*x + A_obs * x_obs = A*x - b
%
% only rows of gradient corresponding to unknown values used:
% min_x,y,z ||Ax - b||_1 + a/2||x-y||^2 + b/2 ||Ay - b||^2
% Get rid of operator A by using dual variable:
% min_x,y,z ||z - b||_1 + a/2||x-y||^2 + b/2 ||Ay - b||^2
% Think of all primal variables as one vector w = [x;y] containing ONLY
% the elements on the unlabeled nodes
% choose parameters for solving the problem:
alpha = .1;
beta = 1;
% compute the Diff matrix of the graph
G = gsp_adj2vec(G);
% operators based on the Diff
A_obs = G.Diff(:, ind_obs);
A_unobs = G.Diff(:, not(ind_obs));
n_unobs = nnz(not(ind_obs));
zeros_x = zeros(n_unobs, 1);
zeros_y = zeros_x;
% graph TV term is the l-1 on the dual:
b = - A_obs * x_obs(ind_obs);
f1_params.y = b;
f1.eval = @(z) norm(z - b, 1);
f1.prox = @(z, gamma) prox_l1(z, gamma, f1_params);
f1.L = @(w) A_unobs * w(1:n_unobs);
f1.Lt = @(z) [A_unobs' * z; zeros_y];
f1.norm_L = normest(A_unobs);
%f1.norm_L = sqrt(G.lmax);
% use prox for the second function
% w = [x;y], w0 = [x0;y0]
% prox_g(x) = ((1+g)x + gy)/(1+2g)
f2.prox = @(w, gamma) [ ((1+alpha*gamma) * w(1:n_unobs) + ...
alpha*gamma * w(n_unobs+1:end))/(1+2*alpha*gamma);...
(alpha*gamma * w(1:n_unobs) + ...
(1+alpha*gamma) * w(n_unobs+1:end))/(1+2*alpha*gamma)];
f2.eval = @(w) alpha/2 * norm(w(1:n_unobs) - w(n_unobs+1:end))^2;
% use gradient for the last function
AtA = A_unobs' * A_unobs;
Atb = A_unobs' * b;
f3.grad = @(w) [zeros_x; beta * (AtA*w(n_unobs+1:end) - Atb)];
f3.eval = @(w) beta/2 * norm(A_unobs * w(n_unobs+1:end) - b)^2;
f3.beta = beta * normest(AtA);
%% solve the problem
param_solver.verbose = 1;
param_solver.algo = 'FBF_PRIMAL_DUAL';
param_solver.normalized_timestep = 0.99;
tic;
[sol_mine, info_sol] = solvep([sol_tv(not(ind_obs)); sol_tik(not(ind_obs))], {f1,f2,f3}, param_solver);
toc
sol_mine_x = x_obs;
sol_mine_x(ind_unobs) = sol_mine(1:n_unobs);
sol_mine_y = x_obs;
sol_mine_y(ind_unobs) = sol_mine(n_unobs+1:end);
%% Method 4: L-1 of Lx ||Lx||_1
tic;
param_solver = rmfield(param_solver, 'algo');
sol_Lx_l1 = gsp_regression_Lx_l1(G, ind_obs, x_obs, 0, param_solver);
toc
%% Compute the errors
ans_tv = round(sol_tv);
ans_tik = round(sol_tik);
ans_mine = round(sol_mine_x);
ans_Lxl1 = round(sol_Lx_l1);
err_tv = nnz(ans_tv ~= x_0) / G.N
err_tik = nnz(ans_tik ~= x_0) / G.N
err_mine = nnz(ans_mine ~= x_0) / G.N
err_Lxl1 = nnz(ans_Lxl1 ~= x_0) / G.N
%% Print the result
paramplot.show_edges = 1;
val_lims = lin_map([-.05, 1.05], [min([x_0; sol_tv; sol_tik]), max([x_0; sol_tv; sol_tik])], [0, 1]);
%% Plot the original graph
if not(strcmp(G.type, 'path'))
figure(1)
paramplot.vertex_highlight = ind_obs;
gsp_plot_signal(G, x_0, paramplot)
caxis(val_lims);
title('Original signal: highlighted known values')
% % Let show depleted graph
% figure(2)
% gsp_plot_signal(G,depleted_graph_value,paramplot)
% caxis([-1 1])
% title('Measurement')
% Let show the reconstructed graph
figure(3)
paramplot.vertex_highlight = (ans_tv ~= x_0);
gsp_plot_signal(G, ans_tv, paramplot)
caxis(val_lims)
title('TV solution: highlighted mistakes')
% Let show the reconstructed graph
figure(4)
paramplot.vertex_highlight = (ans_tik ~= x_0);
gsp_plot_signal(G, ans_tik, paramplot)
caxis(val_lims)
title('Tikhonov solution: highlighted mistakes')
% Let show the reconstructed graph
figure(5)
paramplot.vertex_highlight = (ans_mine ~= x_0);
gsp_plot_signal(G, ans_mine, paramplot)
caxis(val_lims)
title('My solution: highlighted mistakes')
% Let show the reconstructed graph
figure(6)
paramplot.vertex_highlight = (ans_Lxl1 ~= x_0);
gsp_plot_signal(G, ans_Lxl1, paramplot)
caxis(val_lims)
title('||Lx||_1 solution: highlighted mistakes')
%%
else strcmp(G.type, 'path')
figure;
plot(G.coords(:, 1), sol_tv)
hold on;
plot(G.coords(:, 1), sol_tik)
plot(G.coords(:, 1), x_0);
plot(G.coords(:, 1), sol_Lx_l1)
plot(G.coords(ind_obs, 1), x_0(ind_obs), 'o');
legend('TV', 'Tikhonov', '||Lx||_1', 'true');
end
figure; plot(sol_mine_x)
hold on; plot(sol_mine_y)
c4science · Help