% Statistical and Computational Methods in Brain Image Analysis
% Chapter 08. Cortical Surface Analysis
%
% (C) Moo K. Chung 
% University of Wisconsin-Madison
% mkchung@wisc.edu


% hippocampus surface
load hippocampus.mat;

figure;figure_wire(hippoleft,'k','w') 
view([140 10])

% surface normal computation
coord=hippoleft.vertices';
tri=hippoleft.faces;

path(path, './toolbox_graph');
[n nf]=compute_normal(coord',tri);


% normal vector computation at the first vertex
i_point=1;
[nbr, degree] = FINDnbr(tri);
nbr_list=nbr(i_point,1:degree(i_point));
coord_nbr = coord(:,nbr_list)-kron(ones(1,degree(i_point)),coord(:,i_point));

deno= sqrt(n(1)^2+n(2)^2);
if deno < 0.1
    Q=eye(3);
else
    Q1 = [[n(1)/deno,n(2)/deno,0];
        [-n(2)/deno,n(1)/deno,0];
        [0,     0,      1]];
    Q2 = [[n(3),0, -deno];
        [0,      1,      0];
        [deno,0,n(3)]];
    Q=Q2*Q1;
end;
%if Q*n(:,1) = (0,0,1), the algorithm is correct

% rotate the surface coordinates by Q
rot_coord_nbr = Q*coord_nbr;

z = [rot_coord_nbr(3,:)]';
x = [rot_coord_nbr(1,:)]';
y = [rot_coord_nbr(2,:)]';

% least squares estimation
X=[x, y, x.^2, 2*x.*y, y.^2]; 
beta = (X\z)'

