octave_packages/geometry-1.5.0/geom2d/inertiaEllipse.m

   1 %% Copyright (c) 2011, INRA
   2 %% 2008-2011, David Legland <david.legland@grignon.inra.fr>
   3 %% 2011 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
   4 %%
   5 %% All rights reserved.
   6 %% (simplified BSD License)
   7 %%
   8 %% Redistribution and use in source and binary forms, with or without
   9 %% modification, are permitted provided that the following conditions are met:
  10 %%
  11 %% 1. Redistributions of source code must retain the above copyright notice, this
  12 %%    list of conditions and the following disclaimer.
  13 %%
  14 %% 2. Redistributions in binary form must reproduce the above copyright notice,
  15 %%    this list of conditions and the following disclaimer in the documentation
  16 %%    and/or other materials provided with the distribution.
  17 %%
  18 %% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  19 %% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  20 %% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  21 %% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
  22 %% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  23 %% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  24 %% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  25 %% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  26 %% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  27 %% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  28 %% POSSIBILITY OF SUCH DAMAGE.
  29 %%
  30 %% The views and conclusions contained in the software and documentation are
  31 %% those of the authors and should not be interpreted as representing official
  32 %% policies, either expressed or implied, of copyright holder.
  33
  34 %% -*- texinfo -*-
  35 %% @deftypefn {Function File} {@var{ell} = } inertiaEllipse (@var{pts})
  36 %% Inertia ellipse of a set of points
  37 %%
  38 %%   ELL = inertiaEllipse(PTS);
  39 %%   where PTS is a N*2 array containing coordinates of N points, computes
  40 %%   the inertia ellispe of the set of points.
  41 %%
  42 %%   The result has the form:
  43 %%   ELL = [XC YC A B THETA],
  44 %%   with XC and YC being the center of mass of the point set, A and B are
  45 %%   the lengths of the inertia ellipse (see below), and THETA is the angle
  46 %%   of the main inertia axis with the horizontal (counted in degrees
  47 %%   between 0 and 180).
  48 %%   A and B are the standard deviations of the point coordinates when
  49 %%   ellipse is aligned with the inertia axes.
  50 %%
  51 %%   @example
  52 %%   pts = randn(100, 2);
  53 %%   pts = transformPoint(pts, createScaling(5, 2));
  54 %%   pts = transformPoint(pts, createRotation(pi/6));
  55 %%   pts = transformPoint(pts, createTranslation(3, 4));
  56 %%   ell = inertiaEllipse(pts);
  57 %%   figure(1); clf; hold on;
  58 %%   drawPoint(pts);
  59 %%   drawEllipse(ell, 'linewidth', 2, 'color', 'r');
  60 %% @end example
  61 %%
  62 %%   @seealso{ellipses2d, drawEllipse}
  63 %% @end deftypefn
  64
  65 function ell = inertiaEllipse(points)
  66
  67   % ellipse center
  68   xc = mean(points(:,1));
  69   yc = mean(points(:,2));
  70
  71   % recenter points
  72   x = points(:,1) - xc;
  73   y = points(:,2) - yc;
  74
  75   % number of points
  76   n = size(points, 1);
  77
  78   % inertia parameters
  79   Ixx = sum(x.^2) / n;
  80   Iyy = sum(y.^2) / n;
  81   Ixy = sum(x.*y) / n;
  82
  83   % compute ellipse semi-axis lengths
  84   common = sqrt( (Ixx - Iyy)^2 + 4 * Ixy^2);
  85   ra = sqrt(2) * sqrt(Ixx + Iyy + common);
  86   rb = sqrt(2) * sqrt(Ixx + Iyy - common);
  87
  88   % compute ellipse angle in degrees
  89   theta = atan2(2 * Ixy, Ixx - Iyy) / 2;
  90   theta = rad2deg(theta);
  91
  92   % create the resulting inertia ellipse
  93   ell = [xc yc ra rb theta];
  94
  95 endfunction
  96