% Statistical and Computational Methods in Brain Image Analysis
% Chapter 01. Introduction to Brain Images
%
% (C) Moo K. Chung 
% University of Wisconsin-Madison
% mkchung@wisc.edu

%------------------------------------
% BINARY VOLUME DATA
% This dataset is published in the following paper. If you are using this
% data set, please reference the 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));

% plot for amygdala volume
figure;
plot(al,ar,'or', 'MarkerEdgeColor','k', 'MarkerFaceColor','k', 'MarkerSize',7)
hold on
plot(cl,cr,'ob', 'MarkerEdgeColor','k', 'MarkerFaceColor','w', 'MarkerSize',7)
legend('Autism','Control')
xlabel('Left Amygdala')
ylabel('Right Amygdala')
set(gcf,'Color','w');

%-------------------------------------------
% SURFACE MESH DATA
% This dataset is published in the following paper. If you are using this
% data set, please reference the 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 chung.2010.NI.mat

% TEMPLATE

% amygdala surface templates
left_template.vertices = squeeze(mean(left_surf,1));
left_template.faces=sphere.faces;

right_template.vertices = squeeze(mean(right_surf,1));
right_template.faces=sphere.faces;

% left amygdala template
figure; figure_wire(left_template,'yellow','white');
% right amygdala template
figure; figure_wire(right_template,'yellow','white');


% left_disp is a matrix of size 46 (subject number) x 2562 (vertex number)
% x 3 (x,y,z coordinates)

temp=reshape(left_template.vertices, 1, 2562,3);
temp=repmat(temp, [46 1 1]);
left_disp= left_surf - temp;

% right_disp is a matrix of size 46 (subject number) x 2562 (vertex number)
% x 3 (x,y,z coordinates)

temp=reshape(right_template.vertices, 1, 2562,3);
temp=repmat(temp, [46 1 1]);
right_disp= right_surf - temp;

% left_length is of size 46 (subject number ) x 2562 (vertex number). 
% To display the displacement vector field of subject 2, we run
left_length= sqrt(sum(left_disp.^2,3));
figure_orgami(left_template, left_length(2,:))
colormap('hot')


% right_length is of size 46 (subject number ) x 2562 (vertex number). 
% To display the displacement vector field of subject 2, we run
right_length= sqrt(sum(right_disp.^2,3));
figure_orgami(right_template, right_length(2,:))
colormap('hot')


%----------------------------------------------------
% LANDMARK DATA

load mandible-landmarks.mat

%-------------------------
% id(44) = F155-12-08 

figure; subplot(1,2,1);
hold on; figure_patch(F155_12_08, 'y', 0.5);
q=landmarks(:,:,44);
C=20*ones(24,1);
hold on; scatter3(q(:,1),q(:,2),q(:,3), C, C, 'filled', 'r')
view([90 60]);  zoom(1.2); camlight; 
title(id(44))

% id(1) = F203-01-03
subplot(1,2,2);
figure_patch(F203_01_03, [0.74 0.71 0.61], 0.5);
p=landmarks(:,:,1);
hold on; scatter3(p(:,1),p(:,2),p(:,3), C, C, 'filled', 'b')
view([90 60]);  zoom(1.2); camlight; 
title(id(1))

%--------------------------
% Affine transforms between the two sets of landamrks. 

A =affine_transform(p',q')
F203_01_03a= affine_surface(A, F203_01_03);
figure; subplot(1,2,1);
figure_patch(F203_01_03a, [0.74 0.71 0.61], 0.5);

paffine= p*A(1:3,1:3)' + repmat(A(1:3,4), 1, 24)';
hold on; scatter3(paffine(:,1),paffine(:,2),paffine(:,3), C, C, 'filled', 'b')
view([90 60]);  zoom(1.2); camlight; camlight
title('Affine transform of F203-01-03')

%-------------------------
% checking landmarks

subplot(1,2,2);
figure_patch(F155_12_08, 'y', 0.5);
hold on; scatter3(q(:,1),q(:,2),q(:,3), C, C, 'filled', 'r')
hold on; figure_patch(F203_01_03a, [0.74 0.71 0.61], 0.5);
hold on; scatter3(paffine(:,1),paffine(:,2),paffine(:,3), C, C, 'filled', 'b')
view([90 60]);  zoom(1.2); camlight; camlight
title('Landmark alignment')
    
%-------------------------------------------------------------------
% VECTOR DATA: DISPLACEMENT VECTOR FIELDS

load mandible77subjects.mat   

% subject 1
sub1.vertices=voxelsize(44)*vertices(:,:,1);
sub1.faces=faces;
figure;
subplot(2,2,1); figure_patch(sub1,[0.74 0.71 0.61],0.7);
view([90 60]); camlight
title(id(1))

% superimposition of n surfaces
subplot(2,2,2);
for i=1:20
    surf.faces=faces;
    surf.vertices=voxelsize(44)*vertices(:,:,i);
    hold on; figure_patch(surf,[0.74 0.71 0.61],0.1);
end;
view([90 60]); camlight;
title('20 surfaces')

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

% surface displacemnet from template to the subject 1
disp = sub1.vertices - template.vertices;
displength=sqrt(sum(disp.^2,2));
mean(displength)

subplot(2,2,4); ;figure_surf(template,displength);
figure_quiver3surface(template, disp, 10);
view([90 60]); camlight
title('Displacement field')
colormap('hot')

%------------------------------------------------------
% CURVE DATA

% FA map
nii = load_nii('mean_diffeomorphic_initial6_fa.nii');  
d=nii.img;
surf = isosurface(d); 

% surface coordinates has to be flipped to align properly
coord=surf.vertices;
temp=coord;
temp(:,2)=coord(:,1);
temp(:,1)=coord(:,2);
surf.vertices=temp;

figure;
figure_patch(surf,[0.5 0.5 0.5],0.2);
axis off
set(gcf,'Color','w');
hold on;


% White matter fiber tracts
load SL.mat

hold on;
for i=1:20:10000 % subsampling 1 for every 30 tracts
    tract=SL{i}';
    if size(tract,2)> 10
        hold on
        plot3(tract(1,:),tract(2,:),tract(3,:),'b')
    end;
end;






