%GSP_DEMO_WAVELET Introduction to spectral graph wavelet with the GSPBox
% 
%   The wavelets are a special type of filterbank. In this demo, we will
%   show how you can very easily construct a wavelet frame and apply it to
%   a signal. If you want to do find an interactive demo of the wavelet, we
%   encourage you to use the sgwt_demo2 of the sgwt toolbox. It can be
%   downloaded at:
%
%   http://wiki.epfl.ch/sgwt/documents/sgwt_toolbox-1.02.zip
%
%   The sgwt toolbox has the same core as the GSPBox and all his functions
%   have equivalent in the GSPBox ( Except the demos ;-) ).
%
%   In this demo we will show you how to compute the wavelet coefficients
%   of a graph and visualize them. First, let's load a graph :
%
%       G = gsp_bunny();
%
%   This graph is a nearest-neighbor graph of a pointcloud of the Stanford
%   bunny. It will allow us to get interesting visual results using
%   wavelets. 
%
%   At this stage we could compute the full Fourier basis using
%   GSP_COMPUTE_FOURIER_BASIS, but this would take a lot of time, and can
%   be avoided by using Chebychev polynomials approximations. This operation
%   is implemented in most function and is thus completely transparent.
%
%   Simple filtering
%   ----------------
%
%   Before tackling wavelets, we can see the effect of one filter localized
%   on the graph. So we can first design a few heat kernel filters :
%       
%       taus = [1, 10, 25, 50];
%       Hk = gsp_design_heat(G, taus);
%
%   Let's now create a signal as a Kronecker located on one vertex (e.g.
%   the vertex 100) :
%
%       S = zeros(G.N, 1);
%       vertex_delta = 83;
%       S(vertex_delta) = 1;
% 
%       Sf_vec = gsp_filter_analysis(G, Hk, S);
%       Sf = gsp_vec2mat(Sf_vec, length(taus));
%
%   and plot the filtered signal :
%
%       param_plot.cp = [0.1223, -0.3828, 12.3666];
% 
%       figure;
%       subplot(221)
%       gsp_plot_signal(G,Sf(:,1), param_plot);
%       axis square
%       title(sprintf('Heat diffusion tau = %d', taus(1)));
%       subplot(222)
%       gsp_plot_signal(G,Sf(:,2), param_plot);
%       axis square
%       title(sprintf('Heat diffusion tau = %d', taus(2)));
%       subplot(223)
%       gsp_plot_signal(G,Sf(:,3), param_plot);
%       axis square
%       title(sprintf('Heat diffusion tau = %d', taus(3)));
%       subplot(224)
%       gsp_plot_signal(G,Sf(:,4), param_plot);
%       axis square
%       title(sprintf('Heat diffusion tau = %d', taus(4)));
%
%   Figure 1: Heat diffusion at different scales
%
%
%   Visualizing wavelets atoms
%   --------------------------
%   
%   Let's now replace the filtering by the heat kernel by a filter bank of
%   wavelets. We can create a filter bank using one of the design functions
%   such as GSP_DESIGN_MEXICAN_HAT :
%
%         Nf = 6;
%         Wk = gsp_design_mexican_hat(G, Nf);
%
%   We can plot the filter bank spectrum :
%
%         figure;
%         gsp_plot_filter(G,Wk);
%   
%   Figure 2: Wavelets filterbank (Original)
%
%
%   As we can see, the wavelets atoms are stacked on the low frequency part
%   of the spectrum. If we want to get a better coverage of the graph
%   spectrum we can use the function GSP_DESIGN_WARPED_TRANSLATES :
%
%         param_filter.filter = Wk;
%         Wkw = gsp_design_warped_translates(G,Nf,param_filter);
%
%   Now let's plot the new filter bank :
%
%         figure;
%         gsp_plot_filter(G,Wkw);
%   
%   Figure 3: Wavelet filterbank (spectrum adaptated)
%
%
%   We can see that the wavelet atoms are much more spread along the graph
%   spectrum. We can visualize the filtering by one atom as we did with the
%   heat kernel, by placing a Kronecker delta at one specific vertex and
%   filter using the filter bank :
% 
%         S = zeros(G.N*Nf,Nf);
%         S(vertex_delta) = 1;
%         for ii=1:Nf
%             S(vertex_delta+(ii-1)*G.N,ii) = 1;
%         end
% 
%         Sf = gsp_filter_synthesis(G,Wkw,S);
%   
%   We can plot the resulting signal for the different scales :
%
%         figure;
%         subplot(221)
%         gsp_plot_signal(G,Sf(:,1), param_plot);
%         axis square
%         mu = mean(Sf(:,1));
%         sigma = std(Sf(:,1));
%         c_scale = 4;
%         caxis([mu - c_scale*sigma, mu + c_scale*sigma]);
%         title('Wavelet 1');
%
%         subplot(222)
%         gsp_plot_signal(G,Sf(:,2), param_plot);
%         axis square
%         mu = mean(Sf(:,2));
%         sigma = std(Sf(:,2));
%         caxis([mu - c_scale*sigma, mu + c_scale*sigma]);
%         title('Wavelet 2');
%
%         subplot(223)
%         gsp_plot_signal(G,Sf(:,3), param_plot);
%         axis square
%         mu = mean(Sf(:,3));
%         sigma = std(Sf(:,3));
%         caxis([mu - c_scale*sigma, mu + c_scale*sigma]);
%         title('Wavelet 3');
%
%         subplot(224)
%         gsp_plot_signal(G,Sf(:,4), param_plot);
%         axis square
%         mu = mean(Sf(:,4));
%         sigma = std(Sf(:,4));
%         caxis([mu - c_scale*sigma, mu + c_scale*sigma]);
%         title('Wavelet 4');
%
%   Figure 4: A few wavelets atoms
%
%
%   Curvature estimation
%   --------------------
%   
%   As a last and more applied example, let us try to estimate the
%   curvature of the underlying 3D model by only using only spectral
%   filtering on the graph. 
%
%   A simple way to accomplish that is to use the
%   coordinates map [x, y, z] and filter it using the wavelets defined
%   above. We obtain a 3-dimensional signal [hat(x), hat(y), hat(z)]*
%   which describes variation along the 3 coordinates :
%
%         s_map = G.coords;
%         s_map_out = gsp_filter_analysis(G, Wk, s_map);
%         s_map_out = gsp_vec2mat(s_map_out, Nf);
%   
%   Finally we can get the curvature estimation by taking the l_1 or
%   l_2 norm of the filtered signal :
%
%         dd = s_map_out(:,:,1).^2 + s_map_out(:,:,2).^2 + s_map_out(:,:,3).^2;
%         dd = sqrt(dd);
%
%   Let's now plot the result to observe that we indeed have a measure of
%   the curvature :
%
%         figure;
%         subplot(221)
%         gsp_plot_signal(G,dd(:,2), param_plot);
%         axis square
%         title('Curvature estimation scale 1');
%         subplot(222)
%         gsp_plot_signal(G,dd(:,3), param_plot);
%         axis square
%         title('Curvature estimation scale 2');
%         subplot(223)
%         gsp_plot_signal(G,dd(:,4), param_plot);
%         axis square
%         title('Curvature estimation scale 3');
%         subplot(224)
%         gsp_plot_signal(G,dd(:,5), param_plot);
%         axis square
%         title('Curvature estimation scale 4');
%
%   Figure 5: Curvature estimation using wavelet feature
%
%
%
%   Url: https://epfl-lts2.github.io/gspbox-html/doc/demos/gsp_demo_wavelet.html

% Copyright (C) 2013-2016 Nathanael Perraudin, Johan Paratte, David I Shuman.
% This file is part of GSPbox version 0.7.5
%
% 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: Johan Paratte
% Date : 21 August 2014

%% Initialization
clear;
close all;

%% Load the graph of the bunny
G = gsp_bunny();
G = gsp_estimate_lmax(G);

%% Heat kernel
taus = [10, 15, 20, 30];
Hk = gsp_design_heat(G, taus);

S = zeros(G.N, 1);
vertex_delta = 83;
S(vertex_delta) = 1;

Sf_vec = gsp_filter_analysis(G, Hk, S);
Sf = gsp_vec2mat(Sf_vec, length(taus));

param_plot.cp = [0.1223, -0.3828, 12.3666];

subplot(221)
gsp_plot_signal(G,Sf(:,1), param_plot);
axis square
title(sprintf('Heat diffusion tau = %d', taus(1)));
subplot(222)
gsp_plot_signal(G,Sf(:,2), param_plot);
axis square
title(sprintf('Heat diffusion tau = %d', taus(2)));
subplot(223)
gsp_plot_signal(G,Sf(:,3), param_plot);
axis square
title(sprintf('Heat diffusion tau = %d', taus(3)));
subplot(224)
gsp_plot_signal(G,Sf(:,4), param_plot);
axis square
title(sprintf('Heat diffusion tau = %d', taus(4)));


%% Wavelets
Nf = 6;

Wk = gsp_design_mexican_hat(G, Nf);

figure;
gsp_plot_filter(G,Wk);

param_filter.filter = Wk;
Wkw = gsp_design_warped_translates(G,Nf,param_filter);
 
figure;
gsp_plot_filter(G,Wkw);

S = zeros(G.N*Nf,Nf);
S(vertex_delta) = 1;
for ii=1:Nf
    S(vertex_delta+(ii-1)*G.N,ii) = 1;
end

Sf = gsp_filter_synthesis(G,Wk,S);

figure;
subplot(221)
gsp_plot_signal(G,Sf(:,1), param_plot);
axis square
mu = mean(Sf(:,1));
sigma = std(Sf(:,1));
c_scale = 4;
caxis([mu - c_scale*sigma, mu + c_scale*sigma]);
title('Wavelet 1');
subplot(222)
gsp_plot_signal(G,Sf(:,2), param_plot);
axis square
mu = mean(Sf(:,2));
sigma = std(Sf(:,2));
caxis([mu - c_scale*sigma, mu + c_scale*sigma]);
title('Wavelet 2');
subplot(223)
gsp_plot_signal(G,Sf(:,3), param_plot);
axis square
mu = mean(Sf(:,3));
sigma = std(Sf(:,3));
caxis([mu - c_scale*sigma, mu + c_scale*sigma]);
title('Wavelet 3');
subplot(224)
gsp_plot_signal(G,Sf(:,4), param_plot);
axis square
mu = mean(Sf(:,4));
sigma = std(Sf(:,4));
caxis([mu - c_scale*sigma, mu + c_scale*sigma]);
title('Wavelet 4');

%% Curvature estimation

s_map = G.coords;
s_map_out = gsp_filter_analysis(G, Wk, s_map);
s_map_out = gsp_vec2mat(s_map_out, Nf);

dd = s_map_out(:,:,1).^2 + s_map_out(:,:,2).^2 + s_map_out(:,:,3).^2;
dd = sqrt(dd);

figure;
subplot(221)
gsp_plot_signal(G,dd(:,2), param_plot);
axis square
title('Curvature estimation scale 1');
subplot(222)
gsp_plot_signal(G,dd(:,3), param_plot);
axis square
title('Curvature estimation scale 2');
subplot(223)
gsp_plot_signal(G,dd(:,4), param_plot);
axis square
title('Curvature estimation scale 3');
subplot(224)
gsp_plot_signal(G,dd(:,5), param_plot);
axis square
title('Curvature estimation scale 4');


    

