polygonInertiaEllipse

POLYGONINERTIAELLIPSE Compute ellipse with same inertia moments as polygon.

   ELLI = polygonInertiaEllipse(POLY)

   Example
     % convert an ellipse to polygon, and check that inertia ellipse is
     % close to original ellipse
     elli = [50 50 50 30 20];
     poly = ellipseToPolygon(elli, 1000);
     polygonInertiaEllipse(poly)
     ans =
        50.0000   50.0000   49.9998   29.9999   20.0000

     % compute inertia ellipse of more complex figure
     img = imread('circles.png');
     img = imfill(img, 'holes');
     figure; imshow(img); hold on;
     B = bwboundaries(img);
     poly = B{1}(:,[2 1]);
     drawPolygon(poly, 'r');
     elli = polygonInertiaEllipse(poly);
     drawEllipse(elli, 'color', 'g', 'linewidth', 2);


   See also
     polygons2d, polygonSecondAreaMoments, polygonCentroid, inertiaEllipse
     ellipseToPolygon

0001 function elli = polygonInertiaEllipse(poly)
0002 %POLYGONINERTIAELLIPSE Compute ellipse with same inertia moments as polygon.
0003 %
0004 %   ELLI = polygonInertiaEllipse(POLY)
0005 %
0006 %   Example
0007 %     % convert an ellipse to polygon, and check that inertia ellipse is
0008 %     % close to original ellipse
0009 %     elli = [50 50 50 30 20];
0010 %     poly = ellipseToPolygon(elli, 1000);
0011 %     polygonInertiaEllipse(poly)
0012 %     ans =
0013 %        50.0000   50.0000   49.9998   29.9999   20.0000
0014 %
0015 %     % compute inertia ellipse of more complex figure
0016 %     img = imread('circles.png');
0017 %     img = imfill(img, 'holes');
0018 %     figure; imshow(img); hold on;
0019 %     B = bwboundaries(img);
0020 %     poly = B{1}(:,[2 1]);
0021 %     drawPolygon(poly, 'r');
0022 %     elli = polygonInertiaEllipse(poly);
0023 %     drawEllipse(elli, 'color', 'g', 'linewidth', 2);
0024 %
0025 %
0026 %   See also
0027 %     polygons2d, polygonSecondAreaMoments, polygonCentroid, inertiaEllipse
0028 %     ellipseToPolygon
0029 %
0030  
0031 % ------
0032 % Author: David Legland
0033 % e-mail: david.legland@inra.fr
0034 % Created: 2017-09-08,    using Matlab 9.1.0.441655 (R2016b)
0035 % Copyright 2017 INRA - Cepia Software Platform.
0036 
0037 warning('MatGeom:deprecated', ...
0038     'function ''polygonInertiaEllipse'' is obsolete, use ''polygonEquivalentEllipse'' instead');
0039 
0040 % first re-center the polygon
0041 centroid = polygonCentroid(poly);
0042 poly = bsxfun(@minus, poly, centroid);
0043 
0044 % compute non-normalized inertia moments
0045 [Ix, Iy, Ixy] = polygonSecondAreaMoments(poly);
0046 
0047 % noralaize with polygon area
0048 area = polygonArea(poly);
0049 Ix = Ix / area;
0050 Iy = Iy / area;
0051 Ixy = Ixy / area;
0052 
0053 % compute ellipse semi-axis lengths
0054 common = sqrt( (Ix - Iy)^2 + 4 * Ixy^2);
0055 ra = sqrt(2) * sqrt(Ix + Iy + common);
0056 rb = sqrt(2) * sqrt(Ix + Iy - common);
0057 
0058 % compute ellipse angle and convert into degrees
0059 % (different formula from the inertiaEllipse function, as the definition
0060 % for Ix and Iy do not refer to same axes)
0061 theta = atan2(2 * Ixy, Iy - Ix) / 2;
0062 theta = theta * 180 / pi;
0063 
0064 % compute centroid and concatenate results into ellipse format
0065 elli = [centroid ra rb theta];