function [data] = ellipse_fit2d(PC)
% function [data] = ellipse_fit2d(PC)
% 
% This function fits an ellipse to some data given by PC.
%
% Input Params:
%  PC - this is an Mx2 matrix of x, y coordinates that specify a region
%
% Output Params:
%  data - this is a structure that contains the information about the ellipse
%  as follows:
%
%  data.Orientation - the principle direction of the ellipse
%  data.MajorAxisLength - length along the principle direction
%  data.MinorAxisLength - length along the minor axis
%  data.Centroid - this is an 1x2 vector of x and y coordinates
%  data.Area - this is the length of the data 'PC', i.e. the number of
%              pixels for that sector.
%
% Dr. A. I. Hanna (2006).
if nargin<1
    data = [];
    return;
end
if size(PC,1)==1
    data.Orientation = 0;
    data.MajorAxisLength = 1;
    data.MinorAxisLength = 1;
    data.Centroid(1) = PC(1);
    data.Centroid(2) = PC(2);
    data.Area = 1;
    return;
end
mean_PC = mean(PC);
if length(mean_PC)==1
    mean_PC = [mean_PC, mean_PC];
end
PC = PC - ones(size(PC,1),1)*mean(PC);

% find the eigenvalues of this set of data
[V, D] = eig(PC'*PC);
% we want the vector of most variance first, i.e. the principle direction
% of growth.
[Y, I] = sort(diag(D), 'descend');
V = V(:, I);

% Now align the data with i, j, k, so we can calculate the range of the sector along
% the major and minor axes.
ProjC = (V'*PC')';
smax = max(abs(ProjC));


% theta = linspace(0, 2*pi, 100);
% x = cos(theta);
% y = sin(theta);
% if length(smax) == 1
%     C = [x(:)*smax(1) y(:)*smax(1)];
% else
%     C = [x(:)*smax(1) y(:)*smax(2)];
% end
% C = (V*C')';
% x = C(:,1); y = C(:,2);

% calculate the angle of V(:,1) with the x-axis
v_temp = V*[1;0];
theta = atan2(v_temp(2), v_temp(1));
if theta <-pi/2
    theta = theta + pi;
end
if theta>pi/2
    theta = theta - pi;
end
% Uncomment below if you want to see the ellipses being fit.
% figure(1); clf; hold on;
% plot(PC(:,1), PC(:,2), '.');
% plot([0 V(1,1)*10],[0 V(2,1)*10], 'r');
% plot([0 V(1,2)*10],[0 V(2,2)*10], 'g');
% axis image xy;
% title(sprintf('theta = %s, smax = %s, smin = %s', num2str(theta), num2str(max(smax)), num2str(min(smax))));
% drawnow;
% APC = (rotmat(-theta)*PC')';
% smax = max(abs(APC));
% if smax(2)>smax(1)
%     theta = theta + pi/2;
% end
% figure(2); clf;
% plot(PC(:,1), PC(:,2), '.')
data.Orientation = theta;
data.MajorAxisLength = smax(1);
data.MinorAxisLength = smax(2);
data.Centroid(1) = mean_PC(1);
data.Centroid(2) = mean_PC(2);
data.Area = size(PC,1);
