% lecture 04
% Sept 24, 2010
% Moo K. Chung mkchung@wisc.edu
% 
%
% 1D function simulation given in Chung et al. 2009 IPMI and MICCAI
%-----------------------------------------------------

x=[0:0.01:1]';

% unknown signal. we assumet mu is known and
% we try to estimate it.
mu=10*(x - 0.5).^2 + cos(7*pi*x)/2

figure; plot(x, mu, '-')
% noise
e=normrnd(0,0.2,length(x),1);

% functional measurements
Y=mu+e;

hold on;
plot(x,Y,'ro', 'MarkerEdgeColor',[0.5 0.5 0.5],...
    'MarkerFaceColor',[0.5 0.5 0.5], 'MarkerSize',5)
 
%-------------------------------------
% linear regression

r0=Y; % initial residual is taken as the original data.

n=length(x);
X=[ones(n,1),x-0.5];
beta1 =inv(X'*X)*X'*r0
rhat1 =X*beta1

figure;
plot(x,Y,'ro', 'MarkerEdgeColor',[0.5 0.5 0.5],...
    'MarkerFaceColor',[0.5 0.5 0.5], 'MarkerSize',5)
hold on;
plot(x,Yhat1,'r','LineWidth',2);

%-----------------------
% residual of linear fit

r1 = Y-Yhat1;

figure;
plot(x,Y,'ro', 'MarkerEdgeColor',[0.5 0.5 0.5],...
    'MarkerFaceColor',[0.5 0.5 0.5], 'MarkerSize',5)
hold on;
plot(x,r1,'r','LineWidth',2);


% Quadratic fit of residual r1

X=(x-0.5).^2;
beta2 =inv(X'*X)*X'*r1
rhat2 =X*beta2

figure;
plot(x,r1,'ro', 'MarkerEdgeColor',[0.5 0.5 0.5],...
    'MarkerFaceColor',[0.5 0.5 0.5], 'MarkerSize',5)
hold on
plot(x,rhat1,'r','LineWidth',2);
set(gcf,'Color','w') 

% Overall model of linaer + quadratic fit.

Yhat = rhat1 + rhat2;
figure;
plot(x,Y,'ro', 'MarkerEdgeColor',[0.5 0.5 0.5],...
    'MarkerFaceColor',[0.5 0.5 0.5], 'MarkerSize',5)
hold on
plot(x,Yhat,'r','LineWidth',2);
set(gcf,'Color','w') 

%-----------------------------
% residual upto quadratic fit

r2= Y- rhat1 - rhat2;
figure;
plot(x,Y,'ro', 'MarkerEdgeColor',[0.5 0.5 0.5],...
    'MarkerFaceColor',[0.5 0.5 0.5], 'MarkerSize',5)
hold on
plot(x,r2,'r','LineWidth',2);
set(gcf,'Color','w') 


%---------------------------------
% Iterative residual fitting (IRF) algorithm
% Chung et al., 2007. IEEE Transactions on Medical Imaging.
% we repeat the above process iteratively.


% polynomial based
r=Y;
beta=zeros(11,1);
X=ones(n,1);
beta =inv(X'*X)*X'*r;
rhat =X*beta;
Yhat=rhat;


for i=1:10
    r=Y-Yhat;
    
    X=(x-0.5).^i;
    %X=psi(i,x);
    beta =inv(X'*X)*X'*r;
    rhat =X*beta; %i-th order fit of the residual
    
    Yhat=Yhat+rhat;

    figure;
    plot(x,Y,'ro', 'MarkerEdgeColor',[0.5 0.5 0.5], 'MarkerFaceColor',[0.5 0.5 0.5], 'MarkerSize',5)
    hold on
    plot(x,Yhat,'r','LineWidth',2);

end;

% cosine based

psi=inline('sqrt(2)*cos(pi*l.*x)')



%-------------------------------
% compare IRF to LSE.
n=length(x);
X=[ones(n,1)];
for i=1:10
    X=[X, (x-0.5).^i];
end;
beta =inv(X'*X)*X'*Y
Yhat =X*beta

figure;
plot(x,Y,'ro', 'MarkerEdgeColor',[0.5 0.5 0.5], 'MarkerFaceColor',[0.5 0.5 0.5], 'MarkerSize',5)
hold on
plot(x,Yhat,'r','LineWidth',2);
set(gcf,'Color','w')

%--------------------------
% how the sum of squared resiaul changes over 
% the number of iterations

r=Y;
beta=zeros(11,1);
X=ones(n,1);
beta =inv(X'*X)*X'*r;
rhat =X*beta;
Yhat=rhat;

figure;
for i=1:30
    r=Y-Yhat;
    hold on;
    res=sum(r.^2);
    plot(i,res,'ro-', 'MarkerEdgeColor',[0.5 0.5 0.5], 'MarkerFaceColor',[0.5 0.5 0.5], 'MarkerSize',5)
    X=psi(i,x);
    beta =inv(X'*X)*X'*r;
    rhat =X*beta; 
    Yhat=Yhat+rhat;
end;

%-----------------------------------------
% SPHARM representation of simulated binary objects


% Cube
% This simulation is given in Chung et al. 2010. NeuroImage
vol=zeros(70,70,70);
vol(40:50,40:50,40:50)=1;
figure;imagesc(squeeze(vol(:,50,:)));colorbar

% triangulation of a binary object
surf=isosurface(vol)
figure;
figure_wire(surf,'white','white')

% lighting effect
figure;figure_patch(surf,[0.74 0.71 0.61],0.8)
daspect([1 1 1]); view(48,20); axis tight; axis vis3d off; camlight; 
lighting gouraud; 
lighting phong; material shiny; camlight; 


% sphere
vol=zeros(70,70,70);

for i=1:70
    for j=1:70
        for k=1:70
            if round((i-35)^2 + (j-35)^2 + (k-35)^2)< 10^2
                vol(i,j,k)=1;
            end;
        end;
    end;
end;

surf=isosurface(vol)
figure;
figure_wire(surf,'white','white')

% Sphere with a bump
% This simulation is given in Chung et al. 2010 NeuroImage

for i=1:70
    for j=1:70
        for k=1:70
            if round((i-30)^2 + (j-30)^2 + (k-40)^2)< 4^2
                vol(i,j,k)=1;
            end;
        end;
    end;
end;

surf=isosurface(vol)
figure;
figure_wire(surf,'white','white')

% For the SPHARM-representation, we need to flatten the surface 
% Let's flatten the cube
vol=zeros(70,70,70);
vol(40:50,40:50,40:50)=1;
surf=isosurface(vol)

% enclose the cube with a sphere and assign value 1 to the cube
% and -1 to the sphere
[amyg,sphere,amygsphere]=CREATEenclosedamyg(vol,surf);
figure;imagesc(squeeze(amygsphere(:,25,:)));colorbar

% perform heat diffusion
stream=LAPLACE3Dsmooth(amygsphere,amyg,-sphere, 10);
figure;imagesc(squeeze(stream(:,25,:)));colorbar


sphere=isosurface(amyg);
for alpha=1:5
    sphere=LAPLACEcontour(stream,sphere, 1 - 2*alpha/5);
    sphere=REGULARIZEarea(sphere, 0.5);
    figure; figure_wire(sphere,'white','white')    
end;

[surf_smooth, fourier]=SPHARMsmooth2(surf,sphere,42,0.001);
figure_wire(surf_smooth,'white','white')

%-------------------------------------------------
% The limitation of heat diffusion based flattening


vol=zeros(70,70,70);
vol(35:50,35:50,35:50)=1;
vol(40:45,40:45,40:50)=0;

surf=isosurface(vol);
figure_wire(surf,'white', 'white');


[amyg,sphere,amygsphere]=CREATEenclosedamyg(vol,surf);
figure;imagesc(squeeze(amygsphere(:,25,:)));colorbar

stream=LAPLACE3Dsmooth(amygsphere,amyg,-sphere, 10);
figure;imagesc(squeeze(stream(:,25,:)));colorbar

sphere=isosurface(amyg);
for alpha=1:5
    sphere=LAPLACEcontour(stream,sphere, 1 - 2*alpha/5);
    sphere=REGULARIZEarea(sphere, 0.5);
    figure; figure_wire(sphere,'yellow','white')    
end;

% This is a publishable research problem 
% if you can fix the algorithm. Not that hard.


%--------------------------------
% corpus callosum data generated by an undergraduate student
% read Thomas Hoffmann's B.Sc thesis written in 2004.
%
% load corpusboundary.mat
% There are 27 subjects (12 controls, 15 autistic)
% control=0; autistic=1; a(:,3)= age

a=[  1     0    15
     2     0    18
     3     0    18
     4     0    16
     5     1    15
     6     1    17
     7     1    13
     8     0    15
     9     1    12
    10     1    15
    11     0    13
    12     1    25
    13     1    14
    14     1    15
    15     0    18
    16     1    14
    17     1    24
    18     0    21
    19     1    18
    20     0    15
    21     0    17
    22     0    16
    23     1    10
    24     1    12
    25     0    23
    26     1    22
    27     1    12];


% To get the coordinates of the subject 3. 
cc=corpus(3).curveXYpair
x=cc(:,1);y=cc(:,2);
plot(x,y,'-')

figure;
for i=1:27
    cc=corpus(i).curveXYpair;
    x=cc(:,1);y=cc(:,2);
    if a(i,2)==0
        hold on;
        plot(x,y,'-');
    end
     if a(i,2)==1
        hold on;
        plot(x,y,'-r');
    end
end;

% corpus callosum curves need to be aligned (registered) 
% before any type of analysis. 






