% NIMS DTI tutorial
% 2012 July 9, 10
%
%(c) 2012 Moo K. Chung
% University of Wisconsin-Madison
% mkchung@wisc.edu
% http://www.stat.wisc.edu/~mchung
%
% This file/datat have been download from
% tezpur.keck.waisman.wisc.edu/~chung/DTI

input=load_nii('sub004_tensor_T1Space.nii');

size(input.img)
d= squeeze(input.img);
dim=size(d)

d1= d(:,:,:,1);   %xx
d2= d(:,:,:,2);   %xy
d3= d(:,:,:,3);   %yy
d4= d(:,:,:,4);   %xz
d5= d(:,:,:,5);   %yz
d6= d(:,:,:,6);   %zz

sn=85
figure;
subplot(3,3,1); imagesc(d1(:,:,sn)); colorbar
subplot(3,3,2); imagesc(d2(:,:,sn)); colorbar
subplot(3,3,5); imagesc(d3(:,:,sn)); colorbar
subplot(3,3,3); imagesc(d4(:,:,sn)); colorbar
subplot(3,3,6); imagesc(d5(:,:,sn)); colorbar
subplot(3,3,9); imagesc(d6(:,:,sn)); colorbar


% trace normalization
tr= d1 +d3 + d6;

d1=d1./tr;
d2=d2./tr;
d3=d3./tr;
d4=d4./tr;
d5=d5./tr;
d6=d6./tr;

figure;
subplot(3,3,1); imagesc(d1(:,:,sn)); colorbar
subplot(3,3,2); imagesc(d2(:,:,sn)); colorbar
subplot(3,3,5); imagesc(d3(:,:,sn)); colorbar
subplot(3,3,3); imagesc(d4(:,:,sn)); colorbar
subplot(3,3,6); imagesc(d5(:,:,sn)); colorbar
subplot(3,3,9); imagesc(d6(:,:,sn)); colorbar

%Finding the maximum eigenvalue in 2D slice (x-y plane)

d1s = reshape(d1(:,:,sn), dim(1)*dim(2),1);
d2s = reshape(d2(:,:,sn), dim(1)*dim(2),1);
d3s = reshape(d3(:,:,sn), dim(1)*dim(2),1);

for i=1:dim(1)*dim(2)
  D=[d1s(i) d2s(i)
      d2s(i) d3s(i)];
  [v,lambda]=eig(D);
  % i = 1 is not brain 
    V(i,:,:)=v;  % eigenvectors
    L(i,:)=[lambda(1,1), lambda(2,2)];  % eigenvalues
end; 

for i=1:dim(1)*dim(2)
    if (sum(L(i,:)) ~= 0)
        if (L(i,1) >= L(i,2))
            maxL(i)=L(i,1);
            maxV(i,:)=squeeze(V(i,1,:))';
        else
            maxL(i)=L(i,2);
            maxV(i,:)=squeeze(V(i,2,:))';
        end;
    end;
end;

% principal eigenvalue
maxL=reshape(maxL,dim(1), dim(2));
figure;imagesc(maxL); colorbar;


% principal eigenvector
maxVx=reshape(maxV(:,1), dim(1), dim(2));
maxVy=reshape(maxV(:,2), dim(1), dim(2));


scaleVx=maxVx.*maxL;
scaleVy=maxVy.*maxL;

hold on;quiver(scaleVy,scaleVx,'k')
colorbar

%-------------------------------------------------------------------------
% This part of tutorial explains about the basics on DTI manipulation
% based on the following publication:
% Chung, M.K., Adluru, N., Lee, J.E., Lazar, M., Lainhart, J.E., Alexander, A.L., 2010. 
% Cosine series representation of 3D curves and its application to white matter fiber 
% bundles in diffusion tensor imaging. Statistics and Its Interface. 3:69-80.
% http://www.stat.wisc.edu/%7Emchung/papers/chung.2010.SII.pdf

% FA map of a single subject
I=analyze75read('favol');
s=isosurface(I);
figure;figure_patch(s,'yellow', 0.1);hold on;
axis off
set(gcf,'Color','w') 

load bundle_tracts.mat
for i=1:50:5311 %subsample at every 50 tract
    tract=cell2mat(fiber_bundle(i));    
    hold on
    plot3(tract(:,1),tract(:,2),tract(:,3),'g')
end
view([0 90])

% cosine series representation

i=1
tract=cell2mat(fiber_bundle(i));
figure; plot3(tract(:,1),tract(:,2),tract(:,3),'.g')

[arc_length  para]=parameterize_arclength(tract');

[wfs beta] =WFS_tracts(tract',para,19);
hold on; plot3(wfs(:,1),wfs(:,2),wfs(:,3),'-r')


[wfs beta] =WFS_tracts(tract',para,30);
figure; plot3(tract(:,1),tract(:,2),tract(:,3),'.g')
hold on; plot3(wfs(:,1),wfs(:,2),wfs(:,3),'-r')


% optimal degree selection

SSE=zeros(3,50);
for i=1:50
    [wfs beta]=WFS_tracts(tract',para,i);
    SSE(:,i)= L2norm(tract'-wfs');
end;

% plot of SSE
figure;plot(SSE(1,:));
hold on;plot(SSE(2,:),'r');
hold on;plot(SSE(3,:),'g');

% Gibbs phenomenon 
x=1:300;  
y1=zeros(1,100); 
y2=ones(1,100);
y3=zeros(1,100);
y=[y1 y2 y3]
z=1:300;

figure; plot3(x, y, z,'g')

tract=[x' y' z']'
[arc_length  para]=parameterize_arclength(tract);
[wfs beta] =WFS_tracts(tract,para,200)
hold on; plot3(wfs(:,1),wfs(:,2),wfs(:,3),'-r')

%--------------------------------------------
% Rips complex demo is based on the publication:
% Lee, H., Chung, M.K., Kang, H., Kim, B.-N., Lee, D.S. 2011. Computing the
% shape of brain network using graph filtration and Gromov-Haudorff metric. 
% Lecture Notes in Computer Science (LNCS). 6892:302-309.
% http://www.stat.wisc.edu/%7Emchung/papers/lee.2011.MICCAI.pdf


endp = mouse2points(500,500);

% Rips complex construction

figure; imshow('toy-key.tif')
[adj, prob, vertices] = points2graph(endp,50, 'rips'); 
figure_graph_color(adj,vertices, ones(size(adj,1),1))

figure; imshow('toy-key.tif')
[adj, prob, vertices] = points2graph(endp,100, 'rips'); 
figure_graph_color(adj,vertices, ones(size(adj,1),1))

figure; imshow('toy-key.tif')
[adj, prob, vertices] = points2graph(endp,150, 'rips'); 
figure_graph_color(adj,vertices, ones(size(adj,1),1))

%------------------------------------------
% Epsilon-neighbor construction
% Graph construction is based on the publication:
% Chung, M.K., Adluru, N., Dalton, K.M., Alexander, A.L., Davidson, R.J. 2011. 
% Scalable brain network construction on white matter fibers. SPIE Medical Imaging.
% http://www.stat.wisc.edu/%7Emchung/papers/chung.2011.SPIE.pdf

nii=load_nii('mean_diffeomorphic_initial6_fa.nii');
d=nii.img;
surf = isosurface(d); 
coord=surf.vertices;
temp=coord;
temp(:,2)=coord(:,1);
temp(:,1)=coord(:,2);
surf.vertices=temp;
figure_patch(surf,[0.74 0.71 0.61],0.2); hold on;


load SL.mat
for i=1:30: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;

endp=get_endp_tracts(SL,30);
hold on; plot3(endp(:,1,1),endp(:,1,2), endp(:,1,3),'.r')
hold on; plot3(endp(:,2,1),endp(:,2,2), endp(:,2,3),'.r')

[adj, prob, vertices]=tract2graph(SL,10);

figure_patch(surf,[0.74 0.71 0.61],0.2); hold on;
figure_dtigraph(adj,vertices);

figure; imagesc(-adj); colormap('bone');


