orientedBox

ORIENTEDBOX Minimum-width oriented bounding box of a set of points.

   OBOX = orientedBox(PTS)
   Computes the oriented bounding box of a set of points. Oriented box is
   defined by a center, two dimensions (the length and the width), and the
   orientation of the length axis. Orientation is counted in degrees, 
   counter-clockwise.

   Example
     % Draw oriented bounding box of an ellipse
     elli = [30 40 40 20 30];
     pts = ellipseToPolygon(elli, 120);
     obox = orientedBox(pts);
     figure; hold on;
     drawEllipse(elli);
     drawOrientedBox(obox, 'm');

   See also
   drawOrientedBox, orientedBoxToPolygon

0001 function obox = orientedBox(points)
0002 %ORIENTEDBOX Minimum-width oriented bounding box of a set of points.
0003 %
0004 %   OBOX = orientedBox(PTS)
0005 %   Computes the oriented bounding box of a set of points. Oriented box is
0006 %   defined by a center, two dimensions (the length and the width), and the
0007 %   orientation of the length axis. Orientation is counted in degrees,
0008 %   counter-clockwise.
0009 %
0010 %   Example
0011 %     % Draw oriented bounding box of an ellipse
0012 %     elli = [30 40 40 20 30];
0013 %     pts = ellipseToPolygon(elli, 120);
0014 %     obox = orientedBox(pts);
0015 %     figure; hold on;
0016 %     drawEllipse(elli);
0017 %     drawOrientedBox(obox, 'm');
0018 %
0019 %   See also
0020 %   drawOrientedBox, orientedBoxToPolygon
0021 %
0022 
0023 % ------
0024 % Author: David Legland
0025 % e-mail: david.legland@nantes.inra.fr
0026 % Created: 2012-03-29,    using Matlab 7.9.0.529 (R2009b)
0027 % Copyright 2012 INRA - Cepia Software Platform.
0028 
0029 
0030 %% initialisations
0031 
0032 % first, compute convex hull of the polygon
0033 inds = convhull(points(:,1), points(:,2));
0034 hull = points(inds, :);
0035 
0036 % if first and last points are the same, remove the last one
0037 if inds(1) == inds(end)
0038     hull = hull(1:end-1, :);
0039 end
0040 
0041 % compute convex hull centroid, that corresponds to approximate
0042 % location of rectangle center
0043 center = mean(hull, 1);
0044 hull = bsxfun(@minus, hull, center);
0045 
0046 % number of hull vertices
0047 nV = size(hull, 1);
0048 
0049 % default values
0050 rotatedAngle = 0;
0051 minWidth = inf;
0052 minAngle = 0;
0053 
0054 % avoid degenerated cases
0055 if nV < 3
0056     return;
0057 end
0058 
0059 % indices of vertices in extreme y directions
0060 [tmp, indA] = min(hull(:, 2)); %#ok<ASGLU>
0061 [tmp, indB] = max(hull(:, 2)); %#ok<ASGLU>
0062 
0063 caliperA = [ 1 0];    % Caliper A points along the positive x-axis
0064 caliperB = [-1 0];    % Caliper B points along the negative x-axis
0065 
0066 
0067 %% Find the direction with minimum width (rotating caliper algorithm)
0068 
0069 while rotatedAngle < pi
0070     % compute the direction vectors corresponding to each edge
0071     indA2 = mod(indA, nV) + 1;
0072     vectorA = hull(indA2, :) - hull(indA, :);
0073     
0074     indB2 = mod(indB, nV) + 1;
0075     vectorB = hull(indB2, :) - hull(indB, :);
0076     
0077     % Determine the angle between each caliper and the next adjacent edge
0078     % in the polygon
0079     angleA = vectorAngle(caliperA, vectorA);
0080     angleB = vectorAngle(caliperB, vectorB);
0081     
0082     % Determine the smallest of these angles
0083     angleIncrement = min(angleA, angleB);
0084     
0085     % Rotate the calipers by the smallest angle
0086     caliperA = rotateVector(caliperA, angleIncrement);
0087     caliperB = rotateVector(caliperB, angleIncrement);
0088     
0089     rotatedAngle = rotatedAngle + angleIncrement;
0090     
0091     % compute current width, and update opposite vertex
0092     if angleA < angleB
0093         line = createLine(hull(indA, :), hull(indA2, :));
0094         width = distancePointLine(hull(indB, :), line);
0095         indA = mod(indA, nV) + 1;
0096     
0097     else
0098         line = createLine(hull(indB, :), hull(indB2, :));
0099         width = distancePointLine(hull(indA, :), line);
0100         indB = mod(indB, nV) + 1;
0101 
0102     end
0103     
0104     % update minimum width and corresponding angle if needed
0105     if width < minWidth
0106         minWidth = width;
0107         minAngle = rotatedAngle;
0108     end
0109 end
0110 
0111 
0112 %% Compute box dimensions
0113 
0114 % orientation of the main axis
0115 theta = rad2deg(minAngle);
0116 
0117 % pre-compute trigonometric functions
0118 cot = cos(minAngle);
0119 sit = sin(minAngle);
0120 
0121 % elongation in direction of rectangle length
0122 x = hull(:,1);
0123 y = hull(:,2);
0124 x2  =   x * cot + y * sit;
0125 y2  = - x * sit + y * cot;
0126 
0127 % compute extension along main directions
0128 xmin = min(x2);    xmax = max(x2);
0129 ymin = min(y2);    ymax = max(y2);
0130 
0131 % position of the center with respect to the centroid compute before
0132 dl = (xmax + xmin)/2;
0133 dw = (ymax + ymin)/2;
0134 
0135 % change  coordinate from rectangle to user-space
0136 dx  = dl * cot - dw * sit;
0137 dy  = dl * sit + dw * cot;
0138 
0139 % coordinates of oriented box center
0140 center = center + [dx dy];
0141 
0142 % size of the rectangle
0143 rectLength  = xmax - xmin;
0144 rectWidth   = ymax - ymin;
0145 
0146 % concatenate rectangle data
0147 obox = [center rectLength rectWidth theta];
0148