% lecture 06
% Moo K. Chung mkchung@wisc.edu
% 
% Fourier Analysis on white matter tracts

id=textread('REPORT.tracts.sept.2009.txt','%s')
% tracts for the 1st subject
SL=get_streamlines(id{1},[1.5 1.75 2.25]);

tract=SL{1000}'; %1000th tract out of total 10000 tracts
figure; plot3(tract(1,:),tract(2,:),tract(3,:),'.b')

% map x,y,z coordinates onto the unit interval [0,1]
[arc_length  para]=parameterize_arclength(tract);
figure;
subplot(3,1,1); plot(para, tract(1,:));
subplot(3,1,2); plot(para, tract(2,:));
subplot(3,1,3); plot(para, tract(3,:));

%cosine series representation in [0,1]
[wfs beta]=WFS_tracts(tract,para,2);
hold on; plot3(wfs(1,:),wfs(2,:),wfs(3,:),'k')

% Simulation
% Let's simulate white matter fiber tracts
% Deterministic simulation

t=0:0.1:10
tract=[t.*sin(t); t.*cos(t); t];
figure; plot3(tract(1,:),tract(2,:),tract(3,:),'.b')
%
[arc_length  para]=parameterize_arclength(tract);  %computing arclength
[wfs beta]=WFS_tracts(tract,para,19);  %computing cosine series expansion
hold on; plot3(wfs(1,:),wfs(2,:),wfs(3,:),'r')

% Stocastic simulation
% Let's add noise and generate 20 more tracts
% We are tring to generate similarly shaped random tracts
% Let tract1 be a clinical population.
t=0:0.1:10
figure;
view(3);
tract1=zeros(3,101,20);
for i=1:20
    tract1(:,:,i)=[t.*sin(t+normrnd(0,0.1)); t.*cos(t+normrnd(0,0.1)); t+normrnd(0,0.1)];
    hold on; plot3(tract1(1,:,i),tract1(2,:,i),tract1(3,:,i),'--r');
end;

% Additional 20 more random tracts will serve as a normal population.
tract2=zeros(3,101,20);
for i=1:20
    tract2(:,:,i)=[(t+normrnd(0,0.5)).*sin(t+0.1); (t+normrnd(0,0.5)).*cos(t-0.1); t-0.1];
    hold on; plot3(tract2(1,:,i),tract2(2,:,i),tract2(3,:,i),'--b');
end;
set(gcf,'Color','white','InvertHardcopy','off');

% Let's see if Hotelling's T-square test procedure on the Fourier 
% coefficients will be able to discriminate the two groups.

BETA1=zeros(20,3,20); %clinical
BETA2=zeros(20,3,20); %normal

for i=1:20
    [arc_length  para]=parameterize_arclength(tract1(:,:,i));
    [wfs beta] =WFS_tracts(tract1(:,:,i),para,19);
    BETA1(:,:,i)=beta;
    
    [arc_length  para]=parameterize_arclength(tract2(:,:,i));
    [wfs beta] =WFS_tracts(tract2(:,:,i),para,19);
    BETA2(:,:,i)=beta;
end;

figure;errorbar(mean(BETA1(:,1,:),3),std(BETA1(:,1,:),0,3))
figure;errorbar(mean(BETA2(:,1,:),3),std(BETA2(:,1,:),0,3))

[h p]=hotelT2(BETA1,BETA2);
figure; bar (-log10(p))


%HW11. Clearly describe what Hotelling's T-square test procedure is. 
%      What is the underlying statistical assumptions in Hotelling's T-square
%      test procedure? Are these assumptions satisified in the above
%      example? Why we are not using the two-sample t-test in Chung et al., 2010 
%      and the above example. 
%     
%      Sample covariance matrix needs to be computed for Hotelling's
%      T-square statistic. Explain what covariance matrix is. Covariance
%      matrix are usually modelled as Wishart distribution.

%HW12. Literature review. Given the 3D graph model of structural
%      connectivity, discuss how others performed the network complexity
%      analsysis. A certain clinical population may have over- or 
%      under-connectivity of brain compared normal contorls. Discuss how
%      you will test over- and over-connectivity hypothesis.
%      Then using the 3D graph model given in NIP.lecture05.dti.m,
%      determine which subject has more complex brain network.

%      Hint: There can be many solutions. One approach is to use the idea of
%            small-world connectivity. 

%  HW9-16 due on November 6. 1:00pm. You can still happily drop the course.


% For lecture 07, please read marrelec.2006.ni.pdf. This explains about 
% Wishart distribution and how it is used in brain connectivity analysis.
% Please read jian.2007.ipmi.pdf. This will prepare you for Peter Kim's 
% talk on October 23. Note: assigned papers will be covered during in-class 
% exam on November 27 (worth 20%). 

%HW13. Explain clearly what Wishart distribution is. Randomly simulate 1000 
%      Wishart distribution in MATLAB. Explain where Wishart distribution 
%      can be used in brain imaging.

%HW14. Intercorrelation analysis of cortical thickness. Following Lerch et
%      al. (2006. NeuroImage), compute the intercorrelation of cortical 
%      thickness and display the result using our autism data set.
%      You can choose any seed vortex to construct the correlation map.
%      Determine if there is any correlation difference between the groups.

%Hint: Use corr command






