BLI_convexhull_2d: optimize AABB fitting method

Previously the bounding box was calculated for each rotation,
even though the loop would early exit if the area exceeded the best
known area, this meant the areas of the bounds was checked for
each point too. This scaled poorly, close to O(n^2).

Replace this with simple logic that keeps track of the indices
at the bounds, advancing them forward for each rotation.

This is around ~140x faster for a hull with 10k points on my system.

Note that there is potential for other improvements although
the cases where this new method doesn't perform well are quite specific
and faster than checking all points, see code-comments for details.

source/blender/blenlib/intern/convexhull_2d.cc

@@ -13,11 +13,14 @@
 
 #include "MEM_guardedalloc.h"
 
+#include "BLI_bounds.hh"
 #include "BLI_convexhull_2d.h"
 #include "BLI_math_vector.h"
 #include "BLI_strict_flags.h"
 #include "BLI_utildefines.h"
 
+using namespace blender;
+
 /* Copyright 2001, softSurfer (http://www.softsurfer.com)
  * This code may be freely used and modified for any purpose
  * providing that this copyright notice is included with it.
@@ -208,43 +211,130 @@ int BLI_convexhull_2d(const float (*points)[2], const int points_num, int r_poin
 /** \name Utility Convex-Hull Functions
  * \{ */
 
+static float mul_v2_v2_cw_x(const float2 &mat, const float2 &vec)
+{
+  return (mat[0] * vec[0]) + (mat[1] * vec[1]);
+}
+
+static float mul_v2_v2_cw_y(const float2 &mat, const float2 &vec)
+{
+  return (mat[1] * vec[0]) - (mat[0] * vec[1]);
+}
+
+/**
+ * When using the rotating calipers, step one half of the caliper to a new index.
+ *
+ * Note that this relies on `points_hull` being ordered CCW which #BLI_convexhull_2d ensures.
+ */
+template<int Axis, int AxisSign>
+static float convexhull_2d_compute_extent_on_axis(const float (*points_hull)[2],
+                                                  const int points_hull_num,
+                                                  const float2 &dvec,
+                                                  int *index_p)
+{
+  /* NOTE(@ideasman42): This could be optimized to use a search strategy
+   * that computes the upper bounds and narrows down the result instead of
+   * simply checking every point until the new maximum is reached.
+   * From looking into I couldn't find cases where doing significant benefits,
+   * especially when compared with the complexity of using more involved logic for
+   * the common case, where only a few steps are needed.
+   * Typically the number of vertices to scan is small (around [0..8]).
+   * And while a high-detail hull with single outliner points will cause stepping over
+   * many more points, in practice there are rarely more than a few of these in a convex-hull.
+   * Nevertheless, a high-poly hull that has subtle curves containing many points as well as
+   * some sharp-corners wont perform as well with this method. */
+
+  const int index_init = *index_p;
+  int index_best = index_init;
+  float value_init = (Axis == 0) ? mul_v2_v2_cw_x(dvec, points_hull[index_best]) :
+                                   mul_v2_v2_cw_y(dvec, points_hull[index_best]);
+  float value_best = value_init;
+  /* Simply scan up the array. */
+  for (int count = 1; count < points_hull_num; count++) {
+    const int index_test = (index_init + count) % points_hull_num;
+    const float value_test = (Axis == 0) ? mul_v2_v2_cw_x(dvec, points_hull[index_test]) :
+                                           mul_v2_v2_cw_y(dvec, points_hull[index_test]);
+    if ((AxisSign == -1) ? (value_test > value_best) : (value_test < value_best)) {
+      break;
+    }
+    value_best = value_test;
+    index_best = index_test;
+  }
+
+  *index_p = index_best;
+  return value_best;
+}
+
 float BLI_convexhull_aabb_fit_hull_2d(const float (*points_hull)[2], int points_hull_num)
 {
   float area_best = FLT_MAX;
-  float dvec_best[2]; /* best angle, delay atan2 */
+  float2 dvec_best; /* Track the best angle as a unit vector, delaying `atan2`. */
+  bool is_first = true;
+
+  /* Initialize to zero because the first pass uses the first index to set the bounds. */
+  blender::Bounds<int> bounds_index[2] = {{0, 0}, {0, 0}};
 
   for (int i = 0, i_prev = points_hull_num - 1; i < points_hull_num; i_prev = i++) {
-    const float *ev_a = points_hull[i];
-    const float *ev_b = points_hull[i_prev];
-    float dvec[2]; /* 2d rotation matrix */
-
-    sub_v2_v2v2(dvec, ev_a, ev_b);
-    if (UNLIKELY(normalize_v2(dvec) == 0.0f)) {
+    /* 2D rotation matrix. */
+    float dvec_length = 0.0f;
+    const float2 dvec = math::normalize_and_get_length(
+        float2(points_hull[i]) - float2(points_hull[i_prev]), dvec_length);
+    if (UNLIKELY(dvec_length == 0.0f)) {
       continue;
     }
-    /* rotation matrix */
-    float min[2] = {FLT_MAX, FLT_MAX}, max[2] = {-FLT_MAX, -FLT_MAX};
-    float area_test;
 
-    for (int j = 0; j < points_hull_num; j++) {
-      float tvec[2];
-      mul_v2_v2_cw(tvec, dvec, points_hull[j]);
+    if (UNLIKELY(is_first)) {
+      is_first = false;
 
-      min[0] = min_ff(min[0], tvec[0]);
-      min[1] = min_ff(min[1], tvec[1]);
+      blender::Bounds<float> bounds[2];
 
-      max[0] = max_ff(max[0], tvec[0]);
-      max[1] = max_ff(max[1], tvec[1]);
+      bounds[0].min = bounds[0].max = mul_v2_v2_cw_x(dvec, points_hull[0]);
+      bounds[1].min = bounds[1].max = mul_v2_v2_cw_y(dvec, points_hull[0]);
 
-      area_test = (max[0] - min[0]) * (max[1] - min[1]);
-      if (area_test > area_best) {
-        break;
+      bounds_index[0].min = bounds_index[0].max = 0;
+      bounds_index[1].min = bounds_index[1].max = 0;
+
+      for (int j = 1; j < points_hull_num; j++) {
+        const float2 tvec = {
+            mul_v2_v2_cw_x(dvec, points_hull[j]),
+            mul_v2_v2_cw_y(dvec, points_hull[j]),
+        };
+        for (int axis = 0; axis < 2; axis++) {
+          if (tvec[axis] < bounds[axis].min) {
+            bounds[axis].min = tvec[axis];
+            bounds_index[axis].min = j;
+          }
+          if (tvec[axis] > bounds[axis].max) {
+            bounds[axis].max = tvec[axis];
+            bounds_index[axis].max = j;
+          }
+        }
       }
+
+      area_best = (bounds[0].max - bounds[0].min) * (bounds[1].max - bounds[1].min);
+      dvec_best = dvec;
+      continue;
     }
 
+    /* Step the calipers to the new rotation `dvec`, returning the bounds at the same time. */
+    blender::Bounds<float> bounds_test[2] = {
+        {convexhull_2d_compute_extent_on_axis<0, -1>(
+             points_hull, points_hull_num, dvec, &bounds_index[0].min),
+         convexhull_2d_compute_extent_on_axis<0, 1>(
+             points_hull, points_hull_num, dvec, &bounds_index[0].max)},
+        {convexhull_2d_compute_extent_on_axis<1, -1>(
+             points_hull, points_hull_num, dvec, &bounds_index[1].min),
+         convexhull_2d_compute_extent_on_axis<1, 1>(
+             points_hull, points_hull_num, dvec, &bounds_index[1].max)},
+
+    };
+
+    const float area_test = (bounds_test[0].max - bounds_test[0].min) *
+                            (bounds_test[1].max - bounds_test[1].min);
+
     if (area_test < area_best) {
       area_best = area_test;
-      copy_v2_v2(dvec_best, dvec);
+      dvec_best = dvec;
     }
   }