%------------------------------------------
% HIPPOCAMPUS SHAPE MODELING
%
% Lecture 8. October 29, 2010
% (c) Moo K. Chung, mkchung@wisc.edu

% template MRI 
fileinfo = analyze75info('Template_0000');
template = analyze75read(fileinfo);
figure;imagesc(squeeze(template(152,:,:))); colormap('bone')
figure;imagesc(squeeze(template(:,152,:))); colormap('bone')
figure;imagesc(squeeze(template(:,:,152))); colormap('bone')
axis on
set(gcf,'Color','w');

% hippocampus manual masking
% left hippocampus located at  152, 116, 105
fileinfo = analyze75info('Hippo_Masks');
vol = analyze75read(fileinfo);

figure;
imagesc(squeeze(vol(:,152,:)));
load colormap_rwb.mat
colormap(colormap_rwb)


E=squeeze(template(:,152,:));
E= E/10;
BW1 = edge(E); %,'prewitt');
hold on; 
h=imagesc(-1*BW1);
set(h, 'AlphaData', 0.3);
axis off
set(gcf,'Color','w');


figure;
imagesc(squeeze(vol(110,:,:)));
load colormap_rwb.mat
colormap(colormap_rwb)


E=squeeze(template(110,:,:));
E= E/10;
BW1 = edge(E); %,'prewitt');
hold on; 
h=imagesc(-1*BW1);
set(h, 'AlphaData', 0.3);
axis off
set(gcf,'Color','w');

%-------------------------
% hippocampus surface model
load hippocampus.mat   
% contains left and right hippocampi surface


figure;figure_patch(hippoleft,[0.7 0.7 0.6],0.5)
hold on;
figure_patch(hipporight,[0.7 0.7 0.9],0.5)
axis on
% right hippo is blueish.


% image has to be flipped vertically
faces=hippoleft.faces;
vertices=-hippoleft.vertices;
%temp=vertices(:,3);
%vertices(:,3)=vertices(:,1);
%vertices(:,1)=temp;
hippoleft.vertices=vertices;

figure;%figure_wire(hippoleft, 'white', 'white');
h = patch('faces', faces, 'vertices', vertices, ...
	'FaceColor', 'w', 'EdgeColor', 'k', 'FaceAlpha', 1); % transparancy can be adjusted by FaceAlpha
axis image; box on; view(3);
set(gcf,'Color','w');
view(146,24)


%------------------------------
% perform finite element method

[A, C] =FEM(hippoleft);
figure;imagesc(-abs(A)); colorbar, colormap('bone')
set(gcf,'Color','w') 

figure; imagesc(-abs(C)); colorbar, colormap('bone')
set(gcf,'Color','w') 
%print -r300 -dtiff sparseC

[V,D] = eigs(C,A, 100,'sm');
figure; plot (sort(diag(D)));

% visualize first 10 basis

for i=100:-1:90
    figure;
    figure_trimesh(hippoleft,V(:,i),'rwb');
    caxis([-0.04 0.04]);
    view(146,24)
    camlight('headlight'); 
    %print('-dtiff', '-r300',strcat('hippo-eigen', num2str(i),'.tif'))
end;

%------------------------------------
% 2nd basis of elongated objects

figure;
hold on; figure_trimesh(hippoleft,V(:,99),'rwb');
colormap('default')
caxis([-0.04 0.04]);
view(146,24)
camlight('headlight');
%print('-dtiff', '-r300',strcat('hippo-eigen', num2str(199),'.tif'))


% mandible

load ('mandible_volume');
mand_vol = reducevolume(mand_vol_fixed,[10 10 10]);
mand_iso = isosurface(mand_vol);
clear mand_vol_fixed

figure;
hold on; figure_wire(mand_iso, 'white', 'white');

[A, C] =FEM(mand_iso);
[V,D] = eigs(C,A, 2,'sm'); 

figure;
figure_trimesh(mand_iso,V(:,1),'rwb');
colormap('default')
caxis([-0.04 0.04]);
view(50,0)
camlight('headlight');
%print('-dtiff', '-r300',strcat('mand-eigen', num2str(199),'.tif'))


%---------------------------
% skeleton extraction in mandible

psi=V(:,1);

figure;
figure_trimesh(mand_iso,psi,'rwb');
colormap('default')
caxis([-0.04 0.04]);
view(50,0)
camlight('headlight');


coord=mand_iso.vertices;
maxp=max(psi);
minp=min(psi);

interval=maxp-minp;
ring_size=50;
ring=[];
alpha(0.5)
for i=0:ring_size
    h=minp + interval/ring_size*i;
    ind=find(abs(psi - h) < interval/ring_size);
    
    if length(ind)~=0
        c=mean(coord(ind,:),1);
        hold on
        %plot3(c(1),c(2),c(3),'.k', 'MarkerSize',2));
        
              scatter3(c(1),c(2),c(3),15,'k', 'filled')
        ring=[ring; c];
    end;
end;

%print('-dtiff', '-r300',strcat('hippo-eigen', num2str(199),'.tif'))


%-----------------------------------------------
% 3D deformation fields projected on the template


fname_x='Econ_1001_T1_InverseWarpxvec.nii' 
fname_y='Econ_1001_T1_InverseWarpyvec.nii' 
fname_z='Econ_1001_T1_InverseWarpzvec.nii' 

x = load_nii(fname_x);
ximg=x.img;
y = load_nii(fname_y);
yimg=y.img;
z = load_nii(fname_z);
zimg=z.img;

k=5
[x,y]=meshgrid(1:k:256,1:k:256);

figure;imagesc(squeeze(template(:,152,:))); %colormap('bone')
hold on;
quiver(x,y, squeeze(ximg(1:k:256,152,1:k:256)),...
    squeeze(zimg(1:k:256,152,1:k:256)), ...
    'LineWidth',1.5, 'Color','b' );
colormap('bone')
set(gcf,'Color','w');
axis off;
%print('-dtiff', '-r300',strcat('deform', num2str(2),'.tif'))

%----------------------------------------------
% visualization of displacement vector fields on surface

load left_disp.mat
coord=squeeze(left_disp(2,:,:)); coord=0.5*coord';
figure;figure_patch(hippoleft,[0.9 0.9 0.9],0.3)
hold on;
n=1:4:2238
quiver3(vertices(n,1), vertices(n,2), vertices(n,3), ...
    coord( 1, n ), coord( 2, n ), coord( 3, n ), 0, ...
    'LineWidth',1, 'Color','k' );
box on; view(3)
view(146,24)

%print -r300 -dtiff hippo-gray7

%---------------------------------------------
% after interpolating 3D vector fields onto the surface
% subject 1 displacement on the template

faces=hippoleft.faces;
vertices=hippoleft.vertices;
temp=vertices(:,3);
vertices(:,3)=vertices(:,2);
vertices(:,2)=temp;
hippoleft.vertices=vertices;

figure;figure_wire(hippoleft, 'white', 'white');
h = patch('faces', faces, 'vertices', vertices, ...
	'FaceColor', 'w', 'EdgeColor', 'k', 'FaceAlpha', 1); % transparancy can be adjusted by FaceAlpha
axis image; box on; view(3);


disp_h=left_disp(1:86,:,:);


coord=squeeze(mean(left_disp,1)); coord=coord';

hold on;
quiver3(vertices(:,1), vertices(:,2), vertices(:,3), ...
    coord( 1, : ), coord( 2, : ), coord( 3, : ), 0, ...
    'LineWidth',2, 'Color','yellow' );


%----------------------------------------------------
% smooth out surface coordinates

f=hippoleft.vertices(:,2);

figure;
figure_trimesh(hippoleft,f,'rwb');
colormap('default');
%caxis([100 140]);
view(146,24)
camlight('headlight');
%print('-dtiff', '-r300',strcat('hippo-smooth', num2str(1),'.tif'))

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

hippolefts = lb_smooth([],hippoleft,0.01, 200, V, D);

f=hippolefts.vertices(:,2);
figure;
figure_trimesh(hippolefts,f,'rwb');
colormap('default');
%caxis([100 140]);
view(146,24)
camlight('headlight');
%print('-dtiff', '-r300',strcat('hippo-smooth', num2str(6),'.tif'))

