% Statistical and Computational Methods in Brain Image Analysis
% Chapter 11. Weighted Spherical Harmonic Representation
%
% (C) Moo K. Chung 
% University of Wisconsin-Madison
% mkchung@wisc.edu

load autism.surface.mat

% 1st subject, autism, inner surface

surf.vertices=squeeze(autisminner(1,:,:))';
surf.faces=tri;
figure; figure_wire(surf, 'yellow', 'white');

surf.vertices=squeeze(autismouter(1,:,:))';
surf.faces=tri;
figure; figure_wire(surf, 'yellow', 'white');


% two spherical angles

load sphere40962.mat

[theta varphi]=SPHARMangles(sphere40962);
figure; figure_trimesh(surf,theta,'rwb');
figure; figure_trimesh(surf,varphi,'rwb');

% spherical harmonic basis
Ylm=Y_lm(3,1,theta',varphi');
figure; figure_trimesh(sphere40962,Ylm,'rwb');


% weighted-SPHARM 
[surfsmooth, fourier]=SPHARMsmooth2(surf,sphere40962,18,0.01);
figure_wire(surfsmooth,'yellow','white');

% manipulating SPHARM coefficients

figure;
subplot(1,3,1)
imagesc(-abs(fourier.x)); colorbar
subplot(1,3,2)
imagesc(-abs(fourier.y)); colorbar;
subplot(1,3,3)
imagesc(-abs(fourier.z)); colorbar

colormap('gray'); 
figure_bg('white')

fourierv=SPHARMvectorize(fourier,18)

% reconstructing SPHARM from given coefficients

surf=SPHARMrepresent2(sphere40962, fourier, 18,0);
figure;figure_wire(surf,'white','white'); axis off


surf=SPHARMrepresent2(sphere40962, fourier, 18,0.01);
figure;figure_wire(surf,'white','white'); axis off

surf=SPHARMrepresent2(sphere40962, fourier, 18,0.1);
figure;figure_wire(surf,'white','white'); axis off


%-----------------------------------------
%Surface registration via SPHARM-correspondence

load autism.surface.mat

surfout1.vertices=squeeze(autismouter(1,:,:))';
surfout1.faces=tri;
surfin1.vertices=squeeze(autisminner(1,:,:))';
surfin1.faces=tri;


surfout2.vertices=squeeze(autismouter(2,:,:))';
surfout2.faces=tri;
surfin2.vertices=squeeze(autisminner(2,:,:))';
surfin2.faces=tri;

%thickness
thick1 = L2norm(surfout1.vertices-surfin1.vertices);
thick2 = L2norm(surfout2.vertices-surfin2.vertices);

load sphere40962.mat

%smoothing outer surface
[surfout1, fourier1]=SPHARMsmooth2(surfout1,sphere40962,42,0.001);
[surfout2, fourier2]=SPHARMsmooth2(surfout2,sphere40962,42,0.001);

%smoothing thickness
surfthick1.vertices=[thick1 thick1 thick1];
surfthick1.faces=tri;
[surfthick1, fourierthick1]=SPHARMsmooth2(surfthick1,sphere40962,42,0.001);
thick1smooth = surfthick1.vertices(:,1);

surfthick2.vertices=[thick2 thick2 thick2];
surfthick2.faces=tri;
[surfthick2, fourierthick2]=SPHARMsmooth2(surfthick2,sphere40962,42,0.001);
thick2smooth = surfthick2.vertices(:,1);

%visualization
figure;  figure_trimesh(surfout1,thick1smooth);
colormap('hot')
shading interp
view([0 90])
camlight headlight

figure;  figure_trimesh(sphere40962, thick1smooth);
colormap('hot')
shading interp
view([0 90])
camlight headlight

figure;  figure_trimesh(surfout2,thick2smooth);
colormap('hot')
shading interp
view([0 90])
camlight headlight

figure;  figure_trimesh(sphere40962, thick2smooth);
colormap('hot')
shading interp
view([0 90])
camlight headlight

% Let us transform the spherical coordinates 
% to the rectagular coordinates where y= varphi, x= theta

theta_d=0:0.01:pi;
varphi_d=0:0.01:2*pi;
m=length(theta_d);
n=length(varphi_d);

theta=kron(ones(1,n),theta_d');
theta=reshape(theta,1,m*n);

varphi=kron(ones(m,1),varphi_d);
varphi=reshape(varphi,1,m*n);
    
Ylm= Y_lm(3,1,theta,varphi);
square=reshape(Ylm,m,n);
figure;imagesc(square); colormap('hot'); colorbar
set(gcf,'Color','w')

Ylm= Y_lm(10,2,theta,varphi);
square=reshape(Ylm,m,n);
figure;imagesc(square); colormap('hot'); colorbar
set(gcf,'Color','w')


%cortical thickness in rectangular coordinates

square1=SPHARM2square(fourier1.x,42, 0.001);
square2=SPHARM2square(fourier2.x,42, 0.001);

%displacement vector between surface 1 to 2

disp.vertices=surfout2.vertices - surfout1.vertices;
disp.faces=surfout1.faces;

figure; figure_wire(surfout1, 'white', 'white');
hold on; figure_quiver2(surfout1, disp.vertices, 20,'r')













