1## Copyright (C) 2021 David Legland
2## All rights reserved.
3##
4## Redistribution and use in source and binary forms, with or without
5## modification, are permitted provided that the following conditions are met:
6##
7##     1 Redistributions of source code must retain the above copyright notice,
8##       this list of conditions and the following disclaimer.
9##     2 Redistributions in binary form must reproduce the above copyright
10##       notice, this list of conditions and the following disclaimer in the
11##       documentation and/or other materials provided with the distribution.
12##
13## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
14## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
15## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
16## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
17## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
18## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
19## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
20## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
21## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
22## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
23##
24## The views and conclusions contained in the software and documentation are
25## those of the authors and should not be interpreted as representing official
26## policies, either expressed or implied, of the copyright holders.
27
28function obox = orientedBox(points)
29%ORIENTEDBOX Minimum-width oriented bounding box of a set of points.
30%
31%   OBOX = orientedBox(PTS)
32%   Computes the oriented bounding box of a set of points. Oriented box is
33%   defined by a center, two dimensions (the length and the width), and the
34%   orientation of the length axis. Orientation is counted in degrees,
35%   counter-clockwise.
36%
37%   Example
38%     % Draw oriented bounding box of an ellipse
39%     elli = [30 40 40 20 30];
40%     pts = ellipseToPolygon(elli, 120);
41%     obox = orientedBox(pts);
42%     figure; hold on;
43%     drawEllipse(elli);
44%     drawOrientedBox(obox, 'm');
45%
46%   See also
47%   drawOrientedBox, orientedBoxToPolygon
48%
49
50% ------
51% Author: David Legland
52% e-mail: david.legland@nantes.inra.fr
53% Created: 2012-03-29,    using Matlab 7.9.0.529 (R2009b)
54% Copyright 2012 INRA - Cepia Software Platform.
55
56
57%% initialisations
58
59% first, compute convex hull of the polygon
60inds = convhull(points(:,1), points(:,2));
61hull = points(inds, :);
62
63% if first and last points are the same, remove the last one
64if inds(1) == inds(end)
65    hull = hull(1:end-1, :);
66end
67
68% compute convex hull centroid, that corresponds to approximate
69% location of rectangle center
70center = mean(hull, 1);
71hull = bsxfun(@minus, hull, center);
72
73% number of hull vertices
74nV = size(hull, 1);
75
76% default values
77rotatedAngle = 0;
78minWidth = inf;
79minAngle = 0;
80
81% avoid degenerated cases
82if nV < 3
83    return;
84end
85
86% indices of vertices in extreme y directions
87[tmp, indA] = min(hull(:, 2)); %#ok<ASGLU>
88[tmp, indB] = max(hull(:, 2)); %#ok<ASGLU>
89
90caliperA = [ 1 0];    % Caliper A points along the positive x-axis
91caliperB = [-1 0];    % Caliper B points along the negative x-axis
92
93
94%% Find the direction with minimum width (rotating caliper algorithm)
95
96while rotatedAngle < pi
97    % compute the direction vectors corresponding to each edge
98    indA2 = mod(indA, nV) + 1;
99    vectorA = hull(indA2, :) - hull(indA, :);
100
101    indB2 = mod(indB, nV) + 1;
102    vectorB = hull(indB2, :) - hull(indB, :);
103
104    % Determine the angle between each caliper and the next adjacent edge
105    % in the polygon
106    angleA = vectorAngle(caliperA, vectorA);
107    angleB = vectorAngle(caliperB, vectorB);
108
109    % Determine the smallest of these angles
110    angleIncrement = min(angleA, angleB);
111
112    % Rotate the calipers by the smallest angle
113    caliperA = rotateVector(caliperA, angleIncrement);
114    caliperB = rotateVector(caliperB, angleIncrement);
115
116    rotatedAngle = rotatedAngle + angleIncrement;
117
118    % compute current width, and update opposite vertex
119    if angleA < angleB
120        line = createLine(hull(indA, :), hull(indA2, :));
121        width = distancePointLine(hull(indB, :), line);
122        indA = mod(indA, nV) + 1;
123
124    else
125        line = createLine(hull(indB, :), hull(indB2, :));
126        width = distancePointLine(hull(indA, :), line);
127        indB = mod(indB, nV) + 1;
128
129    end
130
131    % update minimum width and corresponding angle if needed
132    if width < minWidth
133        minWidth = width;
134        minAngle = rotatedAngle;
135    end
136end
137
138
139%% Compute box dimensions
140
141% orientation of the main axis
142theta = rad2deg(minAngle);
143
144% pre-compute trigonometric functions
145cot = cos(minAngle);
146sit = sin(minAngle);
147
148% elongation in direction of rectangle length
149x = hull(:,1);
150y = hull(:,2);
151x2  =   x * cot + y * sit;
152y2  = - x * sit + y * cot;
153
154% compute extension along main directions
155xmin = min(x2);    xmax = max(x2);
156ymin = min(y2);    ymax = max(y2);
157
158% position of the center with respect to the centroid compute before
159dl = (xmax + xmin)/2;
160dw = (ymax + ymin)/2;
161
162% change  coordinate from rectangle to user-space
163dx  = dl * cot - dw * sit;
164dy  = dl * sit + dw * cot;
165
166% coordinates of oriented box center
167center = center + [dx dy];
168
169% size of the rectangle
170rectLength  = xmax - xmin;
171rectWidth   = ymax - ymin;
172
173% concatenate rectangle data
174obox = [center rectLength rectWidth theta];
175
176