polygonSymmetryAxis

POLYGONSYMMETRYAXIS Try to identify symmetry axis of polygon.

   LINE = polygonSymmetryAxis(POLY)
   Returns a line that minimize difference between the polygon POLY and
   its reflection with the line.
   The difference metric between the two polygons is the sum of distances
   between each vertex of original polygon to the reflected polygon.

   Example
     % identify symmetry axis of an ellipse
     elli = [50 50 40 20 30];
     poly = ellipseToPolygon(elli, 100);
     line = polygonSymmetryAxis(poly);
     figure; hold on;
     drawEllipse(elli);
     axis equal; axis ([0 100 0 100]);
     drawLine(line);

   See also
   transforms2d, transformPoint, distancePointPolygon

0001 function axis = polygonSymmetryAxis(poly)
0002 %POLYGONSYMMETRYAXIS Try to identify symmetry axis of polygon.
0003 %
0004 %   LINE = polygonSymmetryAxis(POLY)
0005 %   Returns a line that minimize difference between the polygon POLY and
0006 %   its reflection with the line.
0007 %   The difference metric between the two polygons is the sum of distances
0008 %   between each vertex of original polygon to the reflected polygon.
0009 %
0010 %   Example
0011 %     % identify symmetry axis of an ellipse
0012 %     elli = [50 50 40 20 30];
0013 %     poly = ellipseToPolygon(elli, 100);
0014 %     line = polygonSymmetryAxis(poly);
0015 %     figure; hold on;
0016 %     drawEllipse(elli);
0017 %     axis equal; axis ([0 100 0 100]);
0018 %     drawLine(line);
0019 %
0020 %   See also
0021 %   transforms2d, transformPoint, distancePointPolygon
0022  
0023 % ------
0024 % Author: David Legland
0025 % e-mail: david.legland@nantes.inra.fr
0026 % Created: 2015-04-28,    using Matlab 8.4.0.150421 (R2014b)
0027 % Copyright 2015 INRA - Cepia Software Platform.
0028 
0029 % start by centering the polygon
0030 center = polygonCentroid(poly);
0031 poly = bsxfun(@minus, poly, center);
0032 
0033 % first performs a rough search with 8 angles
0034 initAngles = linspace(0, pi, 9);
0035 initAngles(end) = [];
0036 initRes = zeros(8, 1);
0037 for i = 1:8
0038     line = createLine([0 0 cos(initAngles(i)) sin(initAngles(i))]);
0039     rotMat = createLineReflection(line);
0040     polyRot = transformPoint(poly, rotMat);
0041     initRes(i) = sum(distancePointPolygon(poly, polyRot).^2);
0042 end
0043 
0044 % keep the angle that gives best result
0045 [dummy, indMin] = min(initRes); %#ok<ASGLU>
0046 initAngle = initAngles(indMin);
0047 
0048 % Compute angle that best fit between polygon and its symmetric along line
0049 thetaMin = fminbnd(...
0050     @(theta) sum(distancePointPolygon(transformPoint(poly, createLineReflection([0 0 cos(theta) sin(theta)])), poly).^2), ...
0051     initAngle-pi/4, initAngle+pi/4);
0052 
0053 % format as a line
0054 axis = [center cos(thetaMin) sin(thetaMin)];