0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028 function sgwt_demo1
0029 close all
0030 fprintf('Welcome to SGWT demo #1\n');
0031
0032 npoints=500;
0033 fprintf('Creating Swiss Roll point cloud with %g points\n',npoints);
0034 dataparams=struct('n',npoints,'dataset',-1','noise',0,'state',0);
0035 r=create_synthetic_dataset(dataparams);
0036 x=rescale_center(r.x);
0037
0038
0039 fprintf('Computing edge weights and graph Laplacian\n');
0040
0041 d=distanz(x);
0042 s=.1;
0043 A=exp(-d.^2/(2*s^2));
0044 L=full(sgwt_laplacian(A));
0045
0046 fprintf('Measuring largest eigenvalue, lmax = ');
0047
0048
0049 Nscales=4;
0050 lmax=sgwt_rough_lmax(L);
0051 fprintf('%g\n',lmax);
0052
0053 fprintf('Designing transform in spectral domain\n');
0054 [g,gp,t]=sgwt_filter_design(lmax,Nscales);
0055 arange=[0 lmax];
0056
0057 figure
0058 sgwt_view_design(g,t,arange);
0059 ylim([0 3])
0060 set(gcf,'position',[0 780,600,300])
0061
0062 m=50;
0063 fprintf('Computing Chebyshev polynomials of order %g for fast transform \n',m);
0064 for k=1:numel(g)
0065 c{k}=sgwt_cheby_coeff(g{k},m,m+1,arange);
0066 end
0067
0068
0069 jcenter=32;
0070 fprintf('Computing forward transform of delta at vertex %g\n',jcenter);
0071 N=size(L,1);
0072 d=sgwt_delta(N,jcenter);
0073
0074 wpall=sgwt_cheby_op(d,L,c,arange);
0075
0076 fprintf('Displaying wavelets\n');
0077 msize=100;
0078 cp=[-1.4,-16.9,3.4];
0079
0080
0081
0082 ws=300;
0083 figure;
0084 xp=0; yp=ws+100;
0085 set(gcf,'position',[xp,yp,ws-10,ws+10]);
0086 scatter3(x(1,:),x(2,:),x(3,:),msize,d,'.');
0087 set(gcf,'Colormap',[.5 .5 .5;1 0 0]);
0088 clean_axes(cp);
0089 title(sprintf('Vertex %g',jcenter));
0090
0091
0092 for n=1:Nscales+1
0093 wp=wpall{n};
0094 figure
0095 xp=mod(n,3)*(ws+10);
0096 yp=(1-floor((n)/3))*(ws+100);
0097 set(gcf,'position',[xp,yp,ws-10,ws+10]);
0098 scatter3(x(1,:),x(2,:),x(3,:),msize,wp,'.');
0099 colormap jet
0100 caxis([-1 1]*max(abs(wp)));
0101 clean_axes(cp);
0102
0103 hcb=colorbar('location','north');
0104 cxt=get(hcb,'Xtick');
0105 cxt=[cxt(1),0,cxt(end)];
0106 set(hcb,'Xtick',cxt);
0107 cpos=get(hcb,'Position');
0108 cpos(4)=.02;
0109 set(hcb,'Position',cpos);
0110 set(hcb,'Position',[.25 .91 .6 .02]);
0111
0112 if n==1
0113 title('Scaling function');
0114 else
0115 title(sprintf('Wavelet at scale j=%g, t_j = %0.2f',n-1,t(end+1-(n-1))));
0116
0117 end
0118 end
0119
0120
0121 function clean_axes(cp)
0122 xlim([-1 1]);ylim([-1 1]);zlim([-1 1]);
0123 set(gca,'Xtick',[-1 0 1]);
0124 set(gca,'Ytick',[-1 0 1]);
0125 set(gca,'Ztick',[-1 0 1]);
0126 axis square
0127 set(gca,'CameraPosition',cp);
0128
0129
0130
0131
0132
0133
0134 function r=rescale_center(x)
0135 N=size(x,2);
0136 d=size(x,1);
0137 x=x-repmat(mean(x,2),[1,N]);
0138 c=max(abs(x(:)));
0139 r=x/c;