polygonEquivalentEllipse

 Compute equivalent ellipse with same second order moments as polygon.

   ELLI = polygonEquivalentEllipse(POLY)

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

     % compute equivalent 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 = polygonEquivalentEllipse(poly);
     drawEllipse(elli, 'color', 'g', 'linewidth', 2);


   See also
     polygons2d, polygonSecondAreaMoments, polygonCentroid,
     equivalentEllipse, ellipseToPolygon

0001 function elli = polygonEquivalentEllipse(poly)
0002 % Compute equivalent ellipse with same second order moments as polygon.
0003 %
0004 %   ELLI = polygonEquivalentEllipse(POLY)
0005 %
0006 %   Example
0007 %     % convert an ellipse to polygon, and check that it equivalent ellipse
0008 %     is close to original ellipse
0009 %     elli = [50 50 50 30 20];
0010 %     poly = ellipseToPolygon(elli, 1000);
0011 %     polygonEquivalentEllipse(poly)
0012 %     ans =
0013 %        50.0000   50.0000   49.9998   29.9999   20.0000
0014 %
0015 %     % compute equivalent 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 = polygonEquivalentEllipse(poly);
0023 %     drawEllipse(elli, 'color', 'g', 'linewidth', 2);
0024 %
0025 %
0026 %   See also
0027 %     polygons2d, polygonSecondAreaMoments, polygonCentroid,
0028 %     equivalentEllipse, 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 % first re-center the polygon
0038 centroid = polygonCentroid(poly);
0039 poly = bsxfun(@minus, poly, centroid);
0040 
0041 % compute non-normalized inertia moments
0042 [Ix, Iy, Ixy] = polygonSecondAreaMoments(poly);
0043 
0044 % noralaize with polygon area
0045 area = polygonArea(poly);
0046 Ix = Ix / area;
0047 Iy = Iy / area;
0048 Ixy = Ixy / area;
0049 
0050 % compute ellipse semi-axis lengths
0051 common = sqrt( (Ix - Iy)^2 + 4 * Ixy^2);
0052 ra = sqrt(2) * sqrt(Ix + Iy + common);
0053 rb = sqrt(2) * sqrt(Ix + Iy - common);
0054 
0055 % compute ellipse angle and convert into degrees
0056 % (different formula from the equivalentEllipse function, as the definition
0057 % for Ix and Iy do not refer to same axes)
0058 theta = atan2(2 * Ixy, Iy - Ix) / 2;
0059 theta = theta * 180 / pi;
0060 
0061 % compute centroid and concatenate results into ellipse format
0062 elli = [centroid ra rb theta];