% Line_fit.m takes a set of N points in X(N,n).
% It then fits a line of best fit to them.
% The line is of the form X(j,i) = a(i) + b(i)*t(j) 
% The return values are b, a, t, epsilon.
% b is a unit n-vector pointing in the direction of the line.
% a is the n-centre of the line.
% epsilon is the root mean square distance from the line.
% t is a N-vector parameterizing the position of each point on the line

function [a,b,t,epsilon] = line_fit(X)
[N,n] = size(X);

     % Find a, the mean position
     a(n) = 0;
     for j=1:n
       for i=1:N
       a(j) = a(j) + X(i,j);
       end
     end
     a = a /N ;

     % Determine the n*n Correlation Matrix 
     R = -X'*X/N;
     sum_sq = 0;
     for i=1:n
         sum_sq = sum_sq - R(i,i) - a(i)*a(i);
     end
     R = R + a'*a;
     for i=1:n
         R(i,i) = R(i,i) + sum_sq;
     end
  

     % Find epsilon and b.
     [V,D] = eig(R);
     D = diag(D);
     [epsilon,ind] = min(D);
      epsilon = sqrt(epsilon);
     b = V(:,ind);
     b = b' / norm(b);

     % Finally calculate t.
     t = (X*b' - a*b');