function [c,p,flatness] = bestFitPlane( vxs )
%[c,p,flatness] = bestFitPlane( vxs )
%   Find the best-fit plane through the vertexes VXS, an N*D array where D
%   is the number of dimensions.  The result is returned in the form of C,
%   the centroid of the vertexes, and P, the unit normal vector to the plane
%   through C, both as row vectors.  FLATNESS is a measure of how flat the
%   set of points is, and is the ratio of the smallest to largest
%   eigenvalue of the covariance matrix.

    N = size(vxs,1);
    c = sum(vxs,1)/N;
    m = (vxs' * vxs) / N - c'*c;
    [v,d] = eig(m);
    p = v(:,1)';
    p = p/sqrt(sum(p.*p));
    flatness = d(1,1)/d(size(d,1),size(d,2));
    return;
    
    figure(1);
    clf
    hold on
    plotpts( gca, vxs, 'ko' );
    plotvecs( gca, c, c+v(:,1)', 'r' );
    axis equal
    view([35,30]);
    hold off
end
