inertiaEllipse

INERTIAELLIPSE Inertia ellipse of a set of points.

   Note: Deprecated! Use equivalentEllipse instead.

   ELL = inertiaEllipse(PTS);
   where PTS is a N*2 array containing coordinates of N points, computes
   the inertia ellipse of the set of points.

   The result has the form:
   ELL = [XC YC A B THETA],
   with XC and YC being the center of mass of the point set, A and B are
   the lengths of the inertia ellipse (see below), and THETA is the angle
   of the main inertia axis with the horizontal (counted in degrees
   between 0 and 180). 
   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 = inertiaEllipse(pts);
   figure(1); clf; hold on;
   drawPoint(pts);
   drawEllipse(ell, 'linewidth', 2, 'color', 'r');

   See also
     equivalentEllipse

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