1 /*
2  * This program is free software; you can redistribute it and/or
3  * modify it under the terms of the GNU General Public License
4  * as published by the Free Software Foundation; either version 2
5  * of the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program; if not, write to the Free Software Foundation,
14  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15  */
17 /** \file \ingroup bli
18  */
21 #include <stdlib.h>
22 #include <string.h>
24 #include "MEM_guardedalloc.h"
26 #include "BLI_convexhull_2d.h"
27 #include "BLI_math.h"
28 #include "BLI_strict_flags.h"
29 #include "BLI_utildefines.h"
31 /* Copyright 2001, softSurfer (www.softsurfer.com)
32  * This code may be freely used and modified for any purpose
33  * providing that this copyright notice is included with it.
34  * SoftSurfer makes no warranty for this code, and cannot be held
35  * liable for any real or imagined damage resulting from its use.
36  * Users of this code must verify correctness for their application.
37  * http://softsurfer.com/Archive/algorithm_0203/algorithm_0203.htm
38  */
40 /** \name Main Convex-Hull Calculation
41  * \{ */
43 /**
44  * tests if a point is Left|On|Right of an infinite line.
45  *    Input:  three points P0, P1, and P2
46  * \returns > 0.0 for P2 left of the line through P0 and P1.
47  *          = 0.0 for P2 on the line.
48  *          < 0.0 for P2 right of the line.
49  */
50 static float is_left(const float p0, const float p1, const float p2)
51 {
52         return (p1 - p0) * (p2 - p0) - (p2 - p0) * (p1 - p0);
53 }
55 /**
56  * A.M. Andrew's monotone chain 2D convex hull algorithm
57  *
58  * \param points: An array of 2D points presorted by increasing x and y-coords.
59  * \param n: The number of points in points.
60  * \param r_points: An array of the convex hull vertex indices (max is n).
61  * \returns the number of points in r_points.
62  */
63 int BLI_convexhull_2d_sorted(const float (*points), const int n, int r_points[])
64 {
65         /* the output array r_points[] will be used as the stack */
66         int bot = 0;
67         int top = -1;  /* indices for bottom and top of the stack */
68         int i;  /* array scan index */
69         int minmin, minmax;
70         int maxmin, maxmax;
71         float xmax;
73         /* Get the indices of points with min x-coord and min|max y-coord */
74         float xmin = points;
75         for (i = 1; i < n; i++) {
76                 if (points[i] != xmin) {
77                         break;
78                 }
79         }
81         minmin = 0;
82         minmax = i - 1;
83         if (minmax == n - 1) {  /* degenerate case: all x-coords == xmin */
84                 r_points[++top] = minmin;
85                 if (points[minmax] != points[minmin])  /* a nontrivial segment */
86                         r_points[++top] = minmax;
87                 r_points[++top] = minmin;  /* add polygon endpoint */
88                 return top + 1;
89         }
91         /* Get the indices of points with max x-coord and min|max y-coord */
93         maxmax = n - 1;
94         xmax = points[n - 1];
95         for (i = n - 2; i >= 0; i--) {
96                 if (points[i] != xmax) {
97                         break;
98                 }
99         }
100         maxmin = i + 1;
102         /* Compute the lower hull on the stack r_points */
103         r_points[++top] = minmin;  /* push minmin point onto stack */
104         i = minmax;
105         while (++i <= maxmin) {
106                 /* the lower line joins points[minmin] with points[maxmin] */
107                 if (is_left(points[minmin], points[maxmin], points[i]) >= 0 && i < maxmin) {
108                         continue;  /* ignore points[i] above or on the lower line */
109                 }
111                 while (top > 0) {  /* there are at least 2 points on the stack */
112                         /* test if points[i] is left of the line at the stack top */
113                         if (is_left(points[r_points[top - 1]], points[r_points[top]], points[i]) > 0.0f) {
114                                 break;  /* points[i] is a new hull vertex */
115                         }
116                         else {
117                                 top--;  /* pop top point off stack */
118                         }
119                 }
121                 r_points[++top] = i;  /* push points[i] onto stack */
122         }
124         /* Next, compute the upper hull on the stack r_points above the bottom hull */
125         if (maxmax != maxmin) {  /* if distinct xmax points */
126                 r_points[++top] = maxmax;  /* push maxmax point onto stack */
127         }
129         bot = top;  /* the bottom point of the upper hull stack */
130         i = maxmin;
131         while (--i >= minmax) {
132                 /* the upper line joins points[maxmax] with points[minmax] */
133                 if (is_left(points[maxmax], points[minmax], points[i]) >= 0 && i > minmax) {
134                         continue;  /* ignore points[i] below or on the upper line */
135                 }
137                 while (top > bot) {  /* at least 2 points on the upper stack */
138                         /* test if points[i] is left of the line at the stack top */
139                         if (is_left(points[r_points[top - 1]], points[r_points[top]], points[i]) > 0.0f) {
140                                 break;  /* points[i] is a new hull vertex */
141                         }
142                         else {
143                                 top--;  /* pop top point off stack */
144                         }
145                 }
147                 if (points[i] == points[r_points] && points[i] == points[r_points]) {
148                         return top + 1;  /* special case (mgomes) */
149                 }
151                 r_points[++top] = i;  /* push points[i] onto stack */
152         }
154         if (minmax != minmin) {
155                 r_points[++top] = minmin;  /* push joining endpoint onto stack */
156         }
158         return top + 1;
159 }
161 struct PointRef {
162         const float *pt;  /* 2d vector */
163 };
165 static int pointref_cmp_yx(const void *a_, const void *b_)
166 {
167         const struct PointRef *a = a_;
168         const struct PointRef *b = b_;
170         if      (a->pt > b->pt) return  1;
171         else if (a->pt < b->pt) return -1;
173         if      (a->pt > b->pt) return  1;
174         else if (a->pt < b->pt) return -1;
176         else                          return  0;
177 }
179 /**
180  * A.M. Andrew's monotone chain 2D convex hull algorithm
181  *
182  * \param points: An array of 2D points.
183  * \param n: The number of points in points.
184  * \param r_points: An array of the convex hull vertex indices (max is n).
185  * _must_ be allocated as ``n * 2`` because of how its used internally,
186  * even though the final result will be no more than \a n in size.
187  * \returns the number of points in r_points.
188  */
189 int BLI_convexhull_2d(const float (*points), const int n, int r_points[])
190 {
191         struct PointRef *points_ref = MEM_mallocN(sizeof(*points_ref) * (size_t)n, __func__);
192         float (*points_sort) = MEM_mallocN(sizeof(*points_sort) * (size_t)n, __func__);
193         int *points_map;
194         int tot, i;
196         for (i = 0; i < n; i++) {
197                 points_ref[i].pt = points[i];
198         }
200         /* Sort the points by X, then by Y (required by the algorithm) */
201         qsort(points_ref, (size_t)n, sizeof(struct PointRef), pointref_cmp_yx);
203         for (i = 0; i < n; i++) {
204                 memcpy(points_sort[i], points_ref[i].pt, sizeof(float));
205         }
207         tot = BLI_convexhull_2d_sorted(points_sort, n, r_points);
209         /* map back to the original index values */
210         points_map = (int *)points_sort;  /* abuse float array for temp storage */
211         for (i = 0; i < tot; i++) {
212                 points_map[i] = (int)((const float(*))points_ref[r_points[i]].pt - points);
213         }
215         memcpy(r_points, points_map, (size_t)tot * sizeof(*points_map));
217         MEM_freeN(points_ref);
218         MEM_freeN(points_sort);
220         return tot;
221 }
223 /** \} */
226 /* -------------------------------------------------------------------- */
227 /* Helper functions */
229 /** \name Utility Convex-Hull Functions
230  * \{ */
232 /**
233  * \return The best angle for fitting the convex hull to an axis aligned bounding box.
234  *
235  * Intended to be used with #BLI_convexhull_2d
236  *
237  * \param points_hull: Ordered hull points
238  * (result of #BLI_convexhull_2d mapped to a contiguous array).
239  *
240  * \note we could return the index of the best edge too if its needed.
241  */
242 float BLI_convexhull_aabb_fit_hull_2d(const float (*points_hull), unsigned int n)
243 {
244         unsigned int i, i_prev;
245         float area_best = FLT_MAX;
246         float dvec_best;  /* best angle, delay atan2 */
248         i_prev = n - 1;
249         for (i = 0; i < n; i++) {
250                 const float *ev_a = points_hull[i];
251                 const float *ev_b = points_hull[i_prev];
252                 float dvec;  /* 2d rotation matrix */
254                 sub_v2_v2v2(dvec, ev_a, ev_b);
255                 if (normalize_v2(dvec) != 0.0f) {
256                         /* rotation matrix */
257                         float min = {FLT_MAX, FLT_MAX}, max = {-FLT_MAX, -FLT_MAX};
258                         unsigned int j;
259                         float area;
261                         for (j = 0; j < n; j++) {
262                                 float tvec;
263                                 mul_v2_v2_cw(tvec, dvec, points_hull[j]);
265                                 min = min_ff(min, tvec);
266                                 min = min_ff(min, tvec);
268                                 max = max_ff(max, tvec);
269                                 max = max_ff(max, tvec);
271                                 area = (max - min) * (max - min);
272                                 if (area > area_best) {
273                                         break;
274                                 }
275                         }
277                         if (area < area_best) {
278                                 area_best = area;
279                                 copy_v2_v2(dvec_best, dvec);
280                         }
281                 }
283                 i_prev = i;
284         }
286         return (area_best != FLT_MAX) ? atan2f(dvec_best, dvec_best) : 0.0f;
287 }
289 /**
290  * Wrap #BLI_convexhull_aabb_fit_hull_2d and do the convex hull calculation.
291  *
292  * \param points: arbitrary 2d points.
293  */
294 float BLI_convexhull_aabb_fit_points_2d(const float (*points), unsigned int n)
295 {
296         int *index_map;
297         int tot;
299         float angle;
301         index_map = MEM_mallocN(sizeof(*index_map) * n * 2, __func__);
303         tot = BLI_convexhull_2d(points, (int)n, index_map);
305         if (tot) {
306                 float (*points_hull);
307                 int j;
309                 points_hull = MEM_mallocN(sizeof(*points_hull) * (size_t)tot, __func__);
310                 for (j = 0; j < tot; j++) {
311                         copy_v2_v2(points_hull[j], points[index_map[j]]);
312                 }
314                 angle = BLI_convexhull_aabb_fit_hull_2d(points_hull, (unsigned int)tot);
315                 MEM_freeN(points_hull);
316         }
317         else {
318                 angle = 0.0f;
319         }
321         MEM_freeN(index_map);
323         return angle;
324 }
326 /** \} */