Blender  V3.3
convexhull_2d.c
Go to the documentation of this file.
1 /* SPDX-License-Identifier: GPL-2.0-or-later */
2 
7 #include <stdlib.h>
8 #include <string.h>
9 
10 #include "MEM_guardedalloc.h"
11 
12 #include "BLI_convexhull_2d.h"
13 #include "BLI_math.h"
14 #include "BLI_strict_flags.h"
15 #include "BLI_utildefines.h"
16 
17 /* Copyright 2001, softSurfer (http://www.softsurfer.com)
18  * This code may be freely used and modified for any purpose
19  * providing that this copyright notice is included with it.
20  * SoftSurfer makes no warranty for this code, and cannot be held
21  * liable for any real or imagined damage resulting from its use.
22  * Users of this code must verify correctness for their application.
23  * http://softsurfer.com/Archive/algorithm_0203/algorithm_0203.htm
24  */
25 
26 /* -------------------------------------------------------------------- */
37 static float is_left(const float p0[2], const float p1[2], const float p2[2])
38 {
39  return (p1[0] - p0[0]) * (p2[1] - p0[1]) - (p2[0] - p0[0]) * (p1[1] - p0[1]);
40 }
41 
42 int BLI_convexhull_2d_sorted(const float (*points)[2], const int n, int r_points[])
43 {
44  /* the output array r_points[] will be used as the stack */
45  int bot = 0;
46  int top = -1; /* indices for bottom and top of the stack */
47  int i; /* array scan index */
48  int minmin, minmax;
49  int maxmin, maxmax;
50  float xmax;
51 
52  /* Get the indices of points with min x-coord and min|max y-coord */
53  float xmin = points[0][0];
54  for (i = 1; i < n; i++) {
55  if (points[i][0] != xmin) {
56  break;
57  }
58  }
59 
60  minmin = 0;
61  minmax = i - 1;
62  if (minmax == n - 1) { /* degenerate case: all x-coords == xmin */
63  r_points[++top] = minmin;
64  if (points[minmax][1] != points[minmin][1]) {
65  /* a nontrivial segment */
66  r_points[++top] = minmax;
67  }
68  return top + 1;
69  }
70 
71  /* Get the indices of points with max x-coord and min|max y-coord */
72 
73  maxmax = n - 1;
74  xmax = points[n - 1][0];
75  for (i = n - 2; i >= 0; i--) {
76  if (points[i][0] != xmax) {
77  break;
78  }
79  }
80  maxmin = i + 1;
81 
82  /* Compute the lower hull on the stack r_points */
83  r_points[++top] = minmin; /* push minmin point onto stack */
84  i = minmax;
85  while (++i <= maxmin) {
86  /* the lower line joins points[minmin] with points[maxmin] */
87  if (is_left(points[minmin], points[maxmin], points[i]) >= 0 && i < maxmin) {
88  continue; /* ignore points[i] above or on the lower line */
89  }
90 
91  while (top > 0) { /* there are at least 2 points on the stack */
92  /* test if points[i] is left of the line at the stack top */
93  if (is_left(points[r_points[top - 1]], points[r_points[top]], points[i]) > 0.0f) {
94  break; /* points[i] is a new hull vertex */
95  }
96  top--; /* pop top point off stack */
97  }
98 
99  r_points[++top] = i; /* push points[i] onto stack */
100  }
101 
102  /* Next, compute the upper hull on the stack r_points above the bottom hull */
103  if (maxmax != maxmin) { /* if distinct xmax points */
104  r_points[++top] = maxmax; /* push maxmax point onto stack */
105  }
106 
107  bot = top; /* the bottom point of the upper hull stack */
108  i = maxmin;
109  while (--i >= minmax) {
110  /* the upper line joins points[maxmax] with points[minmax] */
111  if (is_left(points[maxmax], points[minmax], points[i]) >= 0 && i > minmax) {
112  continue; /* ignore points[i] below or on the upper line */
113  }
114 
115  while (top > bot) { /* at least 2 points on the upper stack */
116  /* test if points[i] is left of the line at the stack top */
117  if (is_left(points[r_points[top - 1]], points[r_points[top]], points[i]) > 0.0f) {
118  break; /* points[i] is a new hull vertex */
119  }
120  top--; /* pop top point off stack */
121  }
122 
123  if (points[i][0] == points[r_points[0]][0] && points[i][1] == points[r_points[0]][1]) {
124  return top + 1; /* special case (mgomes) */
125  }
126 
127  r_points[++top] = i; /* push points[i] onto stack */
128  }
129 
130  if (minmax != minmin) {
131  r_points[++top] = minmin; /* push joining endpoint onto stack */
132  }
133 
134  return top + 1;
135 }
136 
137 struct PointRef {
138  const float *pt; /* 2d vector */
139 };
140 
141 static int pointref_cmp_yx(const void *a_, const void *b_)
142 {
143  const struct PointRef *a = a_;
144  const struct PointRef *b = b_;
145 
146  if (a->pt[1] > b->pt[1]) {
147  return 1;
148  }
149  if (a->pt[1] < b->pt[1]) {
150  return -1;
151  }
152 
153  if (a->pt[0] > b->pt[0]) {
154  return 1;
155  }
156  if (a->pt[0] < b->pt[0]) {
157  return -1;
158  }
159  return 0;
160 }
161 
162 int BLI_convexhull_2d(const float (*points)[2], const int n, int r_points[])
163 {
164  struct PointRef *points_ref = MEM_mallocN(sizeof(*points_ref) * (size_t)n, __func__);
165  float(*points_sort)[2] = MEM_mallocN(sizeof(*points_sort) * (size_t)n, __func__);
166  int *points_map;
167  int points_hull_num, i;
168 
169  for (i = 0; i < n; i++) {
170  points_ref[i].pt = points[i];
171  }
172 
173  /* Sort the points by X, then by Y (required by the algorithm) */
174  qsort(points_ref, (size_t)n, sizeof(struct PointRef), pointref_cmp_yx);
175 
176  for (i = 0; i < n; i++) {
177  memcpy(points_sort[i], points_ref[i].pt, sizeof(float[2]));
178  }
179 
180  points_hull_num = BLI_convexhull_2d_sorted(points_sort, n, r_points);
181 
182  /* map back to the original index values */
183  points_map = (int *)points_sort; /* abuse float array for temp storage */
184  for (i = 0; i < points_hull_num; i++) {
185  points_map[i] = (int)((const float(*)[2])points_ref[r_points[i]].pt - points);
186  }
187 
188  memcpy(r_points, points_map, (size_t)points_hull_num * sizeof(*points_map));
189 
190  MEM_freeN(points_ref);
191  MEM_freeN(points_sort);
192 
193  return points_hull_num;
194 }
195 
198 /* Helper functions */
199 
200 /* -------------------------------------------------------------------- */
204 float BLI_convexhull_aabb_fit_hull_2d(const float (*points_hull)[2], unsigned int n)
205 {
206  unsigned int i, i_prev;
207  float area_best = FLT_MAX;
208  float dvec_best[2]; /* best angle, delay atan2 */
209 
210  i_prev = n - 1;
211  for (i = 0; i < n; i++) {
212  const float *ev_a = points_hull[i];
213  const float *ev_b = points_hull[i_prev];
214  float dvec[2]; /* 2d rotation matrix */
215 
216  sub_v2_v2v2(dvec, ev_a, ev_b);
217  if (normalize_v2(dvec) != 0.0f) {
218  /* rotation matrix */
219  float min[2] = {FLT_MAX, FLT_MAX}, max[2] = {-FLT_MAX, -FLT_MAX};
220  unsigned int j;
221  float area;
222 
223  for (j = 0; j < n; j++) {
224  float tvec[2];
225  mul_v2_v2_cw(tvec, dvec, points_hull[j]);
226 
227  min[0] = min_ff(min[0], tvec[0]);
228  min[1] = min_ff(min[1], tvec[1]);
229 
230  max[0] = max_ff(max[0], tvec[0]);
231  max[1] = max_ff(max[1], tvec[1]);
232 
233  area = (max[0] - min[0]) * (max[1] - min[1]);
234  if (area > area_best) {
235  break;
236  }
237  }
238 
239  if (area < area_best) {
240  area_best = area;
241  copy_v2_v2(dvec_best, dvec);
242  }
243  }
244 
245  i_prev = i;
246  }
247 
248  return (area_best != FLT_MAX) ? atan2f(dvec_best[0], dvec_best[1]) : 0.0f;
249 }
250 
251 float BLI_convexhull_aabb_fit_points_2d(const float (*points)[2], unsigned int n)
252 {
253  int *index_map;
254  int points_hull_num;
255 
256  float angle;
257 
258  index_map = MEM_mallocN(sizeof(*index_map) * n * 2, __func__);
259 
260  points_hull_num = BLI_convexhull_2d(points, (int)n, index_map);
261 
262  if (points_hull_num) {
263  float(*points_hull)[2];
264  int j;
265 
266  points_hull = MEM_mallocN(sizeof(*points_hull) * (size_t)points_hull_num, __func__);
267  for (j = 0; j < points_hull_num; j++) {
268  copy_v2_v2(points_hull[j], points[index_map[j]]);
269  }
270 
271  angle = BLI_convexhull_aabb_fit_hull_2d(points_hull, (unsigned int)points_hull_num);
272  MEM_freeN(points_hull);
273  }
274  else {
275  angle = 0.0f;
276  }
277 
278  MEM_freeN(index_map);
279 
280  return angle;
281 }
282 
typedef float(TangentPoint)[2]
MINLINE float max_ff(float a, float b)
MINLINE float min_ff(float a, float b)
MINLINE void mul_v2_v2_cw(float r[2], const float mat[2], const float vec[2])
MINLINE void copy_v2_v2(float r[2], const float a[2])
MINLINE void sub_v2_v2v2(float r[2], const float a[2], const float b[2])
MINLINE float normalize_v2(float r[2])
Strict compiler flags for areas of code we want to ensure don't do conversions without us knowing abo...
_GL_VOID GLfloat value _GL_VOID_RET _GL_VOID const GLuint GLboolean *residences _GL_BOOL_RET _GL_VOID GLsizei GLfloat GLfloat GLfloat GLfloat const GLubyte *bitmap _GL_VOID_RET _GL_VOID GLenum const void *lists _GL_VOID_RET _GL_VOID const GLdouble *equation _GL_VOID_RET _GL_VOID GLdouble GLdouble blue _GL_VOID_RET _GL_VOID GLfloat GLfloat blue _GL_VOID_RET _GL_VOID GLint GLint blue _GL_VOID_RET _GL_VOID GLshort GLshort blue _GL_VOID_RET _GL_VOID GLubyte GLubyte blue _GL_VOID_RET _GL_VOID GLuint GLuint blue _GL_VOID_RET _GL_VOID GLushort GLushort blue _GL_VOID_RET _GL_VOID GLbyte GLbyte GLbyte alpha _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble alpha _GL_VOID_RET _GL_VOID GLfloat GLfloat GLfloat alpha _GL_VOID_RET _GL_VOID GLint GLint GLint alpha _GL_VOID_RET _GL_VOID GLshort GLshort GLshort alpha _GL_VOID_RET _GL_VOID GLubyte GLubyte GLubyte alpha _GL_VOID_RET _GL_VOID GLuint GLuint GLuint alpha _GL_VOID_RET _GL_VOID GLushort GLushort GLushort alpha _GL_VOID_RET _GL_VOID GLenum mode _GL_VOID_RET _GL_VOID GLint GLsizei GLsizei GLenum type _GL_VOID_RET _GL_VOID GLsizei GLenum GLenum const void *pixels _GL_VOID_RET _GL_VOID const void *pointer _GL_VOID_RET _GL_VOID GLdouble v _GL_VOID_RET _GL_VOID GLfloat v _GL_VOID_RET _GL_VOID GLint GLint i2 _GL_VOID_RET _GL_VOID GLint j _GL_VOID_RET _GL_VOID GLfloat param _GL_VOID_RET _GL_VOID GLint param _GL_VOID_RET _GL_VOID GLdouble GLdouble GLdouble top
Read Guarded memory(de)allocation.
SIMD_FORCE_INLINE btScalar angle(const btVector3 &v) const
Return the angle between this and another vector.
Definition: btVector3.h:356
static int pointref_cmp_yx(const void *a_, const void *b_)
float BLI_convexhull_aabb_fit_points_2d(const float(*points)[2], unsigned int n)
int BLI_convexhull_2d_sorted(const float(*points)[2], const int n, int r_points[])
Definition: convexhull_2d.c:42
float BLI_convexhull_aabb_fit_hull_2d(const float(*points_hull)[2], unsigned int n)
int BLI_convexhull_2d(const float(*points)[2], const int n, int r_points[])
static float is_left(const float p0[2], const float p1[2], const float p2[2])
Definition: convexhull_2d.c:37
void(* MEM_freeN)(void *vmemh)
Definition: mallocn.c:27
void *(* MEM_mallocN)(size_t len, const char *str)
Definition: mallocn.c:33
#define atan2f(x, y)
Definition: metal/compat.h:227
static unsigned a[3]
Definition: RandGen.cpp:78
static void area(int d1, int d2, int e1, int e2, float weights[2])
static const pxr::TfToken b("b", pxr::TfToken::Immortal)
#define min(a, b)
Definition: sort.c:35
const float * pt
float max