equivalentEllipse

 Equivalent ellipse of a set of points.

   ELL = equivalentEllipse(PTS);
   Computes the ellips with the same moments up to the second order as the
   set of points specified by the N-by-2 array PTS.

   The result has the following form:
   ELL = [XC YC A B THETA],
   with XC and YC being the center of mass of the point set, A and B being
   the lengths of the equivalent ellipse (see below), and THETA being the
   angle of the first principal axis with the horizontal (counted in
   degrees between 0 and 180 in counter-clockwise direction). 
   A and B are the standard deviations of the point coordinates when
   ellipse is aligned with the principal axes.

   Example
     pts = randn(100, 2);
     pts = transformPoint(pts, createScaling(5, 2));
     pts = transformPoint(pts, createRotation(pi/6));
     pts = transformPoint(pts, createTranslation(3, 4));
     ell = equivalentEllipse(pts);
     figure(1); clf; hold on;
     drawPoint(pts);
     drawEllipse(ell, 'linewidth', 2, 'color', 'r');

   See also
     ellipses2d, drawEllipse, equivalentEllipsoid, principalAxes,
     principalAxesTransform

0001 function ell = equivalentEllipse(points)
0002 % Equivalent ellipse of a set of points.
0003 %
0004 %   ELL = equivalentEllipse(PTS);
0005 %   Computes the ellips with the same moments up to the second order as the
0006 %   set of points specified by the N-by-2 array PTS.
0007 %
0008 %   The result has the following form:
0009 %   ELL = [XC YC A B THETA],
0010 %   with XC and YC being the center of mass of the point set, A and B being
0011 %   the lengths of the equivalent ellipse (see below), and THETA being the
0012 %   angle of the first principal axis with the horizontal (counted in
0013 %   degrees between 0 and 180 in counter-clockwise direction).
0014 %   A and B are the standard deviations of the point coordinates when
0015 %   ellipse is aligned with the principal axes.
0016 %
0017 %   Example
0018 %     pts = randn(100, 2);
0019 %     pts = transformPoint(pts, createScaling(5, 2));
0020 %     pts = transformPoint(pts, createRotation(pi/6));
0021 %     pts = transformPoint(pts, createTranslation(3, 4));
0022 %     ell = equivalentEllipse(pts);
0023 %     figure(1); clf; hold on;
0024 %     drawPoint(pts);
0025 %     drawEllipse(ell, 'linewidth', 2, 'color', 'r');
0026 %
0027 %   See also
0028 %     ellipses2d, drawEllipse, equivalentEllipsoid, principalAxes,
0029 %     principalAxesTransform
0030 %
0031 
0032 % ------
0033 % Author: David Legland
0034 % e-mail: david.legland@inra.fr
0035 % Created: 2008-02-21,    using Matlab 7.4.0.287 (R2007a)
0036 % Copyright 2008 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.
0037 
0038 % HISTORY
0039 % 2009-07-29 take into account ellipse orientation
0040 % 2011-03-12 rewrite using equivalent moments
0041 
0042 % ellipse center
0043 xc = mean(points(:,1));
0044 yc = mean(points(:,2));
0045 
0046 % recenter points
0047 x = points(:,1) - xc;
0048 y = points(:,2) - yc;
0049 
0050 % number of points
0051 n = size(points, 1);
0052 
0053 % equivalent parameters
0054 Ixx = sum(x.^2) / n;
0055 Iyy = sum(y.^2) / n;
0056 Ixy = sum(x.*y) / n;
0057 
0058 % compute ellipse semi-axis lengths
0059 common = sqrt( (Ixx - Iyy)^2 + 4 * Ixy^2);
0060 ra = sqrt(2) * sqrt(Ixx + Iyy + common);
0061 rb = sqrt(2) * sqrt(Ixx + Iyy - common);
0062 
0063 % compute ellipse angle in degrees
0064 theta = atan2(2 * Ixy, Ixx - Iyy) / 2;
0065 theta = rad2deg(theta);
0066 
0067 % create the resulting equivalent ellipse
0068 ell = [xc yc ra rb theta];