% Statistical and Computational Methods in Brain Image Analysis
% Chapter 12. Multivariate Surface Shape Analysis
%
% (C) Moo K. Chung 
% University of Wisconsin-Madison
% mkchung@wisc.edu

% This chaper is based on the dataset published in the following paper
%
% 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. 

load amygdala.volume.mat

% amygdala surface visualization
vol=squeeze(leftvol(1,:,:,:));
surf = isosurface(vol);
figure;figure_wire(surf,'yellow', 'white')

vol=squeeze(rightvol(1,:,:,:));
surf = isosurface(vol);
hold on ;figure_wire(surf,'green', 'white')
 
% volume computation
left=zeros(22,1);
right=zeros(22,1);
for i=1:22
    left(i)=sum(sum(sum(squeeze(leftvol(i,:,:,:)))));
    right(i)=sum(sum(sum(squeeze(rightvol(i,:,:,:)))));
end;

al= left(find(group));
ar=right(find(group));

cl= left(find(~group));
cr=right(find(~group));

%-----------------------------------------
% SURFACE FLATTENING: CUBE

% cube simulation
vol=zeros(70,70,70);
vol(40:50,40:50,40:50)=1;
%figure;imagesc(squeeze(vol(:,50,:)));colorbar

surf=isosurface(vol)
figure;figure_wire(surf,'yellow', 'white')
axis off

% enclose cube with a sphere
[amyg,sphere,amygsphere]=CREATEenclosedamyg(vol,surf);
figure;imagesc(squeeze(amygsphere(:,25,:)));
colorbar; colormap('hot')

% solve heat diffusion
n_steps=5;
stream=LAPLACE3Dsmooth(amygsphere,amyg,-sphere, n_steps);
figure;imagesc(squeeze(stream(:,25,:)));
colorbar; colormap('hot')

% follow stremlines
sphere=isosurface(amyg);
for alpha=1:n_steps
    sphere=LAPLACEcontour(stream,sphere, 1 - 2*alpha/n_steps);
    sphere=REGULARIZEarea(sphere, 0.5);
    figure;figure_wire(sphere,'yellow', 'white'); axis off
end;


%--------------------------------------------------
% SURFACE FLATTENING: AMYGDALA
 
 
 vol=squeeze(leftvol(1,:,:,:));
 surf = isosurface(vol);
 figure;figure_wire(surf,'yellow', 'white')
 
 [amyg,sphere,amygsphere]=CREATEenclosedamyg(vol,surf);
 figure;imagesc(squeeze(amygsphere(:,25,:)));colorbar
 
 n_steps=5;
 stream=LAPLACE3Dsmooth(amygsphere,amyg,-sphere, n_steps);
 figure;imagesc(squeeze(stream(:,25,:)));colorbar
 
 sphere=isosurface(amyg);
 for alpha=1:n_steps
     sphere=LAPLACEcontour(stream,sphere, 1 - 2*alpha/n_steps);
     sphere=REGULARIZEarea(sphere, 0.5);
     figure;figure_wire(sphere,'yellow', 'white')
 end;

 