% lecture 07
% Moo K. Chung mkchung@wisc.edu
% 
% Best model selection

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

[arc_length  para]=parameterize_arclength(tract);


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

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');
set(gcf,'Color','white','InvertHardcopy','off');

% F stat
F=zeros(3,49);
pvalue=zeros(3,49);
for i=1:49
    F(:,i)=(SSE(:,i)-SSE(:,i+1))./(SSE(:,i)/(length(para)-i-2));
    pvalue(:,i)= 1-fcdf(F(:,i),1,length(para)-i-2);
end

figure;plot(F(1,:));
hold on;plot(F(2,:),'r');
hold on;plot(F(3,:),'k');

figure;plot(pvalue(1,:));
hold on;plot(pvalue(2,:),'r');
hold on;plot(pvalue(3,:),'k');


% Estimating normal vectors

load autism_cortical_surface.mat

surf.vertices=squeeze(autism_outer(1,:,:))';
surf.faces=tri;
figure; figure_wire(surf, 'yellow');


[nbr, degree] = FINDnbr(tri);

% neighborhood of vertex 1
nbr(1,:)'
degree(1)

nbr(1,1:degree(1)) % this will remove 0 entry at the end.

% neighborhood of vertex 100
nbr(100,:)'
degree(100)

% Neighboring vertices of vertex 1
coord=squeeze(autism_outer(1,:,:))';
center=coord(1,:);
coord_nbr = coord(nbr(1,1:degree(1)),:)

plot3(center(1),center(2),center(3),'*r') %center
hold on; 
plot3(coord_nbr(:,1), coord_nbr(:,2), coord_nbr(:,3),'*b')

C=cov(coord_nbr)
[N,gamma] = eig(C)
% translate to the origin

coord_nbr - center % this won't work. matrix dimensions don't match
%  5x3        1x3
for i=1:5
    temp(i,:)=coord_nbr(i,:)-center
end;

coord_trans= coord_nbr - repmat(center,5,1)

D=cov(coord_trans) % covariance matrix of coordinates
 % why? covariance matrix is translation invariant!

% compute normal vector
[N,gamma] = eig(C)
N(:,1)

plot3(center(1),center(2),center(3),'*r') %center
hold on; 
plot3(coord_nbr(:,1), coord_nbr(:,2), coord_nbr(:,3),'*b')
hold on;
plot3([center(1) center(1)+N(1,1)],[center(2) center(2)+N(2,1)]...
    ,[center(3) center(3)+N(3,1)],'r') % normal vector
hold on;
plot3([center(1) center(1)+N(1,2)],[center(2) center(2)+N(2,2)]...
    ,[center(3) center(3)+N(3,2)],'k')  % tangent vector
hold on;
plot3([center(1) center(1)+N(1,3)],[center(2) center(2)+N(2,3)]...
    ,[center(3) center(3)+N(3,3)],'k')  % tangent vector

% There is no HW15 and HW16. Someone asked me for 
% submitting only 6 problems by next week. 
% Simply submit the solutions to six problems HW9-14 
% by November 6 1:00pm.

%HW17. Estimate normal vectors on the above cortical surface 
%      without using PCA. Compare your method with PCA. 
%      Which one performs beter? Other than fitting a 
%      quadratic surface, why estimating surface noramls is
%      important in brain imaging applications. Discuss. 

%HW18. Either compute surface curvatures or come up with your 
%      own measure of surface bending/folding on the above cortical
%      surface. Why cortical bending/folding is biologically 
%      important? What do you think is causing the cortical folding? Discuss. 









