source/blender/editors/uvedit/uvedit_unwrap_ops.cc

@@ -1390,6 +1390,9 @@ static void uvedit_pack_islands_multi(const Scene *scene,
     selection_center[1] = (selection_min_co[1] + selection_max_co[1]) / 2.0f;
   }
 
+  MemArena *arena = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, __func__);
+  Heap *heap = BLI_heap_new();
+
   float scale[2] = {1.0f, 1.0f};
   blender::Vector<blender::geometry::PackIsland *> pack_island_vector;
   for (int i = 0; i < island_vector.size(); i++) {
@@ -1398,7 +1401,29 @@ static void uvedit_pack_islands_multi(const Scene *scene,
     pack_island->bounds_rect = face_island->bounds_rect;
     pack_island->caller_index = i;
     pack_island_vector.append(pack_island);
+
+    for (int i = 0; i < face_island->faces_len; i++) {
+      BMFace *f = face_island->faces[i];
+
+      /* Storage. */
+      blender::Array<blender::float2> uvs(f->len);
+
+      /* Obtain UVs of polygon. */
+      BMLoop *l;
+      BMIter iter;
+      int j;
+      BM_ITER_ELEM_INDEX (l, &iter, f, BM_LOOPS_OF_FACE, j) {
+        copy_v2_v2(uvs[j], BM_ELEM_CD_GET_FLOAT_P(l, face_island->offsets.uv));
+      }
+
+      pack_island->add_polygon(uvs, arena, heap);
+
+      BLI_memarena_clear(arena);
+    }
+    pack_island->finalize_geometry(*params, arena, heap);
   }
+  BLI_heap_free(heap, nullptr);
+  BLI_memarena_free(arena);
   pack_islands(pack_island_vector, *params, scale);
 
   float base_offset[2] = {0.0f, 0.0f};
@@ -1523,6 +1548,8 @@ static int pack_islands_exec(bContext *C, wmOperator *op)
   pack_island_params.margin_method = eUVPackIsland_MarginMethod(
       RNA_enum_get(op->ptr, "margin_method"));
   pack_island_params.margin = RNA_float_get(op->ptr, "margin");
+  pack_island_params.shape_method = eUVPackIsland_ShapeMethod(
+      RNA_enum_get(op->ptr, "shape_method"));
 
   UVMapUDIM_Params closest_udim_buf;
   UVMapUDIM_Params *closest_udim = nullptr;
@@ -1559,6 +1586,14 @@ static const EnumPropertyItem pack_margin_method_items[] = {
     {0, nullptr, 0, nullptr, nullptr},
 };
 
+static const EnumPropertyItem pack_shape_method_items[] = {
+    {ED_UVPACK_SHAPE_CONCAVE, "CONCAVE", 0, "Exact shape (Concave)", "Uses exact geometry"},
+    {ED_UVPACK_SHAPE_CONVEX, "CONVEX", 0, "Boundary shape (Convex)", "Uses convex hull"},
+    RNA_ENUM_ITEM_SEPR,
+    {ED_UVPACK_SHAPE_AABB, "AABB", 0, "Bounding box", "Uses bounding boxes"},
+    {0, nullptr, 0, nullptr, nullptr},
+};
+
 void UV_OT_pack_islands(wmOperatorType *ot)
 {
   static const EnumPropertyItem pack_target[] = {
@@ -1593,6 +1628,12 @@ void UV_OT_pack_islands(wmOperatorType *ot)
                "");
   RNA_def_float_factor(
       ot->srna, "margin", 0.001f, 0.0f, 1.0f, "Margin", "Space between islands", 0.0f, 1.0f);
+  RNA_def_enum(ot->srna,
+               "shape_method",
+               pack_shape_method_items,
+               ED_UVPACK_SHAPE_CONCAVE,
+               "Shape Method",
+               "");
 }
 
 /** \} */

source/blender/geometry/GEO_uv_pack.hh

@@ -1,7 +1,10 @@
 /* SPDX-License-Identifier: GPL-2.0-or-later */
 
+#include "BLI_heap.h"
 #include "BLI_math_matrix.hh"
+#include "BLI_memarena.h"
 #include "BLI_span.hh"
+#include "BLI_vector.hh"
 
 #include "DNA_space_types.h"
 #include "DNA_vec_types.h"
@@ -20,6 +23,12 @@ enum eUVPackIsland_MarginMethod {
   ED_UVPACK_MARGIN_FRACTION,   /* Specify a precise fraction of final UV output. */
 };
 
+enum eUVPackIsland_ShapeMethod {
+  ED_UVPACK_SHAPE_AABB = 0, /* Use Axis-Aligned Bounding-Boxes. */
+  ED_UVPACK_SHAPE_CONVEX,   /* Use convex hull. */
+  ED_UVPACK_SHAPE_CONCAVE,  /* Use concave hull. */
+};
+
 namespace blender::geometry {
 
 /** See also #UnwrapOptions. */
@@ -51,6 +60,8 @@ class UVPackIsland_Params {
   eUVPackIsland_MarginMethod margin_method;
   /** Additional translation for bottom left corner. */
   float udim_base_offset[2];
+  /** Which shape to use when packing. */
+  eUVPackIsland_ShapeMethod shape_method;
 };
 
 class PackIsland {
@@ -58,6 +69,14 @@ class PackIsland {
   rctf bounds_rect;
   float2 pre_translate; /* Output. */
   int caller_index;     /* Unchanged by #pack_islands, used by caller. */
+
+  void add_triangle(const float2 uv0, const float2 uv1, const float2 uv2);
+  void add_polygon(const blender::Span<float2> uvs, MemArena *arena, Heap *heap);
+  void finalize_geometry(const UVPackIsland_Params &params, MemArena *arena, Heap *heap);
+
+ private:
+  blender::Vector<float2> triangleVertices;
+  friend class Occupancy;
 };
 
 void pack_islands(const Span<PackIsland *> &islands,

source/blender/geometry/intern/uv_pack.cc

@@ -11,6 +11,8 @@
 #include "BLI_convexhull_2d.h"
 #include "BLI_listbase.h"
 #include "BLI_math.h"
+#include "BLI_polyfill_2d.h"
+#include "BLI_polyfill_2d_beautify.h"
 #include "BLI_rect.h"
 #include "BLI_vector.hh"
 
@@ -22,6 +24,119 @@
 
 namespace blender::geometry {
 
+/* Compute signed distance squared to a line passing through `uva` and `uvb`.
+ */
+static float dist_signed_squared_to_edge(float2 probe, float2 uva, float2 uvb)
+{
+  const float2 edge = uvb - uva;
+  const float2 side = probe - uva;
+
+  const float edge_length_squared = blender::math::length_squared(edge);
+  /* Tolerance here is to avoid division by zero later. */
+  if (edge_length_squared < 1e-40f) {
+    return blender::math::length_squared(side);
+  }
+
+  const float numerator = edge.x * side.y - edge.y * side.x; /* c.f. cross product. */
+  const float numerator_ssq = numerator >= 0.0f ? numerator * numerator : -numerator * numerator;
+
+  return numerator_ssq / edge_length_squared;
+}
+
+void PackIsland::add_triangle(const float2 uv0, const float2 uv1, const float2 uv2)
+{
+  /* Be careful with winding. */
+  if (dist_signed_squared_to_edge(uv0, uv1, uv2) < 0.0f) {
+    triangleVertices.append(uv0);
+    triangleVertices.append(uv1);
+    triangleVertices.append(uv2);
+  }
+  else {
+    triangleVertices.append(uv0);
+    triangleVertices.append(uv2);
+    triangleVertices.append(uv1);
+  }
+}
+
+void PackIsland::add_polygon(const blender::Span<float2> uvs, MemArena *arena, Heap *heap)
+{
+  int vert_count = int(uvs.size());
+  BLI_assert(vert_count >= 3);
+  int nfilltri = vert_count - 2;
+  if (nfilltri == 1) {
+    /* Trivial case, just one triangle. */
+    add_triangle(uvs[0], uvs[1], uvs[2]);
+    return;
+  }
+
+  /* Storage. */
+  uint(*tris)[3] = static_cast<uint(*)[3]>(
+      BLI_memarena_alloc(arena, sizeof(*tris) * size_t(nfilltri)));
+  float(*source)[2] = static_cast<float(*)[2]>(
+      BLI_memarena_alloc(arena, sizeof(*source) * size_t(vert_count)));
+
+  /* Copy input. */
+  for (int i = 0; i < vert_count; i++) {
+    copy_v2_v2(source[i], uvs[i]);
+  }
+
+  /* Triangulate. */
+  BLI_polyfill_calc_arena(source, vert_count, 0, tris, arena);
+
+  /* Beautify improves performance of packer. (Optional)
+   * Long thin triangles, especially at 45 degree angles,
+   * can trigger worst-case performance in #trace_triangle.
+   * Using `Beautify` brings more inputs into average-case.
+   */
+  BLI_polyfill_beautify(source, vert_count, tris, arena, heap);
+
+  /* Add as triangles. */
+  for (int j = 0; j < nfilltri; j++) {
+    uint *tri = tris[j];
+    add_triangle(source[tri[0]], source[tri[1]], source[tri[2]]);
+  }
+
+  BLI_heap_clear(heap, nullptr);
+}
+
+void PackIsland::finalize_geometry(const UVPackIsland_Params &params, MemArena *arena, Heap *heap)
+{
+  BLI_assert(triangleVertices.size() >= 3);
+  const eUVPackIsland_ShapeMethod shape_method = params.shape_method;
+  if (shape_method == ED_UVPACK_SHAPE_CONVEX) {
+    /* Compute convex hull of existing triangles. */
+    if (triangleVertices.size() <= 3) {
+      return; /* Trivial case, nothing to do. */
+    }
+
+    int vert_count = int(triangleVertices.size());
+
+    /* Allocate storage. */
+    int *index_map = static_cast<int *>(
+        BLI_memarena_alloc(arena, sizeof(*index_map) * vert_count));
+    float(*source)[2] = static_cast<float(*)[2]>(
+        BLI_memarena_alloc(arena, sizeof(*source) * size_t(vert_count)));
+
+    /* Prepare input for convex hull. */
+    for (int i = 0; i < vert_count; i++) {
+      copy_v2_v2(source[i], triangleVertices[i]);
+    }
+
+    /* Compute convex hull. */
+    int convex_len = BLI_convexhull_2d(source, vert_count, index_map);
+
+    /* Write back. */
+    triangleVertices.clear();
+    blender::Array<float2> convexVertices(convex_len);
+    for (int i = 0; i < convex_len; i++) {
+      convexVertices[i] = source[index_map[i]];
+    }
+    add_polygon(convexVertices, arena, heap);
+
+    BLI_heap_clear(heap, nullptr);
+  }
+}
+
 UVPackIsland_Params::UVPackIsland_Params()
 {
   /* TEMPORARY, set every thing to "zero" for backwards compatibility. */
@@ -36,6 +151,7 @@ UVPackIsland_Params::UVPackIsland_Params()
   margin_method = ED_UVPACK_MARGIN_SCALED;
   udim_base_offset[0] = 0.0f;
   udim_base_offset[1] = 0.0f;
+  shape_method = ED_UVPACK_SHAPE_AABB;
 }
 
 /* Compact representation for AABB packers. */
@@ -110,10 +226,317 @@ static void pack_islands_alpaca_turbo(const Span<UVAABBIsland *> islands,
   *r_max_v = next_v1;
 }
 
+/**
+ * Helper class for the `xatlas` strategy.
+ * Accelerates geometry queries by approximating exact queries with a bitmap.
+ * Includes some book keeping variables to simplify the algorithm.
+ */
+class Occupancy {
+ public:
+  Occupancy(const float initial_scale);
+
+  void increase_scale(); /* Resize the scale of the bitmap and clear it. */
+
+  /* Write or Query a triangle on the bitmap. */
+  float trace_triangle(const float2 &uv0,
+                       const float2 &uv1,
+                       const float2 &uv2,
+                       const float margin,
+                       const bool write) const;
+
+  /* Write or Query an island on the bitmap. */
+  float trace_island(PackIsland *island,
+                     const float scale,
+                     const float margin,
+                     const float2 &uv,
+                     const bool write) const;
+
+  int bitmap_radix;              /* Width and Height of `bitmap`. */
+  float bitmap_scale_reciprocal; /* == 1.0f / `bitmap_scale`. */
+ private:
+  mutable blender::Array<float> bitmap;
+
+  mutable float2 witness;         /* Witness to a previously known occupied pixel. */
+  mutable float witness_distance; /* Signed distance to nearest placed island. */
+  mutable uint triangle_hint;     /* Hint to a previously suspected overlapping triangle. */
+
+  const float terminal = 1048576.0f; /* A "very" large number, much bigger than 4 * bitmap_radix */
+};
+
+Occupancy::Occupancy(const float initial_scale)
+    : bitmap_radix(800), bitmap(bitmap_radix * bitmap_radix, false)
+{
+  increase_scale();
+  bitmap_scale_reciprocal = bitmap_radix / initial_scale;
+}
+
+void Occupancy::increase_scale()
+{
+  bitmap_scale_reciprocal *= 0.5f;
+  for (int i = 0; i < bitmap_radix * bitmap_radix; i++) {
+    bitmap[i] = terminal;
+  }
+  witness.x = -1;
+  witness.y = -1;
+  witness_distance = 0.0f;
+  triangle_hint = 0;
+}
+
+static float signed_distance_fat_triangle(const float2 probe,
+                                          const float2 uv0,
+                                          const float2 uv1,
+                                          const float2 uv2)
+{
+  /* Be careful with ordering, uv0 <- uv1 <- uv2 <- uv0 <- uv1 etc. */
+  const float dist01_ssq = dist_signed_squared_to_edge(probe, uv0, uv1);
+  const float dist12_ssq = dist_signed_squared_to_edge(probe, uv1, uv2);
+  const float dist20_ssq = dist_signed_squared_to_edge(probe, uv2, uv0);
+  float result_ssq = max_fff(dist01_ssq, dist12_ssq, dist20_ssq);
+  if (result_ssq < 0.0f) {
+    return -sqrtf(-result_ssq);
+  }
+  BLI_assert(result_ssq >= 0.0f);
+  result_ssq = std::min(result_ssq, blender::math::length_squared(probe - uv0));
+  result_ssq = std::min(result_ssq, blender::math::length_squared(probe - uv1));
+  result_ssq = std::min(result_ssq, blender::math::length_squared(probe - uv2));
+  BLI_assert(result_ssq >= 0.0f);
+  return sqrtf(result_ssq);
+}
+
+float Occupancy::trace_triangle(const float2 &uv0,
+                                const float2 &uv1,
+                                const float2 &uv2,
+                                const float margin,
+                                const bool write) const
+{
+  const float x0 = min_fff(uv0.x, uv1.x, uv2.x);
+  const float y0 = min_fff(uv0.y, uv1.y, uv2.y);
+  const float x1 = max_fff(uv0.x, uv1.x, uv2.x);
+  const float y1 = max_fff(uv0.y, uv1.y, uv2.y);
+  float spread = write ? margin * 2 : 0.0f;
+  int ix0 = std::max(int(floorf((x0 - spread) * bitmap_scale_reciprocal)), 0);
+  int iy0 = std::max(int(floorf((y0 - spread) * bitmap_scale_reciprocal)), 0);
+  int ix1 = std::min(int(floorf((x1 + spread) * bitmap_scale_reciprocal + 2)), bitmap_radix);
+  int iy1 = std::min(int(floorf((y1 + spread) * bitmap_scale_reciprocal + 2)), bitmap_radix);
+
+  const float2 uv0s = uv0 * bitmap_scale_reciprocal;
+  const float2 uv1s = uv1 * bitmap_scale_reciprocal;
+  const float2 uv2s = uv2 * bitmap_scale_reciprocal;
+
+  /* TODO: Better epsilon handling here could reduce search size. */
+  float epsilon = 0.7071f; /* == sqrt(0.5f), rounded up by 0.00002f. */
+  epsilon = std::max(epsilon, 2 * margin * bitmap_scale_reciprocal);
+
+  if (!write) {
+    if (ix0 <= witness.x && witness.x < ix1) {
+      if (iy0 <= witness.y && witness.y < iy1) {
+        const float distance = signed_distance_fat_triangle(witness, uv0s, uv1s, uv2s);
+        const float extent = epsilon - distance - witness_distance;
+        if (extent > 0.0f) {
+          return extent; /* Witness observes occupied. */
+        }
+      }
+    }
+  }
+
+  /* Iterate in opposite direction to outer search to improve witness effectiveness. */
+  for (int y = iy1 - 1; y >= iy0; y--) {
+    for (int x = ix1 - 1; x >= ix0; x--) {
+      float *hotspot = &bitmap[y * bitmap_radix + x];
+      if (!write && *hotspot > epsilon) {
+        continue;
+      }
+      const float2 probe(x, y);
+      const float distance = signed_distance_fat_triangle(probe, uv0s, uv1s, uv2s);
+      if (write) {
+        *hotspot = min_ff(distance, *hotspot);
+        continue;
+      }
+      const float extent = epsilon - distance - *hotspot;
+      if (extent > 0.0f) {
+        witness = probe;
+        witness_distance = *hotspot;
+        return extent; /* Occupied. */
+      }
+    }
+  }
+  return -1.0f; /* Available. */
+}
+
+float Occupancy::trace_island(PackIsland *island,
+                              const float scale,
+                              const float margin,
+                              const float2 &uv,
+                              const bool write) const
+{
+
+  if (!write) {
+    if (uv.x <= 0.0f || uv.y <= 0.0f) {
+      return std::max(-uv.x, -uv.y); /* Occupied. */
+    }
+  }
+  const float2 origin(island->bounds_rect.xmin, island->bounds_rect.ymin);
+  const float2 delta = uv - origin * scale;
+  uint vert_count = uint(island->triangleVertices.size());
+  for (uint i = 0; i < vert_count; i += 3) {
+    uint j = (i + triangle_hint) % vert_count;
+    float extent = trace_triangle(delta + island->triangleVertices[j] * scale,
+                                  delta + island->triangleVertices[j + 1] * scale,
+                                  delta + island->triangleVertices[j + 2] * scale,
+                                  margin,
+                                  write);
+
+    if (!write && extent >= 0.0f) {
+      triangle_hint = j;
+      return extent; /* Occupied. */
+    }
+  }
+  return -1.0f; /* Available. */
+}
+
+static float2 find_best_fit_for_island(
+    PackIsland *island, int scan_line, Occupancy &occupancy, const float scale, const float margin)
+{
+  const float scan_line_bscaled = scan_line / occupancy.bitmap_scale_reciprocal;
+  const float size_x_scaled = BLI_rctf_size_x(&island->bounds_rect) * scale;
+  const float size_y_scaled = BLI_rctf_size_y(&island->bounds_rect) * scale;
+
+  int t = 0;
+  while (t <= scan_line) {
+    const float t_bscaled = t / occupancy.bitmap_scale_reciprocal;
+
+    /* Scan using an "Alpaca"-style search, both horizontally and vertically at the same time. */
+
+    const float2 horiz(scan_line_bscaled - size_x_scaled, t_bscaled - size_y_scaled);
+    const float extentH = occupancy.trace_island(island, scale, margin, horiz, false);
+    if (extentH < 0.0f) {
+      return horiz;
+    }
+
+    const float2 vert(t_bscaled - size_x_scaled, scan_line_bscaled - size_y_scaled);
+    const float extentV = occupancy.trace_island(island, scale, margin, vert, false);
+    if (extentV < 0.0f) {
+      return vert;
+    }
+
+    const float min_extent = std::min(extentH, extentV);
+
+    t = t + std::max(1, int(min_extent));
+  }
+  return float2(-1, -1);
+}
+
+static float guess_initial_scale(const Span<PackIsland *> islands,
+                                 const float scale,
+                                 const float margin)
+{
+  float sum = 1e-40f;
+  for (int64_t i : islands.index_range()) {
+    PackIsland *island = islands[i];
+    sum += BLI_rctf_size_x(&island->bounds_rect) * scale + 2 * margin;
+    sum += BLI_rctf_size_y(&island->bounds_rect) * scale + 2 * margin;
+  }
+  return sqrtf(sum) / 3.0f;
+}
+
+/**
+ * Pack irregular islands using the "xatlas" strategy, with no rotation.
+ *
+ * Loosely based on the 'xatlas' code by Jonathan Young
+ * from https://github.com/jpcy/xatlas
+ *
+ * A brute force packer (BFPacker) with accelerators:
+ * - Uses a Bitmap Occupancy class.
+ * - Uses a "Witness Pixel" and a "Triangle Hint".
+ * - Write with `margin * 2`, read with `margin == 0`.
+ * - Lazy resetting of BF search.
+ *
+ * Performance would normally be `O(n^4)`, however the occupancy
+ * bitmap_radix is fixed, which gives a reduced time complexity of `O(n^3)`.
+ */
+static void pack_island_xatlas(const Span<UVAABBIsland *> island_indices,
+                               const Span<PackIsland *> islands,
+                               BoxPack *box_array,
+                               const float scale,
+                               const float margin,
+                               float *r_max_u,
+                               float *r_max_v)
+{
+  Occupancy occupancy(guess_initial_scale(islands, scale, margin));
+  float max_u = 0.0f;
+  float max_v = 0.0f;
+
+  int scan_line = 0;
+  int i = 0;
+
+  /* The following `while` loop is setting up a three-way race:
+   * for (scan_line=0; scan_line<bitmap_radix; scan_line++)
+   * for (i : island_indices.index_range())
+   * while (bitmap_scale_reciprocal > 0) { bitmap_scale_reciprocal *= 0.5f; }
+   */
+
+  while (i < island_indices.size()) {
+    PackIsland *island = islands[island_indices[i]->index];
+    const float2 best = find_best_fit_for_island(island, scan_line, occupancy, scale, margin);
+
+    if (best.x <= -1.0f) {
+      /* Unable to find a fit on this scan_line. */
+
+      if (i < 10) {
+        scan_line++;
+      }
+      else {
+        /* Increasing by 2 here has the effect of changing the sampling pattern.
+         * The parameter '2' is not "free" in the sense that changing it requires
+         * a change to `bitmap_radix` and then retuning `alpaca_cutoff`.
+         * Possible values here *could* be 1, 2 or 3, however the only *reasonable*
+         * choice is 2.
+         */
+        scan_line += 2;
+      }
+      if (scan_line < occupancy.bitmap_radix) {
+        continue; /* Try again on next scan_line. */
+      }
+
+      /* Enlarge search parameters. */
+      scan_line = 0;
+      occupancy.increase_scale();
+
+      /* Redraw already placed islands. (Greedy.) */
+      for (int j = 0; j < i; j++) {
+        BoxPack *box = box_array + j;
+        occupancy.trace_island(
+            islands[island_indices[j]->index], scale, margin, float2(box->x, box->y), true);
+      }
+      continue;
+    }
+
+    /* Place island. */
+    BoxPack *box = box_array + i;
+    box->x = best.x;
+    box->y = best.y;
+    max_u = std::max(box->x + BLI_rctf_size_x(&island->bounds_rect) * scale + 2 * margin, max_u);
+    max_v = std::max(box->y + BLI_rctf_size_y(&island->bounds_rect) * scale + 2 * margin, max_v);
+    occupancy.trace_island(island, scale, margin, float2(box->x, box->y), true);
+    i++; /* Next island. */
+
+    if (i < 128 || (i & 31) == 16) {
+      scan_line = 0; /* Restart completely. */
+    }
+    else {
+      scan_line = std::max(0, scan_line - 25); /* `-25` must by odd. */
+    }
+  }
+
+  *r_max_u = max_u;
+  *r_max_v = max_v;
+}
+
 static float pack_islands_scale_margin(const Span<PackIsland *> islands,
                                        BoxPack *box_array,
                                        const float scale,
-                                       const float margin)
+                                       const float margin,
+                                       const UVPackIsland_Params &params)
 {
   /* #BLI_box_pack_2d produces layouts with high packing efficiency, but has `O(n^3)`
    * time complexity, causing poor performance if there are lots of islands. See: #102843.
@@ -148,10 +571,18 @@ static float pack_islands_scale_margin(const Span<PackIsland *> islands,
     return b->uv_diagonal.x * b->uv_diagonal.y < a->uv_diagonal.x * a->uv_diagonal.y;
   });
 
-  /* Partition island_vector, largest will go to box_pack, the rest alpaca_turbo.
+  /* Partition `islands`, largest will go to a slow packer, the rest alpaca_turbo.
    * See discussion above for details. */
-  const int64_t alpaca_cutoff = int64_t(1024); /* TODO: Tune constant. */
-  int64_t max_box_pack = std::min(alpaca_cutoff, islands.size());
+  int64_t alpaca_cutoff = int64_t(
+      1024); /* Regular situation, pack 1024 islands with slow packer. */
+  int64_t alpaca_cutoff_fast = int64_t(
+      80); /* Reduced problem size, only 80 islands with slow packer. */
+  if (params.margin_method == ED_UVPACK_MARGIN_FRACTION) {
+    if (margin > 0.0f) {
+      alpaca_cutoff = alpaca_cutoff_fast;
+    }
+  }
+  const int64_t max_box_pack = std::min(alpaca_cutoff, islands.size());
 
   /* Prepare for box_pack_2d. */
   for (const int64_t i : islands.index_range()) {
@@ -165,7 +596,21 @@ static float pack_islands_scale_margin(const Span<PackIsland *> islands,
   /* Call box_pack_2d (slow for large N.) */
   float max_u = 0.0f;
   float max_v = 0.0f;
-  BLI_box_pack_2d(box_array, int(max_box_pack), &max_u, &max_v);
+  switch (params.shape_method) {
+    case ED_UVPACK_SHAPE_CONVEX:
+    case ED_UVPACK_SHAPE_CONCAVE:
+      pack_island_xatlas(aabbs.as_span().take_front(max_box_pack),
+                         islands,
+                         box_array,
+                         scale,
+                         margin,
+                         &max_u,
+                         &max_v);
+      break;
+    default:
+      BLI_box_pack_2d(box_array, int(max_box_pack), &max_u, &max_v);
+      break;
+  }
 
   /* At this stage, `max_u` and `max_v` contain the box_pack UVs. */
 
@@ -192,7 +637,8 @@ static float pack_islands_scale_margin(const Span<PackIsland *> islands,
 
 static float pack_islands_margin_fraction(const Span<PackIsland *> &island_vector,
                                           BoxPack *box_array,
-                                          const float margin_fraction)
+                                          const float margin_fraction,
+                                          const UVPackIsland_Params &params)
 {
   /*
    * Root finding using a combined search / modified-secant method.
@@ -260,7 +706,7 @@ static float pack_islands_margin_fraction(const Span<PackIsland *> &island_vecto
     /* Evaluate our `f`. */
     scale_last = scale;
     float max_uv = pack_islands_scale_margin(
-        island_vector, box_array, scale_last, margin_fraction);
+        island_vector, box_array, scale_last, margin_fraction, params);
     float value = sqrtf(max_uv) - 1.0f;
 
     if (value <= 0.0f) {
@@ -284,7 +730,7 @@ static float pack_islands_margin_fraction(const Span<PackIsland *> &island_vecto
     if (scale_last != scale_low) {
       scale_last = scale_low;
       float max_uv = pack_islands_scale_margin(
-          island_vector, box_array, scale_last, margin_fraction);
+          island_vector, box_array, scale_last, margin_fraction, params);
       UNUSED_VARS(max_uv);
       /* TODO (?): `if (max_uv < 1.0f) { scale_last /= max_uv; }` */
     }
@@ -326,7 +772,7 @@ static BoxPack *pack_islands_box_array(const Span<PackIsland *> &island_vector,
 
   if (params.margin == 0.0f) {
     /* Special case for zero margin. Margin_method is ignored as all formulas give same result. */
-    const float max_uv = pack_islands_scale_margin(island_vector, box_array, 1.0f, 0.0f);
+    const float max_uv = pack_islands_scale_margin(island_vector, box_array, 1.0f, 0.0f, params);
     r_scale[0] = 1.0f / max_uv;
     r_scale[1] = r_scale[0];
     return box_array;
@@ -334,7 +780,8 @@ static BoxPack *pack_islands_box_array(const Span<PackIsland *> &island_vector,
 
   if (params.margin_method == ED_UVPACK_MARGIN_FRACTION) {
     /* Uses a line search on scale. ~10x slower than other method. */
-    const float scale = pack_islands_margin_fraction(island_vector, box_array, params.margin);
+    const float scale = pack_islands_margin_fraction(
+        island_vector, box_array, params.margin, params);
     r_scale[0] = scale;
     r_scale[1] = scale;
     /* pack_islands_margin_fraction will pad FaceIslands, return early. */
@@ -355,7 +802,7 @@ static BoxPack *pack_islands_box_array(const Span<PackIsland *> &island_vector,
       BLI_assert_unreachable();
   }
 
-  const float max_uv = pack_islands_scale_margin(island_vector, box_array, 1.0f, margin);
+  const float max_uv = pack_islands_scale_margin(island_vector, box_array, 1.0f, margin, params);
   r_scale[0] = 1.0f / max_uv;
   r_scale[1] = r_scale[0];