% MIDUS DATA SET
% (C) Moo K. Chung 2013
% University of Wisconsin-Madison
% mkchung@wisc.edu
% Image registration and template construction was done by S.-G. Kim of SNU
%
% Version 1.0 created August 9, 2013

%----------------------
% Visualizing amygdala and hippocampus template for MIDUS data set

load midus.mat

amygl=template{1,1}
amygr=template{1,2}
hippol=template{2,1}
hippor=template{2,2}

figure;
hold on; figure_wire(amygl,'black','red');
hold on; figure_wire(amygr,'black','blue');
hold on; figure_wire(hippol,'black','yellow');
hold on; figure_wire(hippor,'black','green');
view([-120 15])

%--------------------------------------------------------
% Getting individual surface out of the template
% Template + displacement gives the surface.

% displaying the 10th subject

n=10
amygl_ind.vertices = amygl.vertices + displacement{n,1,1};
amygl_ind.faces=amygl.faces;

amygr_ind.vertices = amygr.vertices + displacement{n,1,2};
amygr_ind.faces=amygr.faces;

hippol_ind.vertices = hippol.vertices + displacement{n,2,1};
hippol_ind.faces= hippol.faces;

hippor_ind.vertices = hippor.vertices + displacement{n,2,2};
hippor_ind.faces=hippor.faces;

figure;
hold on; figure_wire(amygl_ind,'black','red');
hold on; figure_wire(amygr_ind,'black','blue');
hold on; figure_wire(hippol_ind,'black','yellow');
hold on; figure_wire(hippor_ind,'black','green');
view([-120 15])

%-----------------------------------------------
% SAVING DATA INTO TEXT FORMAT

temp=amygl.vertices
save('amygl.vertices','temp', '-ascii')

temp=amygl.faces
save('amygl.faces','temp', '-ascii')

temp=amygr.vertices
save('amygr.vertices','temp', '-ascii')

temp=amygr.faces
save('amygr.faces','temp', '-ascii')

temp=hippol.vertices
save('hippol.vertices','temp', '-ascii')

temp=hippol.faces
save('hippol.faces','temp', '-ascii')

temp=hippor.vertices
save('hippor.vertices','temp', '-ascii')

temp=hippor.faces
save('hippor.faces','temp', '-ascii')


for i=1:69
    temp=displacement{i,1,1};
    strname=strcat('amygl.disp', num2str(i));
    save(strname, 'temp', '-ascii')
end;

for i=1:69
    temp=displacement{i,1,2};
    strname=strcat('amygr.disp', num2str(i));
    save(strname, 'temp', '-ascii')
end;

for i=1:69
    temp=displacement{i,2,1};
    strname=strcat('hippol.disp', num2str(i));
    save(strname, 'temp', '-ascii')
end;

for i=1:69
    temp=displacement{i,2,2};
    strname=strcat('hippor.disp', num2str(i));
    save(strname, 'temp', '-ascii')
end;

temp='age'
save('age.txt', 'temp', '-ascii')

temp='brain'
save('brain.txt', 'temp', '-ascii')

temp='gender'
save('gender.txt', 'temp', '-ascii')

temp='isfemale'
save('isfemale.txt', 'temp', '-ascii')


%-----------------------------------
% Area preserving flattening

figure; subplot(1,2,1); figure_wire(amygl,'r','w');
sphal=mesh_areapreserve(amygl,10000);
subplot(1,2,2); figure_wire(sphal,'r','w');

figure; subplot(1,2,1); figure_wire(amygr,'b','w');
sphar=mesh_areapreserve(amygr,10000);
subplot(1,2,2); figure_wire(sphar,'b','w');

figure; subplot(1,2,1); figure_wire(hippol,'y','w');
sphhl=mesh_areapreserve(hippol,10000);
subplot(1,2,2); figure_wire(sphhl,'y','w');

figure; subplot(1,2,1); figure_wire(hippor,'g','w');
sphhr=mesh_areapreserve(hippor,10000);
subplot(1,2,2); figure_wire(sphar,'g','w');


%--------------------------------------
% SPHARM representation

figure; subplot(1,2,1); figure_wire(amygl,'r','w');
[SPHARMal, fourier]=SPHARMsmooth2(amygl,sphal,20,0);
subplot(1,2,2); figure_wire(SPHARMal,'r','w');


figure; subplot(1,2,1); figure_wire(amygr,'b','w');
[SPHARMar, fourier]=SPHARMsmooth2(amygr,sphar,20,0);
subplot(1,2,2); figure_wire(SPHARMar,'b','w');


figure; subplot(1,2,1); figure_wire(hippol,'y','w');
[SPHARMhl, fourier]=SPHARMsmooth2(hippol,sphhl,20,0);
subplot(1,2,2); figure_wire(SPHARMhl,'y','w');

figure; subplot(1,2,1); figure_wire(hippor,'g','w');
[SPHARMhr, fourier]=SPHARMsmooth2(hippor,sphhr,20,0);
subplot(1,2,2); figure_wire(SPHARMhr,'g','w');

%------------------------------------------------------------
% A better SPAHRM representation with better area 
% preserving mapping and resampling


figure; subplot(1,3,1); figure_wire(hippol,'y','w');
sphhl=mesh_areapreserve2(hippol,1000);
subplot(1,3,2); figure_wire(sphhl,'y','w');
[SPHARMhl, fourier]=SPHARMsmooth2(hippol,sphhl,40,0);
subplot(1,3,3); figure_wire(SPHARMhl,'y','w');

% resampling spherical mesh
sphere6= sphere_tri('ico',6,1);
figure; subplot(1,3,1);  figure_wire(sphere6, 'y', 'w');

[S2hlres, ind] = sphere_resample(sphhl, sphere6);
hippol_res.faces=S2hlres.faces;
hippol_res.vertices=hippol.vertices(ind,:);
subplot(1,3,2);  figure_wire(hippol_res, 'y', 'w');

[hippolSP, fourier]=SPHARMsmooth2(hippol_res,S2hlres,40,0);
subplot(1,3,3); figure_wire(hippolSP,'y','w');


% resampling hippocampus surface
sphere4= sphere_tri('ico',4,1)

SPHARMsurf1=SPHARMrepresent2(sphere4,fourier,40,0)
figure; subplot(1,2,1); figure_wire(SPHARMsurf1,'y','w');

SPHARMsurf2=SPHARMrepresent2(sphhl,fourier,40,0)
subplot(1,2,2); figure_wire(SPHARMsurf2,'y','w');

%--------------------------------------
% Reconstruction Error Plot via MSE
% Here we plot the errors between 


residual=hippol.vertices - SPHARMsurf2.vertices;
MSE=sqrt(sum(residual.^2, 2));

figure; figure_surf(hippol , MSE)
view([-120 15])
colorbar; camlight; alpha(0.8)
title('Reconstruction Error')


%----------------------------------------
% Statistical Analysis with SurfStat
% It is based on Chung, M.K., Worsley, K.J., Nacewicz, B.M., Dalton, K.M., Davidson, R.J. 2010. 
% General multivariate linear modeling of surface shapes using SurfStat. NeuroImage. 53:491-505.  


% On the template surface, we need to test for gender and age effect of total 69 surfaces. This is
% equivalent to testing the corresponding 69 displacement vector fields. 
% The displacement of the left hippocampus are stored in displacement{:,2,1};

% Testing gender effect
% There is no gender effect on the left hippocampus

Gender=term(gender);
Age=term(age);
Brain=term(brain);

disphl=zeros(69,2444,3); % displacement vector fields for 69 subjects at 2444 mesh vertices 
for i=1:69
    disphl(i,:,:)=displacement{i,2,1};
end

% This is SurfStat mesh format
hippol.coord= hippol.vertices';
hippol.tri=hippol.faces';

slm1 = SurfStatLinMod( disphl, 1+Brain + Gender+Age, hippol);
slm2 = SurfStatT( slm1, gender);
figure; figure_surf( hippol ,slm2.t);
colorbar; view([-120 20]); camlight

% t-stat threshold of 4.73 corresponds to the multiple comparisons corrected pvalue of 0.05
thresholdhl=randomfield_threshold(slm, 0.05) 

%-------------------------
% Testing age effect
% There is no gender effect on the left hippocampus


slm3 = SurfStatT( slm1, age);
figure; figure_surf( hippol ,slm3.t);
colorbar; view([-120 20]); camlight




