Home > sgwt_toolbox > demo > demo1.m

demo1

PURPOSE ^

demo1 : SGWT for swiss roll data set

SYNOPSIS ^

function demo1

DESCRIPTION ^

 demo1 : SGWT for swiss roll data set

 This demo builds the SGWT for the swiss roll synthetic data set. It
 computes a set of scales adapted to the computed upper bound on the
 spectrum of the graph Laplacian, and displays the scaling function and
 the scaled wavlet kernels, as well as the corresponding frame bounds. It
 then computes the wavelets centered on a single vertex, and displays
 them. This essentally reproduces figure 3 from 
 Hammond,Vangergheynst, Gribonval 2010.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % demo1 : SGWT for swiss roll data set
0002 %
0003 % This demo builds the SGWT for the swiss roll synthetic data set. It
0004 % computes a set of scales adapted to the computed upper bound on the
0005 % spectrum of the graph Laplacian, and displays the scaling function and
0006 % the scaled wavlet kernels, as well as the corresponding frame bounds. It
0007 % then computes the wavelets centered on a single vertex, and displays
0008 % them. This essentally reproduces figure 3 from
0009 % Hammond,Vangergheynst, Gribonval 2010.
0010 
0011 % Create swiss roll point cloud
0012 function demo1 
0013 close all
0014 fprintf('Welcome to SGWT demo #1\n');
0015 
0016 npoints=500;
0017 fprintf('Creating Swiss Roll point cloud with %g points\n',npoints);
0018 dataparams=struct('n',npoints,'dataset',-1','noise',0,'state',0);
0019 r=create_synthetic_dataset(dataparams);
0020 x=rescale_center(r.x);
0021 
0022 
0023 fprintf('Computing edge weights and graph Laplacian\n');
0024 % Compute Weighted graph adjacency matrix, and graph Laplacian
0025 d=distanz(x);
0026 s=.1;
0027 A=exp(-d.^2/(2*s^2)); 
0028 L=full(sgwt_laplacian(A));
0029 
0030 fprintf('Measuring largest eigenvalue, lmax = ');
0031 
0032 %% Design filters for transform
0033 Nscales=4;
0034 lmax=sgwt_rough_lmax(L);
0035 fprintf('%g\n',lmax);
0036 
0037 fprintf('Designing transform in spectral domain\n'); 
0038 [g,gp,t]=sgwt_filter_design(lmax,Nscales);
0039 arange=[0 lmax];
0040 %% Display filter design in spectral domain
0041 figure
0042 sgwt_view_design(g,t,arange);
0043 ylim([0 3])
0044 set(gcf,'position',[0 780,600,300])
0045 %% Chebyshev polynomial approximation
0046 m=50; % order of polynomial approximation
0047 fprintf('Computing Chebyshev polynomials of order %g for fast transform \n',m);
0048 for k=1:numel(g)
0049   c{k}=sgwt_cheby_coeff(g{k},m,m+1,arange);
0050 end
0051 
0052 %% compute transform of delta at one vertex
0053 jcenter=32; % vertex to center wavelets to be shown
0054 fprintf('Computing forward transform of delta at vertex %g\n',jcenter);
0055 N=size(L,1);
0056 d=sgwt_delta(N,jcenter);
0057 % forward transform, using chebyshev approximation
0058 wpall=sgwt_cheby_op(d,L,c,arange);
0059 
0060 fprintf('Displaying wavelets\n');
0061 msize=100;
0062 cp=[-1.4,-16.9,3.4]; % camera position
0063 %% Visualize result
0064 
0065 % show original point
0066 ws=300;
0067 figure;
0068 xp=0; yp=ws+100;
0069 set(gcf,'position',[xp,yp,ws-10,ws+10]);
0070 scatter3(x(1,:),x(2,:),x(3,:),msize,d,'.');
0071 set(gcf,'Colormap',[.5 .5 .5;1 0 0]);
0072 clean_axes(cp);
0073 title(sprintf('Vertex %g',jcenter));
0074 
0075 % show wavelets
0076 for n=1:Nscales+1
0077     wp=wpall{n};
0078     figure
0079     xp=mod(n,3)*(ws+10);
0080     yp=(1-floor((n)/3))*(ws+100);
0081     set(gcf,'position',[xp,yp,ws-10,ws+10]);
0082     scatter3(x(1,:),x(2,:),x(3,:),msize,wp,'.');
0083     colormap jet
0084     caxis([-1 1]*max(abs(wp)));
0085     clean_axes(cp);
0086 
0087     hcb=colorbar('location','north');
0088     cxt=get(hcb,'Xtick');
0089     cxt=[cxt(1),0,cxt(end)];
0090     set(hcb,'Xtick',cxt);
0091     cpos=get(hcb,'Position');
0092     cpos(4)=.02; % make colorbar thinner
0093     set(hcb,'Position',cpos);
0094     set(hcb,'Position',[.25 .91 .6 .02]);
0095     
0096     if n==1
0097       title('Scaling function');
0098     else      
0099       title(sprintf('Wavelet at scale j=%g, t_j = %0.2f',n-1,t(end+1-(n-1))));
0100 
0101     end
0102 end
0103 
0104 
0105 function clean_axes(cp)
0106 xlim([-1 1]);ylim([-1 1]);zlim([-1 1]);
0107 set(gca,'Xtick',[-1 0 1]);
0108 set(gca,'Ytick',[-1 0 1]);
0109 set(gca,'Ztick',[-1 0 1]);
0110 axis square
0111 set(gca,'CameraPosition',cp);

Generated on Fri 30-Apr-2010 17:49:57 by m2html © 2003