minimumCaliperDiameter

MINIMUMCALIPERDIAMETER Minimum caliper diameter of a set of points.

   WIDTH = minimumCaliperDiameter(POINTS)
   Computes the minimum width of a set of points. As polygons and
   polylines are represented as point lists, this function works also for
   polygons and polylines.

   [WIDTH THETA] = minimumCaliperDiameter(POINTS)
   Also returns the direction of minimum width. The direction corresponds
   to the horizontal angle of the edge that minimizes the width. THETA is
   given in radians, between 0 and PI.


   Example
   % Compute minimal caliper diameter, and check coords of rotated points
   % have expected extent
     points = randn(30, 2);
     [width theta] = minimumCaliperDiameter(points);
     points2 = transformPoint(points, createRotation(-theta));
     diff = max(points2) - min(points2);
     abs(width - diff(2)) < 1e-10
     ans =
         1

   References
   Algorithms use rotating caliper. Implementation was based on that of
   Wikipedia:
   http://en.wikipedia.org/wiki/Rotating_calipers

   See also
   polygons2d, convexHull, orientedBox

0001 function [min_width, min_angle] = minimumCaliperDiameter(points)
0002 %MINIMUMCALIPERDIAMETER Minimum caliper diameter of a set of points.
0003 %
0004 %   WIDTH = minimumCaliperDiameter(POINTS)
0005 %   Computes the minimum width of a set of points. As polygons and
0006 %   polylines are represented as point lists, this function works also for
0007 %   polygons and polylines.
0008 %
0009 %   [WIDTH THETA] = minimumCaliperDiameter(POINTS)
0010 %   Also returns the direction of minimum width. The direction corresponds
0011 %   to the horizontal angle of the edge that minimizes the width. THETA is
0012 %   given in radians, between 0 and PI.
0013 %
0014 %
0015 %   Example
0016 %   % Compute minimal caliper diameter, and check coords of rotated points
0017 %   % have expected extent
0018 %     points = randn(30, 2);
0019 %     [width theta] = minimumCaliperDiameter(points);
0020 %     points2 = transformPoint(points, createRotation(-theta));
0021 %     diff = max(points2) - min(points2);
0022 %     abs(width - diff(2)) < 1e-10
0023 %     ans =
0024 %         1
0025 %
0026 %   References
0027 %   Algorithms use rotating caliper. Implementation was based on that of
0028 %   Wikipedia:
0029 %   http://en.wikipedia.org/wiki/Rotating_calipers
0030 %
0031 %   See also
0032 %   polygons2d, convexHull, orientedBox
0033 %
0034 
0035 % ------
0036 % Author: David Legland
0037 % e-mail: david.legland@grignon.inra.fr
0038 % Created: 2011-04-08,    using Matlab 7.9.0.529 (R2009b)
0039 % Copyright 2011 INRA - Cepia Software Platform.
0040 
0041 % first, compute convex hull of the polygon
0042 inds = convhull(points(:,1), points(:,2));
0043 hull = points(inds, :);
0044 
0045 % if first and last points are the same, remove the last one
0046 if inds(1) == inds(end)
0047     hull = hull(1:end-1, :);
0048 end
0049 
0050 % number of hull vertices
0051 nV = size(hull, 1);
0052 
0053 % default values
0054 rotated_angle = 0;
0055 min_width = inf;
0056 min_angle = 0;
0057 
0058 % avoid degenerated cases
0059 if nV < 3
0060     return;
0061 end
0062 
0063 [tmp, p_a] = min(hull(:, 2)); %#ok<ASGLU>
0064 [tmp, p_b] = max(hull(:, 2)); %#ok<ASGLU>
0065 
0066 caliper_a = [ 1 0];    % Caliper A points along the positive x-axis
0067 caliper_b = [-1 0];    % Caliper B points along the negative x-axis
0068 
0069 while rotated_angle < pi
0070     % compute the direction vectors corresponding to each edge
0071     ind_a2 = mod(p_a, nV) + 1;
0072     vector_a = hull(ind_a2, :) - hull(p_a, :);
0073     
0074     ind_b2 = mod(p_b, nV) + 1;
0075     vector_b = hull(ind_b2, :) - hull(p_b, :);
0076     
0077     % Determine the angle between each caliper and the next adjacent edge
0078     % in the polygon
0079     angle_a = vectorAngle(caliper_a, vector_a);
0080     angle_b = vectorAngle(caliper_b, vector_b);
0081     
0082     % Determine the smallest of these angles
0083     minAngle = min(angle_a, angle_b);
0084     
0085     % Rotate the calipers by the smallest angle
0086     caliper_a = rotateVector(caliper_a, minAngle);
0087     caliper_b = rotateVector(caliper_b, minAngle);
0088     
0089     rotated_angle = rotated_angle + minAngle;
0090     
0091     % compute current width, and update opposite vertex
0092     if angle_a < angle_b
0093         line = createLine(hull(p_a, :), hull(ind_a2, :));
0094         width = distancePointLine(hull(p_b, :), line);
0095         p_a = mod(p_a, nV) + 1;
0096     
0097     else
0098         line = createLine(hull(p_b, :), hull(ind_b2, :));
0099         width = distancePointLine(hull(p_a, :), line);
0100         p_b = mod(p_b, nV) + 1;
0101 
0102     end
0103     
0104     % update minimum width and corresponding angle if needed
0105     if width < min_width
0106         min_width = width;
0107         min_angle = rotated_angle;
0108     end
0109 
0110 end
0111