function [av_ellip] = average_ellipse_vec_add(M)

maj_a = M(:,4);
min_a = M(:,5);
theta = M(:,6);
[pa] = add_lengths(maj_a, theta);
[ma] = add_lengths(min_a, theta+pi/2);
[ma] = proj_min_axes(pa, ma);

%figure(1); clf; hold on;
%plot([0 pa(1)],[0 pa(2)], 'r');
%plot([0 ma(1)],[0 ma(2)], 'b');
%axis 'equal'

theta = atan2(pa(2), pa(1));
if theta<0
    theta = theta+pi;
end
av_ellip = [M(1,1:2), mean(M(:,3)), norm(pa), norm(ma), theta];
return;

function [v] = add_lengths(a, theta)
v = [0;0];
for i=1:length(a)
    l = a(i);
    phi = theta(i);
    dv = [l*cos(phi); l*sin(phi)];
    % plot([v(1) v(1) + dv(1)],[v(2) v(2) + dv(2)]);
    % drawnow;
    v = v + dv;
end
v = v/length(a);
return;

function [ma] = proj_min_axes(pa, ma)

na = flipud(pa);
na(1) = -1*na(1);
na = na./norm(na);
P = (na*na')/(na'*na);
%plot([0 ma(1)],[0 ma(2)], 'r');
ma = P*ma;
%plot([0 ma(1)],[0 ma(2)], 'g');



return;
