equivalentEllipse

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 %EQUIVALENTELLIPSE 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@inrae.fr
0035 % Created: 2008-02-21, using Matlab 7.4.0.287 (R2007a)
0036 % Copyright 2008-2024 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas
0037 
0038 % ellipse center
0039 xc = mean(points(:,1));
0040 yc = mean(points(:,2));
0041 
0042 % recenter points
0043 x = points(:,1) - xc;
0044 y = points(:,2) - yc;
0045 
0046 % number of points
0047 n = size(points, 1);
0048 
0049 % equivalent parameters
0050 Ixx = sum(x.^2) / n;
0051 Iyy = sum(y.^2) / n;
0052 Ixy = sum(x.*y) / n;
0053 
0054 % compute ellipse semi-axis lengths
0055 common = sqrt( (Ixx - Iyy)^2 + 4 * Ixy^2);
0056 ra = sqrt(2) * sqrt(Ixx + Iyy + common);
0057 rb = sqrt(2) * sqrt(Ixx + Iyy - common);
0058 
0059 % compute ellipse angle in degrees
0060 theta = atan2(2 * Ixy, Ixx - Iyy) / 2;
0061 theta = rad2deg(theta);
0062 
0063 % create the resulting equivalent ellipse
0064 ell = [xc yc ra rb theta];