% Statistical and Computational Methods in Brain Image Analysis
% Chapter 13 Laplace-Beltrami Eigenfunctions for Surface Data
%
% (C) Moo K. Chung
% University of Wisconsin-Madison
% mkchung@wisc.edu


% heat kernel smoothing in 2D image slice
L=[300 300];
for j=0:5
    l=[j j];
    psi=eigenfunction_box(l,L);
    figure; imagesc(psi); colorbar;
end;

A = imread('deform.tif');
A=double(A);
imagesc(A);colorbar

% Fourier expansion
expansion=zeros(300,300);
fourier=zeros(100,100);
L=[300 300]
degree=100;

for j=0:degree-1
    for k=0:degree-1
        l=[j k]; 
        psi=eigenfunction_box(l,L);
        coeff=sum(sum(psi.*A))/300^2; 
        fourier(j+1,k+1) = coeff;
        expansion=expansion + coeff*psi; 
    end;
end;
figure; imagesc(expansion); colorbar


% heat kernel expansion
sigma=10;
degree=100;
expansion=zeros(300,300);

for j=0:degree-1
    for k=0:degree-1
        l=[j k];
        psi=eigenfunction_box(l,L);
        lambda= (j*pi/L(1))^2 + (k*pi/L(2))^2;
        coeff = fourier(j+1,k+1);
        expansion=expansion + exp(-sigma*lambda)*coeff*psi; 
    end;
end;
figure; imagesc(expansion); colorbar


%-----------------------------------------------------
% heat kernel smoothing on hippocampus

load hippocampus.mat;

figure;figure_wire(hippoleft,'k','w') 
view([140 10])

[A, C] =FEM(hippoleft);
[V, D] = eigs(C,A,1000,'sm'); 

%heat kernel smoothing on hippocampus
hippolefts = lb_smooth([],hippoleft, 0.2, 1000, V, D);
figure;figure_wire(hippolefts,'k','w') 
view([140 10])


%heat kernel smoothing on function defined on hippocampus
signal=hippoleft.vertices(:,2)
signal = signal + normrnd(0,10,2338,1);
figure; figure_trimesh(hippoleft,signal)
caxis([90 150]); colormap('hot')
view([140 10]); camlight

smoothed = lb_smooth(signal,hippoleft, 0.2, 1000, V, D);
figure; figure_trimesh(hippolefts,smoothed);
caxis([90 150]); colormap('hot')
view([140 10]); camlight

%------------------------------------------------
% validation against spherical harmonics

sphere = sphere_tri('ico',5,1);
figure;figure_wire(sphere, 'black', 'white');
[A, C]= FEM(sphere);
[V, D]= eigs(C,A,30,'sm');
sort(diag(D))






