Cleanup: improve 2D convex hull

Improve correctness, API, comments, memory usage and performance
of the 2D convex hull calculation.

Pre-requisite for UV packing improvements.

Differential Revision: https://developer.blender.org/D16055

3 files changed, 46 insertions, 67 deletions

diff --git a/source/blender/blenlib/BLI_convexhull_2d.h b/source/blender/blenlib/BLI_convexhull_2d.h
index 0b4c3d486fb..f4e4e4d66f1 100644
--- a/source/blender/blenlib/BLI_convexhull_2d.h
+++ b/source/blender/blenlib/BLI_convexhull_2d.h
@@ -11,43 +11,26 @@ extern "C" {
 #endif
 
 /**
- * A.M. Andrew's monotone chain 2D convex hull algorithm.
- *
- * \param points: An array of 2D points presorted by increasing x and y-coords.
- * \param n: The number of points in points.
- * \param r_points: An array of the convex hull vertex indices (max is n).
- * \returns the number of points in r_points.
- */
-int BLI_convexhull_2d_sorted(const float (*points)[2], int n, int r_points[]);
-/**
- * A.M. Andrew's monotone chain 2D convex hull algorithm.
+ * Extract 2D convex hull.
  *
  * \param points: An array of 2D points.
  * \param n: The number of points in points.
  * \param r_points: An array of the convex hull vertex indices (max is n).
- * _must_ be allocated as `n * 2` because of how its used internally,
- * even though the final result will be no more than \a n in size.
- * \returns the number of points in r_points.
- */
-int BLI_convexhull_2d(const float (*points)[2], int n, int r_points[]);
-
-/**
- * \return The best angle for fitting the convex hull to an axis aligned bounding box.
- *
- * Intended to be used with #BLI_convexhull_2d
+ * \return The number of indices in r_points.
  *
- * \param points_hull: Ordered hull points
- * (result of #BLI_convexhull_2d mapped to a contiguous array).
+ * \note Performance is O(n.log(n)), same as qsort.
  *
- * \note we could return the index of the best edge too if its needed.
  */
-float BLI_convexhull_aabb_fit_hull_2d(const float (*points_hull)[2], unsigned int n);
+int BLI_convexhull_2d(const float (*points)[2], int n, int r_points[/* n */]);
+
 /**
- * Wrap #BLI_convexhull_aabb_fit_hull_2d and do the convex hull calculation.
+ * \return The best angle for fitting the points to an axis aligned bounding box.
+ *
+ * \note We could return the index of the best edge too if its needed.
  *
- * \param points: arbitrary 2d points.
+ * \param points: Arbitrary 2d points.
  */
-float BLI_convexhull_aabb_fit_points_2d(const float (*points)[2], unsigned int n);
+float BLI_convexhull_aabb_fit_points_2d(const float (*points)[2], int n);
 
 #ifdef __cplusplus
 }
diff --git a/source/blender/blenlib/intern/convexhull_2d.c b/source/blender/blenlib/intern/convexhull_2d.c
index ee5d000b72f..fee6241a817 100644
--- a/source/blender/blenlib/intern/convexhull_2d.c
+++ b/source/blender/blenlib/intern/convexhull_2d.c
@@ -39,8 +39,9 @@ static float is_left(const float p0[2], const float p1[2], const float p2[2])
   return (p1[0] - p0[0]) * (p2[1] - p0[1]) - (p2[0] - p0[0]) * (p1[1] - p0[1]);
 }
 
-int BLI_convexhull_2d_sorted(const float (*points)[2], const int n, int r_points[])
+static int BLI_convexhull_2d_sorted(const float (*points)[2], const int n, int r_points[])
 {
+  BLI_assert(n >= 2); /* Doesn't handle trivial cases. */
   /* the output array r_points[] will be used as the stack */
   int bot = 0;
   int top = -1; /* indices for bottom and top of the stack */
@@ -66,6 +67,7 @@ int BLI_convexhull_2d_sorted(const float (*points)[2], const int n, int r_points
       r_points[++top] = minmax;
     }
     r_points[++top] = minmin; /* add polygon endpoint */
+    BLI_assert(top + 1 <= n);
     return top + 1;
   }
 
@@ -122,16 +124,18 @@ int BLI_convexhull_2d_sorted(const float (*points)[2], const int n, int r_points
     }
 
     if (points[i][0] == points[r_points[0]][0] && points[i][1] == points[r_points[0]][1]) {
+      BLI_assert(top + 1 <= n);
       return top + 1; /* special case (mgomes) */
     }
 
     r_points[++top] = i; /* push points[i] onto stack */
   }
 
-  if (minmax != minmin) {
+  if (minmax != minmin && r_points[0] != minmin) {
     r_points[++top] = minmin; /* push joining endpoint onto stack */
   }
 
+  BLI_assert(top + 1 <= n);
   return top + 1;
 }
 
@@ -162,35 +166,38 @@ static int pointref_cmp_yx(const void *a_, const void *b_)
 
 int BLI_convexhull_2d(const float (*points)[2], const int n, int r_points[])
 {
-  struct PointRef *points_ref = MEM_mallocN(sizeof(*points_ref) * (size_t)n, __func__);
-  float(*points_sort)[2] = MEM_mallocN(sizeof(*points_sort) * (size_t)n, __func__);
-  int *points_map;
-  int points_hull_num, i;
+  BLI_assert(n >= 0);
+  if (n < 2) {
+    if (n == 1) {
+      r_points[0] = 0;
+    }
+    return n;
+  }
+  struct PointRef *points_ref = MEM_mallocN(sizeof(*points_ref) * n, __func__);
+  float(*points_sort)[2] = MEM_mallocN(sizeof(*points_sort) * n, __func__);
 
-  for (i = 0; i < n; i++) {
+  for (int i = 0; i < n; i++) {
     points_ref[i].pt = points[i];
   }
 
-  /* Sort the points by X, then by Y (required by the algorithm) */
+  /* Sort the points by X, then by Y. */
   qsort(points_ref, (size_t)n, sizeof(struct PointRef), pointref_cmp_yx);
 
-  for (i = 0; i < n; i++) {
+  for (int i = 0; i < n; i++) {
     memcpy(points_sort[i], points_ref[i].pt, sizeof(float[2]));
   }
 
-  points_hull_num = BLI_convexhull_2d_sorted(points_sort, n, r_points);
+  int points_hull_num = BLI_convexhull_2d_sorted(points_sort, n, r_points);
 
-  /* map back to the original index values */
-  points_map = (int *)points_sort; /* abuse float array for temp storage */
-  for (i = 0; i < points_hull_num; i++) {
-    points_map[i] = (int)((const float(*)[2])points_ref[r_points[i]].pt - points);
+  /* Map back to the unsorted index values. */
+  for (int i = 0; i < points_hull_num; i++) {
+    r_points[i] = (int)((const float(*)[2])points_ref[r_points[i]].pt - points);
   }
 
-  memcpy(r_points, points_map, (size_t)points_hull_num * sizeof(*points_map));
-
   MEM_freeN(points_ref);
   MEM_freeN(points_sort);
 
+  BLI_assert(points_hull_num <= n);
   return points_hull_num;
 }
 
@@ -202,14 +209,13 @@ int BLI_convexhull_2d(const float (*points)[2], const int n, int r_points[])
 /** \name Utility Convex-Hull Functions
  * \{ */
 
-float BLI_convexhull_aabb_fit_hull_2d(const float (*points_hull)[2], uint n)
+static float BLI_convexhull_aabb_fit_hull_2d(const float (*points_hull)[2], int n)
 {
-  uint i, i_prev;
   float area_best = FLT_MAX;
   float dvec_best[2]; /* best angle, delay atan2 */
 
-  i_prev = n - 1;
-  for (i = 0; i < n; i++) {
+  int i_prev = n - 1;
+  for (int i = 0; i < n; i++) {
     const float *ev_a = points_hull[i];
     const float *ev_b = points_hull[i_prev];
     float dvec[2]; /* 2d rotation matrix */
@@ -218,10 +224,9 @@ float BLI_convexhull_aabb_fit_hull_2d(const float (*points_hull)[2], uint n)
     if (normalize_v2(dvec) != 0.0f) {
       /* rotation matrix */
       float min[2] = {FLT_MAX, FLT_MAX}, max[2] = {-FLT_MAX, -FLT_MAX};
-      uint j;
       float area;
 
-      for (j = 0; j < n; j++) {
+      for (int j = 0; j < n; j++) {
         float tvec[2];
         mul_v2_v2_cw(tvec, dvec, points_hull[j]);
 
@@ -249,32 +254,23 @@ float BLI_convexhull_aabb_fit_hull_2d(const float (*points_hull)[2], uint n)
   return (area_best != FLT_MAX) ? atan2f(dvec_best[0], dvec_best[1]) : 0.0f;
 }
 
-float BLI_convexhull_aabb_fit_points_2d(const float (*points)[2], uint n)
+float BLI_convexhull_aabb_fit_points_2d(const float (*points)[2], int n)
 {
-  int *index_map;
-  int points_hull_num;
+  float angle = 0.0f;
 
-  float angle;
+  int *index_map = MEM_mallocN(sizeof(*index_map) * n, __func__);
 
-  index_map = MEM_mallocN(sizeof(*index_map) * n * 2, __func__);
+  int points_hull_num = BLI_convexhull_2d(points, n, index_map);
 
-  points_hull_num = BLI_convexhull_2d(points, (int)n, index_map);
-
-  if (points_hull_num) {
-    float(*points_hull)[2];
-    int j;
-
-    points_hull = MEM_mallocN(sizeof(*points_hull) * (size_t)points_hull_num, __func__);
-    for (j = 0; j < points_hull_num; j++) {
+  if (points_hull_num > 1) {
+    float(*points_hull)[2] = MEM_mallocN(sizeof(*points_hull) * points_hull_num, __func__);
+    for (int j = 0; j < points_hull_num; j++) {
       copy_v2_v2(points_hull[j], points[index_map[j]]);
     }
 
-    angle = BLI_convexhull_aabb_fit_hull_2d(points_hull, (uint)points_hull_num);
+    angle = BLI_convexhull_aabb_fit_hull_2d(points_hull, points_hull_num);
     MEM_freeN(points_hull);
   }
-  else {
-    angle = 0.0f;
-  }
 
   MEM_freeN(index_map);
 
diff --git a/source/blender/python/mathutils/mathutils_geometry.c b/source/blender/python/mathutils/mathutils_geometry.c
index 28deebcf5ac..c8f9bfb0ed7 100644
--- a/source/blender/python/mathutils/mathutils_geometry.c
+++ b/source/blender/python/mathutils/mathutils_geometry.c
@@ -1465,7 +1465,7 @@ static PyObject *M_Geometry_convex_hull_2d(PyObject *UNUSED(self), PyObject *poi
     int *index_map;
     Py_ssize_t len_ret, i;
 
-    index_map = MEM_mallocN(sizeof(*index_map) * len * 2, __func__);
+    index_map = MEM_mallocN(sizeof(*index_map) * len, __func__);
 
     /* Non Python function */
     len_ret = BLI_convexhull_2d(points, len, index_map);