% 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))


%------------------------------------------
% MANDIBLE GROWTH ANALYSIS

load mandible77subjects.mat   

% template surface
template.faces=faces;
template.vertices=voxelsize(44)*mean(vertices,3);
subplot(2,2,3); figure_patch(template,[0.74 0.71 0.61],0.7);
view([90 60]); camlight;
title('Template')

%-----------------------------------
% Extracting gender and age info from id

str=cell2mat(id);
group=str(:,1)
group=[group=='M']  % M=1, F=0
year = str2num(str(:,6:7));
month=str2num(str(:,9:10));
age = (12*year+month)/12;

%-----------------------------------------
% 21639 vertices are subsampled to spped up computation 

tt=reducepatch(template,0.1);
ind=mesh_commonvertex(template, tt);
p= vertices(ind,:,:);

%-------------
% age bins
age1= age<7;
age2=(7<=age)&(age<13);
age3=13<=age;

%---------------
% displacement lengths in each age bins

p1= p(:,:,age1);
p2= p(:,:,age2);
p3= p(:,:,age3);

p1L=squeeze(sqrt(sum(p1.^2,2)));
p2L=squeeze(sqrt(sum(p2.^2,2)));
p3L=squeeze(sqrt(sum(p3.^2,2)));

p12=mean(p2,3)-mean(p1,3);
p23=mean(p3,3)-mean(p2,3);

p12L=sqrt(sum(p12.^2,2));
p23L=sqrt(sum(p23.^2,2));

%-----------------------
% visualizing age effects 
s=voxelsize(44);

figure; figure_surf(tt, s*p12L);
figure_quiver3surface(tt, s*5*p12, 5);
caxis([0 10*s]);
colorbar
title('Between age group I and II')
colormap('hot')
view([90 60]); 
camlight;

figure; figure_surf(tt, s*p23L);
figure_quiver3surface(tt, s*5*p23, 5);
caxis([0 10*s]);
colorbar
title('Between age group II and III')
colormap('hot'); 
view([90 60]); 
camlight; 


%-----------------------------------------------------
%Obtain LB-eigenfunctions
% eigenvalues D is sorted in increasing order. 
% V is the eigenfunctions of the LB-operator.

[A, C]= FEM(tt);
[V, D]= eigs(C,A, 1000,'sa');  

% Displaying 5th eigenfunction
figure; figure_surf(tt,V(:,5))
view([90 60]); camlight;

%----------------------------------------
% Heat kernel smoothing on displacement lengths

p1Lsmooth=lb_smooth(p1L,tt, 20, 200, V, D);
p2Lsmooth=lb_smooth(p2L,tt, 20, 200, V, D);
p3Lsmooth=lb_smooth(p3L,tt, 20, 200, V, D);

%-------------------------
% Two sample T-test on the length difference
% Note T^2 = F test.

tfield=[];
for i=1:length(p1Lsmooth)
[h p, ci, t] = ttest2(p2Lsmooth(i,:), p1Lsmooth(i,:),[],[],'equal');
tfield=[tfield t.tstat];
end;

s=voxelsize(44); %voxel size
figure; figure_surf(tt, tfield);
figure_quiver3surface(tt, s*5*p12, 4);
view([90 60]); colormap('hot'); colorbar
title('T-stat between age group I and II')
caxis([-6 6]); colorbar; camlight

tfield=[];
for i=1:length(p1Lsmooth)
[h p, ci, t] = ttest2(p3Lsmooth(i,:), p2Lsmooth(i,:), [],[], 'equal');
tfield=[tfield t.tstat];
end;

figure; figure_surf(tt, tfield);
figure_quiver3surface(tt, s*5*p12, 4);
view([90 60]); colormap('hot'); colorbar
title('T-stat between age group II and III')
caxis([-6 6]); colorbar; camlight


%----------------------------------------------
% Hotelling's T2 result
% The resulting F-stat is the square of T-stat. 

tfield=[];
for i=1:length(p1Lsmooth)
    h= hotelT2(p2Lsmooth(i,:)', p1Lsmooth(i,:)');
    tfield=[tfield  h.t];
end;

figure; figure_surf(tt, tfield);
figure_quiver3surface(tt, s*5*p12, 4);
view([90 60]); 

caxis([0 16]); colorbar
camlight; colormap('hot')

c=colormap;
colormap(c(end:-1:1,:));
title('Hotelling T-square between age group I and II')


%----------------------------------
% RANDOM FIELD THEORY: corrected p-value for F-map

t=[9 12 16]
a=1; b=sum(age1) + sum(age2) -2;
sigma=20;

mu0=2;

% surface has to be scaled properly by the voxel size
ss.faces = tt.faces;
ss.vertices= s*tt.vertices;
area=GETarea(ss);
mu2=sum(area.faces)/2;

rho0 = 1 - fcdf(t,a,b);

rho2 = 1/(4*sigma^2*pi)*gamma((a+b-2)/2)/ gamma(a/2)/gamma(b/2)...
    *(a*t/b).^((a-2)/2).*(1+a*t/b).^(-(a+b-2)/2).*((b-1)*a*t/b-(a-1));

p=mu0*rho0+rho2*mu2

% F-stat
F.tstat=[9 12 16]
F.df= [a b];
stat=stat_Fcorrected(F, sigma,ss)


%----------------------
% p-value plot and thresholding

F.tstat=[8:0.01:20]
a=1; b=sum(age1) + sum(age2) -2;
F.df= [a b];

stat=stat_Fcorrected(F,20,ss)
figure; plot(F.tstat, stat.pcorrected, '-k', 'LIneWidth', 2)
title('Corrected p-value')
xlabel('F-statistic')
ylabel('p-value')

ind=find(stat.pcorrected<0.01);
F.tstat(ind(1))



