% lecture 03 
% Sept 17, 2010
% Moo K. Chung mkchung@wisc.edu
% 
%
% Simulation given in Chung et al. 2009 IPMI and MICCAI
%-----------------------------------------------------

x=[0:0.01:1]';

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

% noise
e=normrnd(0,0.2,length(x),1);

% functional measurements
Y=mu+e;

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

%-------------------------------------
% linear regression

n=length(x);

X=[ones(n,1),x];

% We need to solve Y=X*beta by inverting the matrix.

beta =inv(X'*X)*X'*Y

Yhat =X*beta

plot(x,Yhat,'r','LineWidth',2);
set(gcf,'Color','w') 

% Quadratic fit

n=length(x);
X=[ones(n,1),x, x.^2];
beta =inv(X'*X)*X'*Y
Yhat =X*beta

hold on
plot(x,Yhat,'r','LineWidth',2);
set(gcf,'Color','w') 

% 6th order polynomial fit

n=length(x);
X=[ones(n,1),x, x.^2, x.^3, x.^4, x.^5, x.^6];
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')

% 20th order polynomial fit

n=length(x);
X=[ones(n,1)];
for i=1:20
    X=[X, x.^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')

%you are seeing the limitation of polynomial basis. polynomials
%have problem near the boundary of data. How to fix this easily?


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

%-------------------------------
% Let's change the basis
% Cosine series expansion: 
% Chung et al., 2010. Statistics and Its Interface. 3:69-80.

% define cosine basis:
psi=inline('sqrt(2)*cos(pi*l.*x)')
plot(x,psi(5,x),'ro', 'MarkerEdgeColor',[0.5 0.5 0.5], 'MarkerFaceColor',[0.5 0.5 0.5], 'MarkerSize',5)


X=[ones(n,1)];
for i=0:10
    X=[X, psi(i,x)];
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')

% The polynomial fit tends to wiggle more near the boundaries. 
% This is the limitation of polynomial basis. 
% How to fix this then? This is a possible research topic.

%--------------------------------
% spherical harmonic basis

% download lh.pial lh.sphere

% outer cortical surface
[vertices, faces] = read_surf('lh.pial');
surf.vertices=vertices;
surf.faces=faces+1; 
figure_wire(surf,'white','white');

% corresponding sphere
[vertices, faces] = read_surf('lh.sphere');
sphere.vertices=vertices;
sphere.faces=faces+1;
figure_wire(sphere,'white','white');

% two spherical angles
[theta varphi]=EULERangles(sphere);
figure_trimesh(surf,theta,'rwb');
figure_trimesh(surf,varphi,'rwb');

% spherical harmonic basis
Ylm=Y_lm(3,1,theta',varphi');
figure_trimesh(sphere,Ylm,'rwb');

% Let us transform the spherical coordinates 
% into the rectagular coordinates where y= varphi, x= theta

theta_d=0:0.01:pi;
varphi_d=0:0.01:2*pi;
m=length(theta_d);
n=length(varphi_d);

theta=kron(ones(1,n),theta_d');
figure;imagesc(theta)
% vectorzie theta
theta=reshape(theta,1,m*n);

varphi=kron(ones(m,1),varphi_d);
figure;imagesc(varphi)
% vectorize varphi
varphi=reshape(varphi,1,m*n);
    
Ylm= Y_lm(3,1,theta,varphi);
square=reshape(Ylm,m,n);
figure;imagesc(square); colormap('bone'); colorbar

figure; [C, h] = contour(square);
set(h,'ShowText','on','TextStep',get(h,'LevelStep'))
colormap('bone'); colorbar

%---------------------------------------------
% Spherical harmonic representation
% based on the iterative residual fitting algorithm.
% The detail will be studied in Lecture 4. 
[surf_smooth, fourier]=SPHARMsmooth2(surf,sphere,1,0);
figure;figure_wire(surf_smooth,'white','white');

[surf_smooth, fourier]=SPHARMsmooth2(surf,sphere,2,0);
figure;figure_wire(surf_smooth,'white','white');



