Cleanup: Move UV edit parameterize code to geometry module

This will allow reusing it elsewhere, such as in a geometry node.

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

5 files changed, 57 insertions, 5133 deletions

diff --git a/source/blender/editors/uvedit/CMakeLists.txt b/source/blender/editors/uvedit/CMakeLists.txt
index f8a192e3254..761e7cd091e 100644
--- a/source/blender/editors/uvedit/CMakeLists.txt
+++ b/source/blender/editors/uvedit/CMakeLists.txt
@@ -7,6 +7,7 @@ set(INC
   ../../blentranslation
   ../../bmesh
   ../../depsgraph
+  ../../geometry
   ../../gpu
   ../../makesdna
   ../../makesrna
@@ -24,7 +25,6 @@ set(SRC
   uvedit_draw.c
   uvedit_islands.c
   uvedit_ops.c
-  uvedit_parametrizer.c
   uvedit_path.c
   uvedit_rip.c
   uvedit_select.c
@@ -32,7 +32,6 @@ set(SRC
   uvedit_unwrap_ops.c
 
   uvedit_intern.h
-  uvedit_parametrizer.h
 )
 
 set(LIB
diff --git a/source/blender/editors/uvedit/uvedit_islands.c b/source/blender/editors/uvedit/uvedit_islands.c
index 59992d23e2e..e1752ae5a29 100644
--- a/source/blender/editors/uvedit/uvedit_islands.c
+++ b/source/blender/editors/uvedit/uvedit_islands.c
@@ -5,7 +5,7 @@
  *
  * Utilities for manipulating UV islands.
  *
- * \note This is similar to `uvedit_parametrizer.c`,
+ * \note This is similar to `GEO_uv_parametrizer.h`,
  * however the data structures there don't support arbitrary topology
  * such as an edge with 3 or more faces using it.
  * This API uses #BMesh data structures and doesn't have limitations for manifold meshes.
diff --git a/source/blender/editors/uvedit/uvedit_parametrizer.c b/source/blender/editors/uvedit/uvedit_parametrizer.c
deleted file mode 100644
index 5ad326c19e5..00000000000
--- a/source/blender/editors/uvedit/uvedit_parametrizer.c
+++ /dev/null
@@ -1,4978 +0,0 @@
-/* SPDX-License-Identifier: GPL-2.0-or-later */
-
-/** \file
- * \ingroup eduv
- */
-
-#include "MEM_guardedalloc.h"
-
-#include "BLI_boxpack_2d.h"
-#include "BLI_convexhull_2d.h"
-#include "BLI_heap.h"
-#include "BLI_math.h"
-#include "BLI_memarena.h"
-#include "BLI_polyfill_2d.h"
-#include "BLI_polyfill_2d_beautify.h"
-#include "BLI_rand.h"
-#include "BLI_utildefines.h"
-
-#include "uvedit_parametrizer.h"
-
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#include "BLI_sys_types.h" /* for intptr_t support */
-
-#include "eigen_capi.h"
-
-/* Utils */
-
-#define param_assert(condition) \
-  if (!(condition)) { /*printf("Assertion %s:%d\n", __FILE__, __LINE__); abort();*/ \
-  } \
-  (void)0
-#define param_warning(message) \
-  {/*printf("Warning %s:%d: %s\n", __FILE__, __LINE__, message);*/}(void)0
-
-typedef enum PBool {
-  P_TRUE = 1,
-  P_FALSE = 0,
-} PBool;
-
-/* Special Purpose Hash */
-
-typedef intptr_t PHashKey;
-
-typedef struct PHashLink {
-  struct PHashLink *next;
-  PHashKey key;
-} PHashLink;
-
-typedef struct PHash {
-  PHashLink **list;
-  PHashLink **buckets;
-  int size, cursize, cursize_id;
-} PHash;
-
-struct PChart;
-struct PEdge;
-struct PFace;
-struct PHandle;
-struct PVert;
-
-/* Simplices */
-
-typedef struct PVert {
-  struct PVert *nextlink;
-
-  union PVertUnion {
-    PHashKey key;       /* Construct. */
-    int id;             /* ABF/LSCM matrix index. */
-    float distortion;   /* Area smoothing. */
-    HeapNode *heaplink; /* Edge collapsing. */
-  } u;
-
-  struct PEdge *edge;
-  float co[3];
-  float uv[2];
-  uchar flag;
-
-} PVert;
-
-typedef struct PEdge {
-  struct PEdge *nextlink;
-
-  union PEdgeUnion {
-    PHashKey key;               /* Construct. */
-    int id;                     /* ABF matrix index. */
-    HeapNode *heaplink;         /* Fill holes. */
-    struct PEdge *nextcollapse; /* Simplification. */
-  } u;
-
-  struct PVert *vert;
-  struct PEdge *pair;
-  struct PEdge *next;
-  struct PFace *face;
-  float *orig_uv, old_uv[2];
-  ushort flag;
-
-} PEdge;
-
-typedef struct PFace {
-  struct PFace *nextlink;
-
-  union PFaceUnion {
-    PHashKey key; /* Construct. */
-    int chart;    /* Construct splitting. */
-    float area3d; /* Stretch. */
-    int id;       /* ABF matrix index. */
-  } u;
-
-  struct PEdge *edge;
-  uchar flag;
-} PFace;
-
-enum PVertFlag {
-  PVERT_PIN = 1,
-  PVERT_SELECT = 2,
-  PVERT_INTERIOR = 4,
-  PVERT_COLLAPSE = 8,
-  PVERT_SPLIT = 16,
-};
-
-enum PEdgeFlag {
-  PEDGE_SEAM = 1,
-  PEDGE_VERTEX_SPLIT = 2,
-  PEDGE_PIN = 4,
-  PEDGE_SELECT = 8,
-  PEDGE_DONE = 16,
-  PEDGE_FILLED = 32,
-  PEDGE_COLLAPSE = 64,
-  PEDGE_COLLAPSE_EDGE = 128,
-  PEDGE_COLLAPSE_PAIR = 256,
-};
-
-/* for flipping faces */
-#define PEDGE_VERTEX_FLAGS (PEDGE_PIN)
-
-enum PFaceFlag {
-  PFACE_CONNECTED = 1,
-  PFACE_FILLED = 2,
-  PFACE_COLLAPSE = 4,
-};
-
-/* Chart */
-
-typedef struct PChart {
-  PVert *verts;
-  PEdge *edges;
-  PFace *faces;
-  int nverts, nedges, nfaces;
-
-  PVert *collapsed_verts;
-  PEdge *collapsed_edges;
-  PFace *collapsed_faces;
-
-  union PChartUnion {
-    struct PChartLscm {
-      LinearSolver *context;
-      float *abf_alpha;
-      PVert *pin1, *pin2;
-      PVert *single_pin;
-      float single_pin_area;
-      float single_pin_uv[2];
-    } lscm;
-    struct PChartPack {
-      float rescale, area;
-      float size[2] /* , trans[2] */;
-    } pack;
-  } u;
-
-  uchar flag;
-  struct PHandle *handle;
-} PChart;
-
-enum PChartFlag {
-  PCHART_HAS_PINS = 1,
-};
-
-enum PHandleState {
-  PHANDLE_STATE_ALLOCATED,
-  PHANDLE_STATE_CONSTRUCTED,
-  PHANDLE_STATE_LSCM,
-  PHANDLE_STATE_STRETCH,
-};
-
-typedef struct PHandle {
-  enum PHandleState state;
-  MemArena *arena;
-  MemArena *polyfill_arena;
-  Heap *polyfill_heap;
-
-  PChart *construction_chart;
-  PHash *hash_verts;
-  PHash *hash_edges;
-  PHash *hash_faces;
-
-  PChart **charts;
-  int ncharts;
-
-  float aspx, aspy;
-
-  RNG *rng;
-  float blend;
-  char do_aspect;
-} PHandle;
-
-/* PHash
- * - special purpose hash that keeps all its elements in a single linked list.
- * - after construction, this hash is thrown away, and the list remains.
- * - removing elements is not possible efficiently.
- */
-
-static int PHashSizes[] = {
-    1,       3,       5,       11,      17,       37,       67,       131,       257,       521,
-    1031,    2053,    4099,    8209,    16411,    32771,    65537,    131101,    262147,    524309,
-    1048583, 2097169, 4194319, 8388617, 16777259, 33554467, 67108879, 134217757, 268435459,
-};
-
-#define PHASH_hash(ph, item) (((uintptr_t)(item)) % ((uint)(ph)->cursize))
-#define PHASH_edge(v1, v2) (((v1) < (v2)) ? ((v1)*39) ^ ((v2)*31) : ((v1)*31) ^ ((v2)*39))
-
-static PHash *phash_new(PHashLink **list, int sizehint)
-{
-  PHash *ph = (PHash *)MEM_callocN(sizeof(PHash), "PHash");
-  ph->size = 0;
-  ph->cursize_id = 0;
-  ph->list = list;
-
-  while (PHashSizes[ph->cursize_id] < sizehint) {
-    ph->cursize_id++;
-  }
-
-  ph->cursize = PHashSizes[ph->cursize_id];
-  ph->buckets = (PHashLink **)MEM_callocN(ph->cursize * sizeof(*ph->buckets), "PHashBuckets");
-
-  return ph;
-}
-
-static void phash_delete(PHash *ph)
-{
-  MEM_freeN(ph->buckets);
-  MEM_freeN(ph);
-}
-
-static int phash_size(PHash *ph)
-{
-  return ph->size;
-}
-
-static void phash_insert(PHash *ph, PHashLink *link)
-{
-  int size = ph->cursize;
-  uintptr_t hash = PHASH_hash(ph, link->key);
-  PHashLink *lookup = ph->buckets[hash];
-
-  if (lookup == NULL) {
-    /* insert in front of the list */
-    ph->buckets[hash] = link;
-    link->next = *(ph->list);
-    *(ph->list) = link;
-  }
-  else {
-    /* insert after existing element */
-    link->next = lookup->next;
-    lookup->next = link;
-  }
-
-  ph->size++;
-
-  if (ph->size > (size * 3)) {
-    PHashLink *next = NULL, *first = *(ph->list);
-
-    ph->cursize = PHashSizes[++ph->cursize_id];
-    MEM_freeN(ph->buckets);
-    ph->buckets = (PHashLink **)MEM_callocN(ph->cursize * sizeof(*ph->buckets), "PHashBuckets");
-    ph->size = 0;
-    *(ph->list) = NULL;
-
-    for (link = first; link; link = next) {
-      next = link->next;
-      phash_insert(ph, link);
-    }
-  }
-}
-
-static PHashLink *phash_lookup(PHash *ph, PHashKey key)
-{
-  PHashLink *link;
-  uintptr_t hash = PHASH_hash(ph, key);
-
-  for (link = ph->buckets[hash]; link; link = link->next) {
-    if (link->key == key) {
-      return link;
-    }
-    if (PHASH_hash(ph, link->key) != hash) {
-      return NULL;
-    }
-  }
-
-  return link;
-}
-
-static PHashLink *phash_next(PHash *ph, PHashKey key, PHashLink *link)
-{
-  uintptr_t hash = PHASH_hash(ph, key);
-
-  for (link = link->next; link; link = link->next) {
-    if (link->key == key) {
-      return link;
-    }
-    if (PHASH_hash(ph, link->key) != hash) {
-      return NULL;
-    }
-  }
-
-  return link;
-}
-
-/* Geometry */
-
-static float p_vec_angle_cos(const float v1[3], const float v2[3], const float v3[3])
-{
-  float d1[3], d2[3];
-
-  d1[0] = v1[0] - v2[0];
-  d1[1] = v1[1] - v2[1];
-  d1[2] = v1[2] - v2[2];
-
-  d2[0] = v3[0] - v2[0];
-  d2[1] = v3[1] - v2[1];
-  d2[2] = v3[2] - v2[2];
-
-  normalize_v3(d1);
-  normalize_v3(d2);
-
-  return d1[0] * d2[0] + d1[1] * d2[1] + d1[2] * d2[2];
-}
-
-static float p_vec_angle(const float v1[3], const float v2[3], const float v3[3])
-{
-  float dot = p_vec_angle_cos(v1, v2, v3);
-
-  if (dot <= -1.0f) {
-    return (float)M_PI;
-  }
-  if (dot >= 1.0f) {
-    return 0.0f;
-  }
-  return acosf(dot);
-}
-
-static float p_vec2_angle(const float v1[2], const float v2[2], const float v3[2])
-{
-  float u1[3], u2[3], u3[3];
-
-  u1[0] = v1[0];
-  u1[1] = v1[1];
-  u1[2] = 0.0f;
-  u2[0] = v2[0];
-  u2[1] = v2[1];
-  u2[2] = 0.0f;
-  u3[0] = v3[0];
-  u3[1] = v3[1];
-  u3[2] = 0.0f;
-
-  return p_vec_angle(u1, u2, u3);
-}
-
-static void p_triangle_angles(
-    const float v1[3], const float v2[3], const float v3[3], float *r_a1, float *r_a2, float *r_a3)
-{
-  *r_a1 = p_vec_angle(v3, v1, v2);
-  *r_a2 = p_vec_angle(v1, v2, v3);
-  *r_a3 = (float)M_PI - *r_a2 - *r_a1;
-}
-
-static void p_face_angles(PFace *f, float *r_a1, float *r_a2, float *r_a3)
-{
-  PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-  PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
-
-  p_triangle_angles(v1->co, v2->co, v3->co, r_a1, r_a2, r_a3);
-}
-
-static float p_face_area(PFace *f)
-{
-  PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-  PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
-
-  return area_tri_v3(v1->co, v2->co, v3->co);
-}
-
-static float p_area_signed(const float v1[2], const float v2[2], const float v3[2])
-{
-  return 0.5f * (((v2[0] - v1[0]) * (v3[1] - v1[1])) - ((v3[0] - v1[0]) * (v2[1] - v1[1])));
-}
-
-static float p_face_uv_area_signed(PFace *f)
-{
-  PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-  PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
-
-  return 0.5f * (((v2->uv[0] - v1->uv[0]) * (v3->uv[1] - v1->uv[1])) -
-                 ((v3->uv[0] - v1->uv[0]) * (v2->uv[1] - v1->uv[1])));
-}
-
-static float p_edge_length(PEdge *e)
-{
-  PVert *v1 = e->vert, *v2 = e->next->vert;
-  float d[3];
-
-  d[0] = v2->co[0] - v1->co[0];
-  d[1] = v2->co[1] - v1->co[1];
-  d[2] = v2->co[2] - v1->co[2];
-
-  return sqrtf(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
-}
-
-static float p_edge_uv_length(PEdge *e)
-{
-  PVert *v1 = e->vert, *v2 = e->next->vert;
-  float d[3];
-
-  d[0] = v2->uv[0] - v1->uv[0];
-  d[1] = v2->uv[1] - v1->uv[1];
-
-  return sqrtf(d[0] * d[0] + d[1] * d[1]);
-}
-
-static void p_chart_uv_bbox(PChart *chart, float minv[2], float maxv[2])
-{
-  PVert *v;
-
-  INIT_MINMAX2(minv, maxv);
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    minmax_v2v2_v2(minv, maxv, v->uv);
-  }
-}
-
-static float p_chart_uv_area(PChart *chart)
-{
-  float area = 0.0f;
-
-  for (PFace *f = chart->faces; f; f = f->nextlink) {
-    area += fabsf(p_face_uv_area_signed(f));
-  }
-
-  return area;
-}
-
-static void p_chart_uv_scale(PChart *chart, float scale)
-{
-  PVert *v;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    v->uv[0] *= scale;
-    v->uv[1] *= scale;
-  }
-}
-
-static void p_chart_uv_scale_xy(PChart *chart, float x, float y)
-{
-  PVert *v;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    v->uv[0] *= x;
-    v->uv[1] *= y;
-  }
-}
-
-static void p_chart_uv_translate(PChart *chart, const float trans[2])
-{
-  PVert *v;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    v->uv[0] += trans[0];
-    v->uv[1] += trans[1];
-  }
-}
-
-static void p_chart_uv_transform(PChart *chart, const float mat[2][2])
-{
-  PVert *v;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    mul_m2_v2(mat, v->uv);
-  }
-}
-
-static void p_chart_uv_to_array(PChart *chart, float (*points)[2])
-{
-  PVert *v;
-  uint i = 0;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    copy_v2_v2(points[i++], v->uv);
-  }
-}
-
-static void UNUSED_FUNCTION(p_chart_uv_from_array)(PChart *chart, float (*points)[2])
-{
-  PVert *v;
-  uint i = 0;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    copy_v2_v2(v->uv, points[i++]);
-  }
-}
-
-static PBool p_intersect_line_2d_dir(const float v1[2],
-                                     const float dir1[2],
-                                     const float v2[2],
-                                     const float dir2[2],
-                                     float r_isect[2])
-{
-  float lmbda, div;
-
-  div = dir2[0] * dir1[1] - dir2[1] * dir1[0];
-
-  if (div == 0.0f) {
-    return P_FALSE;
-  }
-
-  lmbda = ((v1[1] - v2[1]) * dir1[0] - (v1[0] - v2[0]) * dir1[1]) / div;
-  r_isect[0] = v1[0] + lmbda * dir2[0];
-  r_isect[1] = v1[1] + lmbda * dir2[1];
-
-  return P_TRUE;
-}
-
-#if 0
-static PBool p_intersect_line_2d(const float v1[2],
-                                 const float v2[2],
-                                 const float v3[2],
-                                 const float v4[2],
-                                 const float r_isect[2])
-{
-  float dir1[2], dir2[2];
-
-  dir1[0] = v4[0] - v3[0];
-  dir1[1] = v4[1] - v3[1];
-
-  dir2[0] = v2[0] - v1[0];
-  dir2[1] = v2[1] - v1[1];
-
-  if (!p_intersect_line_2d_dir(v1, dir1, v2, dir2, isect)) {
-    /* parallel - should never happen in theory for polygon kernel, but
-     * let's give a point nearby in case things go wrong */
-    isect[0] = (v1[0] + v2[0]) * 0.5f;
-    isect[1] = (v1[1] + v2[1]) * 0.5f;
-    return P_FALSE;
-  }
-
-  return P_TRUE;
-}
-#endif
-
-/* Topological Utilities */
-
-static PEdge *p_wheel_edge_next(PEdge *e)
-{
-  return e->next->next->pair;
-}
-
-static PEdge *p_wheel_edge_prev(PEdge *e)
-{
-  return (e->pair) ? e->pair->next : NULL;
-}
-
-static PEdge *p_boundary_edge_next(PEdge *e)
-{
-  return e->next->vert->edge;
-}
-
-static PEdge *p_boundary_edge_prev(PEdge *e)
-{
-  PEdge *we = e, *last;
-
-  do {
-    last = we;
-    we = p_wheel_edge_next(we);
-  } while (we && (we != e));
-
-  return last->next->next;
-}
-
-static PBool p_vert_interior(PVert *v)
-{
-  return (v->edge->pair != NULL);
-}
-
-static void p_face_flip(PFace *f)
-{
-  PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-  PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
-  int f1 = e1->flag, f2 = e2->flag, f3 = e3->flag;
-  float *orig_uv1 = e1->orig_uv, *orig_uv2 = e2->orig_uv, *orig_uv3 = e3->orig_uv;
-
-  e1->vert = v2;
-  e1->next = e3;
-  e1->orig_uv = orig_uv2;
-  e1->flag = (f1 & ~PEDGE_VERTEX_FLAGS) | (f2 & PEDGE_VERTEX_FLAGS);
-
-  e2->vert = v3;
-  e2->next = e1;
-  e2->orig_uv = orig_uv3;
-  e2->flag = (f2 & ~PEDGE_VERTEX_FLAGS) | (f3 & PEDGE_VERTEX_FLAGS);
-
-  e3->vert = v1;
-  e3->next = e2;
-  e3->orig_uv = orig_uv1;
-  e3->flag = (f3 & ~PEDGE_VERTEX_FLAGS) | (f1 & PEDGE_VERTEX_FLAGS);
-}
-
-#if 0
-static void p_chart_topological_sanity_check(PChart *chart)
-{
-  PVert *v;
-  PEdge *e;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    param_test_equals_ptr("v->edge->vert", v, v->edge->vert);
-  }
-
-  for (e = chart->edges; e; e = e->nextlink) {
-    if (e->pair) {
-      param_test_equals_ptr("e->pair->pair", e, e->pair->pair);
-      param_test_equals_ptr("pair->vert", e->vert, e->pair->next->vert);
-      param_test_equals_ptr("pair->next->vert", e->next->vert, e->pair->vert);
-    }
-  }
-}
-#endif
-
-/* Loading / Flushing */
-
-static void p_vert_load_pin_select_uvs(PHandle *handle, PVert *v)
-{
-  PEdge *e;
-  int nedges = 0, npins = 0;
-  float pinuv[2];
-
-  v->uv[0] = v->uv[1] = 0.0f;
-  pinuv[0] = pinuv[1] = 0.0f;
-  e = v->edge;
-  do {
-    if (e->orig_uv) {
-      if (e->flag & PEDGE_SELECT) {
-        v->flag |= PVERT_SELECT;
-      }
-
-      if (e->flag & PEDGE_PIN) {
-        pinuv[0] += e->orig_uv[0] * handle->aspx;
-        pinuv[1] += e->orig_uv[1] * handle->aspy;
-        npins++;
-      }
-      else {
-        v->uv[0] += e->orig_uv[0] * handle->aspx;
-        v->uv[1] += e->orig_uv[1] * handle->aspy;
-      }
-
-      nedges++;
-    }
-
-    e = p_wheel_edge_next(e);
-  } while (e && e != (v->edge));
-
-  if (npins > 0) {
-    v->uv[0] = pinuv[0] / npins;
-    v->uv[1] = pinuv[1] / npins;
-    v->flag |= PVERT_PIN;
-  }
-  else if (nedges > 0) {
-    v->uv[0] /= nedges;
-    v->uv[1] /= nedges;
-  }
-}
-
-static void p_flush_uvs(PHandle *handle, PChart *chart)
-{
-  PEdge *e;
-
-  for (e = chart->edges; e; e = e->nextlink) {
-    if (e->orig_uv) {
-      e->orig_uv[0] = e->vert->uv[0] / handle->aspx;
-      e->orig_uv[1] = e->vert->uv[1] / handle->aspy;
-    }
-  }
-}
-
-static void p_flush_uvs_blend(PHandle *handle, PChart *chart, float blend)
-{
-  PEdge *e;
-  float invblend = 1.0f - blend;
-
-  for (e = chart->edges; e; e = e->nextlink) {
-    if (e->orig_uv) {
-      e->orig_uv[0] = blend * e->old_uv[0] + invblend * e->vert->uv[0] / handle->aspx;
-      e->orig_uv[1] = blend * e->old_uv[1] + invblend * e->vert->uv[1] / handle->aspy;
-    }
-  }
-}
-
-static void p_face_backup_uvs(PFace *f)
-{
-  PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-
-  if (e1->orig_uv) {
-    e1->old_uv[0] = e1->orig_uv[0];
-    e1->old_uv[1] = e1->orig_uv[1];
-  }
-  if (e2->orig_uv) {
-    e2->old_uv[0] = e2->orig_uv[0];
-    e2->old_uv[1] = e2->orig_uv[1];
-  }
-  if (e3->orig_uv) {
-    e3->old_uv[0] = e3->orig_uv[0];
-    e3->old_uv[1] = e3->orig_uv[1];
-  }
-}
-
-static void p_face_restore_uvs(PFace *f)
-{
-  PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-
-  if (e1->orig_uv) {
-    e1->orig_uv[0] = e1->old_uv[0];
-    e1->orig_uv[1] = e1->old_uv[1];
-  }
-  if (e2->orig_uv) {
-    e2->orig_uv[0] = e2->old_uv[0];
-    e2->orig_uv[1] = e2->old_uv[1];
-  }
-  if (e3->orig_uv) {
-    e3->orig_uv[0] = e3->old_uv[0];
-    e3->orig_uv[1] = e3->old_uv[1];
-  }
-}
-
-/* Construction (use only during construction, relies on u.key being set */
-
-static PVert *p_vert_add(PHandle *handle, PHashKey key, const float co[3], PEdge *e)
-{
-  PVert *v = (PVert *)BLI_memarena_alloc(handle->arena, sizeof(*v));
-  copy_v3_v3(v->co, co);
-
-  /* Sanity check, a single nan/inf point causes the entire result to be invalid.
-   * Note that values within the calculation may _become_ non-finite,
-   * so the rest of the code still needs to take this possibility into account. */
-  for (int i = 0; i < 3; i++) {
-    if (UNLIKELY(!isfinite(v->co[i]))) {
-      v->co[i] = 0.0f;
-    }
-  }
-
-  v->u.key = key;
-  v->edge = e;
-  v->flag = 0;
-
-  phash_insert(handle->hash_verts, (PHashLink *)v);
-
-  return v;
-}
-
-static PVert *p_vert_lookup(PHandle *handle, PHashKey key, const float co[3], PEdge *e)
-{
-  PVert *v = (PVert *)phash_lookup(handle->hash_verts, key);
-
-  if (v) {
-    return v;
-  }
-  return p_vert_add(handle, key, co, e);
-}
-
-static PVert *p_vert_copy(PChart *chart, PVert *v)
-{
-  PVert *nv = (PVert *)BLI_memarena_alloc(chart->handle->arena, sizeof(*nv));
-
-  copy_v3_v3(nv->co, v->co);
-  nv->uv[0] = v->uv[0];
-  nv->uv[1] = v->uv[1];
-  nv->u.key = v->u.key;
-  nv->edge = v->edge;
-  nv->flag = v->flag;
-
-  return nv;
-}
-
-static PEdge *p_edge_lookup(PHandle *handle, const PHashKey *vkeys)
-{
-  PHashKey key = PHASH_edge(vkeys[0], vkeys[1]);
-  PEdge *e = (PEdge *)phash_lookup(handle->hash_edges, key);
-
-  while (e) {
-    if ((e->vert->u.key == vkeys[0]) && (e->next->vert->u.key == vkeys[1])) {
-      return e;
-    }
-    if ((e->vert->u.key == vkeys[1]) && (e->next->vert->u.key == vkeys[0])) {
-      return e;
-    }
-
-    e = (PEdge *)phash_next(handle->hash_edges, key, (PHashLink *)e);
-  }
-
-  return NULL;
-}
-
-static int p_face_exists(ParamHandle *phandle, ParamKey *pvkeys, int i1, int i2, int i3)
-{
-  PHandle *handle = (PHandle *)phandle;
-  PHashKey *vkeys = (PHashKey *)pvkeys;
-  PHashKey key = PHASH_edge(vkeys[i1], vkeys[i2]);
-  PEdge *e = (PEdge *)phash_lookup(handle->hash_edges, key);
-
-  while (e) {
-    if ((e->vert->u.key == vkeys[i1]) && (e->next->vert->u.key == vkeys[i2])) {
-      if (e->next->next->vert->u.key == vkeys[i3]) {
-        return P_TRUE;
-      }
-    }
-    else if ((e->vert->u.key == vkeys[i2]) && (e->next->vert->u.key == vkeys[i1])) {
-      if (e->next->next->vert->u.key == vkeys[i3]) {
-        return P_TRUE;
-      }
-    }
-
-    e = (PEdge *)phash_next(handle->hash_edges, key, (PHashLink *)e);
-  }
-
-  return P_FALSE;
-}
-
-static PChart *p_chart_new(PHandle *handle)
-{
-  PChart *chart = (PChart *)MEM_callocN(sizeof(*chart), "PChart");
-  chart->handle = handle;
-
-  return chart;
-}
-
-static void p_chart_delete(PChart *chart)
-{
-  /* the actual links are free by memarena */
-  MEM_freeN(chart);
-}
-
-static PBool p_edge_implicit_seam(PEdge *e, PEdge *ep)
-{
-  float *uv1, *uv2, *uvp1, *uvp2;
-  float limit[2];
-
-  limit[0] = 0.00001;
-  limit[1] = 0.00001;
-
-  uv1 = e->orig_uv;
-  uv2 = e->next->orig_uv;
-
-  if (e->vert->u.key == ep->vert->u.key) {
-    uvp1 = ep->orig_uv;
-    uvp2 = ep->next->orig_uv;
-  }
-  else {
-    uvp1 = ep->next->orig_uv;
-    uvp2 = ep->orig_uv;
-  }
-
-  if ((fabsf(uv1[0] - uvp1[0]) > limit[0]) || (fabsf(uv1[1] - uvp1[1]) > limit[1])) {
-    e->flag |= PEDGE_SEAM;
-    ep->flag |= PEDGE_SEAM;
-    return P_TRUE;
-  }
-  if ((fabsf(uv2[0] - uvp2[0]) > limit[0]) || (fabsf(uv2[1] - uvp2[1]) > limit[1])) {
-    e->flag |= PEDGE_SEAM;
-    ep->flag |= PEDGE_SEAM;
-    return P_TRUE;
-  }
-
-  return P_FALSE;
-}
-
-static PBool p_edge_has_pair(PHandle *handle, PEdge *e, PBool topology_from_uvs, PEdge **r_pair)
-{
-  PHashKey key;
-  PEdge *pe;
-  PVert *v1, *v2;
-  PHashKey key1 = e->vert->u.key;
-  PHashKey key2 = e->next->vert->u.key;
-
-  if (e->flag & PEDGE_SEAM) {
-    return P_FALSE;
-  }
-
-  key = PHASH_edge(key1, key2);
-  pe = (PEdge *)phash_lookup(handle->hash_edges, key);
-  *r_pair = NULL;
-
-  while (pe) {
-    if (pe != e) {
-      v1 = pe->vert;
-      v2 = pe->next->vert;
-
-      if (((v1->u.key == key1) && (v2->u.key == key2)) ||
-          ((v1->u.key == key2) && (v2->u.key == key1))) {
-
-        /* don't connect seams and t-junctions */
-        if ((pe->flag & PEDGE_SEAM) || *r_pair ||
-            (topology_from_uvs && p_edge_implicit_seam(e, pe))) {
-          *r_pair = NULL;
-          return P_FALSE;
-        }
-
-        *r_pair = pe;
-      }
-    }
-
-    pe = (PEdge *)phash_next(handle->hash_edges, key, (PHashLink *)pe);
-  }
-
-  if (*r_pair && (e->vert == (*r_pair)->vert)) {
-    if ((*r_pair)->next->pair || (*r_pair)->next->next->pair) {
-      /* non unfoldable, maybe mobius ring or klein bottle */
-      *r_pair = NULL;
-      return P_FALSE;
-    }
-  }
-
-  return (*r_pair != NULL);
-}
-
-static PBool p_edge_connect_pair(PHandle *handle,
-                                 PEdge *e,
-                                 PBool topology_from_uvs,
-                                 PEdge ***stack)
-{
-  PEdge *pair = NULL;
-
-  if (!e->pair && p_edge_has_pair(handle, e, topology_from_uvs, &pair)) {
-    if (e->vert == pair->vert) {
-      p_face_flip(pair->face);
-    }
-
-    e->pair = pair;
-    pair->pair = e;
-
-    if (!(pair->face->flag & PFACE_CONNECTED)) {
-      **stack = pair;
-      (*stack)++;
-    }
-  }
-
-  return (e->pair != NULL);
-}
-
-static int p_connect_pairs(PHandle *handle, PBool topology_from_uvs)
-{
-  PEdge **stackbase = MEM_mallocN(sizeof(*stackbase) * phash_size(handle->hash_faces),
-                                  "Pstackbase");
-  PEdge **stack = stackbase;
-  PFace *f, *first;
-  PEdge *e, *e1, *e2;
-  PChart *chart = handle->construction_chart;
-  int ncharts = 0;
-
-  /* Connect pairs, count edges, set vertex-edge pointer to a pair-less edge. */
-  for (first = chart->faces; first; first = first->nextlink) {
-    if (first->flag & PFACE_CONNECTED) {
-      continue;
-    }
-
-    *stack = first->edge;
-    stack++;
-
-    while (stack != stackbase) {
-      stack--;
-      e = *stack;
-      e1 = e->next;
-      e2 = e1->next;
-
-      f = e->face;
-      f->flag |= PFACE_CONNECTED;
-
-      /* assign verts to charts so we can sort them later */
-      f->u.chart = ncharts;
-
-      if (!p_edge_connect_pair(handle, e, topology_from_uvs, &stack)) {
-        e->vert->edge = e;
-      }
-      if (!p_edge_connect_pair(handle, e1, topology_from_uvs, &stack)) {
-        e1->vert->edge = e1;
-      }
-      if (!p_edge_connect_pair(handle, e2, topology_from_uvs, &stack)) {
-        e2->vert->edge = e2;
-      }
-    }
-
-    ncharts++;
-  }
-
-  MEM_freeN(stackbase);
-
-  return ncharts;
-}
-
-static void p_split_vert(PChart *chart, PEdge *e)
-{
-  PEdge *we, *lastwe = NULL;
-  PVert *v = e->vert;
-  PBool copy = P_TRUE;
-
-  if (e->flag & PEDGE_PIN) {
-    chart->flag |= PCHART_HAS_PINS;
-  }
-
-  if (e->flag & PEDGE_VERTEX_SPLIT) {
-    return;
-  }
-
-  /* rewind to start */
-  lastwe = e;
-  for (we = p_wheel_edge_prev(e); we && (we != e); we = p_wheel_edge_prev(we)) {
-    lastwe = we;
-  }
-
-  /* go over all edges in wheel */
-  for (we = lastwe; we; we = p_wheel_edge_next(we)) {
-    if (we->flag & PEDGE_VERTEX_SPLIT) {
-      break;
-    }
-
-    we->flag |= PEDGE_VERTEX_SPLIT;
-
-    if (we == v->edge) {
-      /* found it, no need to copy */
-      copy = P_FALSE;
-      v->nextlink = chart->verts;
-      chart->verts = v;
-      chart->nverts++;
-    }
-  }
-
-  if (copy) {
-    /* not found, copying */
-    v->flag |= PVERT_SPLIT;
-    v = p_vert_copy(chart, v);
-    v->flag |= PVERT_SPLIT;
-
-    v->nextlink = chart->verts;
-    chart->verts = v;
-    chart->nverts++;
-
-    v->edge = lastwe;
-
-    we = lastwe;
-    do {
-      we->vert = v;
-      we = p_wheel_edge_next(we);
-    } while (we && (we != lastwe));
-  }
-}
-
-static PChart **p_split_charts(PHandle *handle, PChart *chart, int ncharts)
-{
-  PChart **charts = MEM_mallocN(sizeof(*charts) * ncharts, "PCharts"), *nchart;
-  PFace *f, *nextf;
-  int i;
-
-  for (i = 0; i < ncharts; i++) {
-    charts[i] = p_chart_new(handle);
-  }
-
-  f = chart->faces;
-  while (f) {
-    PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-    nextf = f->nextlink;
-
-    nchart = charts[f->u.chart];
-
-    f->nextlink = nchart->faces;
-    nchart->faces = f;
-    e1->nextlink = nchart->edges;
-    nchart->edges = e1;
-    e2->nextlink = nchart->edges;
-    nchart->edges = e2;
-    e3->nextlink = nchart->edges;
-    nchart->edges = e3;
-
-    nchart->nfaces++;
-    nchart->nedges += 3;
-
-    p_split_vert(nchart, e1);
-    p_split_vert(nchart, e2);
-    p_split_vert(nchart, e3);
-
-    f = nextf;
-  }
-
-  return charts;
-}
-
-static PFace *p_face_add(PHandle *handle)
-{
-  PFace *f;
-  PEdge *e1, *e2, *e3;
-
-  /* allocate */
-  f = (PFace *)BLI_memarena_alloc(handle->arena, sizeof(*f));
-  f->flag = 0; /* init ! */
-
-  e1 = (PEdge *)BLI_memarena_alloc(handle->arena, sizeof(*e1));
-  e2 = (PEdge *)BLI_memarena_alloc(handle->arena, sizeof(*e2));
-  e3 = (PEdge *)BLI_memarena_alloc(handle->arena, sizeof(*e3));
-
-  /* set up edges */
-  f->edge = e1;
-  e1->face = e2->face = e3->face = f;
-
-  e1->next = e2;
-  e2->next = e3;
-  e3->next = e1;
-
-  e1->pair = NULL;
-  e2->pair = NULL;
-  e3->pair = NULL;
-
-  e1->flag = 0;
-  e2->flag = 0;
-  e3->flag = 0;
-
-  return f;
-}
-
-static PFace *p_face_add_construct(PHandle *handle,
-                                   ParamKey key,
-                                   const ParamKey *vkeys,
-                                   float *co[4],
-                                   float *uv[4],
-                                   int i1,
-                                   int i2,
-                                   int i3,
-                                   const ParamBool *pin,
-                                   const ParamBool *select)
-{
-  PFace *f = p_face_add(handle);
-  PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-
-  e1->vert = p_vert_lookup(handle, vkeys[i1], co[i1], e1);
-  e2->vert = p_vert_lookup(handle, vkeys[i2], co[i2], e2);
-  e3->vert = p_vert_lookup(handle, vkeys[i3], co[i3], e3);
-
-  e1->orig_uv = uv[i1];
-  e2->orig_uv = uv[i2];
-  e3->orig_uv = uv[i3];
-
-  if (pin) {
-    if (pin[i1]) {
-      e1->flag |= PEDGE_PIN;
-    }
-    if (pin[i2]) {
-      e2->flag |= PEDGE_PIN;
-    }
-    if (pin[i3]) {
-      e3->flag |= PEDGE_PIN;
-    }
-  }
-
-  if (select) {
-    if (select[i1]) {
-      e1->flag |= PEDGE_SELECT;
-    }
-    if (select[i2]) {
-      e2->flag |= PEDGE_SELECT;
-    }
-    if (select[i3]) {
-      e3->flag |= PEDGE_SELECT;
-    }
-  }
-
-  /* insert into hash */
-  f->u.key = key;
-  phash_insert(handle->hash_faces, (PHashLink *)f);
-
-  e1->u.key = PHASH_edge(vkeys[i1], vkeys[i2]);
-  e2->u.key = PHASH_edge(vkeys[i2], vkeys[i3]);
-  e3->u.key = PHASH_edge(vkeys[i3], vkeys[i1]);
-
-  phash_insert(handle->hash_edges, (PHashLink *)e1);
-  phash_insert(handle->hash_edges, (PHashLink *)e2);
-  phash_insert(handle->hash_edges, (PHashLink *)e3);
-
-  return f;
-}
-
-static PFace *p_face_add_fill(PChart *chart, PVert *v1, PVert *v2, PVert *v3)
-{
-  PFace *f = p_face_add(chart->handle);
-  PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-
-  e1->vert = v1;
-  e2->vert = v2;
-  e3->vert = v3;
-
-  e1->orig_uv = e2->orig_uv = e3->orig_uv = NULL;
-
-  f->nextlink = chart->faces;
-  chart->faces = f;
-  e1->nextlink = chart->edges;
-  chart->edges = e1;
-  e2->nextlink = chart->edges;
-  chart->edges = e2;
-  e3->nextlink = chart->edges;
-  chart->edges = e3;
-
-  chart->nfaces++;
-  chart->nedges += 3;
-
-  return f;
-}
-
-static PBool p_quad_split_direction(PHandle *handle, float **co, PHashKey *vkeys)
-{
-  /* Slight bias to prefer one edge over the other in case they are equal, so
-   * that in symmetric models we choose the same split direction instead of
-   * depending on floating point errors to decide. */
-  float bias = 1.0f + 1e-6f;
-  float fac = len_v3v3(co[0], co[2]) * bias - len_v3v3(co[1], co[3]);
-  PBool dir = (fac <= 0.0f);
-
-  /* The face exists check is there because of a special case:
-   * when two quads share three vertices, they can each be split into two triangles,
-   * resulting in two identical triangles. For example in Suzanne's nose. */
-  if (dir) {
-    if (p_face_exists(handle, vkeys, 0, 1, 2) || p_face_exists(handle, vkeys, 0, 2, 3)) {
-      return !dir;
-    }
-  }
-  else {
-    if (p_face_exists(handle, vkeys, 0, 1, 3) || p_face_exists(handle, vkeys, 1, 2, 3)) {
-      return !dir;
-    }
-  }
-
-  return dir;
-}
-
-/* Construction: boundary filling */
-
-static void p_chart_boundaries(PChart *chart, int *r_nboundaries, PEdge **r_outer)
-{
-  PEdge *e, *be;
-  float len, maxlen = -1.0;
-
-  if (r_nboundaries) {
-    *r_nboundaries = 0;
-  }
-  if (r_outer) {
-    *r_outer = NULL;
-  }
-
-  for (e = chart->edges; e; e = e->nextlink) {
-    if (e->pair || (e->flag & PEDGE_DONE)) {
-      continue;
-    }
-
-    if (r_nboundaries) {
-      (*r_nboundaries)++;
-    }
-
-    len = 0.0f;
-
-    be = e;
-    do {
-      be->flag |= PEDGE_DONE;
-      len += p_edge_length(be);
-      be = be->next->vert->edge;
-    } while (be != e);
-
-    if (r_outer && (len > maxlen)) {
-      *r_outer = e;
-      maxlen = len;
-    }
-  }
-
-  for (e = chart->edges; e; e = e->nextlink) {
-    e->flag &= ~PEDGE_DONE;
-  }
-}
-
-static float p_edge_boundary_angle(PEdge *e)
-{
-  PEdge *we;
-  PVert *v, *v1, *v2;
-  float angle;
-  int n = 0;
-
-  v = e->vert;
-
-  /* concave angle check -- could be better */
-  angle = M_PI;
-
-  we = v->edge;
-  do {
-    v1 = we->next->vert;
-    v2 = we->next->next->vert;
-    angle -= p_vec_angle(v1->co, v->co, v2->co);
-
-    we = we->next->next->pair;
-    n++;
-  } while (we && (we != v->edge));
-
-  return angle;
-}
-
-static void p_chart_fill_boundary(PChart *chart, PEdge *be, int nedges)
-{
-  PEdge *e, *e1, *e2;
-
-  PFace *f;
-  struct Heap *heap = BLI_heap_new();
-  float angle;
-
-  e = be;
-  do {
-    angle = p_edge_boundary_angle(e);
-    e->u.heaplink = BLI_heap_insert(heap, angle, e);
-
-    e = p_boundary_edge_next(e);
-  } while (e != be);
-
-  if (nedges == 2) {
-    /* no real boundary, but an isolated seam */
-    e = be->next->vert->edge;
-    e->pair = be;
-    be->pair = e;
-
-    BLI_heap_remove(heap, e->u.heaplink);
-    BLI_heap_remove(heap, be->u.heaplink);
-  }
-  else {
-    while (nedges > 2) {
-      PEdge *ne, *ne1, *ne2;
-
-      e = (PEdge *)BLI_heap_pop_min(heap);
-
-      e1 = p_boundary_edge_prev(e);
-      e2 = p_boundary_edge_next(e);
-
-      BLI_heap_remove(heap, e1->u.heaplink);
-      BLI_heap_remove(heap, e2->u.heaplink);
-      e->u.heaplink = e1->u.heaplink = e2->u.heaplink = NULL;
-
-      e->flag |= PEDGE_FILLED;
-      e1->flag |= PEDGE_FILLED;
-
-      f = p_face_add_fill(chart, e->vert, e1->vert, e2->vert);
-      f->flag |= PFACE_FILLED;
-
-      ne = f->edge->next->next;
-      ne1 = f->edge;
-      ne2 = f->edge->next;
-
-      ne->flag = ne1->flag = ne2->flag = PEDGE_FILLED;
-
-      e->pair = ne;
-      ne->pair = e;
-      e1->pair = ne1;
-      ne1->pair = e1;
-
-      ne->vert = e2->vert;
-      ne1->vert = e->vert;
-      ne2->vert = e1->vert;
-
-      if (nedges == 3) {
-        e2->pair = ne2;
-        ne2->pair = e2;
-      }
-      else {
-        ne2->vert->edge = ne2;
-
-        ne2->u.heaplink = BLI_heap_insert(heap, p_edge_boundary_angle(ne2), ne2);
-        e2->u.heaplink = BLI_heap_insert(heap, p_edge_boundary_angle(e2), e2);
-      }
-
-      nedges--;
-    }
-  }
-
-  BLI_heap_free(heap, NULL);
-}
-
-static void p_chart_fill_boundaries(PChart *chart, PEdge *outer)
-{
-  PEdge *e, *be; /* *enext - as yet unused */
-  int nedges;
-
-  for (e = chart->edges; e; e = e->nextlink) {
-    /* enext = e->nextlink; - as yet unused */
-
-    if (e->pair || (e->flag & PEDGE_FILLED)) {
-      continue;
-    }
-
-    nedges = 0;
-    be = e;
-    do {
-      be->flag |= PEDGE_FILLED;
-      be = be->next->vert->edge;
-      nedges++;
-    } while (be != e);
-
-    if (e != outer) {
-      p_chart_fill_boundary(chart, e, nedges);
-    }
-  }
-}
-
-#if 0
-/* Polygon kernel for inserting uv's non overlapping */
-
-static int p_polygon_point_in(const float cp1[2], const float cp2[2], const float p[2])
-{
-  if ((cp1[0] == p[0]) && (cp1[1] == p[1])) {
-    return 2;
-  }
-  else if ((cp2[0] == p[0]) && (cp2[1] == p[1])) {
-    return 3;
-  }
-  else {
-    return (p_area_signed(cp1, cp2, p) >= 0.0f);
-  }
-}
-
-static void p_polygon_kernel_clip(float (*oldpoints)[2],
-                                  int noldpoints,
-                                  float (*newpoints)[2],
-                                  int *r_nnewpoints,
-                                  const float cp1[2],
-                                  const float cp2[2])
-{
-  float *p2, *p1, isect[2];
-  int i, p2in, p1in;
-
-  p1 = oldpoints[noldpoints - 1];
-  p1in = p_polygon_point_in(cp1, cp2, p1);
-  *r_nnewpoints = 0;
-
-  for (i = 0; i < noldpoints; i++) {
-    p2 = oldpoints[i];
-    p2in = p_polygon_point_in(cp1, cp2, p2);
-
-    if ((p2in >= 2) || (p1in && p2in)) {
-      newpoints[*r_nnewpoints][0] = p2[0];
-      newpoints[*r_nnewpoints][1] = p2[1];
-      (*r_nnewpoints)++;
-    }
-    else if (p1in && !p2in) {
-      if (p1in != 3) {
-        p_intersect_line_2d(p1, p2, cp1, cp2, isect);
-        newpoints[*r_nnewpoints][0] = isect[0];
-        newpoints[*r_nnewpoints][1] = isect[1];
-        (*r_nnewpoints)++;
-      }
-    }
-    else if (!p1in && p2in) {
-      p_intersect_line_2d(p1, p2, cp1, cp2, isect);
-      newpoints[*r_nnewpoints][0] = isect[0];
-      newpoints[*r_nnewpoints][1] = isect[1];
-      (*r_nnewpoints)++;
-
-      newpoints[*r_nnewpoints][0] = p2[0];
-      newpoints[*r_nnewpoints][1] = p2[1];
-      (*r_nnewpoints)++;
-    }
-
-    p1in = p2in;
-    p1 = p2;
-  }
-}
-
-static void p_polygon_kernel_center(float (*points)[2], int npoints, float *center)
-{
-  int i, size, nnewpoints = npoints;
-  float(*oldpoints)[2], (*newpoints)[2], *p1, *p2;
-
-  size = npoints * 3;
-  oldpoints = MEM_mallocN(sizeof(float[2]) * size, "PPolygonOldPoints");
-  newpoints = MEM_mallocN(sizeof(float[2]) * size, "PPolygonNewPoints");
-
-  memcpy(oldpoints, points, sizeof(float[2]) * npoints);
-
-  for (i = 0; i < npoints; i++) {
-    p1 = points[i];
-    p2 = points[(i + 1) % npoints];
-    p_polygon_kernel_clip(oldpoints, nnewpoints, newpoints, &nnewpoints, p1, p2);
-
-    if (nnewpoints == 0) {
-      /* degenerate case, use center of original polygon */
-      memcpy(oldpoints, points, sizeof(float[2]) * npoints);
-      nnewpoints = npoints;
-      break;
-    }
-    else if (nnewpoints == 1) {
-      /* degenerate case, use remaining point */
-      center[0] = newpoints[0][0];
-      center[1] = newpoints[0][1];
-
-      MEM_freeN(oldpoints);
-      MEM_freeN(newpoints);
-
-      return;
-    }
-
-    if (nnewpoints * 2 > size) {
-      size *= 2;
-      MEM_freeN(oldpoints);
-      oldpoints = MEM_mallocN(sizeof(float[2]) * size, "oldpoints");
-      memcpy(oldpoints, newpoints, sizeof(float[2]) * nnewpoints);
-      MEM_freeN(newpoints);
-      newpoints = MEM_mallocN(sizeof(float[2]) * size, "newpoints");
-    }
-    else {
-      float(*sw_points)[2] = oldpoints;
-      oldpoints = newpoints;
-      newpoints = sw_points;
-    }
-  }
-
-  center[0] = center[1] = 0.0f;
-
-  for (i = 0; i < nnewpoints; i++) {
-    center[0] += oldpoints[i][0];
-    center[1] += oldpoints[i][1];
-  }
-
-  center[0] /= nnewpoints;
-  center[1] /= nnewpoints;
-
-  MEM_freeN(oldpoints);
-  MEM_freeN(newpoints);
-}
-#endif
-
-#if 0
-/* Edge Collapser */
-
-int NCOLLAPSE = 1;
-int NCOLLAPSEX = 0;
-
-static float p_vert_cotan(const float v1[3], const float v2[3], const float v3[3])
-{
-  float a[3], b[3], c[3], clen;
-
-  sub_v3_v3v3(a, v2, v1);
-  sub_v3_v3v3(b, v3, v1);
-  cross_v3_v3v3(c, a, b);
-
-  clen = len_v3(c);
-
-  if (clen == 0.0f) {
-    return 0.0f;
-  }
-
-  return dot_v3v3(a, b) / clen;
-}
-
-static PBool p_vert_flipped_wheel_triangle(PVert *v)
-{
-  PEdge *e = v->edge;
-
-  do {
-    if (p_face_uv_area_signed(e->face) < 0.0f) {
-      return P_TRUE;
-    }
-
-    e = p_wheel_edge_next(e);
-  } while (e && (e != v->edge));
-
-  return P_FALSE;
-}
-
-static PBool p_vert_map_harmonic_weights(PVert *v)
-{
-  float weightsum, positionsum[2], olduv[2];
-
-  weightsum = 0.0f;
-  positionsum[0] = positionsum[1] = 0.0f;
-
-  if (p_vert_interior(v)) {
-    PEdge *e = v->edge;
-
-    do {
-      float t1, t2, weight;
-      PVert *v1, *v2;
-
-      v1 = e->next->vert;
-      v2 = e->next->next->vert;
-      t1 = p_vert_cotan(v2->co, e->vert->co, v1->co);
-
-      v1 = e->pair->next->vert;
-      v2 = e->pair->next->next->vert;
-      t2 = p_vert_cotan(v2->co, e->pair->vert->co, v1->co);
-
-      weight = 0.5f * (t1 + t2);
-      weightsum += weight;
-      positionsum[0] += weight * e->pair->vert->uv[0];
-      positionsum[1] += weight * e->pair->vert->uv[1];
-
-      e = p_wheel_edge_next(e);
-    } while (e && (e != v->edge));
-  }
-  else {
-    PEdge *e = v->edge;
-
-    do {
-      float t1, t2;
-      PVert *v1, *v2;
-
-      v2 = e->next->vert;
-      v1 = e->next->next->vert;
-
-      t1 = p_vert_cotan(v1->co, v->co, v2->co);
-      t2 = p_vert_cotan(v2->co, v->co, v1->co);
-
-      weightsum += t1 + t2;
-      positionsum[0] += (v2->uv[1] - v1->uv[1]) + (t1 * v2->uv[0] + t2 * v1->uv[0]);
-      positionsum[1] += (v1->uv[0] - v2->uv[0]) + (t1 * v2->uv[1] + t2 * v1->uv[1]);
-
-      e = p_wheel_edge_next(e);
-    } while (e && (e != v->edge));
-  }
-
-  if (weightsum != 0.0f) {
-    weightsum = 1.0f / weightsum;
-    positionsum[0] *= weightsum;
-    positionsum[1] *= weightsum;
-  }
-
-  olduv[0] = v->uv[0];
-  olduv[1] = v->uv[1];
-  v->uv[0] = positionsum[0];
-  v->uv[1] = positionsum[1];
-
-  if (p_vert_flipped_wheel_triangle(v)) {
-    v->uv[0] = olduv[0];
-    v->uv[1] = olduv[1];
-
-    return P_FALSE;
-  }
-
-  return P_TRUE;
-}
-
-static void p_vert_harmonic_insert(PVert *v)
-{
-  PEdge *e;
-
-  if (!p_vert_map_harmonic_weights(v)) {
-    /* do polygon kernel center insertion: this is quite slow, but should
-     * only be needed for 0.01 % of verts or so, when insert with harmonic
-     * weights fails */
-
-    int npoints = 0, i;
-    float(*points)[2];
-
-    e = v->edge;
-    do {
-      npoints++;
-      e = p_wheel_edge_next(e);
-    } while (e && (e != v->edge));
-
-    if (e == NULL) {
-      npoints++;
-    }
-
-    points = MEM_mallocN(sizeof(float[2]) * npoints, "PHarmonicPoints");
-
-    e = v->edge;
-    i = 0;
-    do {
-      PEdge *nexte = p_wheel_edge_next(e);
-
-      points[i][0] = e->next->vert->uv[0];
-      points[i][1] = e->next->vert->uv[1];
-
-      if (nexte == NULL) {
-        i++;
-        points[i][0] = e->next->next->vert->uv[0];
-        points[i][1] = e->next->next->vert->uv[1];
-        break;
-      }
-
-      e = nexte;
-      i++;
-    } while (e != v->edge);
-
-    p_polygon_kernel_center(points, npoints, v->uv);
-
-    MEM_freeN(points);
-  }
-
-  e = v->edge;
-  do {
-    if (!(e->next->vert->flag & PVERT_PIN)) {
-      p_vert_map_harmonic_weights(e->next->vert);
-    }
-    e = p_wheel_edge_next(e);
-  } while (e && (e != v->edge));
-
-  p_vert_map_harmonic_weights(v);
-}
-
-static void p_vert_fix_edge_pointer(PVert *v)
-{
-  PEdge *start = v->edge;
-
-  /* set v->edge pointer to the edge with no pair, if there is one */
-  while (v->edge->pair) {
-    v->edge = p_wheel_edge_prev(v->edge);
-
-    if (v->edge == start) {
-      break;
-    }
-  }
-}
-
-static void p_collapsing_verts(PEdge *edge, PEdge *pair, PVert **r_newv, PVert **r_keepv)
-{
-  /* the two vertices that are involved in the collapse */
-  if (edge) {
-    *r_newv = edge->vert;
-    *r_keepv = edge->next->vert;
-  }
-  else {
-    *r_newv = pair->next->vert;
-    *r_keepv = pair->vert;
-  }
-}
-
-static void p_collapse_edge(PEdge *edge, PEdge *pair)
-{
-  PVert *oldv, *keepv;
-  PEdge *e;
-
-  p_collapsing_verts(edge, pair, &oldv, &keepv);
-
-  /* change e->vert pointers from old vertex to the target vertex */
-  e = oldv->edge;
-  do {
-    if ((e != edge) && !(pair && pair->next == e)) {
-      e->vert = keepv;
-    }
-
-    e = p_wheel_edge_next(e);
-  } while (e && (e != oldv->edge));
-
-  /* set keepv->edge pointer */
-  if ((edge && (keepv->edge == edge->next)) || (keepv->edge == pair)) {
-    if (edge && edge->next->pair) {
-      keepv->edge = edge->next->pair->next;
-    }
-    else if (pair && pair->next->next->pair) {
-      keepv->edge = pair->next->next->pair;
-    }
-    else if (edge && edge->next->next->pair) {
-      keepv->edge = edge->next->next->pair;
-    }
-    else {
-      keepv->edge = pair->next->pair->next;
-    }
-  }
-
-  /* update pairs and v->edge pointers */
-  if (edge) {
-    PEdge *e1 = edge->next, *e2 = e1->next;
-
-    if (e1->pair) {
-      e1->pair->pair = e2->pair;
-    }
-
-    if (e2->pair) {
-      e2->pair->pair = e1->pair;
-      e2->vert->edge = p_wheel_edge_prev(e2);
-    }
-    else {
-      e2->vert->edge = p_wheel_edge_next(e2);
-    }
-
-    p_vert_fix_edge_pointer(e2->vert);
-  }
-
-  if (pair) {
-    PEdge *e1 = pair->next, *e2 = e1->next;
-
-    if (e1->pair) {
-      e1->pair->pair = e2->pair;
-    }
-
-    if (e2->pair) {
-      e2->pair->pair = e1->pair;
-      e2->vert->edge = p_wheel_edge_prev(e2);
-    }
-    else {
-      e2->vert->edge = p_wheel_edge_next(e2);
-    }
-
-    p_vert_fix_edge_pointer(e2->vert);
-  }
-
-  p_vert_fix_edge_pointer(keepv);
-
-  /* mark for move to collapsed list later */
-  oldv->flag |= PVERT_COLLAPSE;
-
-  if (edge) {
-    PFace *f = edge->face;
-    PEdge *e1 = edge->next, *e2 = e1->next;
-
-    f->flag |= PFACE_COLLAPSE;
-    edge->flag |= PEDGE_COLLAPSE;
-    e1->flag |= PEDGE_COLLAPSE;
-    e2->flag |= PEDGE_COLLAPSE;
-  }
-
-  if (pair) {
-    PFace *f = pair->face;
-    PEdge *e1 = pair->next, *e2 = e1->next;
-
-    f->flag |= PFACE_COLLAPSE;
-    pair->flag |= PEDGE_COLLAPSE;
-    e1->flag |= PEDGE_COLLAPSE;
-    e2->flag |= PEDGE_COLLAPSE;
-  }
-}
-
-static void p_split_vertex(PEdge *edge, PEdge *pair)
-{
-  PVert *newv, *keepv;
-  PEdge *e;
-
-  p_collapsing_verts(edge, pair, &newv, &keepv);
-
-  /* update edge pairs */
-  if (edge) {
-    PEdge *e1 = edge->next, *e2 = e1->next;
-
-    if (e1->pair) {
-      e1->pair->pair = e1;
-    }
-    if (e2->pair) {
-      e2->pair->pair = e2;
-    }
-
-    e2->vert->edge = e2;
-    p_vert_fix_edge_pointer(e2->vert);
-    keepv->edge = e1;
-  }
-
-  if (pair) {
-    PEdge *e1 = pair->next, *e2 = e1->next;
-
-    if (e1->pair) {
-      e1->pair->pair = e1;
-    }
-    if (e2->pair) {
-      e2->pair->pair = e2;
-    }
-
-    e2->vert->edge = e2;
-    p_vert_fix_edge_pointer(e2->vert);
-    keepv->edge = pair;
-  }
-
-  p_vert_fix_edge_pointer(keepv);
-
-  /* set e->vert pointers to restored vertex */
-  e = newv->edge;
-  do {
-    e->vert = newv;
-    e = p_wheel_edge_next(e);
-  } while (e && (e != newv->edge));
-}
-
-static PBool p_collapse_allowed_topologic(PEdge *edge, PEdge *pair)
-{
-  PVert *oldv, *keepv;
-
-  p_collapsing_verts(edge, pair, &oldv, &keepv);
-
-  /* boundary edges */
-  if (!edge || !pair) {
-    /* avoid collapsing chart into an edge */
-    if (edge && !edge->next->pair && !edge->next->next->pair) {
-      return P_FALSE;
-    }
-    else if (pair && !pair->next->pair && !pair->next->next->pair) {
-      return P_FALSE;
-    }
-  }
-  /* avoid merging two boundaries (oldv and keepv are on the 'other side' of
-   * the chart) */
-  else if (!p_vert_interior(oldv) && !p_vert_interior(keepv)) {
-    return P_FALSE;
-  }
-
-  return P_TRUE;
-}
-
-static PBool p_collapse_normal_flipped(float *v1, float *v2, float *vold, float *vnew)
-{
-  float nold[3], nnew[3], sub1[3], sub2[3];
-
-  sub_v3_v3v3(sub1, vold, v1);
-  sub_v3_v3v3(sub2, vold, v2);
-  cross_v3_v3v3(nold, sub1, sub2);
-
-  sub_v3_v3v3(sub1, vnew, v1);
-  sub_v3_v3v3(sub2, vnew, v2);
-  cross_v3_v3v3(nnew, sub1, sub2);
-
-  return (dot_v3v3(nold, nnew) <= 0.0f);
-}
-
-static PBool p_collapse_allowed_geometric(PEdge *edge, PEdge *pair)
-{
-  PVert *oldv, *keepv;
-  PEdge *e;
-  float angulardefect, angle;
-
-  p_collapsing_verts(edge, pair, &oldv, &keepv);
-
-  angulardefect = 2 * M_PI;
-
-  e = oldv->edge;
-  do {
-    float a[3], b[3], minangle, maxangle;
-    PEdge *e1 = e->next, *e2 = e1->next;
-    PVert *v1 = e1->vert, *v2 = e2->vert;
-    int i;
-
-    angle = p_vec_angle(v1->co, oldv->co, v2->co);
-    angulardefect -= angle;
-
-    /* skip collapsing faces */
-    if (v1 == keepv || v2 == keepv) {
-      e = p_wheel_edge_next(e);
-      continue;
-    }
-
-    if (p_collapse_normal_flipped(v1->co, v2->co, oldv->co, keepv->co)) {
-      return P_FALSE;
-    }
-
-    a[0] = angle;
-    a[1] = p_vec_angle(v2->co, v1->co, oldv->co);
-    a[2] = M_PI - a[0] - a[1];
-
-    b[0] = p_vec_angle(v1->co, keepv->co, v2->co);
-    b[1] = p_vec_angle(v2->co, v1->co, keepv->co);
-    b[2] = M_PI - b[0] - b[1];
-
-    /* ABF criterion 1: avoid sharp and obtuse angles. */
-    minangle = 15.0f * M_PI / 180.0f;
-    maxangle = M_PI - minangle;
-
-    for (i = 0; i < 3; i++) {
-      if ((b[i] < a[i]) && (b[i] < minangle)) {
-        return P_FALSE;
-      }
-      else if ((b[i] > a[i]) && (b[i] > maxangle)) {
-        return P_FALSE;
-      }
-    }
-
-    e = p_wheel_edge_next(e);
-  } while (e && (e != oldv->edge));
-
-  if (p_vert_interior(oldv)) {
-    /* HLSCM criterion: angular defect smaller than threshold. */
-    if (fabsf(angulardefect) > (float)(M_PI * 30.0 / 180.0)) {
-      return P_FALSE;
-    }
-  }
-  else {
-    PVert *v1 = p_boundary_edge_next(oldv->edge)->vert;
-    PVert *v2 = p_boundary_edge_prev(oldv->edge)->vert;
-
-    /* ABF++ criterion 2: avoid collapsing verts inwards. */
-    if (p_vert_interior(keepv)) {
-      return P_FALSE;
-    }
-
-    /* Don't collapse significant boundary changes. */
-    angle = p_vec_angle(v1->co, oldv->co, v2->co);
-    if (angle < (M_PI * 160.0 / 180.0)) {
-      return P_FALSE;
-    }
-  }
-
-  return P_TRUE;
-}
-
-static PBool p_collapse_allowed(PEdge *edge, PEdge *pair)
-{
-  PVert *oldv, *keepv;
-
-  p_collapsing_verts(edge, pair, &oldv, &keepv);
-
-  if (oldv->flag & PVERT_PIN) {
-    return P_FALSE;
-  }
-
-  return (p_collapse_allowed_topologic(edge, pair) && p_collapse_allowed_geometric(edge, pair));
-}
-
-static float p_collapse_cost(PEdge *edge, PEdge *pair)
-{
-  /* based on volume and boundary optimization from:
-   * "Fast and Memory Efficient Polygonal Simplification" P. Lindstrom, G. Turk */
-
-  PVert *oldv, *keepv;
-  PEdge *e;
-  PFace *oldf1, *oldf2;
-  float volumecost = 0.0f, areacost = 0.0f, edgevec[3], cost, weight, elen;
-  float shapecost = 0.0f;
-  float shapeold = 0.0f, shapenew = 0.0f;
-  int nshapeold = 0, nshapenew = 0;
-
-  p_collapsing_verts(edge, pair, &oldv, &keepv);
-  oldf1 = (edge) ? edge->face : NULL;
-  oldf2 = (pair) ? pair->face : NULL;
-
-  sub_v3_v3v3(edgevec, keepv->co, oldv->co);
-
-  e = oldv->edge;
-  do {
-    float a1, a2, a3;
-    float *co1 = e->next->vert->co;
-    float *co2 = e->next->next->vert->co;
-
-    if (!ELEM(e->face, oldf1, oldf2)) {
-      float tetrav2[3], tetrav3[3];
-
-      /* tetrahedron volume = (1/3!)*|a.(b x c)| */
-      sub_v3_v3v3(tetrav2, co1, oldv->co);
-      sub_v3_v3v3(tetrav3, co2, oldv->co);
-      volumecost += fabsf(volume_tri_tetrahedron_signed_v3(tetrav2, tetrav3, edgevec));
-
-#  if 0
-      shapecost += dot_v3v3(co1, keepv->co);
-
-      if (p_wheel_edge_next(e) == NULL) {
-        shapecost += dot_v3v3(co2, keepv->co);
-      }
-#  endif
-
-      p_triangle_angles(oldv->co, co1, co2, &a1, &a2, &a3);
-      a1 = a1 - M_PI / 3.0;
-      a2 = a2 - M_PI / 3.0;
-      a3 = a3 - M_PI / 3.0;
-      shapeold = (a1 * a1 + a2 * a2 + a3 * a3) / (M_PI_2 * M_PI_2);
-
-      nshapeold++;
-    }
-    else {
-      p_triangle_angles(keepv->co, co1, co2, &a1, &a2, &a3);
-      a1 = a1 - M_PI / 3.0;
-      a2 = a2 - M_PI / 3.0;
-      a3 = a3 - M_PI / 3.0;
-      shapenew = (a1 * a1 + a2 * a2 + a3 * a3) / (M_PI_2 * M_PI_2);
-
-      nshapenew++;
-    }
-
-    e = p_wheel_edge_next(e);
-  } while (e && (e != oldv->edge));
-
-  if (!p_vert_interior(oldv)) {
-    PVert *v1 = p_boundary_edge_prev(oldv->edge)->vert;
-    PVert *v2 = p_boundary_edge_next(oldv->edge)->vert;
-
-    areacost = area_tri_v3(oldv->co, v1->co, v2->co);
-  }
-
-  elen = len_v3(edgevec);
-  weight = 1.0f; /* 0.2f */
-  cost = weight * volumecost * volumecost + elen * elen * areacost * areacost;
-#  if 0
-  cost += shapecost;
-#  else
-  shapeold /= nshapeold;
-  shapenew /= nshapenew;
-  shapecost = (shapeold + 0.00001) / (shapenew + 0.00001);
-
-  cost *= shapecost;
-#  endif
-
-  return cost;
-}
-
-static void p_collapse_cost_vertex(PVert *vert, float *r_mincost, PEdge **r_mine)
-{
-  PEdge *e, *enext, *pair;
-
-  *r_mine = NULL;
-  *r_mincost = 0.0f;
-  e = vert->edge;
-  do {
-    if (p_collapse_allowed(e, e->pair)) {
-      float cost = p_collapse_cost(e, e->pair);
-
-      if ((*r_mine == NULL) || (cost < *r_mincost)) {
-        *r_mincost = cost;
-        *r_mine = e;
-      }
-    }
-
-    enext = p_wheel_edge_next(e);
-
-    if (enext == NULL) {
-      /* the other boundary edge, where we only have the pair halfedge */
-      pair = e->next->next;
-
-      if (p_collapse_allowed(NULL, pair)) {
-        float cost = p_collapse_cost(NULL, pair);
-
-        if ((*r_mine == NULL) || (cost < *r_mincost)) {
-          *r_mincost = cost;
-          *r_mine = pair;
-        }
-      }
-
-      break;
-    }
-
-    e = enext;
-  } while (e != vert->edge);
-}
-
-static void p_chart_post_collapse_flush(PChart *chart, PEdge *collapsed)
-{
-  /* Move to `collapsed_*`. */
-
-  PVert *v, *nextv = NULL, *verts = chart->verts;
-  PEdge *e, *nexte = NULL, *edges = chart->edges, *laste = NULL;
-  PFace *f, *nextf = NULL, *faces = chart->faces;
-
-  chart->verts = chart->collapsed_verts = NULL;
-  chart->edges = chart->collapsed_edges = NULL;
-  chart->faces = chart->collapsed_faces = NULL;
-
-  chart->nverts = chart->nedges = chart->nfaces = 0;
-
-  for (v = verts; v; v = nextv) {
-    nextv = v->nextlink;
-
-    if (v->flag & PVERT_COLLAPSE) {
-      v->nextlink = chart->collapsed_verts;
-      chart->collapsed_verts = v;
-    }
-    else {
-      v->nextlink = chart->verts;
-      chart->verts = v;
-      chart->nverts++;
-    }
-  }
-
-  for (e = edges; e; e = nexte) {
-    nexte = e->nextlink;
-
-    if (!collapsed || !(e->flag & PEDGE_COLLAPSE_EDGE)) {
-      if (e->flag & PEDGE_COLLAPSE) {
-        e->nextlink = chart->collapsed_edges;
-        chart->collapsed_edges = e;
-      }
-      else {
-        e->nextlink = chart->edges;
-        chart->edges = e;
-        chart->nedges++;
-      }
-    }
-  }
-
-  /* these are added last so they can be popped of in the right order
-   * for splitting */
-  for (e = collapsed; e; e = e->nextlink) {
-    e->nextlink = e->u.nextcollapse;
-    laste = e;
-  }
-  if (laste) {
-    laste->nextlink = chart->collapsed_edges;
-    chart->collapsed_edges = collapsed;
-  }
-
-  for (f = faces; f; f = nextf) {
-    nextf = f->nextlink;
-
-    if (f->flag & PFACE_COLLAPSE) {
-      f->nextlink = chart->collapsed_faces;
-      chart->collapsed_faces = f;
-    }
-    else {
-      f->nextlink = chart->faces;
-      chart->faces = f;
-      chart->nfaces++;
-    }
-  }
-}
-
-static void p_chart_post_split_flush(PChart *chart)
-{
-  /* Move from `collapsed_*`. */
-
-  PVert *v, *nextv = NULL;
-  PEdge *e, *nexte = NULL;
-  PFace *f, *nextf = NULL;
-
-  for (v = chart->collapsed_verts; v; v = nextv) {
-    nextv = v->nextlink;
-    v->nextlink = chart->verts;
-    chart->verts = v;
-    chart->nverts++;
-  }
-
-  for (e = chart->collapsed_edges; e; e = nexte) {
-    nexte = e->nextlink;
-    e->nextlink = chart->edges;
-    chart->edges = e;
-    chart->nedges++;
-  }
-
-  for (f = chart->collapsed_faces; f; f = nextf) {
-    nextf = f->nextlink;
-    f->nextlink = chart->faces;
-    chart->faces = f;
-    chart->nfaces++;
-  }
-
-  chart->collapsed_verts = NULL;
-  chart->collapsed_edges = NULL;
-  chart->collapsed_faces = NULL;
-}
-
-static void p_chart_simplify_compute(PChart *chart)
-{
-  /* Computes a list of edge collapses / vertex splits. The collapsed
-   * simplices go in the `chart->collapsed_*` lists, The original and
-   * collapsed may then be view as stacks, where the next collapse/split
-   * is at the top of the respective lists. */
-
-  Heap *heap = BLI_heap_new();
-  PVert *v, **wheelverts;
-  PEdge *collapsededges = NULL, *e;
-  int nwheelverts, i, ncollapsed = 0;
-
-  wheelverts = MEM_mallocN(sizeof(PVert *) * chart->nverts, "PChartWheelVerts");
-
-  /* insert all potential collapses into heap */
-  for (v = chart->verts; v; v = v->nextlink) {
-    float cost;
-    PEdge *e = NULL;
-
-    p_collapse_cost_vertex(v, &cost, &e);
-
-    if (e) {
-      v->u.heaplink = BLI_heap_insert(heap, cost, e);
-    }
-    else {
-      v->u.heaplink = NULL;
-    }
-  }
-
-  for (e = chart->edges; e; e = e->nextlink) {
-    e->u.nextcollapse = NULL;
-  }
-
-  /* pop edge collapse out of heap one by one */
-  while (!BLI_heap_is_empty(heap)) {
-    if (ncollapsed == NCOLLAPSE) {
-      break;
-    }
-
-    HeapNode *link = BLI_heap_top(heap);
-    PEdge *edge = (PEdge *)BLI_heap_pop_min(heap), *pair = edge->pair;
-    PVert *oldv, *keepv;
-    PEdge *wheele, *nexte;
-
-    /* remember the edges we collapsed */
-    edge->u.nextcollapse = collapsededges;
-    collapsededges = edge;
-
-    if (edge->vert->u.heaplink != link) {
-      edge->flag |= (PEDGE_COLLAPSE_EDGE | PEDGE_COLLAPSE_PAIR);
-      edge->next->vert->u.heaplink = NULL;
-      SWAP(PEdge *, edge, pair);
-    }
-    else {
-      edge->flag |= PEDGE_COLLAPSE_EDGE;
-      edge->vert->u.heaplink = NULL;
-    }
-
-    p_collapsing_verts(edge, pair, &oldv, &keepv);
-
-    /* gather all wheel verts and remember them before collapse */
-    nwheelverts = 0;
-    wheele = oldv->edge;
-
-    do {
-      wheelverts[nwheelverts++] = wheele->next->vert;
-      nexte = p_wheel_edge_next(wheele);
-
-      if (nexte == NULL) {
-        wheelverts[nwheelverts++] = wheele->next->next->vert;
-      }
-
-      wheele = nexte;
-    } while (wheele && (wheele != oldv->edge));
-
-    /* collapse */
-    p_collapse_edge(edge, pair);
-
-    for (i = 0; i < nwheelverts; i++) {
-      float cost;
-      PEdge *collapse = NULL;
-
-      v = wheelverts[i];
-
-      if (v->u.heaplink) {
-        BLI_heap_remove(heap, v->u.heaplink);
-        v->u.heaplink = NULL;
-      }
-
-      p_collapse_cost_vertex(v, &cost, &collapse);
-
-      if (collapse) {
-        v->u.heaplink = BLI_heap_insert(heap, cost, collapse);
-      }
-    }
-
-    ncollapsed++;
-  }
-
-  MEM_freeN(wheelverts);
-  BLI_heap_free(heap, NULL);
-
-  p_chart_post_collapse_flush(chart, collapsededges);
-}
-
-static void p_chart_complexify(PChart *chart)
-{
-  PEdge *e, *pair, *edge;
-  PVert *newv, *keepv;
-  int x = 0;
-
-  for (e = chart->collapsed_edges; e; e = e->nextlink) {
-    if (!(e->flag & PEDGE_COLLAPSE_EDGE)) {
-      break;
-    }
-
-    edge = e;
-    pair = e->pair;
-
-    if (edge->flag & PEDGE_COLLAPSE_PAIR) {
-      SWAP(PEdge *, edge, pair);
-    }
-
-    p_split_vertex(edge, pair);
-    p_collapsing_verts(edge, pair, &newv, &keepv);
-
-    if (x >= NCOLLAPSEX) {
-      newv->uv[0] = keepv->uv[0];
-      newv->uv[1] = keepv->uv[1];
-    }
-    else {
-      p_vert_harmonic_insert(newv);
-      x++;
-    }
-  }
-
-  p_chart_post_split_flush(chart);
-}
-
-#  if 0
-static void p_chart_simplify(PChart *chart)
-{
-  /* Not implemented, needs proper reordering in split_flush. */
-}
-#  endif
-#endif
-
-/* ABF */
-
-#define ABF_MAX_ITER 20
-
-typedef struct PAbfSystem {
-  int ninterior, nfaces, nangles;
-  float *alpha, *beta, *sine, *cosine, *weight;
-  float *bAlpha, *bTriangle, *bInterior;
-  float *lambdaTriangle, *lambdaPlanar, *lambdaLength;
-  float (*J2dt)[3], *bstar, *dstar;
-  float minangle, maxangle;
-} PAbfSystem;
-
-static void p_abf_setup_system(PAbfSystem *sys)
-{
-  int i;
-
-  sys->alpha = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFalpha");
-  sys->beta = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFbeta");
-  sys->sine = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFsine");
-  sys->cosine = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFcosine");
-  sys->weight = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFweight");
-
-  sys->bAlpha = (float *)MEM_mallocN(sizeof(float) * sys->nangles, "ABFbalpha");
-  sys->bTriangle = (float *)MEM_mallocN(sizeof(float) * sys->nfaces, "ABFbtriangle");
-  sys->bInterior = (float *)MEM_mallocN(sizeof(float[2]) * sys->ninterior, "ABFbinterior");
-
-  sys->lambdaTriangle = (float *)MEM_callocN(sizeof(float) * sys->nfaces, "ABFlambdatri");
-  sys->lambdaPlanar = (float *)MEM_callocN(sizeof(float) * sys->ninterior, "ABFlamdaplane");
-  sys->lambdaLength = (float *)MEM_mallocN(sizeof(float) * sys->ninterior, "ABFlambdalen");
-
-  sys->J2dt = MEM_mallocN(sizeof(float) * sys->nangles * 3, "ABFj2dt");
-  sys->bstar = (float *)MEM_mallocN(sizeof(float) * sys->nfaces, "ABFbstar");
-  sys->dstar = (float *)MEM_mallocN(sizeof(float) * sys->nfaces, "ABFdstar");
-
-  for (i = 0; i < sys->ninterior; i++) {
-    sys->lambdaLength[i] = 1.0;
-  }
-
-  sys->minangle = 1.0 * M_PI / 180.0;
-  sys->maxangle = (float)M_PI - sys->minangle;
-}
-
-static void p_abf_free_system(PAbfSystem *sys)
-{
-  MEM_freeN(sys->alpha);
-  MEM_freeN(sys->beta);
-  MEM_freeN(sys->sine);
-  MEM_freeN(sys->cosine);
-  MEM_freeN(sys->weight);
-  MEM_freeN(sys->bAlpha);
-  MEM_freeN(sys->bTriangle);
-  MEM_freeN(sys->bInterior);
-  MEM_freeN(sys->lambdaTriangle);
-  MEM_freeN(sys->lambdaPlanar);
-  MEM_freeN(sys->lambdaLength);
-  MEM_freeN(sys->J2dt);
-  MEM_freeN(sys->bstar);
-  MEM_freeN(sys->dstar);
-}
-
-static void p_abf_compute_sines(PAbfSystem *sys)
-{
-  int i;
-  float *sine = sys->sine, *cosine = sys->cosine, *alpha = sys->alpha;
-
-  for (i = 0; i < sys->nangles; i++, sine++, cosine++, alpha++) {
-    *sine = sinf(*alpha);
-    *cosine = cosf(*alpha);
-  }
-}
-
-static float p_abf_compute_sin_product(PAbfSystem *sys, PVert *v, int aid)
-{
-  PEdge *e, *e1, *e2;
-  float sin1, sin2;
-
-  sin1 = sin2 = 1.0;
-
-  e = v->edge;
-  do {
-    e1 = e->next;
-    e2 = e->next->next;
-
-    if (aid == e1->u.id) {
-      /* we are computing a derivative for this angle,
-       * so we use cos and drop the other part */
-      sin1 *= sys->cosine[e1->u.id];
-      sin2 = 0.0;
-    }
-    else {
-      sin1 *= sys->sine[e1->u.id];
-    }
-
-    if (aid == e2->u.id) {
-      /* see above */
-      sin1 = 0.0;
-      sin2 *= sys->cosine[e2->u.id];
-    }
-    else {
-      sin2 *= sys->sine[e2->u.id];
-    }
-
-    e = e->next->next->pair;
-  } while (e && (e != v->edge));
-
-  return (sin1 - sin2);
-}
-
-static float p_abf_compute_grad_alpha(PAbfSystem *sys, PFace *f, PEdge *e)
-{
-  PVert *v = e->vert, *v1 = e->next->vert, *v2 = e->next->next->vert;
-  float deriv;
-
-  deriv = (sys->alpha[e->u.id] - sys->beta[e->u.id]) * sys->weight[e->u.id];
-  deriv += sys->lambdaTriangle[f->u.id];
-
-  if (v->flag & PVERT_INTERIOR) {
-    deriv += sys->lambdaPlanar[v->u.id];
-  }
-
-  if (v1->flag & PVERT_INTERIOR) {
-    float product = p_abf_compute_sin_product(sys, v1, e->u.id);
-    deriv += sys->lambdaLength[v1->u.id] * product;
-  }
-
-  if (v2->flag & PVERT_INTERIOR) {
-    float product = p_abf_compute_sin_product(sys, v2, e->u.id);
-    deriv += sys->lambdaLength[v2->u.id] * product;
-  }
-
-  return deriv;
-}
-
-static float p_abf_compute_gradient(PAbfSystem *sys, PChart *chart)
-{
-  PFace *f;
-  PEdge *e;
-  PVert *v;
-  float norm = 0.0;
-
-  for (f = chart->faces; f; f = f->nextlink) {
-    PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-    float gtriangle, galpha1, galpha2, galpha3;
-
-    galpha1 = p_abf_compute_grad_alpha(sys, f, e1);
-    galpha2 = p_abf_compute_grad_alpha(sys, f, e2);
-    galpha3 = p_abf_compute_grad_alpha(sys, f, e3);
-
-    sys->bAlpha[e1->u.id] = -galpha1;
-    sys->bAlpha[e2->u.id] = -galpha2;
-    sys->bAlpha[e3->u.id] = -galpha3;
-
-    norm += galpha1 * galpha1 + galpha2 * galpha2 + galpha3 * galpha3;
-
-    gtriangle = sys->alpha[e1->u.id] + sys->alpha[e2->u.id] + sys->alpha[e3->u.id] - (float)M_PI;
-    sys->bTriangle[f->u.id] = -gtriangle;
-    norm += gtriangle * gtriangle;
-  }
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    if (v->flag & PVERT_INTERIOR) {
-      float gplanar = -2 * M_PI, glength;
-
-      e = v->edge;
-      do {
-        gplanar += sys->alpha[e->u.id];
-        e = e->next->next->pair;
-      } while (e && (e != v->edge));
-
-      sys->bInterior[v->u.id] = -gplanar;
-      norm += gplanar * gplanar;
-
-      glength = p_abf_compute_sin_product(sys, v, -1);
-      sys->bInterior[sys->ninterior + v->u.id] = -glength;
-      norm += glength * glength;
-    }
-  }
-
-  return norm;
-}
-
-static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
-{
-  PFace *f;
-  PEdge *e;
-  int i, j, ninterior = sys->ninterior, nvar = 2 * sys->ninterior;
-  PBool success;
-  LinearSolver *context;
-
-  context = EIG_linear_solver_new(0, nvar, 1);
-
-  for (i = 0; i < nvar; i++) {
-    EIG_linear_solver_right_hand_side_add(context, 0, i, sys->bInterior[i]);
-  }
-
-  for (f = chart->faces; f; f = f->nextlink) {
-    float wi1, wi2, wi3, b, si, beta[3], j2[3][3], W[3][3];
-    float row1[6], row2[6], row3[6];
-    int vid[6];
-    PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-    PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
-
-    wi1 = 1.0f / sys->weight[e1->u.id];
-    wi2 = 1.0f / sys->weight[e2->u.id];
-    wi3 = 1.0f / sys->weight[e3->u.id];
-
-    /* bstar1 = (J1*dInv*bAlpha - bTriangle) */
-    b = sys->bAlpha[e1->u.id] * wi1;
-    b += sys->bAlpha[e2->u.id] * wi2;
-    b += sys->bAlpha[e3->u.id] * wi3;
-    b -= sys->bTriangle[f->u.id];
-
-    /* si = J1*d*J1t */
-    si = 1.0f / (wi1 + wi2 + wi3);
-
-    /* J1t*si*bstar1 - bAlpha */
-    beta[0] = b * si - sys->bAlpha[e1->u.id];
-    beta[1] = b * si - sys->bAlpha[e2->u.id];
-    beta[2] = b * si - sys->bAlpha[e3->u.id];
-
-    /* use this later for computing other lambda's */
-    sys->bstar[f->u.id] = b;
-    sys->dstar[f->u.id] = si;
-
-    /* set matrix */
-    W[0][0] = si - sys->weight[e1->u.id];
-    W[0][1] = si;
-    W[0][2] = si;
-    W[1][0] = si;
-    W[1][1] = si - sys->weight[e2->u.id];
-    W[1][2] = si;
-    W[2][0] = si;
-    W[2][1] = si;
-    W[2][2] = si - sys->weight[e3->u.id];
-
-    vid[0] = vid[1] = vid[2] = vid[3] = vid[4] = vid[5] = -1;
-
-    if (v1->flag & PVERT_INTERIOR) {
-      vid[0] = v1->u.id;
-      vid[3] = ninterior + v1->u.id;
-
-      sys->J2dt[e1->u.id][0] = j2[0][0] = 1.0f * wi1;
-      sys->J2dt[e2->u.id][0] = j2[1][0] = p_abf_compute_sin_product(sys, v1, e2->u.id) * wi2;
-      sys->J2dt[e3->u.id][0] = j2[2][0] = p_abf_compute_sin_product(sys, v1, e3->u.id) * wi3;
-
-      EIG_linear_solver_right_hand_side_add(context, 0, v1->u.id, j2[0][0] * beta[0]);
-      EIG_linear_solver_right_hand_side_add(
-          context, 0, ninterior + v1->u.id, j2[1][0] * beta[1] + j2[2][0] * beta[2]);
-
-      row1[0] = j2[0][0] * W[0][0];
-      row2[0] = j2[0][0] * W[1][0];
-      row3[0] = j2[0][0] * W[2][0];
-
-      row1[3] = j2[1][0] * W[0][1] + j2[2][0] * W[0][2];
-      row2[3] = j2[1][0] * W[1][1] + j2[2][0] * W[1][2];
-      row3[3] = j2[1][0] * W[2][1] + j2[2][0] * W[2][2];
-    }
-
-    if (v2->flag & PVERT_INTERIOR) {
-      vid[1] = v2->u.id;
-      vid[4] = ninterior + v2->u.id;
-
-      sys->J2dt[e1->u.id][1] = j2[0][1] = p_abf_compute_sin_product(sys, v2, e1->u.id) * wi1;
-      sys->J2dt[e2->u.id][1] = j2[1][1] = 1.0f * wi2;
-      sys->J2dt[e3->u.id][1] = j2[2][1] = p_abf_compute_sin_product(sys, v2, e3->u.id) * wi3;
-
-      EIG_linear_solver_right_hand_side_add(context, 0, v2->u.id, j2[1][1] * beta[1]);
-      EIG_linear_solver_right_hand_side_add(
-          context, 0, ninterior + v2->u.id, j2[0][1] * beta[0] + j2[2][1] * beta[2]);
-
-      row1[1] = j2[1][1] * W[0][1];
-      row2[1] = j2[1][1] * W[1][1];
-      row3[1] = j2[1][1] * W[2][1];
-
-      row1[4] = j2[0][1] * W[0][0] + j2[2][1] * W[0][2];
-      row2[4] = j2[0][1] * W[1][0] + j2[2][1] * W[1][2];
-      row3[4] = j2[0][1] * W[2][0] + j2[2][1] * W[2][2];
-    }
-
-    if (v3->flag & PVERT_INTERIOR) {
-      vid[2] = v3->u.id;
-      vid[5] = ninterior + v3->u.id;
-
-      sys->J2dt[e1->u.id][2] = j2[0][2] = p_abf_compute_sin_product(sys, v3, e1->u.id) * wi1;
-      sys->J2dt[e2->u.id][2] = j2[1][2] = p_abf_compute_sin_product(sys, v3, e2->u.id) * wi2;
-      sys->J2dt[e3->u.id][2] = j2[2][2] = 1.0f * wi3;
-
-      EIG_linear_solver_right_hand_side_add(context, 0, v3->u.id, j2[2][2] * beta[2]);
-      EIG_linear_solver_right_hand_side_add(
-          context, 0, ninterior + v3->u.id, j2[0][2] * beta[0] + j2[1][2] * beta[1]);
-
-      row1[2] = j2[2][2] * W[0][2];
-      row2[2] = j2[2][2] * W[1][2];
-      row3[2] = j2[2][2] * W[2][2];
-
-      row1[5] = j2[0][2] * W[0][0] + j2[1][2] * W[0][1];
-      row2[5] = j2[0][2] * W[1][0] + j2[1][2] * W[1][1];
-      row3[5] = j2[0][2] * W[2][0] + j2[1][2] * W[2][1];
-    }
-
-    for (i = 0; i < 3; i++) {
-      int r = vid[i];
-
-      if (r == -1) {
-        continue;
-      }
-
-      for (j = 0; j < 6; j++) {
-        int c = vid[j];
-
-        if (c == -1) {
-          continue;
-        }
-
-        if (i == 0) {
-          EIG_linear_solver_matrix_add(context, r, c, j2[0][i] * row1[j]);
-        }
-        else {
-          EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[0][i] * row1[j]);
-        }
-
-        if (i == 1) {
-          EIG_linear_solver_matrix_add(context, r, c, j2[1][i] * row2[j]);
-        }
-        else {
-          EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[1][i] * row2[j]);
-        }
-
-        if (i == 2) {
-          EIG_linear_solver_matrix_add(context, r, c, j2[2][i] * row3[j]);
-        }
-        else {
-          EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[2][i] * row3[j]);
-        }
-      }
-    }
-  }
-
-  success = EIG_linear_solver_solve(context);
-
-  if (success) {
-    for (f = chart->faces; f; f = f->nextlink) {
-      float dlambda1, pre[3], dalpha;
-      PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-      PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
-
-      pre[0] = pre[1] = pre[2] = 0.0;
-
-      if (v1->flag & PVERT_INTERIOR) {
-        float x = EIG_linear_solver_variable_get(context, 0, v1->u.id);
-        float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v1->u.id);
-        pre[0] += sys->J2dt[e1->u.id][0] * x;
-        pre[1] += sys->J2dt[e2->u.id][0] * x2;
-        pre[2] += sys->J2dt[e3->u.id][0] * x2;
-      }
-
-      if (v2->flag & PVERT_INTERIOR) {
-        float x = EIG_linear_solver_variable_get(context, 0, v2->u.id);
-        float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v2->u.id);
-        pre[0] += sys->J2dt[e1->u.id][1] * x2;
-        pre[1] += sys->J2dt[e2->u.id][1] * x;
-        pre[2] += sys->J2dt[e3->u.id][1] * x2;
-      }
-
-      if (v3->flag & PVERT_INTERIOR) {
-        float x = EIG_linear_solver_variable_get(context, 0, v3->u.id);
-        float x2 = EIG_linear_solver_variable_get(context, 0, ninterior + v3->u.id);
-        pre[0] += sys->J2dt[e1->u.id][2] * x2;
-        pre[1] += sys->J2dt[e2->u.id][2] * x2;
-        pre[2] += sys->J2dt[e3->u.id][2] * x;
-      }
-
-      dlambda1 = pre[0] + pre[1] + pre[2];
-      dlambda1 = sys->dstar[f->u.id] * (sys->bstar[f->u.id] - dlambda1);
-
-      sys->lambdaTriangle[f->u.id] += dlambda1;
-
-      dalpha = (sys->bAlpha[e1->u.id] - dlambda1);
-      sys->alpha[e1->u.id] += dalpha / sys->weight[e1->u.id] - pre[0];
-
-      dalpha = (sys->bAlpha[e2->u.id] - dlambda1);
-      sys->alpha[e2->u.id] += dalpha / sys->weight[e2->u.id] - pre[1];
-
-      dalpha = (sys->bAlpha[e3->u.id] - dlambda1);
-      sys->alpha[e3->u.id] += dalpha / sys->weight[e3->u.id] - pre[2];
-
-      /* clamp */
-      e = f->edge;
-      do {
-        if (sys->alpha[e->u.id] > (float)M_PI) {
-          sys->alpha[e->u.id] = (float)M_PI;
-        }
-        else if (sys->alpha[e->u.id] < 0.0f) {
-          sys->alpha[e->u.id] = 0.0f;
-        }
-      } while (e != f->edge);
-    }
-
-    for (i = 0; i < ninterior; i++) {
-      sys->lambdaPlanar[i] += (float)EIG_linear_solver_variable_get(context, 0, i);
-      sys->lambdaLength[i] += (float)EIG_linear_solver_variable_get(context, 0, ninterior + i);
-    }
-  }
-
-  EIG_linear_solver_delete(context);
-
-  return success;
-}
-
-static PBool p_chart_abf_solve(PChart *chart)
-{
-  PVert *v;
-  PFace *f;
-  PEdge *e, *e1, *e2, *e3;
-  PAbfSystem sys;
-  int i;
-  float /* lastnorm, */ /* UNUSED */ limit = (chart->nfaces > 100) ? 1.0f : 0.001f;
-
-  /* setup id's */
-  sys.ninterior = sys.nfaces = sys.nangles = 0;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    if (p_vert_interior(v)) {
-      v->flag |= PVERT_INTERIOR;
-      v->u.id = sys.ninterior++;
-    }
-    else {
-      v->flag &= ~PVERT_INTERIOR;
-    }
-  }
-
-  for (f = chart->faces; f; f = f->nextlink) {
-    e1 = f->edge;
-    e2 = e1->next;
-    e3 = e2->next;
-    f->u.id = sys.nfaces++;
-
-    /* angle id's are conveniently stored in half edges */
-    e1->u.id = sys.nangles++;
-    e2->u.id = sys.nangles++;
-    e3->u.id = sys.nangles++;
-  }
-
-  p_abf_setup_system(&sys);
-
-  /* compute initial angles */
-  for (f = chart->faces; f; f = f->nextlink) {
-    float a1, a2, a3;
-
-    e1 = f->edge;
-    e2 = e1->next;
-    e3 = e2->next;
-    p_face_angles(f, &a1, &a2, &a3);
-
-    if (a1 < sys.minangle) {
-      a1 = sys.minangle;
-    }
-    else if (a1 > sys.maxangle) {
-      a1 = sys.maxangle;
-    }
-    if (a2 < sys.minangle) {
-      a2 = sys.minangle;
-    }
-    else if (a2 > sys.maxangle) {
-      a2 = sys.maxangle;
-    }
-    if (a3 < sys.minangle) {
-      a3 = sys.minangle;
-    }
-    else if (a3 > sys.maxangle) {
-      a3 = sys.maxangle;
-    }
-
-    sys.alpha[e1->u.id] = sys.beta[e1->u.id] = a1;
-    sys.alpha[e2->u.id] = sys.beta[e2->u.id] = a2;
-    sys.alpha[e3->u.id] = sys.beta[e3->u.id] = a3;
-
-    sys.weight[e1->u.id] = 2.0f / (a1 * a1);
-    sys.weight[e2->u.id] = 2.0f / (a2 * a2);
-    sys.weight[e3->u.id] = 2.0f / (a3 * a3);
-  }
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    if (v->flag & PVERT_INTERIOR) {
-      float anglesum = 0.0, scale;
-
-      e = v->edge;
-      do {
-        anglesum += sys.beta[e->u.id];
-        e = e->next->next->pair;
-      } while (e && (e != v->edge));
-
-      scale = (anglesum == 0.0f) ? 0.0f : 2.0f * (float)M_PI / anglesum;
-
-      e = v->edge;
-      do {
-        sys.beta[e->u.id] = sys.alpha[e->u.id] = sys.beta[e->u.id] * scale;
-        e = e->next->next->pair;
-      } while (e && (e != v->edge));
-    }
-  }
-
-  if (sys.ninterior > 0) {
-    p_abf_compute_sines(&sys);
-
-    /* iteration */
-    /* lastnorm = 1e10; */ /* UNUSED */
-
-    for (i = 0; i < ABF_MAX_ITER; i++) {
-      float norm = p_abf_compute_gradient(&sys, chart);
-
-      /* lastnorm = norm; */ /* UNUSED */
-
-      if (norm < limit) {
-        break;
-      }
-
-      if (!p_abf_matrix_invert(&sys, chart)) {
-        param_warning("ABF failed to invert matrix");
-        p_abf_free_system(&sys);
-        return P_FALSE;
-      }
-
-      p_abf_compute_sines(&sys);
-    }
-
-    if (i == ABF_MAX_ITER) {
-      param_warning("ABF maximum iterations reached");
-      p_abf_free_system(&sys);
-      return P_FALSE;
-    }
-  }
-
-  chart->u.lscm.abf_alpha = MEM_dupallocN(sys.alpha);
-  p_abf_free_system(&sys);
-
-  return P_TRUE;
-}
-
-/* Least Squares Conformal Maps */
-
-static void p_chart_pin_positions(PChart *chart, PVert **pin1, PVert **pin2)
-{
-  if (!*pin1 || !*pin2 || *pin1 == *pin2) {
-    /* degenerate case */
-    PFace *f = chart->faces;
-    *pin1 = f->edge->vert;
-    *pin2 = f->edge->next->vert;
-
-    (*pin1)->uv[0] = 0.0f;
-    (*pin1)->uv[1] = 0.5f;
-    (*pin2)->uv[0] = 1.0f;
-    (*pin2)->uv[1] = 0.5f;
-  }
-  else {
-    int diru, dirv, dirx, diry;
-    float sub[3];
-
-    sub_v3_v3v3(sub, (*pin1)->co, (*pin2)->co);
-    sub[0] = fabsf(sub[0]);
-    sub[1] = fabsf(sub[1]);
-    sub[2] = fabsf(sub[2]);
-
-    if ((sub[0] > sub[1]) && (sub[0] > sub[2])) {
-      dirx = 0;
-      diry = (sub[1] > sub[2]) ? 1 : 2;
-    }
-    else if ((sub[1] > sub[0]) && (sub[1] > sub[2])) {
-      dirx = 1;
-      diry = (sub[0] > sub[2]) ? 0 : 2;
-    }
-    else {
-      dirx = 2;
-      diry = (sub[0] > sub[1]) ? 0 : 1;
-    }
-
-    if (dirx == 2) {
-      diru = 1;
-      dirv = 0;
-    }
-    else {
-      diru = 0;
-      dirv = 1;
-    }
-
-    (*pin1)->uv[diru] = (*pin1)->co[dirx];
-    (*pin1)->uv[dirv] = (*pin1)->co[diry];
-    (*pin2)->uv[diru] = (*pin2)->co[dirx];
-    (*pin2)->uv[dirv] = (*pin2)->co[diry];
-  }
-}
-
-static PBool p_chart_symmetry_pins(PChart *chart, PEdge *outer, PVert **pin1, PVert **pin2)
-{
-  PEdge *be, *lastbe = NULL, *maxe1 = NULL, *maxe2 = NULL, *be1, *be2;
-  PEdge *cure = NULL, *firste1 = NULL, *firste2 = NULL, *nextbe;
-  float maxlen = 0.0f, curlen = 0.0f, totlen = 0.0f, firstlen = 0.0f;
-  float len1, len2;
-
-  /* find longest series of verts split in the chart itself, these are
-   * marked during construction */
-  be = outer;
-  lastbe = p_boundary_edge_prev(be);
-  do {
-    float len = p_edge_length(be);
-    totlen += len;
-
-    nextbe = p_boundary_edge_next(be);
-
-    if ((be->vert->flag & PVERT_SPLIT) ||
-        (lastbe->vert->flag & nextbe->vert->flag & PVERT_SPLIT)) {
-      if (!cure) {
-        if (be == outer) {
-          firste1 = be;
-        }
-        cure = be;
-      }
-      else {
-        curlen += p_edge_length(lastbe);
-      }
-    }
-    else if (cure) {
-      if (curlen > maxlen) {
-        maxlen = curlen;
-        maxe1 = cure;
-        maxe2 = lastbe;
-      }
-
-      if (firste1 == cure) {
-        firstlen = curlen;
-        firste2 = lastbe;
-      }
-
-      curlen = 0.0f;
-      cure = NULL;
-    }
-
-    lastbe = be;
-    be = nextbe;
-  } while (be != outer);
-
-  /* make sure we also count a series of splits over the starting point */
-  if (cure && (cure != outer)) {
-    firstlen += curlen + p_edge_length(be);
-
-    if (firstlen > maxlen) {
-      maxlen = firstlen;
-      maxe1 = cure;
-      maxe2 = firste2;
-    }
-  }
-
-  if (!maxe1 || !maxe2 || (maxlen < 0.5f * totlen)) {
-    return P_FALSE;
-  }
-
-  /* find pin1 in the split vertices */
-  be1 = maxe1;
-  be2 = maxe2;
-  len1 = 0.0f;
-  len2 = 0.0f;
-
-  do {
-    if (len1 < len2) {
-      len1 += p_edge_length(be1);
-      be1 = p_boundary_edge_next(be1);
-    }
-    else {
-      be2 = p_boundary_edge_prev(be2);
-      len2 += p_edge_length(be2);
-    }
-  } while (be1 != be2);
-
-  *pin1 = be1->vert;
-
-  /* find pin2 outside the split vertices */
-  be1 = maxe1;
-  be2 = maxe2;
-  len1 = 0.0f;
-  len2 = 0.0f;
-
-  do {
-    if (len1 < len2) {
-      be1 = p_boundary_edge_prev(be1);
-      len1 += p_edge_length(be1);
-    }
-    else {
-      len2 += p_edge_length(be2);
-      be2 = p_boundary_edge_next(be2);
-    }
-  } while (be1 != be2);
-
-  *pin2 = be1->vert;
-
-  p_chart_pin_positions(chart, pin1, pin2);
-
-  return !equals_v3v3((*pin1)->co, (*pin2)->co);
-}
-
-static void p_chart_extrema_verts(PChart *chart, PVert **pin1, PVert **pin2)
-{
-  float minv[3], maxv[3], dirlen;
-  PVert *v, *minvert[3], *maxvert[3];
-  int i, dir;
-
-  /* find minimum and maximum verts over x/y/z axes */
-  minv[0] = minv[1] = minv[2] = 1e20;
-  maxv[0] = maxv[1] = maxv[2] = -1e20;
-
-  minvert[0] = minvert[1] = minvert[2] = NULL;
-  maxvert[0] = maxvert[1] = maxvert[2] = NULL;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    for (i = 0; i < 3; i++) {
-      if (v->co[i] < minv[i]) {
-        minv[i] = v->co[i];
-        minvert[i] = v;
-      }
-      if (v->co[i] > maxv[i]) {
-        maxv[i] = v->co[i];
-        maxvert[i] = v;
-      }
-    }
-  }
-
-  /* find axes with longest distance */
-  dir = 0;
-  dirlen = -1.0;
-
-  for (i = 0; i < 3; i++) {
-    if (maxv[i] - minv[i] > dirlen) {
-      dir = i;
-      dirlen = maxv[i] - minv[i];
-    }
-  }
-
-  *pin1 = minvert[dir];
-  *pin2 = maxvert[dir];
-
-  p_chart_pin_positions(chart, pin1, pin2);
-}
-
-static void p_chart_lscm_load_solution(PChart *chart)
-{
-  LinearSolver *context = chart->u.lscm.context;
-  PVert *v;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    v->uv[0] = EIG_linear_solver_variable_get(context, 0, 2 * v->u.id);
-    v->uv[1] = EIG_linear_solver_variable_get(context, 0, 2 * v->u.id + 1);
-  }
-}
-
-static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
-{
-  PVert *v, *pin1, *pin2;
-  PBool select = P_FALSE, deselect = P_FALSE;
-  int npins = 0, id = 0;
-
-  /* give vertices matrix indices and count pins */
-  for (v = chart->verts; v; v = v->nextlink) {
-    if (v->flag & PVERT_PIN) {
-      npins++;
-      if (v->flag & PVERT_SELECT) {
-        select = P_TRUE;
-      }
-    }
-
-    if (!(v->flag & PVERT_SELECT)) {
-      deselect = P_TRUE;
-    }
-  }
-
-  if ((live && (!select || !deselect))) {
-    chart->u.lscm.context = NULL;
-  }
-  else {
-#if 0
-    p_chart_simplify_compute(chart);
-    p_chart_topological_sanity_check(chart);
-#endif
-
-    if (npins == 1) {
-      chart->u.lscm.single_pin_area = p_chart_uv_area(chart);
-      for (v = chart->verts; v; v = v->nextlink) {
-        if (v->flag & PVERT_PIN) {
-          chart->u.lscm.single_pin = v;
-          break;
-        }
-      }
-    }
-
-    if (abf) {
-      if (!p_chart_abf_solve(chart)) {
-        param_warning("ABF solving failed: falling back to LSCM.\n");
-      }
-    }
-
-    if (npins <= 1) {
-      /* No pins, let's find some ourself. */
-      PEdge *outer;
-
-      p_chart_boundaries(chart, NULL, &outer);
-
-      /* Outer can be NULL with non-finite coords. */
-      if (!(outer && p_chart_symmetry_pins(chart, outer, &pin1, &pin2))) {
-        p_chart_extrema_verts(chart, &pin1, &pin2);
-      }
-
-      chart->u.lscm.pin1 = pin1;
-      chart->u.lscm.pin2 = pin2;
-    }
-
-    for (v = chart->verts; v; v = v->nextlink) {
-      v->u.id = id++;
-    }
-
-    chart->u.lscm.context = EIG_linear_least_squares_solver_new(
-        2 * chart->nfaces, 2 * chart->nverts, 1);
-  }
-}
-
-static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
-{
-  LinearSolver *context = chart->u.lscm.context;
-  PVert *v, *pin1 = chart->u.lscm.pin1, *pin2 = chart->u.lscm.pin2;
-  PFace *f;
-  const float *alpha = chart->u.lscm.abf_alpha;
-  float area_pinned_up, area_pinned_down;
-  bool flip_faces;
-  int row;
-
-#if 0
-  /* TODO: make loading pins work for simplify/complexify. */
-#endif
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    if (v->flag & PVERT_PIN) {
-      p_vert_load_pin_select_uvs(handle, v); /* reload for live */
-    }
-  }
-
-  if (chart->u.lscm.single_pin) {
-    /* If only one pin, save area and pin for transform later. */
-    copy_v2_v2(chart->u.lscm.single_pin_uv, chart->u.lscm.single_pin->uv);
-  }
-
-  if (chart->u.lscm.pin1) {
-    EIG_linear_solver_variable_lock(context, 2 * pin1->u.id);
-    EIG_linear_solver_variable_lock(context, 2 * pin1->u.id + 1);
-    EIG_linear_solver_variable_lock(context, 2 * pin2->u.id);
-    EIG_linear_solver_variable_lock(context, 2 * pin2->u.id + 1);
-
-    EIG_linear_solver_variable_set(context, 0, 2 * pin1->u.id, pin1->uv[0]);
-    EIG_linear_solver_variable_set(context, 0, 2 * pin1->u.id + 1, pin1->uv[1]);
-    EIG_linear_solver_variable_set(context, 0, 2 * pin2->u.id, pin2->uv[0]);
-    EIG_linear_solver_variable_set(context, 0, 2 * pin2->u.id + 1, pin2->uv[1]);
-  }
-  else {
-    /* set and lock the pins */
-    for (v = chart->verts; v; v = v->nextlink) {
-      if (v->flag & PVERT_PIN) {
-        EIG_linear_solver_variable_lock(context, 2 * v->u.id);
-        EIG_linear_solver_variable_lock(context, 2 * v->u.id + 1);
-
-        EIG_linear_solver_variable_set(context, 0, 2 * v->u.id, v->uv[0]);
-        EIG_linear_solver_variable_set(context, 0, 2 * v->u.id + 1, v->uv[1]);
-      }
-    }
-  }
-
-  /* detect up direction based on pinned vertices */
-  area_pinned_up = 0.0f;
-  area_pinned_down = 0.0f;
-
-  for (f = chart->faces; f; f = f->nextlink) {
-    PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-    PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
-
-    if ((v1->flag & PVERT_PIN) && (v2->flag & PVERT_PIN) && (v3->flag & PVERT_PIN)) {
-      float area = p_face_uv_area_signed(f);
-
-      if (area > 0.0f) {
-        area_pinned_up += area;
-      }
-      else {
-        area_pinned_down -= area;
-      }
-    }
-  }
-
-  flip_faces = (area_pinned_down > area_pinned_up);
-
-  /* construct matrix */
-
-  row = 0;
-  for (f = chart->faces; f; f = f->nextlink) {
-    PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-    PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
-    float a1, a2, a3, ratio, cosine, sine;
-    float sina1, sina2, sina3, sinmax;
-
-    if (alpha) {
-      /* use abf angles if passed on */
-      a1 = *(alpha++);
-      a2 = *(alpha++);
-      a3 = *(alpha++);
-    }
-    else {
-      p_face_angles(f, &a1, &a2, &a3);
-    }
-
-    if (flip_faces) {
-      SWAP(float, a2, a3);
-      SWAP(PEdge *, e2, e3);
-      SWAP(PVert *, v2, v3);
-    }
-
-    sina1 = sinf(a1);
-    sina2 = sinf(a2);
-    sina3 = sinf(a3);
-
-    sinmax = max_fff(sina1, sina2, sina3);
-
-    /* shift vertices to find most stable order */
-    if (sina3 != sinmax) {
-      SHIFT3(PVert *, v1, v2, v3);
-      SHIFT3(float, a1, a2, a3);
-      SHIFT3(float, sina1, sina2, sina3);
-
-      if (sina2 == sinmax) {
-        SHIFT3(PVert *, v1, v2, v3);
-        SHIFT3(float, a1, a2, a3);
-        SHIFT3(float, sina1, sina2, sina3);
-      }
-    }
-
-    /* angle based lscm formulation */
-    ratio = (sina3 == 0.0f) ? 1.0f : sina2 / sina3;
-    cosine = cosf(a1) * ratio;
-    sine = sina1 * ratio;
-
-    EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id, cosine - 1.0f);
-    EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id + 1, -sine);
-    EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id, -cosine);
-    EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id + 1, sine);
-    EIG_linear_solver_matrix_add(context, row, 2 * v3->u.id, 1.0);
-    row++;
-
-    EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id, sine);
-    EIG_linear_solver_matrix_add(context, row, 2 * v1->u.id + 1, cosine - 1.0f);
-    EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id, -sine);
-    EIG_linear_solver_matrix_add(context, row, 2 * v2->u.id + 1, -cosine);
-    EIG_linear_solver_matrix_add(context, row, 2 * v3->u.id + 1, 1.0);
-    row++;
-  }
-
-  if (EIG_linear_solver_solve(context)) {
-    p_chart_lscm_load_solution(chart);
-    return P_TRUE;
-  }
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    v->uv[0] = 0.0f;
-    v->uv[1] = 0.0f;
-  }
-
-  return P_FALSE;
-}
-
-static void p_chart_lscm_transform_single_pin(PChart *chart)
-{
-  PVert *pin = chart->u.lscm.single_pin;
-
-  /* If only one pin, keep UV area the same. */
-  const float new_area = p_chart_uv_area(chart);
-  if (new_area > 0.0f) {
-    const float scale = chart->u.lscm.single_pin_area / new_area;
-    if (scale > 0.0f) {
-      p_chart_uv_scale(chart, sqrtf(scale));
-    }
-  }
-
-  /* Translate to keep the pinned vertex in place. */
-  float offset[2];
-  sub_v2_v2v2(offset, chart->u.lscm.single_pin_uv, pin->uv);
-  p_chart_uv_translate(chart, offset);
-}
-
-static void p_chart_lscm_end(PChart *chart)
-{
-  if (chart->u.lscm.context) {
-    EIG_linear_solver_delete(chart->u.lscm.context);
-  }
-
-  MEM_SAFE_FREE(chart->u.lscm.abf_alpha);
-
-  chart->u.lscm.context = NULL;
-  chart->u.lscm.pin1 = NULL;
-  chart->u.lscm.pin2 = NULL;
-  chart->u.lscm.single_pin = NULL;
-  chart->u.lscm.single_pin_area = 0.0f;
-}
-
-/* Stretch */
-
-#define P_STRETCH_ITER 20
-
-static void p_stretch_pin_boundary(PChart *chart)
-{
-  PVert *v;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    if (v->edge->pair == NULL) {
-      v->flag |= PVERT_PIN;
-    }
-    else {
-      v->flag &= ~PVERT_PIN;
-    }
-  }
-}
-
-static float p_face_stretch(PFace *f)
-{
-  float T, w, tmp[3];
-  float Ps[3], Pt[3];
-  float a, c, area;
-  PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-  PVert *v1 = e1->vert, *v2 = e2->vert, *v3 = e3->vert;
-
-  area = p_face_uv_area_signed(f);
-
-  if (area <= 0.0f) { /* flipped face -> infinite stretch */
-    return 1e10f;
-  }
-
-  w = 1.0f / (2.0f * area);
-
-  /* compute derivatives */
-  copy_v3_v3(Ps, v1->co);
-  mul_v3_fl(Ps, (v2->uv[1] - v3->uv[1]));
-
-  copy_v3_v3(tmp, v2->co);
-  mul_v3_fl(tmp, (v3->uv[1] - v1->uv[1]));
-  add_v3_v3(Ps, tmp);
-
-  copy_v3_v3(tmp, v3->co);
-  mul_v3_fl(tmp, (v1->uv[1] - v2->uv[1]));
-  add_v3_v3(Ps, tmp);
-
-  mul_v3_fl(Ps, w);
-
-  copy_v3_v3(Pt, v1->co);
-  mul_v3_fl(Pt, (v3->uv[0] - v2->uv[0]));
-
-  copy_v3_v3(tmp, v2->co);
-  mul_v3_fl(tmp, (v1->uv[0] - v3->uv[0]));
-  add_v3_v3(Pt, tmp);
-
-  copy_v3_v3(tmp, v3->co);
-  mul_v3_fl(tmp, (v2->uv[0] - v1->uv[0]));
-  add_v3_v3(Pt, tmp);
-
-  mul_v3_fl(Pt, w);
-
-  /* Sander Tensor */
-  a = dot_v3v3(Ps, Ps);
-  c = dot_v3v3(Pt, Pt);
-
-  T = sqrtf(0.5f * (a + c));
-  if (f->flag & PFACE_FILLED) {
-    T *= 0.2f;
-  }
-
-  return T;
-}
-
-static float p_stretch_compute_vertex(PVert *v)
-{
-  PEdge *e = v->edge;
-  float sum = 0.0f;
-
-  do {
-    sum += p_face_stretch(e->face);
-    e = p_wheel_edge_next(e);
-  } while (e && e != (v->edge));
-
-  return sum;
-}
-
-static void p_chart_stretch_minimize(PChart *chart, RNG *rng)
-{
-  PVert *v;
-  PEdge *e;
-  int j, nedges;
-  float orig_stretch, low, stretch_low, high, stretch_high, mid, stretch;
-  float orig_uv[2], dir[2], random_angle, trusted_radius;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    if ((v->flag & PVERT_PIN) || !(v->flag & PVERT_SELECT)) {
-      continue;
-    }
-
-    orig_stretch = p_stretch_compute_vertex(v);
-    orig_uv[0] = v->uv[0];
-    orig_uv[1] = v->uv[1];
-
-    /* move vertex in a random direction */
-    trusted_radius = 0.0f;
-    nedges = 0;
-    e = v->edge;
-
-    do {
-      trusted_radius += p_edge_uv_length(e);
-      nedges++;
-
-      e = p_wheel_edge_next(e);
-    } while (e && e != (v->edge));
-
-    trusted_radius /= 2 * nedges;
-
-    random_angle = BLI_rng_get_float(rng) * 2.0f * (float)M_PI;
-    dir[0] = trusted_radius * cosf(random_angle);
-    dir[1] = trusted_radius * sinf(random_angle);
-
-    /* calculate old and new stretch */
-    low = 0;
-    stretch_low = orig_stretch;
-
-    add_v2_v2v2(v->uv, orig_uv, dir);
-    high = 1;
-    stretch = stretch_high = p_stretch_compute_vertex(v);
-
-    /* binary search for lowest stretch position */
-    for (j = 0; j < P_STRETCH_ITER; j++) {
-      mid = 0.5f * (low + high);
-      v->uv[0] = orig_uv[0] + mid * dir[0];
-      v->uv[1] = orig_uv[1] + mid * dir[1];
-      stretch = p_stretch_compute_vertex(v);
-
-      if (stretch_low < stretch_high) {
-        high = mid;
-        stretch_high = stretch;
-      }
-      else {
-        low = mid;
-        stretch_low = stretch;
-      }
-    }
-
-    /* no luck, stretch has increased, reset to old values */
-    if (stretch >= orig_stretch) {
-      copy_v2_v2(v->uv, orig_uv);
-    }
-  }
-}
-
-/* Minimum area enclosing rectangle for packing */
-
-static int p_compare_geometric_uv(const void *a, const void *b)
-{
-  const PVert *v1 = *(const PVert *const *)a;
-  const PVert *v2 = *(const PVert *const *)b;
-
-  if (v1->uv[0] < v2->uv[0]) {
-    return -1;
-  }
-  if (v1->uv[0] == v2->uv[0]) {
-    if (v1->uv[1] < v2->uv[1]) {
-      return -1;
-    }
-    if (v1->uv[1] == v2->uv[1]) {
-      return 0;
-    }
-    return 1;
-  }
-  return 1;
-}
-
-static PBool p_chart_convex_hull(PChart *chart, PVert ***r_verts, int *r_nverts, int *r_right)
-{
-  /* Graham algorithm, taken from:
-   * http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/117225 */
-
-  PEdge *be, *e;
-  int npoints = 0, i, ulen, llen;
-  PVert **U, **L, **points, **p;
-
-  p_chart_boundaries(chart, NULL, &be);
-
-  if (!be) {
-    return P_FALSE;
-  }
-
-  e = be;
-  do {
-    npoints++;
-    e = p_boundary_edge_next(e);
-  } while (e != be);
-
-  p = points = (PVert **)MEM_mallocN(sizeof(PVert *) * npoints * 2, "PCHullpoints");
-  U = (PVert **)MEM_mallocN(sizeof(PVert *) * npoints, "PCHullU");
-  L = (PVert **)MEM_mallocN(sizeof(PVert *) * npoints, "PCHullL");
-
-  e = be;
-  do {
-    *p = e->vert;
-    p++;
-    e = p_boundary_edge_next(e);
-  } while (e != be);
-
-  qsort(points, npoints, sizeof(PVert *), p_compare_geometric_uv);
-
-  ulen = llen = 0;
-  for (p = points, i = 0; i < npoints; i++, p++) {
-    while ((ulen > 1) && (p_area_signed(U[ulen - 2]->uv, (*p)->uv, U[ulen - 1]->uv) <= 0)) {
-      ulen--;
-    }
-    while ((llen > 1) && (p_area_signed(L[llen - 2]->uv, (*p)->uv, L[llen - 1]->uv) >= 0)) {
-      llen--;
-    }
-
-    U[ulen] = *p;
-    ulen++;
-    L[llen] = *p;
-    llen++;
-  }
-
-  npoints = 0;
-  for (p = points, i = 0; i < ulen; i++, p++, npoints++) {
-    *p = U[i];
-  }
-
-  /* the first and last point in L are left out, since they are also in U */
-  for (i = llen - 2; i > 0; i--, p++, npoints++) {
-    *p = L[i];
-  }
-
-  *r_verts = points;
-  *r_nverts = npoints;
-  *r_right = ulen - 1;
-
-  MEM_freeN(U);
-  MEM_freeN(L);
-
-  return P_TRUE;
-}
-
-static float p_rectangle_area(float *p1, float *dir, float *p2, float *p3, float *p4)
-{
-  /* given 4 points on the rectangle edges and the direction of on edge,
-   * compute the area of the rectangle */
-
-  float orthodir[2], corner1[2], corner2[2], corner3[2];
-
-  orthodir[0] = dir[1];
-  orthodir[1] = -dir[0];
-
-  if (!p_intersect_line_2d_dir(p1, dir, p2, orthodir, corner1)) {
-    return 1e10;
-  }
-
-  if (!p_intersect_line_2d_dir(p1, dir, p4, orthodir, corner2)) {
-    return 1e10;
-  }
-
-  if (!p_intersect_line_2d_dir(p3, dir, p4, orthodir, corner3)) {
-    return 1e10;
-  }
-
-  return len_v2v2(corner1, corner2) * len_v2v2(corner2, corner3);
-}
-
-static float p_chart_minimum_area_angle(PChart *chart)
-{
-  /* minimum area enclosing rectangle with rotating calipers, info:
-   * http://cgm.cs.mcgill.ca/~orm/maer.html */
-
-  float rotated, minarea, minangle, area, len;
-  float *angles, miny, maxy, v[2], a[4], mina;
-  int npoints, right, i_min, i_max, i, idx[4], nextidx;
-  PVert **points, *p1, *p2, *p3, *p4, *p1n;
-
-  /* compute convex hull */
-  if (!p_chart_convex_hull(chart, &points, &npoints, &right)) {
-    return 0.0;
-  }
-
-  /* find left/top/right/bottom points, and compute angle for each point */
-  angles = MEM_mallocN(sizeof(float) * npoints, "PMinAreaAngles");
-
-  i_min = i_max = 0;
-  miny = 1e10;
-  maxy = -1e10;
-
-  for (i = 0; i < npoints; i++) {
-    p1 = (i == 0) ? points[npoints - 1] : points[i - 1];
-    p2 = points[i];
-    p3 = (i == npoints - 1) ? points[0] : points[i + 1];
-
-    angles[i] = (float)M_PI - p_vec2_angle(p1->uv, p2->uv, p3->uv);
-
-    if (points[i]->uv[1] < miny) {
-      miny = points[i]->uv[1];
-      i_min = i;
-    }
-    if (points[i]->uv[1] > maxy) {
-      maxy = points[i]->uv[1];
-      i_max = i;
-    }
-  }
-
-  /* left, top, right, bottom */
-  idx[0] = 0;
-  idx[1] = i_max;
-  idx[2] = right;
-  idx[3] = i_min;
-
-  v[0] = points[idx[0]]->uv[0];
-  v[1] = points[idx[0]]->uv[1] + 1.0f;
-  a[0] = p_vec2_angle(points[(idx[0] + 1) % npoints]->uv, points[idx[0]]->uv, v);
-
-  v[0] = points[idx[1]]->uv[0] + 1.0f;
-  v[1] = points[idx[1]]->uv[1];
-  a[1] = p_vec2_angle(points[(idx[1] + 1) % npoints]->uv, points[idx[1]]->uv, v);
-
-  v[0] = points[idx[2]]->uv[0];
-  v[1] = points[idx[2]]->uv[1] - 1.0f;
-  a[2] = p_vec2_angle(points[(idx[2] + 1) % npoints]->uv, points[idx[2]]->uv, v);
-
-  v[0] = points[idx[3]]->uv[0] - 1.0f;
-  v[1] = points[idx[3]]->uv[1];
-  a[3] = p_vec2_angle(points[(idx[3] + 1) % npoints]->uv, points[idx[3]]->uv, v);
-
-  /* 4 rotating calipers */
-
-  rotated = 0.0;
-  minarea = 1e10;
-  minangle = 0.0;
-
-  while (rotated <= (float)M_PI_2) { /* INVESTIGATE: how far to rotate? */
-    /* rotate with the smallest angle */
-    i_min = 0;
-    mina = 1e10;
-
-    for (i = 0; i < 4; i++) {
-      if (a[i] < mina) {
-        mina = a[i];
-        i_min = i;
-      }
-    }
-
-    rotated += mina;
-    nextidx = (idx[i_min] + 1) % npoints;
-
-    a[i_min] = angles[nextidx];
-    a[(i_min + 1) % 4] = a[(i_min + 1) % 4] - mina;
-    a[(i_min + 2) % 4] = a[(i_min + 2) % 4] - mina;
-    a[(i_min + 3) % 4] = a[(i_min + 3) % 4] - mina;
-
-    /* compute area */
-    p1 = points[idx[i_min]];
-    p1n = points[nextidx];
-    p2 = points[idx[(i_min + 1) % 4]];
-    p3 = points[idx[(i_min + 2) % 4]];
-    p4 = points[idx[(i_min + 3) % 4]];
-
-    len = len_v2v2(p1->uv, p1n->uv);
-
-    if (len > 0.0f) {
-      len = 1.0f / len;
-      v[0] = (p1n->uv[0] - p1->uv[0]) * len;
-      v[1] = (p1n->uv[1] - p1->uv[1]) * len;
-
-      area = p_rectangle_area(p1->uv, v, p2->uv, p3->uv, p4->uv);
-
-      /* remember smallest area */
-      if (area < minarea) {
-        minarea = area;
-        minangle = rotated;
-      }
-    }
-
-    idx[i_min] = nextidx;
-  }
-
-  /* try keeping rotation as small as possible */
-  if (minangle > (float)M_PI_4) {
-    minangle -= (float)M_PI_2;
-  }
-
-  MEM_freeN(angles);
-  MEM_freeN(points);
-
-  return minangle;
-}
-
-static void p_chart_rotate_minimum_area(PChart *chart)
-{
-  float angle = p_chart_minimum_area_angle(chart);
-  float sine = sinf(angle);
-  float cosine = cosf(angle);
-  PVert *v;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    float oldu = v->uv[0], oldv = v->uv[1];
-    v->uv[0] = cosine * oldu - sine * oldv;
-    v->uv[1] = sine * oldu + cosine * oldv;
-  }
-}
-
-static void p_chart_rotate_fit_aabb(PChart *chart)
-{
-  float(*points)[2] = MEM_mallocN(sizeof(*points) * chart->nverts, __func__);
-
-  p_chart_uv_to_array(chart, points);
-
-  float angle = BLI_convexhull_aabb_fit_points_2d(points, chart->nverts);
-
-  MEM_freeN(points);
-
-  if (angle != 0.0f) {
-    float mat[2][2];
-    angle_to_mat2(mat, angle);
-    p_chart_uv_transform(chart, mat);
-  }
-}
-
-/* Area Smoothing */
-
-/* 2d BSP tree for inverse mapping - that's a bit silly. */
-
-typedef struct SmoothTriangle {
-  float co1[2], co2[2], co3[2];
-  float oco1[2], oco2[2], oco3[2];
-} SmoothTriangle;
-
-typedef struct SmoothNode {
-  struct SmoothNode *c1, *c2;
-  SmoothTriangle **tri;
-  float split;
-  int axis, ntri;
-} SmoothNode;
-
-static void p_barycentric_2d(
-    const float v1[2], const float v2[2], const float v3[2], const float p[2], float b[3])
-{
-  float a[2], c[2], h[2], div;
-
-  a[0] = v2[0] - v1[0];
-  a[1] = v2[1] - v1[1];
-  c[0] = v3[0] - v1[0];
-  c[1] = v3[1] - v1[1];
-
-  div = a[0] * c[1] - a[1] * c[0];
-
-  if (div == 0.0f) {
-    b[0] = 1.0f / 3.0f;
-    b[1] = 1.0f / 3.0f;
-    b[2] = 1.0f / 3.0f;
-  }
-  else {
-    h[0] = p[0] - v1[0];
-    h[1] = p[1] - v1[1];
-
-    div = 1.0f / div;
-
-    b[1] = (h[0] * c[1] - h[1] * c[0]) * div;
-    b[2] = (a[0] * h[1] - a[1] * h[0]) * div;
-    b[0] = 1.0f - b[1] - b[2];
-  }
-}
-
-static PBool p_triangle_inside(SmoothTriangle *t, float co[2])
-{
-  float b[3];
-
-  p_barycentric_2d(t->co1, t->co2, t->co3, co, b);
-
-  if ((b[0] >= 0.0f) && (b[1] >= 0.0f) && (b[2] >= 0.0f)) {
-    co[0] = t->oco1[0] * b[0] + t->oco2[0] * b[1] + t->oco3[0] * b[2];
-    co[1] = t->oco1[1] * b[0] + t->oco2[1] * b[1] + t->oco3[1] * b[2];
-    return P_TRUE;
-  }
-
-  return P_FALSE;
-}
-
-static SmoothNode *p_node_new(
-    MemArena *arena, SmoothTriangle **tri, int ntri, float *bmin, float *bmax, int depth)
-{
-  SmoothNode *node = BLI_memarena_alloc(arena, sizeof(*node));
-  int axis, i, t1size = 0, t2size = 0;
-  float split, /* mi, */ /* UNUSED */ mx;
-  SmoothTriangle **t1, **t2, *t;
-
-  node->tri = tri;
-  node->ntri = ntri;
-
-  if (ntri <= 10 || depth >= 15) {
-    return node;
-  }
-
-  t1 = MEM_mallocN(sizeof(*t1) * ntri, "PNodeTri1");
-  t2 = MEM_mallocN(sizeof(*t2) * ntri, "PNodeTri1");
-
-  axis = (bmax[0] - bmin[0] > bmax[1] - bmin[1]) ? 0 : 1;
-  split = 0.5f * (bmin[axis] + bmax[axis]);
-
-  for (i = 0; i < ntri; i++) {
-    t = tri[i];
-
-    if ((t->co1[axis] <= split) || (t->co2[axis] <= split) || (t->co3[axis] <= split)) {
-      t1[t1size] = t;
-      t1size++;
-    }
-    if ((t->co1[axis] >= split) || (t->co2[axis] >= split) || (t->co3[axis] >= split)) {
-      t2[t2size] = t;
-      t2size++;
-    }
-  }
-
-  if ((t1size == t2size) && (t1size == ntri)) {
-    MEM_freeN(t1);
-    MEM_freeN(t2);
-    return node;
-  }
-
-  node->tri = NULL;
-  node->ntri = 0;
-  MEM_freeN(tri);
-
-  node->axis = axis;
-  node->split = split;
-
-  /* mi = bmin[axis]; */ /* UNUSED */
-  mx = bmax[axis];
-  bmax[axis] = split;
-  node->c1 = p_node_new(arena, t1, t1size, bmin, bmax, depth + 1);
-
-  bmin[axis] = bmax[axis];
-  bmax[axis] = mx;
-  node->c2 = p_node_new(arena, t2, t2size, bmin, bmax, depth + 1);
-
-  return node;
-}
-
-static void p_node_delete(SmoothNode *node)
-{
-  if (node->c1) {
-    p_node_delete(node->c1);
-  }
-  if (node->c2) {
-    p_node_delete(node->c2);
-  }
-  if (node->tri) {
-    MEM_freeN(node->tri);
-  }
-}
-
-static PBool p_node_intersect(SmoothNode *node, float co[2])
-{
-  int i;
-
-  if (node->tri) {
-    for (i = 0; i < node->ntri; i++) {
-      if (p_triangle_inside(node->tri[i], co)) {
-        return P_TRUE;
-      }
-    }
-
-    return P_FALSE;
-  }
-
-  if (co[node->axis] < node->split) {
-    return p_node_intersect(node->c1, co);
-  }
-  return p_node_intersect(node->c2, co);
-}
-
-/* smoothing */
-
-static int p_compare_float(const void *a_, const void *b_)
-{
-  const float a = *(const float *)a_;
-  const float b = *(const float *)b_;
-
-  if (a < b) {
-    return -1;
-  }
-  if (a == b) {
-    return 0;
-  }
-  return 1;
-}
-
-static float p_smooth_median_edge_length(PChart *chart)
-{
-  PEdge *e;
-  float *lengths = MEM_mallocN(sizeof(chart->edges) * chart->nedges, "PMedianLength");
-  float median;
-  int i;
-
-  /* ok, so I'm lazy */
-  for (i = 0, e = chart->edges; e; e = e->nextlink, i++) {
-    lengths[i] = p_edge_length(e);
-  }
-
-  qsort(lengths, i, sizeof(float), p_compare_float);
-
-  median = lengths[i / 2];
-  MEM_freeN(lengths);
-
-  return median;
-}
-
-static float p_smooth_distortion(PEdge *e, float avg2d, float avg3d)
-{
-  float len2d = p_edge_uv_length(e) * avg3d;
-  float len3d = p_edge_length(e) * avg2d;
-
-  return (len3d == 0.0f) ? 0.0f : len2d / len3d;
-}
-
-static void p_smooth(PChart *chart)
-{
-  PEdge *e;
-  PVert *v;
-  PFace *f;
-  int j, it2, maxiter2, it;
-  int nedges = chart->nedges, nwheel, gridx, gridy;
-  int edgesx, edgesy, nsize, esize, i, x, y, maxiter;
-  float minv[2], maxv[2], median, invmedian, avglen2d, avglen3d;
-  float center[2], dx, dy, *nodes, dlimit, d, *oldnodesx, *oldnodesy;
-  float *nodesx, *nodesy, *hedges, *vedges, climit, moved, padding;
-  SmoothTriangle *triangles, *t, *t2, **tri, **trip;
-  SmoothNode *root;
-  MemArena *arena;
-
-  if (nedges == 0) {
-    return;
-  }
-
-  p_chart_uv_bbox(chart, minv, maxv);
-  median = p_smooth_median_edge_length(chart) * 0.10f;
-
-  if (median == 0.0f) {
-    return;
-  }
-
-  invmedian = 1.0f / median;
-
-  /* compute edge distortion */
-  avglen2d = avglen3d = 0.0;
-
-  for (e = chart->edges; e; e = e->nextlink) {
-    avglen2d += p_edge_uv_length(e);
-    avglen3d += p_edge_length(e);
-  }
-
-  avglen2d /= nedges;
-  avglen3d /= nedges;
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    v->u.distortion = 0.0;
-    nwheel = 0;
-
-    e = v->edge;
-    do {
-      v->u.distortion += p_smooth_distortion(e, avglen2d, avglen3d);
-      nwheel++;
-
-      e = e->next->next->pair;
-    } while (e && (e != v->edge));
-
-    v->u.distortion /= nwheel;
-  }
-
-  /* need to do excessive grid size checking still */
-  center[0] = 0.5f * (minv[0] + maxv[0]);
-  center[1] = 0.5f * (minv[1] + maxv[1]);
-
-  dx = 0.5f * (maxv[0] - minv[0]);
-  dy = 0.5f * (maxv[1] - minv[1]);
-
-  padding = 0.15f;
-  dx += padding * dx + 2.0f * median;
-  dy += padding * dy + 2.0f * median;
-
-  gridx = (int)(dx * invmedian);
-  gridy = (int)(dy * invmedian);
-
-  minv[0] = center[0] - median * gridx;
-  minv[1] = center[1] - median * gridy;
-  maxv[0] = center[0] + median * gridx;
-  maxv[1] = center[1] + median * gridy;
-
-  /* create grid */
-  gridx = gridx * 2 + 1;
-  gridy = gridy * 2 + 1;
-
-  if ((gridx <= 2) || (gridy <= 2)) {
-    return;
-  }
-
-  edgesx = gridx - 1;
-  edgesy = gridy - 1;
-  nsize = gridx * gridy;
-  esize = edgesx * edgesy;
-
-  nodes = MEM_mallocN(sizeof(float) * nsize, "PSmoothNodes");
-  nodesx = MEM_mallocN(sizeof(float) * nsize, "PSmoothNodesX");
-  nodesy = MEM_mallocN(sizeof(float) * nsize, "PSmoothNodesY");
-  oldnodesx = MEM_mallocN(sizeof(float) * nsize, "PSmoothOldNodesX");
-  oldnodesy = MEM_mallocN(sizeof(float) * nsize, "PSmoothOldNodesY");
-  hedges = MEM_mallocN(sizeof(float) * esize, "PSmoothHEdges");
-  vedges = MEM_mallocN(sizeof(float) * esize, "PSmoothVEdges");
-
-  if (!nodes || !nodesx || !nodesy || !oldnodesx || !oldnodesy || !hedges || !vedges) {
-    if (nodes) {
-      MEM_freeN(nodes);
-    }
-    if (nodesx) {
-      MEM_freeN(nodesx);
-    }
-    if (nodesy) {
-      MEM_freeN(nodesy);
-    }
-    if (oldnodesx) {
-      MEM_freeN(oldnodesx);
-    }
-    if (oldnodesy) {
-      MEM_freeN(oldnodesy);
-    }
-    if (hedges) {
-      MEM_freeN(hedges);
-    }
-    if (vedges) {
-      MEM_freeN(vedges);
-    }
-
-    // printf("Not enough memory for area smoothing grid");
-    return;
-  }
-
-  for (x = 0; x < gridx; x++) {
-    for (y = 0; y < gridy; y++) {
-      i = x + y * gridx;
-
-      nodesx[i] = minv[0] + median * x;
-      nodesy[i] = minv[1] + median * y;
-
-      nodes[i] = 1.0f;
-    }
-  }
-
-  /* embed in grid */
-  for (f = chart->faces; f; f = f->nextlink) {
-    PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
-    float fmin[2], fmax[2];
-    int bx1, by1, bx2, by2;
-
-    INIT_MINMAX2(fmin, fmax);
-
-    minmax_v2v2_v2(fmin, fmax, e1->vert->uv);
-    minmax_v2v2_v2(fmin, fmax, e2->vert->uv);
-    minmax_v2v2_v2(fmin, fmax, e3->vert->uv);
-
-    bx1 = (int)((fmin[0] - minv[0]) * invmedian);
-    by1 = (int)((fmin[1] - minv[1]) * invmedian);
-    bx2 = (int)((fmax[0] - minv[0]) * invmedian + 2);
-    by2 = (int)((fmax[1] - minv[1]) * invmedian + 2);
-
-    for (x = bx1; x < bx2; x++) {
-      for (y = by1; y < by2; y++) {
-        float p[2], b[3];
-
-        i = x + y * gridx;
-
-        p[0] = nodesx[i];
-        p[1] = nodesy[i];
-
-        p_barycentric_2d(e1->vert->uv, e2->vert->uv, e3->vert->uv, p, b);
-
-        if ((b[0] > 0.0f) && (b[1] > 0.0f) && (b[2] > 0.0f)) {
-          nodes[i] = e1->vert->u.distortion * b[0];
-          nodes[i] += e2->vert->u.distortion * b[1];
-          nodes[i] += e3->vert->u.distortion * b[2];
-        }
-      }
-    }
-  }
-
-  /* smooth the grid */
-  maxiter = 10;
-  climit = 0.00001f * nsize;
-
-  for (it = 0; it < maxiter; it++) {
-    moved = 0.0f;
-
-    for (x = 0; x < edgesx; x++) {
-      for (y = 0; y < edgesy; y++) {
-        i = x + y * gridx;
-        j = x + y * edgesx;
-
-        hedges[j] = (nodes[i] + nodes[i + 1]) * 0.5f;
-        vedges[j] = (nodes[i] + nodes[i + gridx]) * 0.5f;
-
-        /* we do *inverse* mapping */
-        hedges[j] = 1.0f / hedges[j];
-        vedges[j] = 1.0f / vedges[j];
-      }
-    }
-
-    maxiter2 = 50;
-    dlimit = 0.0001f;
-
-    for (it2 = 0; it2 < maxiter2; it2++) {
-      d = 0.0f;
-
-      memcpy(oldnodesx, nodesx, sizeof(float) * nsize);
-      memcpy(oldnodesy, nodesy, sizeof(float) * nsize);
-
-      for (x = 1; x < gridx - 1; x++) {
-        for (y = 1; y < gridy - 1; y++) {
-          float p[2], oldp[2], sum1, sum2, diff[2], length;
-
-          i = x + gridx * y;
-          j = x + edgesx * y;
-
-          oldp[0] = oldnodesx[i];
-          oldp[1] = oldnodesy[i];
-
-          sum1 = hedges[j - 1] * oldnodesx[i - 1];
-          sum1 += hedges[j] * oldnodesx[i + 1];
-          sum1 += vedges[j - edgesx] * oldnodesx[i - gridx];
-          sum1 += vedges[j] * oldnodesx[i + gridx];
-
-          sum2 = hedges[j - 1];
-          sum2 += hedges[j];
-          sum2 += vedges[j - edgesx];
-          sum2 += vedges[j];
-
-          nodesx[i] = sum1 / sum2;
-
-          sum1 = hedges[j - 1] * oldnodesy[i - 1];
-          sum1 += hedges[j] * oldnodesy[i + 1];
-          sum1 += vedges[j - edgesx] * oldnodesy[i - gridx];
-          sum1 += vedges[j] * oldnodesy[i + gridx];
-
-          nodesy[i] = sum1 / sum2;
-
-          p[0] = nodesx[i];
-          p[1] = nodesy[i];
-
-          diff[0] = p[0] - oldp[0];
-          diff[1] = p[1] - oldp[1];
-
-          length = len_v2(diff);
-          d = max_ff(d, length);
-          moved += length;
-        }
-      }
-
-      if (d < dlimit) {
-        break;
-      }
-    }
-
-    if (moved < climit) {
-      break;
-    }
-  }
-
-  MEM_freeN(oldnodesx);
-  MEM_freeN(oldnodesy);
-  MEM_freeN(hedges);
-  MEM_freeN(vedges);
-
-  /* Create BSP. */
-  t = triangles = MEM_mallocN(sizeof(SmoothTriangle) * esize * 2, "PSmoothTris");
-  trip = tri = MEM_mallocN(sizeof(SmoothTriangle *) * esize * 2, "PSmoothTriP");
-
-  if (!triangles || !tri) {
-    MEM_freeN(nodes);
-    MEM_freeN(nodesx);
-    MEM_freeN(nodesy);
-
-    if (triangles) {
-      MEM_freeN(triangles);
-    }
-    if (tri) {
-      MEM_freeN(tri);
-    }
-
-    // printf("Not enough memory for area smoothing grid");
-    return;
-  }
-
-  for (x = 0; x < edgesx; x++) {
-    for (y = 0; y < edgesy; y++) {
-      i = x + y * gridx;
-
-      t->co1[0] = nodesx[i];
-      t->co1[1] = nodesy[i];
-
-      t->co2[0] = nodesx[i + 1];
-      t->co2[1] = nodesy[i + 1];
-
-      t->co3[0] = nodesx[i + gridx];
-      t->co3[1] = nodesy[i + gridx];
-
-      t->oco1[0] = minv[0] + x * median;
-      t->oco1[1] = minv[1] + y * median;
-
-      t->oco2[0] = minv[0] + (x + 1) * median;
-      t->oco2[1] = minv[1] + y * median;
-
-      t->oco3[0] = minv[0] + x * median;
-      t->oco3[1] = minv[1] + (y + 1) * median;
-
-      t2 = t + 1;
-
-      t2->co1[0] = nodesx[i + gridx + 1];
-      t2->co1[1] = nodesy[i + gridx + 1];
-
-      t2->oco1[0] = minv[0] + (x + 1) * median;
-      t2->oco1[1] = minv[1] + (y + 1) * median;
-
-      t2->co2[0] = t->co2[0];
-      t2->co2[1] = t->co2[1];
-      t2->oco2[0] = t->oco2[0];
-      t2->oco2[1] = t->oco2[1];
-
-      t2->co3[0] = t->co3[0];
-      t2->co3[1] = t->co3[1];
-      t2->oco3[0] = t->oco3[0];
-      t2->oco3[1] = t->oco3[1];
-
-      *trip = t;
-      trip++;
-      t++;
-      *trip = t;
-      trip++;
-      t++;
-    }
-  }
-
-  MEM_freeN(nodes);
-  MEM_freeN(nodesx);
-  MEM_freeN(nodesy);
-
-  arena = BLI_memarena_new(MEM_SIZE_OPTIMAL(1 << 16), "param smooth arena");
-  root = p_node_new(arena, tri, esize * 2, minv, maxv, 0);
-
-  for (v = chart->verts; v; v = v->nextlink) {
-    if (!p_node_intersect(root, v->uv)) {
-      param_warning("area smoothing error: couldn't find mapping triangle\n");
-    }
-  }
-
-  p_node_delete(root);
-  BLI_memarena_free(arena);
-
-  MEM_freeN(triangles);
-}
-
-/* Exported */
-
-ParamHandle *param_construct_begin(void)
-{
-  PHandle *handle = MEM_callocN(sizeof(*handle), "PHandle");
-  handle->construction_chart = p_chart_new(handle);
-  handle->state = PHANDLE_STATE_ALLOCATED;
-  handle->arena = BLI_memarena_new(MEM_SIZE_OPTIMAL(1 << 16), "param construct arena");
-  handle->polyfill_arena = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "param polyfill arena");
-  handle->polyfill_heap = BLI_heap_new_ex(BLI_POLYFILL_ALLOC_NGON_RESERVE);
-  handle->aspx = 1.0f;
-  handle->aspy = 1.0f;
-  handle->do_aspect = false;
-
-  handle->hash_verts = phash_new((PHashLink **)&handle->construction_chart->verts, 1);
-  handle->hash_edges = phash_new((PHashLink **)&handle->construction_chart->edges, 1);
-  handle->hash_faces = phash_new((PHashLink **)&handle->construction_chart->faces, 1);
-
-  return (ParamHandle *)handle;
-}
-
-void param_aspect_ratio(ParamHandle *handle, float aspx, float aspy)
-{
-  PHandle *phandle = (PHandle *)handle;
-
-  phandle->aspx = aspx;
-  phandle->aspy = aspy;
-  phandle->do_aspect = true;
-}
-
-void param_delete(ParamHandle *handle)
-{
-  PHandle *phandle = (PHandle *)handle;
-  int i;
-
-  param_assert(ELEM(phandle->state, PHANDLE_STATE_ALLOCATED, PHANDLE_STATE_CONSTRUCTED));
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    p_chart_delete(phandle->charts[i]);
-  }
-
-  if (phandle->charts) {
-    MEM_freeN(phandle->charts);
-  }
-
-  if (phandle->construction_chart) {
-    p_chart_delete(phandle->construction_chart);
-
-    phash_delete(phandle->hash_verts);
-    phash_delete(phandle->hash_edges);
-    phash_delete(phandle->hash_faces);
-  }
-
-  BLI_memarena_free(phandle->arena);
-  BLI_memarena_free(phandle->polyfill_arena);
-  BLI_heap_free(phandle->polyfill_heap, NULL);
-  MEM_freeN(phandle);
-}
-
-static void p_add_ngon(ParamHandle *handle,
-                       ParamKey key,
-                       int nverts,
-                       ParamKey *vkeys,
-                       float **co,
-                       float **uv,
-                       ParamBool *pin,
-                       ParamBool *select)
-{
-  /* Allocate memory for polyfill. */
-  PHandle *phandle = (PHandle *)handle;
-  MemArena *arena = phandle->polyfill_arena;
-  Heap *heap = phandle->polyfill_heap;
-  uint nfilltri = nverts - 2;
-  uint(*tris)[3] = BLI_memarena_alloc(arena, sizeof(*tris) * (size_t)nfilltri);
-  float(*projverts)[2] = BLI_memarena_alloc(arena, sizeof(*projverts) * (size_t)nverts);
-
-  /* Calc normal, flipped: to get a positive 2d cross product. */
-  float normal[3];
-  zero_v3(normal);
-
-  const float *co_curr, *co_prev = co[nverts - 1];
-  for (int j = 0; j < nverts; j++) {
-    co_curr = co[j];
-    add_newell_cross_v3_v3v3(normal, co_prev, co_curr);
-    co_prev = co_curr;
-  }
-  if (UNLIKELY(normalize_v3(normal) == 0.0f)) {
-    normal[2] = 1.0f;
-  }
-
-  /* Project verts to 2d. */
-  float axis_mat[3][3];
-  axis_dominant_v3_to_m3_negate(axis_mat, normal);
-  for (int j = 0; j < nverts; j++) {
-    mul_v2_m3v3(projverts[j], axis_mat, co[j]);
-  }
-
-  BLI_polyfill_calc_arena(projverts, nverts, 1, tris, arena);
-
-  /* Beautify helps avoid thin triangles that give numerical problems. */
-  BLI_polyfill_beautify(projverts, nverts, tris, arena, heap);
-
-  /* Add triangles. */
-  for (int j = 0; j < nfilltri; j++) {
-    uint *tri = tris[j];
-    uint v0 = tri[0];
-    uint v1 = tri[1];
-    uint v2 = tri[2];
-
-    ParamKey tri_vkeys[3] = {vkeys[v0], vkeys[v1], vkeys[v2]};
-    float *tri_co[3] = {co[v0], co[v1], co[v2]};
-    float *tri_uv[3] = {uv[v0], uv[v1], uv[v2]};
-    ParamBool tri_pin[3] = {pin[v0], pin[v1], pin[v2]};
-    ParamBool tri_select[3] = {select[v0], select[v1], select[v2]};
-
-    param_face_add(handle, key, 3, tri_vkeys, tri_co, tri_uv, tri_pin, tri_select);
-  }
-
-  BLI_memarena_clear(arena);
-}
-
-void param_face_add(ParamHandle *handle,
-                    ParamKey key,
-                    int nverts,
-                    ParamKey *vkeys,
-                    float *co[4],
-                    float *uv[4],
-                    ParamBool *pin,
-                    ParamBool *select)
-{
-  PHandle *phandle = (PHandle *)handle;
-
-  param_assert(phash_lookup(phandle->hash_faces, key) == NULL);
-  param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
-  param_assert(ELEM(nverts, 3, 4));
-
-  if (nverts > 4) {
-    /* ngon */
-    p_add_ngon(handle, key, nverts, vkeys, co, uv, pin, select);
-  }
-  else if (nverts == 4) {
-    /* quad */
-    if (p_quad_split_direction(phandle, co, vkeys)) {
-      p_face_add_construct(phandle, key, vkeys, co, uv, 0, 1, 2, pin, select);
-      p_face_add_construct(phandle, key, vkeys, co, uv, 0, 2, 3, pin, select);
-    }
-    else {
-      p_face_add_construct(phandle, key, vkeys, co, uv, 0, 1, 3, pin, select);
-      p_face_add_construct(phandle, key, vkeys, co, uv, 1, 2, 3, pin, select);
-    }
-  }
-  else if (!p_face_exists(phandle, vkeys, 0, 1, 2)) {
-    /* triangle */
-    p_face_add_construct(phandle, key, vkeys, co, uv, 0, 1, 2, pin, select);
-  }
-}
-
-void param_edge_set_seam(ParamHandle *handle, ParamKey *vkeys)
-{
-  PHandle *phandle = (PHandle *)handle;
-  PEdge *e;
-
-  param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
-
-  e = p_edge_lookup(phandle, vkeys);
-  if (e) {
-    e->flag |= PEDGE_SEAM;
-  }
-}
-
-void param_construct_end(ParamHandle *handle,
-                         ParamBool fill,
-                         ParamBool topology_from_uvs,
-                         int *count_fail)
-{
-  PHandle *phandle = (PHandle *)handle;
-  PChart *chart = phandle->construction_chart;
-  int i, j, nboundaries = 0;
-  PEdge *outer;
-
-  param_assert(phandle->state == PHANDLE_STATE_ALLOCATED);
-
-  phandle->ncharts = p_connect_pairs(phandle, (PBool)topology_from_uvs);
-  phandle->charts = p_split_charts(phandle, chart, phandle->ncharts);
-
-  p_chart_delete(phandle->construction_chart);
-  phandle->construction_chart = NULL;
-
-  phash_delete(phandle->hash_verts);
-  phash_delete(phandle->hash_edges);
-  phash_delete(phandle->hash_faces);
-  phandle->hash_verts = phandle->hash_edges = phandle->hash_faces = NULL;
-
-  for (i = j = 0; i < phandle->ncharts; i++) {
-    PVert *v;
-    chart = phandle->charts[i];
-
-    p_chart_boundaries(chart, &nboundaries, &outer);
-
-    if (!topology_from_uvs && nboundaries == 0) {
-      p_chart_delete(chart);
-      if (count_fail != NULL) {
-        *count_fail += 1;
-      }
-      continue;
-    }
-
-    phandle->charts[j] = chart;
-    j++;
-
-    if (fill && (nboundaries > 1)) {
-      p_chart_fill_boundaries(chart, outer);
-    }
-
-    for (v = chart->verts; v; v = v->nextlink) {
-      p_vert_load_pin_select_uvs(handle, v);
-    }
-  }
-
-  phandle->ncharts = j;
-
-  phandle->state = PHANDLE_STATE_CONSTRUCTED;
-}
-
-void param_lscm_begin(ParamHandle *handle, ParamBool live, ParamBool abf)
-{
-  PHandle *phandle = (PHandle *)handle;
-  PFace *f;
-  int i;
-
-  param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
-  phandle->state = PHANDLE_STATE_LSCM;
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    for (f = phandle->charts[i]->faces; f; f = f->nextlink) {
-      p_face_backup_uvs(f);
-    }
-    p_chart_lscm_begin(phandle->charts[i], (PBool)live, (PBool)abf);
-  }
-}
-
-void param_lscm_solve(ParamHandle *handle, int *count_changed, int *count_failed)
-{
-  PHandle *phandle = (PHandle *)handle;
-  PChart *chart;
-  int i;
-
-  param_assert(phandle->state == PHANDLE_STATE_LSCM);
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    chart = phandle->charts[i];
-
-    if (chart->u.lscm.context) {
-      const PBool result = p_chart_lscm_solve(phandle, chart);
-
-      if (result && !(chart->flag & PCHART_HAS_PINS)) {
-        p_chart_rotate_minimum_area(chart);
-      }
-      else if (result && chart->u.lscm.single_pin) {
-        p_chart_rotate_fit_aabb(chart);
-        p_chart_lscm_transform_single_pin(chart);
-      }
-
-      if (!result || !(chart->flag & PCHART_HAS_PINS)) {
-        p_chart_lscm_end(chart);
-      }
-
-      if (result) {
-        if (count_changed != NULL) {
-          *count_changed += 1;
-        }
-      }
-      else {
-        if (count_failed != NULL) {
-          *count_failed += 1;
-        }
-      }
-    }
-  }
-}
-
-void param_lscm_end(ParamHandle *handle)
-{
-  PHandle *phandle = (PHandle *)handle;
-  int i;
-
-  param_assert(phandle->state == PHANDLE_STATE_LSCM);
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    p_chart_lscm_end(phandle->charts[i]);
-#if 0
-    p_chart_complexify(phandle->charts[i]);
-#endif
-  }
-
-  phandle->state = PHANDLE_STATE_CONSTRUCTED;
-}
-
-void param_stretch_begin(ParamHandle *handle)
-{
-  PHandle *phandle = (PHandle *)handle;
-  PChart *chart;
-  PVert *v;
-  PFace *f;
-  int i;
-
-  param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
-  phandle->state = PHANDLE_STATE_STRETCH;
-
-  phandle->rng = BLI_rng_new(31415926);
-  phandle->blend = 0.0f;
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    chart = phandle->charts[i];
-
-    for (v = chart->verts; v; v = v->nextlink) {
-      v->flag &= ~PVERT_PIN; /* don't use user-defined pins */
-    }
-
-    p_stretch_pin_boundary(chart);
-
-    for (f = chart->faces; f; f = f->nextlink) {
-      p_face_backup_uvs(f);
-      f->u.area3d = p_face_area(f);
-    }
-  }
-}
-
-void param_stretch_blend(ParamHandle *handle, float blend)
-{
-  PHandle *phandle = (PHandle *)handle;
-
-  param_assert(phandle->state == PHANDLE_STATE_STRETCH);
-  phandle->blend = blend;
-}
-
-void param_stretch_iter(ParamHandle *handle)
-{
-  PHandle *phandle = (PHandle *)handle;
-  PChart *chart;
-  int i;
-
-  param_assert(phandle->state == PHANDLE_STATE_STRETCH);
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    chart = phandle->charts[i];
-    p_chart_stretch_minimize(chart, phandle->rng);
-  }
-}
-
-void param_stretch_end(ParamHandle *handle)
-{
-  PHandle *phandle = (PHandle *)handle;
-
-  param_assert(phandle->state == PHANDLE_STATE_STRETCH);
-  phandle->state = PHANDLE_STATE_CONSTRUCTED;
-
-  BLI_rng_free(phandle->rng);
-  phandle->rng = NULL;
-}
-
-void param_smooth_area(ParamHandle *handle)
-{
-  PHandle *phandle = (PHandle *)handle;
-  int i;
-
-  param_assert(phandle->state == PHANDLE_STATE_CONSTRUCTED);
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    PChart *chart = phandle->charts[i];
-    PVert *v;
-
-    for (v = chart->verts; v; v = v->nextlink) {
-      v->flag &= ~PVERT_PIN;
-    }
-
-    p_smooth(chart);
-  }
-}
-
-/* don't pack, just rotate (used for better packing) */
-static void param_pack_rotate(ParamHandle *handle, bool ignore_pinned)
-{
-  PChart *chart;
-  int i;
-
-  PHandle *phandle = (PHandle *)handle;
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    chart = phandle->charts[i];
-
-    if (ignore_pinned && (chart->flag & PCHART_HAS_PINS)) {
-      continue;
-    }
-
-    p_chart_rotate_fit_aabb(chart);
-  }
-}
-
-void param_pack(ParamHandle *handle, float margin, bool do_rotate, bool ignore_pinned)
-{
-  /* box packing variables */
-  BoxPack *boxarray, *box;
-  float tot_width, tot_height, scale;
-
-  PChart *chart;
-  int i, unpacked = 0;
-  float trans[2];
-  double area = 0.0;
-
-  PHandle *phandle = (PHandle *)handle;
-
-  if (phandle->ncharts == 0) {
-    return;
-  }
-
-  /* this could be its own function */
-  if (do_rotate) {
-    param_pack_rotate(handle, ignore_pinned);
-  }
-
-  if (phandle->aspx != phandle->aspy) {
-    param_scale(handle, 1.0f / phandle->aspx, 1.0f / phandle->aspy);
-  }
-
-  /* we may not use all these boxes */
-  boxarray = MEM_mallocN(phandle->ncharts * sizeof(BoxPack), "BoxPack box");
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    chart = phandle->charts[i];
-
-    if (ignore_pinned && (chart->flag & PCHART_HAS_PINS)) {
-      unpacked++;
-      continue;
-    }
-
-    box = boxarray + (i - unpacked);
-
-    p_chart_uv_bbox(chart, trans, chart->u.pack.size);
-
-    trans[0] = -trans[0];
-    trans[1] = -trans[1];
-
-    p_chart_uv_translate(chart, trans);
-
-    box->w = chart->u.pack.size[0] + trans[0];
-    box->h = chart->u.pack.size[1] + trans[1];
-    box->index = i; /* warning this index skips PCHART_HAS_PINS boxes */
-
-    if (margin > 0.0f) {
-      area += (double)sqrtf(box->w * box->h);
-    }
-  }
-
-  if (margin > 0.0f) {
-    /* multiply the margin by the area to give predictable results not dependent on UV scale,
-     * ...Without using the area running pack multiple times also gives a bad feedback loop.
-     * multiply by 0.1 so the margin value from the UI can be from
-     * 0.0 to 1.0 but not give a massive margin */
-    margin = (margin * (float)area) * 0.1f;
-    unpacked = 0;
-    for (i = 0; i < phandle->ncharts; i++) {
-      chart = phandle->charts[i];
-
-      if (ignore_pinned && (chart->flag & PCHART_HAS_PINS)) {
-        unpacked++;
-        continue;
-      }
-
-      box = boxarray + (i - unpacked);
-      trans[0] = margin;
-      trans[1] = margin;
-      p_chart_uv_translate(chart, trans);
-      box->w += margin * 2;
-      box->h += margin * 2;
-    }
-  }
-
-  BLI_box_pack_2d(boxarray, phandle->ncharts - unpacked, &tot_width, &tot_height);
-
-  if (tot_height > tot_width) {
-    scale = 1.0f / tot_height;
-  }
-  else {
-    scale = 1.0f / tot_width;
-  }
-
-  for (i = 0; i < phandle->ncharts - unpacked; i++) {
-    box = boxarray + i;
-    trans[0] = box->x;
-    trans[1] = box->y;
-
-    chart = phandle->charts[box->index];
-    p_chart_uv_translate(chart, trans);
-    p_chart_uv_scale(chart, scale);
-  }
-  MEM_freeN(boxarray);
-
-  if (phandle->aspx != phandle->aspy) {
-    param_scale(handle, phandle->aspx, phandle->aspy);
-  }
-}
-
-void param_average(ParamHandle *handle, bool ignore_pinned)
-{
-  PChart *chart;
-  int i;
-  float tot_uvarea = 0.0f, tot_facearea = 0.0f;
-  float tot_fac, fac;
-  float minv[2], maxv[2], trans[2];
-  PHandle *phandle = (PHandle *)handle;
-
-  if (phandle->ncharts == 0) {
-    return;
-  }
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    PFace *f;
-    chart = phandle->charts[i];
-
-    if (ignore_pinned && (chart->flag & PCHART_HAS_PINS)) {
-      continue;
-    }
-
-    chart->u.pack.area = 0.0f;    /* 3d area */
-    chart->u.pack.rescale = 0.0f; /* UV area, abusing rescale for tmp storage, oh well :/ */
-
-    for (f = chart->faces; f; f = f->nextlink) {
-      chart->u.pack.area += p_face_area(f);
-      chart->u.pack.rescale += fabsf(p_face_uv_area_signed(f));
-    }
-
-    tot_facearea += chart->u.pack.area;
-    tot_uvarea += chart->u.pack.rescale;
-  }
-
-  if (tot_facearea == tot_uvarea || tot_facearea == 0.0f || tot_uvarea == 0.0f) {
-    /* nothing to do */
-    return;
-  }
-
-  tot_fac = tot_facearea / tot_uvarea;
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    chart = phandle->charts[i];
-
-    if (ignore_pinned && (chart->flag & PCHART_HAS_PINS)) {
-      continue;
-    }
-
-    if (chart->u.pack.area != 0.0f && chart->u.pack.rescale != 0.0f) {
-      fac = chart->u.pack.area / chart->u.pack.rescale;
-
-      /* Get the island center */
-      p_chart_uv_bbox(chart, minv, maxv);
-      trans[0] = (minv[0] + maxv[0]) / -2.0f;
-      trans[1] = (minv[1] + maxv[1]) / -2.0f;
-
-      /* Move center to 0,0 */
-      p_chart_uv_translate(chart, trans);
-      p_chart_uv_scale(chart, sqrtf(fac / tot_fac));
-
-      /* Move to original center */
-      trans[0] = -trans[0];
-      trans[1] = -trans[1];
-      p_chart_uv_translate(chart, trans);
-    }
-  }
-}
-
-void param_scale(ParamHandle *handle, float x, float y)
-{
-  PHandle *phandle = (PHandle *)handle;
-  PChart *chart;
-  int i;
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    chart = phandle->charts[i];
-    p_chart_uv_scale_xy(chart, x, y);
-  }
-}
-
-void param_flush(ParamHandle *handle)
-{
-  PHandle *phandle = (PHandle *)handle;
-  PChart *chart;
-  int i;
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    chart = phandle->charts[i];
-
-    if ((phandle->state == PHANDLE_STATE_LSCM) && !chart->u.lscm.context) {
-      continue;
-    }
-
-    if (phandle->blend == 0.0f) {
-      p_flush_uvs(phandle, chart);
-    }
-    else {
-      p_flush_uvs_blend(phandle, chart, phandle->blend);
-    }
-  }
-}
-
-void param_flush_restore(ParamHandle *handle)
-{
-  PHandle *phandle = (PHandle *)handle;
-  PChart *chart;
-  PFace *f;
-  int i;
-
-  for (i = 0; i < phandle->ncharts; i++) {
-    chart = phandle->charts[i];
-
-    for (f = chart->faces; f; f = f->nextlink) {
-      p_face_restore_uvs(f);
-    }
-  }
-}
diff --git a/source/blender/editors/uvedit/uvedit_parametrizer.h b/source/blender/editors/uvedit/uvedit_parametrizer.h
deleted file mode 100644
index f234fbe2ace..00000000000
--- a/source/blender/editors/uvedit/uvedit_parametrizer.h
+++ /dev/null
@@ -1,99 +0,0 @@
-/* SPDX-License-Identifier: GPL-2.0-or-later */
-
-#pragma once
-
-/** \file
- * \ingroup eduv
- */
-
-#include "BLI_sys_types.h" /* for intptr_t support */
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-typedef void ParamHandle;  /* handle to a set of charts */
-typedef intptr_t ParamKey; /* (hash) key for identifying verts and faces */
-typedef enum ParamBool {
-  PARAM_TRUE = 1,
-  PARAM_FALSE = 0,
-} ParamBool;
-
-/* Chart construction:
- * -------------------
- * - faces and seams may only be added between construct_{begin|end}
- * - the pointers to co and uv are stored, rather than being copied
- * - vertices are implicitly created
- * - in construct_end the mesh will be split up according to the seams
- * - the resulting charts must be:
- * - manifold, connected, open (at least one boundary loop)
- * - output will be written to the uv pointers
- */
-
-ParamHandle *param_construct_begin(void);
-
-void param_aspect_ratio(ParamHandle *handle, float aspx, float aspy);
-
-void param_face_add(ParamHandle *handle,
-                    ParamKey key,
-                    int nverts,
-                    ParamKey *vkeys,
-                    float *co[4],
-                    float *uv[4],
-                    ParamBool *pin,
-                    ParamBool *select);
-
-void param_edge_set_seam(ParamHandle *handle, ParamKey *vkeys);
-
-void param_construct_end(ParamHandle *handle,
-                         ParamBool fill,
-                         ParamBool topology_from_uvs,
-                         int *count_fail);
-void param_delete(ParamHandle *handle);
-
-/* Least Squares Conformal Maps:
- * -----------------------------
- * - charts with less than two pinned vertices are assigned 2 pins
- * - lscm is divided in three steps:
- * - begin: compute matrix and its factorization (expensive)
- * - solve using pinned coordinates (cheap)
- * - end: clean up
- * - uv coordinates are allowed to change within begin/end, for
- *   quick re-solving
- */
-
-void param_lscm_begin(ParamHandle *handle, ParamBool live, ParamBool abf);
-void param_lscm_solve(ParamHandle *handle, int *count_changed, int *count_failed);
-void param_lscm_end(ParamHandle *handle);
-
-/* Stretch */
-
-void param_stretch_begin(ParamHandle *handle);
-void param_stretch_blend(ParamHandle *handle, float blend);
-void param_stretch_iter(ParamHandle *handle);
-void param_stretch_end(ParamHandle *handle);
-
-/* Area Smooth */
-
-void param_smooth_area(ParamHandle *handle);
-
-/* Packing */
-
-void param_pack(ParamHandle *handle, float margin, bool do_rotate, bool ignore_pinned);
-
-/* Average area for all charts */
-
-void param_average(ParamHandle *handle, bool ignore_pinned);
-
-/* Simple x,y scale */
-
-void param_scale(ParamHandle *handle, float x, float y);
-
-/* Flushing */
-
-void param_flush(ParamHandle *handle);
-void param_flush_restore(ParamHandle *handle);
-
-#ifdef __cplusplus
-}
-#endif
diff --git a/source/blender/editors/uvedit/uvedit_unwrap_ops.c b/source/blender/editors/uvedit/uvedit_unwrap_ops.c
index 63300656fda..609fa72d56b 100644
--- a/source/blender/editors/uvedit/uvedit_unwrap_ops.c
+++ b/source/blender/editors/uvedit/uvedit_unwrap_ops.c
@@ -46,6 +46,8 @@
 
 #include "DEG_depsgraph.h"
 
+#include "GEO_uv_parametrizer.h"
+
 #include "PIL_time.h"
 
 #include "UI_interface.h"
@@ -64,7 +66,6 @@
 #include "WM_types.h"
 
 #include "uvedit_intern.h"
-#include "uvedit_parametrizer.h"
 
 /* -------------------------------------------------------------------- */
 /** \name Utility Functions
@@ -319,7 +320,7 @@ static void construct_param_handle_face_add(ParamHandle *handle,
     select[i] = uvedit_uv_select_test(scene, l, cd_loop_uv_offset);
   }
 
-  param_face_add(handle, key, i, vkeys, co, uv, pin, select);
+  GEO_uv_parametrizer_face_add(handle, key, i, vkeys, co, uv, pin, select);
 }
 
 /* See: construct_param_handle_multi to handle multiple objects at once. */
@@ -338,7 +339,7 @@ static ParamHandle *construct_param_handle(const Scene *scene,
 
   const int cd_loop_uv_offset = CustomData_get_offset(&bm->ldata, CD_MLOOPUV);
 
-  handle = param_construct_begin();
+  handle = GEO_uv_parametrizer_construct_begin();
 
   if (options->correct_aspect) {
     float aspx, aspy;
@@ -346,7 +347,7 @@ static ParamHandle *construct_param_handle(const Scene *scene,
     ED_uvedit_get_aspect(ob, &aspx, &aspy);
 
     if (aspx != aspy) {
-      param_aspect_ratio(handle, aspx, aspy);
+      GEO_uv_parametrizer_aspect_ratio(handle, aspx, aspy);
     }
   }
 
@@ -385,15 +386,15 @@ static ParamHandle *construct_param_handle(const Scene *scene,
         ParamKey vkeys[2];
         vkeys[0] = (ParamKey)BM_elem_index_get(eed->v1);
         vkeys[1] = (ParamKey)BM_elem_index_get(eed->v2);
-        param_edge_set_seam(handle, vkeys);
+        GEO_uv_parametrizer_edge_set_seam(handle, vkeys);
       }
     }
   }
 
-  param_construct_end(handle,
-                      options->fill_holes,
-                      options->topology_from_uvs,
-                      result_info ? &result_info->count_failed : NULL);
+  GEO_uv_parametrizer_construct_end(handle,
+                                    options->fill_holes,
+                                    options->topology_from_uvs,
+                                    result_info ? &result_info->count_failed : NULL);
 
   return handle;
 }
@@ -414,7 +415,7 @@ static ParamHandle *construct_param_handle_multi(const Scene *scene,
   BMIter iter, liter;
   int i;
 
-  handle = param_construct_begin();
+  handle = GEO_uv_parametrizer_construct_begin();
 
   if (options->correct_aspect) {
     Object *ob = objects[0];
@@ -422,7 +423,7 @@ static ParamHandle *construct_param_handle_multi(const Scene *scene,
 
     ED_uvedit_get_aspect(ob, &aspx, &aspy);
     if (aspx != aspy) {
-      param_aspect_ratio(handle, aspx, aspy);
+      GEO_uv_parametrizer_aspect_ratio(handle, aspx, aspy);
     }
   }
 
@@ -474,14 +475,15 @@ static ParamHandle *construct_param_handle_multi(const Scene *scene,
           ParamKey vkeys[2];
           vkeys[0] = (ParamKey)BM_elem_index_get(eed->v1);
           vkeys[1] = (ParamKey)BM_elem_index_get(eed->v2);
-          param_edge_set_seam(handle, vkeys);
+          GEO_uv_parametrizer_edge_set_seam(handle, vkeys);
         }
       }
     }
     offset += bm->totface;
   }
 
-  param_construct_end(handle, options->fill_holes, options->topology_from_uvs, count_fail);
+  GEO_uv_parametrizer_construct_end(
+      handle, options->fill_holes, options->topology_from_uvs, count_fail);
 
   return handle;
 }
@@ -560,7 +562,7 @@ static ParamHandle *construct_param_handle_subsurfed(const Scene *scene,
 
   const int cd_loop_uv_offset = CustomData_get_offset(&em->bm->ldata, CD_MLOOPUV);
 
-  handle = param_construct_begin();
+  handle = GEO_uv_parametrizer_construct_begin();
 
   if (options->correct_aspect) {
     float aspx, aspy;
@@ -568,7 +570,7 @@ static ParamHandle *construct_param_handle_subsurfed(const Scene *scene,
     ED_uvedit_get_aspect(ob, &aspx, &aspy);
 
     if (aspx != aspy) {
-      param_aspect_ratio(handle, aspx, aspy);
+      GEO_uv_parametrizer_aspect_ratio(handle, aspx, aspy);
     }
   }
 
@@ -689,7 +691,7 @@ static ParamHandle *construct_param_handle_subsurfed(const Scene *scene,
                                 &pin[3],
                                 &select[3]);
 
-    param_face_add(handle, key, 4, vkeys, co, uv, pin, select);
+    GEO_uv_parametrizer_face_add(handle, key, 4, vkeys, co, uv, pin, select);
   }
 
   /* these are calculated from original mesh too */
@@ -698,14 +700,14 @@ static ParamHandle *construct_param_handle_subsurfed(const Scene *scene,
       ParamKey vkeys[2];
       vkeys[0] = (ParamKey)edge->v1;
       vkeys[1] = (ParamKey)edge->v2;
-      param_edge_set_seam(handle, vkeys);
+      GEO_uv_parametrizer_edge_set_seam(handle, vkeys);
     }
   }
 
-  param_construct_end(handle,
-                      options->fill_holes,
-                      options->topology_from_uvs,
-                      result_info ? &result_info->count_failed : NULL);
+  GEO_uv_parametrizer_construct_end(handle,
+                                    options->fill_holes,
+                                    options->topology_from_uvs,
+                                    result_info ? &result_info->count_failed : NULL);
 
   /* cleanup */
   MEM_freeN(faceMap);
@@ -764,9 +766,9 @@ static bool minimize_stretch_init(bContext *C, wmOperator *op)
   ms->handle = construct_param_handle_multi(scene, objects, objects_len, &options, NULL);
   ms->lasttime = PIL_check_seconds_timer();
 
-  param_stretch_begin(ms->handle);
+  GEO_uv_parametrizer_stretch_begin(ms->handle);
   if (ms->blend != 0.0f) {
-    param_stretch_blend(ms->handle, ms->blend);
+    GEO_uv_parametrizer_stretch_blend(ms->handle, ms->blend);
   }
 
   op->customdata = ms;
@@ -782,8 +784,8 @@ static void minimize_stretch_iteration(bContext *C, wmOperator *op, bool interac
   ToolSettings *ts = scene->toolsettings;
   const bool synced_selection = (ts->uv_flag & UV_SYNC_SELECTION) != 0;
 
-  param_stretch_blend(ms->handle, ms->blend);
-  param_stretch_iter(ms->handle);
+  GEO_uv_parametrizer_stretch_blend(ms->handle, ms->blend);
+  GEO_uv_parametrizer_stretch_iter(ms->handle);
 
   ms->i++;
   RNA_int_set(op->ptr, "iterations", ms->i);
@@ -791,7 +793,7 @@ static void minimize_stretch_iteration(bContext *C, wmOperator *op, bool interac
   if (interactive && (PIL_check_seconds_timer() - ms->lasttime > 0.5)) {
     char str[UI_MAX_DRAW_STR];
 
-    param_flush(ms->handle);
+    GEO_uv_parametrizer_flush(ms->handle);
 
     if (area) {
       BLI_snprintf(str, sizeof(str), TIP_("Minimize Stretch. Blend %.2f"), ms->blend);
@@ -831,14 +833,14 @@ static void minimize_stretch_exit(bContext *C, wmOperator *op, bool cancel)
   }
 
   if (cancel) {
-    param_flush_restore(ms->handle);
+    GEO_uv_parametrizer_flush_restore(ms->handle);
   }
   else {
-    param_flush(ms->handle);
+    GEO_uv_parametrizer_flush(ms->handle);
   }
 
-  param_stretch_end(ms->handle);
-  param_delete(ms->handle);
+  GEO_uv_parametrizer_stretch_end(ms->handle);
+  GEO_uv_parametrizer_delete(ms->handle);
 
   for (uint ob_index = 0; ob_index < ms->objects_len; ob_index++) {
     Object *obedit = ms->objects_edit[ob_index];
@@ -1014,9 +1016,9 @@ static void uvedit_pack_islands(const Scene *scene, Object *ob, BMesh *bm)
 
   ParamHandle *handle;
   handle = construct_param_handle(scene, ob, bm, &options, NULL);
-  param_pack(handle, scene->toolsettings->uvcalc_margin, rotate, ignore_pinned);
-  param_flush(handle);
-  param_delete(handle);
+  GEO_uv_parametrizer_pack(handle, scene->toolsettings->uvcalc_margin, rotate, ignore_pinned);
+  GEO_uv_parametrizer_flush(handle);
+  GEO_uv_parametrizer_delete(handle);
 }
 
 /**
@@ -1034,9 +1036,9 @@ static void uvedit_pack_islands_multi(const Scene *scene,
 {
   ParamHandle *handle;
   handle = construct_param_handle_multi(scene, objects, objects_len, options, NULL);
-  param_pack(handle, scene->toolsettings->uvcalc_margin, rotate, ignore_pinned);
-  param_flush(handle);
-  param_delete(handle);
+  GEO_uv_parametrizer_pack(handle, scene->toolsettings->uvcalc_margin, rotate, ignore_pinned);
+  GEO_uv_parametrizer_flush(handle);
+  GEO_uv_parametrizer_delete(handle);
 
   for (uint ob_index = 0; ob_index < objects_len; ob_index++) {
     Object *obedit = objects[ob_index];
@@ -1167,9 +1169,9 @@ static int average_islands_scale_exec(bContext *C, wmOperator *UNUSED(op))
   }
 
   ParamHandle *handle = construct_param_handle_multi(scene, objects, objects_len, &options, NULL);
-  param_average(handle, false);
-  param_flush(handle);
-  param_delete(handle);
+  GEO_uv_parametrizer_average(handle, false);
+  GEO_uv_parametrizer_flush(handle);
+  GEO_uv_parametrizer_delete(handle);
 
   for (uint ob_index = 0; ob_index < objects_len; ob_index++) {
     Object *obedit = objects[ob_index];
@@ -1239,7 +1241,7 @@ void ED_uvedit_live_unwrap_begin(Scene *scene, Object *obedit)
     handle = construct_param_handle(scene, obedit, em->bm, &options, NULL);
   }
 
-  param_lscm_begin(handle, PARAM_TRUE, abf);
+  GEO_uv_parametrizer_lscm_begin(handle, PARAM_TRUE, abf);
 
   /* Create or increase size of g_live_unwrap.handles array */
   if (g_live_unwrap.handles == NULL) {
@@ -1261,8 +1263,8 @@ void ED_uvedit_live_unwrap_re_solve(void)
 {
   if (g_live_unwrap.handles) {
     for (int i = 0; i < g_live_unwrap.len; i++) {
-      param_lscm_solve(g_live_unwrap.handles[i], NULL, NULL);
-      param_flush(g_live_unwrap.handles[i]);
+      GEO_uv_parametrizer_lscm_solve(g_live_unwrap.handles[i], NULL, NULL);
+      GEO_uv_parametrizer_flush(g_live_unwrap.handles[i]);
     }
   }
 }
@@ -1271,11 +1273,11 @@ void ED_uvedit_live_unwrap_end(short cancel)
 {
   if (g_live_unwrap.handles) {
     for (int i = 0; i < g_live_unwrap.len; i++) {
-      param_lscm_end(g_live_unwrap.handles[i]);
+      GEO_uv_parametrizer_lscm_end(g_live_unwrap.handles[i]);
       if (cancel) {
-        param_flush_restore(g_live_unwrap.handles[i]);
+        GEO_uv_parametrizer_flush_restore(g_live_unwrap.handles[i]);
       }
-      param_delete(g_live_unwrap.handles[i]);
+      GEO_uv_parametrizer_delete(g_live_unwrap.handles[i]);
     }
     MEM_freeN(g_live_unwrap.handles);
     g_live_unwrap.handles = NULL;
@@ -1731,17 +1733,17 @@ static void uvedit_unwrap(const Scene *scene,
     handle = construct_param_handle(scene, obedit, em->bm, options, result_info);
   }
 
-  param_lscm_begin(handle, PARAM_FALSE, scene->toolsettings->unwrapper == 0);
-  param_lscm_solve(handle,
-                   result_info ? &result_info->count_changed : NULL,
-                   result_info ? &result_info->count_failed : NULL);
-  param_lscm_end(handle);
+  GEO_uv_parametrizer_lscm_begin(handle, PARAM_FALSE, scene->toolsettings->unwrapper == 0);
+  GEO_uv_parametrizer_lscm_solve(handle,
+                                 result_info ? &result_info->count_changed : NULL,
+                                 result_info ? &result_info->count_failed : NULL);
+  GEO_uv_parametrizer_lscm_end(handle);
 
-  param_average(handle, true);
+  GEO_uv_parametrizer_average(handle, true);
 
-  param_flush(handle);
+  GEO_uv_parametrizer_flush(handle);
 
-  param_delete(handle);
+  GEO_uv_parametrizer_delete(handle);
 }
 
 static void uvedit_unwrap_multi(const Scene *scene,