Cleanup: use '_len' instead of '_size' w/ BLI API
[blender.git] / source / blender / editors / uvedit / uvedit_parametrizer.c
index 0d2f467767a5770af19dabc42bf7dcf4ec061a2c..f8a93bd2ce3b292ba587bb9da06b1f77e2adecdb 100644 (file)
 
 #include "MEM_guardedalloc.h"
 
-#include "BLI_array.h"
+#include "BLI_utildefines.h"
+#include "BLI_alloca.h"
 #include "BLI_memarena.h"
 #include "BLI_math.h"
 #include "BLI_rand.h"
 #include "BLI_heap.h"
 #include "BLI_boxpack2d.h"
-#include "BLI_utildefines.h"
-
+#include "BLI_convexhull2d.h"
 
-
-#include "ONL_opennl.h"
-
-#include "uvedit_intern.h"
 #include "uvedit_parametrizer.h"
 
 #include <math.h>
 #include <stdio.h>
 #include <string.h>
 
-#include "BLO_sys_types.h"  /* for intptr_t support */
+#include "BLI_sys_types.h"  /* for intptr_t support */
+
+#include "eigen_capi.h"
 
 /* Utils */
 
 #if 0
-       #define param_assert(condition)
-       #define param_warning(message)
-       #define param_test_equals_ptr(condition)
-       #define param_test_equals_int(condition)
+#  define param_assert(condition)
+#  define param_warning(message)
+#  define param_test_equals_ptr(condition)
+#  define param_test_equals_int(condition)
 #else
-       #define param_assert(condition) \
+#  define param_assert(condition) \
                if (!(condition)) \
                        { /*printf("Assertion %s:%d\n", __FILE__, __LINE__); abort();*/ } (void)0
-       #define param_warning(message) \
+#  define param_warning(message) \
                { /*printf("Warning %s:%d: %s\n", __FILE__, __LINE__, message);*/ } (void)0
-#if 0
-       #define param_test_equals_ptr(str, a, b) \
+#  if 0
+#    define param_test_equals_ptr(str, a, b) \
                if (a != b) \
                        { /*printf("Equals %s => %p != %p\n", str, a, b);*/ } (void)0
-       #define param_test_equals_int(str, a, b) \
+#    define param_test_equals_int(str, a, b) \
                if (a != b) \
                        { /*printf("Equals %s => %d != %d\n", str, a, b);*/ } (void)0
-#endif
+#  endif
 #endif
 typedef enum PBool {
        P_TRUE = 1,
@@ -193,13 +191,13 @@ typedef struct PChart {
 
        union PChartUnion {
                struct PChartLscm {
-                       NLContext context;
+                       LinearSolver *context;
                        float *abf_alpha;
                        PVert *pin1, *pin2;
                } lscm;
                struct PChartPack {
                        float rescale, area;
-                       float size[2], trans[2];
+                       float size[2] /* , trans[2] */;
                } pack;
        } u;
 
@@ -250,7 +248,7 @@ static int PHashSizes[] = {
 };
 
 #define PHASH_hash(ph, item) (((uintptr_t) (item)) % ((unsigned int) (ph)->cursize))
-#define PHASH_edge(v1, v2)   ((v1) ^ (v2))
+#define PHASH_edge(v1, v2)   (((v1) < (v2)) ? ((v1) * 39) ^ ((v2) * 31) : ((v1) * 31) ^ ((v2) * 39))
 
 static PHash *phash_new(PHashLink **list, int sizehint)
 {
@@ -371,7 +369,7 @@ static float p_vec_angle(float *v1, float *v2, float *v3)
        else if (dot >= 1.0f)
                return 0.0f;
        else
-               return (float)acos(dot);
+               return acosf(dot);
 }
 
 static float p_vec2_angle(float *v1, float *v2, float *v3)
@@ -432,7 +430,7 @@ static float p_edge_length(PEdge *e)
        d[1] = v2->co[1] - v1->co[1];
        d[2] = v2->co[2] - v1->co[2];
 
-       return sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
+       return sqrtf(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
 }
 
 static float p_edge_uv_length(PEdge *e)
@@ -443,7 +441,7 @@ static float p_edge_uv_length(PEdge *e)
        d[0] = v2->uv[0] - v1->uv[0];
        d[1] = v2->uv[1] - v1->uv[1];
 
-       return sqrt(d[0] * d[0] + d[1] * d[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])
@@ -487,6 +485,36 @@ static void p_chart_uv_translate(PChart *chart, float trans[2])
        }
 }
 
+static void p_chart_uv_transform(PChart *chart, float mat[2][2])
+{
+       PVert *v;
+
+       for (v = chart->verts; v; v = v->nextlink) {
+               mul_m2v2(mat, v->uv);
+       }
+}
+
+static void p_chart_uv_to_array(PChart *chart, float (*points)[2])
+{
+       PVert *v;
+       unsigned int 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;
+       unsigned int i = 0;
+
+       for (v = chart->verts; v; v = v->nextlink) {
+               copy_v2_v2(v->uv, points[i++]);
+       }
+}
+
+
 static PBool p_intersect_line_2d_dir(float *v1, float *dir1, float *v2, float *dir2, float *isect)
 {
        float lmbda, div;
@@ -712,6 +740,16 @@ static PVert *p_vert_add(PHandle *handle, PHashKey key, const float co[3], PEdge
 {
        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;
@@ -1069,7 +1107,7 @@ static PFace *p_face_add(PHandle *handle)
 }
 
 static PFace *p_face_add_construct(PHandle *handle, ParamKey key, ParamKey *vkeys,
-                                   float *co[3], float *uv[3], int i1, int i2, int i3,
+                                   float *co[4], float *uv[4], int i1, int i2, int i3,
                                    ParamBool *pin, ParamBool *select)
 {
        PFace *f = p_face_add(handle);
@@ -1253,7 +1291,7 @@ static void p_chart_fill_boundary(PChart *chart, PEdge *be, int nedges)
                while (nedges > 2) {
                        PEdge *ne, *ne1, *ne2;
 
-                       e = (PEdge *)BLI_heap_popmin(heap);
+                       e = (PEdge *)BLI_heap_pop_min(heap);
 
                        e1 = p_boundary_edge_prev(e);
                        e2 = p_boundary_edge_next(e);
@@ -1855,7 +1893,7 @@ static PBool p_collapse_allowed_geometric(PEdge *edge, PEdge *pair)
 
        if (p_vert_interior(oldv)) {
                /* hlscm criterion: angular defect smaller than threshold */
-               if (fabs(angulardefect) > (M_PI * 30.0 / 180.0))
+               if (fabsf(angulardefect) > (float)(M_PI * 30.0 / 180.0))
                        return P_FALSE;
        }
        else {
@@ -1921,7 +1959,7 @@ static float p_collapse_cost(PEdge *edge, PEdge *pair)
                        sub_v3_v3v3(tetrav3, co2, oldv->co);
                        cross_v3_v3v3(c, tetrav2, tetrav3);
 
-                       volumecost += fabs(dot_v3v3(edgevec, c) / 6.0f);
+                       volumecost += fabsf(dot_v3v3(edgevec, c) / 6.0f);
 #if 0
                        shapecost += dot_v3v3(co1, keepv->co);
 
@@ -2147,12 +2185,12 @@ static void p_chart_simplify_compute(PChart *chart)
                e->u.nextcollapse = NULL;
 
        /* pop edge collapse out of heap one by one */
-       while (!BLI_heap_empty(heap)) {
+       while (!BLI_heap_is_empty(heap)) {
                if (ncollapsed == NCOLLAPSE)
                        break;
 
                HeapNode *link = BLI_heap_top(heap);
-               PEdge *edge = (PEdge *)BLI_heap_popmin(heap), *pair = edge->pair;
+               PEdge *edge = (PEdge *)BLI_heap_pop_min(heap), *pair = edge->pair;
                PVert *oldv, *keepv;
                PEdge *wheele, *nexte;
 
@@ -2322,8 +2360,8 @@ static void p_abf_compute_sines(PAbfSystem *sys)
        float *sine = sys->sine, *cosine = sys->cosine, *alpha = sys->alpha;
 
        for (i = 0; i < sys->nangles; i++, sine++, cosine++, alpha++) {
-               *sine = sin(*alpha);
-               *cosine = cos(*alpha);
+               *sine = sinf(*alpha);
+               *cosine = cosf(*alpha);
        }
 }
 
@@ -2441,16 +2479,12 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
        PEdge *e;
        int i, j, ninterior = sys->ninterior, nvar = 2 * sys->ninterior;
        PBool success;
+       LinearSolver *context;
 
-       nlNewContext();
-       nlSolverParameteri(NL_NB_VARIABLES, nvar);
-
-       nlBegin(NL_SYSTEM);
-
-       nlBegin(NL_MATRIX);
+       context = EIG_linear_solver_new(0, nvar, 1);
 
        for (i = 0; i < nvar; i++)
-               nlRightHandSideAdd(0, i, sys->bInterior[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];
@@ -2496,8 +2530,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        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;
 
-                       nlRightHandSideAdd(0, v1->u.id, j2[0][0] * beta[0]);
-                       nlRightHandSideAdd(0, ninterior + v1->u.id, j2[1][0] * beta[1] + j2[2][0] * beta[2]);
+                       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];
@@ -2516,8 +2550,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        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;
 
-                       nlRightHandSideAdd(0, v2->u.id, j2[1][1] * beta[1]);
-                       nlRightHandSideAdd(0, ninterior + v2->u.id, j2[0][1] * beta[0] + j2[2][1] * beta[2]);
+                       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];
@@ -2536,8 +2570,8 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        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;
 
-                       nlRightHandSideAdd(0, v3->u.id, j2[2][2] * beta[2]);
-                       nlRightHandSideAdd(0, ninterior + v3->u.id, j2[0][2] * beta[0] + j2[1][2] * beta[1]);
+                       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];
@@ -2561,29 +2595,25 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                                        continue;
 
                                if (i == 0)
-                                       nlMatrixAdd(r, c, j2[0][i] * row1[j]);
+                                       EIG_linear_solver_matrix_add(context, r, c, j2[0][i] * row1[j]);
                                else
-                                       nlMatrixAdd(r + ninterior, c, j2[0][i] * row1[j]);
+                                       EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[0][i] * row1[j]);
 
                                if (i == 1)
-                                       nlMatrixAdd(r, c, j2[1][i] * row2[j]);
+                                       EIG_linear_solver_matrix_add(context, r, c, j2[1][i] * row2[j]);
                                else
-                                       nlMatrixAdd(r + ninterior, c, j2[1][i] * row2[j]);
+                                       EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[1][i] * row2[j]);
 
 
                                if (i == 2)
-                                       nlMatrixAdd(r, c, j2[2][i] * row3[j]);
+                                       EIG_linear_solver_matrix_add(context, r, c, j2[2][i] * row3[j]);
                                else
-                                       nlMatrixAdd(r + ninterior, c, j2[2][i] * row3[j]);
+                                       EIG_linear_solver_matrix_add(context, r + ninterior, c, j2[2][i] * row3[j]);
                        }
                }
        }
 
-       nlEnd(NL_MATRIX);
-
-       nlEnd(NL_SYSTEM);
-
-       success = nlSolve();
+       success = EIG_linear_solver_solve(context);
 
        if (success) {
                for (f = chart->faces; f; f = f->nextlink) {
@@ -2594,24 +2624,24 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                        pre[0] = pre[1] = pre[2] = 0.0;
 
                        if (v1->flag & PVERT_INTERIOR) {
-                               float x = nlGetVariable(0, v1->u.id);
-                               float x2 = nlGetVariable(0, ninterior + v1->u.id);
+                               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 = nlGetVariable(0, v2->u.id);
-                               float x2 = nlGetVariable(0, ninterior + v2->u.id);
+                               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 = nlGetVariable(0, v3->u.id);
-                               float x2 = nlGetVariable(0, ninterior + v3->u.id);
+                               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;
@@ -2642,12 +2672,12 @@ static PBool p_abf_matrix_invert(PAbfSystem *sys, PChart *chart)
                }
 
                for (i = 0; i < ninterior; i++) {
-                       sys->lambdaPlanar[i] += nlGetVariable(0, i);
-                       sys->lambdaLength[i] += nlGetVariable(0, 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);
                }
        }
 
-       nlDeleteContext(nlGetCurrent());
+       EIG_linear_solver_delete(context);
 
        return success;
 }
@@ -2774,7 +2804,7 @@ static PBool p_chart_abf_solve(PChart *chart)
 
 static void p_chart_pin_positions(PChart *chart, PVert **pin1, PVert **pin2)
 {
-       if (pin1 == pin2) {
+       if (!*pin1 || !*pin2 || *pin1 == *pin2) {
                /* degenerate case */
                PFace *f = chart->faces;
                *pin1 = f->edge->vert;
@@ -2790,9 +2820,9 @@ static void p_chart_pin_positions(PChart *chart, PVert **pin1, PVert **pin2)
                float sub[3];
 
                sub_v3_v3v3(sub, (*pin1)->co, (*pin2)->co);
-               sub[0] = fabs(sub[0]);
-               sub[1] = fabs(sub[1]);
-               sub[2] = fabs(sub[2]);
+               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;
@@ -2829,7 +2859,7 @@ static PBool p_chart_symmetry_pins(PChart *chart, PEdge *outer, PVert **pin1, PV
        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;
@@ -2925,7 +2955,7 @@ static PBool p_chart_symmetry_pins(PChart *chart, PEdge *outer, PVert **pin1, PV
 
        p_chart_pin_positions(chart, pin1, pin2);
 
-       return P_TRUE;
+       return !equals_v3v3((*pin1)->co, (*pin2)->co);
 }
 
 static void p_chart_extrema_verts(PChart *chart, PVert **pin1, PVert **pin2)
@@ -2973,11 +3003,12 @@ static void p_chart_extrema_verts(PChart *chart, PVert **pin1, PVert **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] = nlGetVariable(0, 2 * v->u.id);
-               v->uv[1] = nlGetVariable(0, 2 * v->u.id + 1);
+               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);
        }
 }
 
@@ -3019,8 +3050,10 @@ static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
 
                        p_chart_boundaries(chart, NULL, &outer);
 
-                       if (!p_chart_symmetry_pins(chart, outer, &pin1, &pin2))
+                       /* 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;
@@ -3032,17 +3065,13 @@ static void p_chart_lscm_begin(PChart *chart, PBool live, PBool abf)
                for (v = chart->verts; v; v = v->nextlink)
                        v->u.id = id++;
 
-               nlNewContext();
-               nlSolverParameteri(NL_NB_VARIABLES, 2 * chart->nverts);
-               nlSolverParameteri(NL_NB_ROWS, 2 * chart->nfaces);
-               nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
-
-               chart->u.lscm.context = nlGetCurrent();
+               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;
        float *alpha = chart->u.lscm.abf_alpha;
@@ -3050,10 +3079,6 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
        bool flip_faces;
        int row;
 
-       nlMakeCurrent(chart->u.lscm.context);
-
-       nlBegin(NL_SYSTEM);
-
 #if 0
        /* TODO: make loading pins work for simplify/complexify. */
 #endif
@@ -3063,25 +3088,25 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
                        p_vert_load_pin_select_uvs(handle, v);  /* reload for live */
 
        if (chart->u.lscm.pin1) {
-               nlLockVariable(2 * pin1->u.id);
-               nlLockVariable(2 * pin1->u.id + 1);
-               nlLockVariable(2 * pin2->u.id);
-               nlLockVariable(2 * pin2->u.id + 1);
+               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);
 
-               nlSetVariable(0, 2 * pin1->u.id, pin1->uv[0]);
-               nlSetVariable(0, 2 * pin1->u.id + 1, pin1->uv[1]);
-               nlSetVariable(0, 2 * pin2->u.id, pin2->uv[0]);
-               nlSetVariable(0, 2 * pin2->u.id + 1, pin2->uv[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) {
-                               nlLockVariable(2 * v->u.id);
-                               nlLockVariable(2 * v->u.id + 1);
+                               EIG_linear_solver_variable_lock(context, 2 * v->u.id);
+                               EIG_linear_solver_variable_lock(context, 2 * v->u.id + 1);
 
-                               nlSetVariable(0, 2 * v->u.id, v->uv[0]);
-                               nlSetVariable(0, 2 * v->u.id + 1, v->uv[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]);
                        }
                }
        }
@@ -3108,8 +3133,6 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 
        /* construct matrix */
 
-       nlBegin(NL_MATRIX);
-
        row = 0;
        for (f = chart->faces; f; f = f->nextlink) {
                PEdge *e1 = f->edge, *e2 = e1->next, *e3 = e2->next;
@@ -3132,9 +3155,9 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
                        SWAP(PVert *, v2, v3);
                }
 
-               sina1 = sin(a1);
-               sina2 = sin(a2);
-               sina3 = sin(a3);
+               sina1 = sinf(a1);
+               sina2 = sinf(a2);
+               sina3 = sinf(a3);
 
                sinmax = max_fff(sina1, sina2, sina3);
 
@@ -3156,44 +3179,22 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
                cosine = cosf(a1) * ratio;
                sine = sina1 * ratio;
 
-#if 0
-               nlBegin(NL_ROW);
-               nlCoefficient(2 * v1->u.id,   cosine - 1.0);
-               nlCoefficient(2 * v1->u.id + 1, -sine);
-               nlCoefficient(2 * v2->u.id,   -cosine);
-               nlCoefficient(2 * v2->u.id + 1, sine);
-               nlCoefficient(2 * v3->u.id,   1.0);
-               nlEnd(NL_ROW);
-
-               nlBegin(NL_ROW);
-               nlCoefficient(2 * v1->u.id,   sine);
-               nlCoefficient(2 * v1->u.id + 1, cosine - 1.0);
-               nlCoefficient(2 * v2->u.id,   -sine);
-               nlCoefficient(2 * v2->u.id + 1, -cosine);
-               nlCoefficient(2 * v3->u.id + 1, 1.0);
-               nlEnd(NL_ROW);
-#else
-               nlMatrixAdd(row, 2 * v1->u.id,   cosine - 1.0f);
-               nlMatrixAdd(row, 2 * v1->u.id + 1, -sine);
-               nlMatrixAdd(row, 2 * v2->u.id,   -cosine);
-               nlMatrixAdd(row, 2 * v2->u.id + 1, sine);
-               nlMatrixAdd(row, 2 * v3->u.id,   1.0);
+               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++;
 
-               nlMatrixAdd(row, 2 * v1->u.id,   sine);
-               nlMatrixAdd(row, 2 * v1->u.id + 1, cosine - 1.0f);
-               nlMatrixAdd(row, 2 * v2->u.id,   -sine);
-               nlMatrixAdd(row, 2 * v2->u.id + 1, -cosine);
-               nlMatrixAdd(row, 2 * v3->u.id + 1, 1.0);
+               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++;
-#endif
        }
 
-       nlEnd(NL_MATRIX);
-
-       nlEnd(NL_SYSTEM);
-
-       if (nlSolveAdvanced(NULL, NL_TRUE)) {
+       if (EIG_linear_solver_solve(context)) {
                p_chart_lscm_load_solution(chart);
                return P_TRUE;
        }
@@ -3210,7 +3211,7 @@ static PBool p_chart_lscm_solve(PHandle *handle, PChart *chart)
 static void p_chart_lscm_end(PChart *chart)
 {
        if (chart->u.lscm.context)
-               nlDeleteContext(chart->u.lscm.context);
+               EIG_linear_solver_delete(chart->u.lscm.context);
        
        if (chart->u.lscm.abf_alpha) {
                MEM_freeN(chart->u.lscm.abf_alpha);
@@ -3283,7 +3284,7 @@ static float p_face_stretch(PFace *f)
        a = dot_v3v3(Ps, Ps);
        c = dot_v3v3(Pt, Pt);
 
-       T =  sqrt(0.5f * (a + c));
+       T =  sqrtf(0.5f * (a + c));
        if (f->flag & PFACE_FILLED)
                T *= 0.2f;
 
@@ -3372,8 +3373,8 @@ static void p_chart_stretch_minimize(PChart *chart, RNG *rng)
 
 static int p_compare_geometric_uv(const void *a, const void *b)
 {
-       PVert *v1 = *(PVert **)a;
-       PVert *v2 = *(PVert **)b;
+       const PVert *v1 = *(const PVert * const *)a;
+       const PVert *v2 = *(const PVert * const *)b;
 
        if (v1->uv[0] < v2->uv[0])
                return -1;
@@ -3599,8 +3600,8 @@ static float p_chart_minimum_area_angle(PChart *chart)
 static void p_chart_rotate_minimum_area(PChart *chart)
 {
        float angle = p_chart_minimum_area_angle(chart);
-       float sine = sin(angle);
-       float cosine = cos(angle);
+       float sine = sinf(angle);
+       float cosine = cosf(angle);
        PVert *v;
 
        for (v = chart->verts; v; v = v->nextlink) {
@@ -3682,8 +3683,8 @@ static SmoothNode *p_node_new(MemArena *arena, SmoothTriangle **tri, int ntri, f
        if (ntri <= 10 || depth >= 15)
                return node;
        
-       t1 = MEM_mallocN(sizeof(SmoothTriangle) * ntri, "PNodeTri1");
-       t2 = MEM_mallocN(sizeof(SmoothTriangle) * ntri, "PNodeTri1");
+       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]);
@@ -3758,11 +3759,14 @@ static PBool p_node_intersect(SmoothNode *node, float co[2])
 
 /* smoothing */
 
-static int p_compare_float(const void *a, const void *b)
+static int p_compare_float(const void *a_, const void *b_)
 {
-       if (*((float *)a) < *((float *)b))
+       const float a = *(const float *)a_;
+       const float b = *(const float *)b_;
+
+       if (a < b)
                return -1;
-       else if (*((float *)a) == *((float *)b))
+       else if (a == b)
                return 0;
        else
                return 1;
@@ -4014,7 +4018,7 @@ static void p_smooth(PChart *chart)
                                        diff[0] = p[0] - oldp[0];
                                        diff[1] = p[1] - oldp[1];
 
-                                       length = sqrt(diff[0] * diff[0] + diff[1] * diff[1]);
+                                       length = len_v2(diff);
                                        d = max_ff(d, length);
                                        moved += length;
                                }
@@ -4094,7 +4098,7 @@ static void p_smooth(PChart *chart)
        MEM_freeN(nodesx);
        MEM_freeN(nodesy);
 
-       arena = BLI_memarena_new(1 << 16, "param smooth arena");
+       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)
@@ -4114,10 +4118,10 @@ 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((1 << 16), "param construct arena");
+       handle->arena = BLI_memarena_new(MEM_SIZE_OPTIMAL(1 << 16), "param construct arena");
        handle->aspx = 1.0f;
        handle->aspy = 1.0f;
-       handle->do_aspect = FALSE;
+       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);
@@ -4132,7 +4136,7 @@ void param_aspect_ratio(ParamHandle *handle, float aspx, float aspy)
 
        phandle->aspx = aspx;
        phandle->aspy = aspy;
-       phandle->do_aspect = TRUE;
+       phandle->do_aspect = true;
 }
 
 void param_delete(ParamHandle *handle)
@@ -4163,17 +4167,18 @@ void param_delete(ParamHandle *handle)
 
 static void p_add_ngon(ParamHandle *handle, ParamKey key, int nverts,
                        ParamKey *vkeys, float **co, float **uv,
-                       ParamBool *pin, ParamBool *select, float normal[3])
+                       ParamBool *pin, ParamBool *select, const float normal[3])
 {
        int *boundary = BLI_array_alloca(boundary, nverts);
-       int i;
 
        /* boundary vertex indexes */
-       for (i = 0; i < nverts; i++)
+       for (int i = 0; i < nverts; i++) {
                boundary[i] = i;
+       }
 
        while (nverts > 2) {
                float minangle = FLT_MAX;
+               float minshape = FLT_MAX;
                int i, mini = 0;
 
                /* find corner with smallest angle */
@@ -4189,8 +4194,20 @@ static void p_add_ngon(ParamHandle *handle, ParamKey key, int nverts,
                        if (normal && (dot_v3v3(n, normal) < 0.0f))
                                angle = (float)(2.0 * M_PI) - angle;
 
-                       if (angle < minangle) {
+                       float other_angle = p_vec_angle(co[v2], co[v0], co[v1]);
+                       float shape = fabsf((float)M_PI - angle - 2.0f * other_angle);
+
+                       if (fabsf(angle - minangle) < 0.01f) {
+                               /* for nearly equal angles, try to get well shaped triangles */
+                               if (shape < minshape) {
+                                       minangle = angle;
+                                       minshape = shape;
+                                       mini = i;
+                               }
+                       }
+                       else if (angle < minangle) {
                                minangle = angle;
+                               minshape = shape;
                                mini = i;
                        }
                }
@@ -4219,7 +4236,7 @@ static void p_add_ngon(ParamHandle *handle, ParamKey key, int nverts,
 }
 
 void param_face_add(ParamHandle *handle, ParamKey key, int nverts,
-                    ParamKey *vkeys, float **co, float **uv,
+                    ParamKey *vkeys, float *co[4], float *uv[4],
                     ParamBool *pin, ParamBool *select, float normal[3])
 {
        PHandle *phandle = (PHandle *)handle;
@@ -4443,8 +4460,42 @@ void param_smooth_area(ParamHandle *handle)
                p_smooth(chart);
        }
 }
-void param_pack(ParamHandle *handle, float margin)
+
+/* don't pack, just rotate (used for better packing) */
+static void param_pack_rotate(ParamHandle *handle)
+{
+       PChart *chart;
+       int i;
+
+       PHandle *phandle = (PHandle *)handle;
+
+       for (i = 0; i < phandle->ncharts; i++) {
+               float (*points)[2];
+               float angle;
+
+               chart = phandle->charts[i];
+
+               if (chart->flag & PCHART_NOPACK) {
+                       continue;
+               }
+
+               points = MEM_mallocN(sizeof(*points) * chart->nverts, __func__);
+
+               p_chart_uv_to_array(chart, points);
+
+               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);
+               }
+       }
+}
+
+void param_pack(ParamHandle *handle, float margin, bool do_rotate)
 {      
        /* box packing variables */
        BoxPack *boxarray, *box;
@@ -4463,6 +4514,11 @@ void param_pack(ParamHandle *handle, float margin)
        if (phandle->aspx != phandle->aspy)
                param_scale(handle, 1.0f / phandle->aspx, 1.0f / phandle->aspy);
        
+       /* this could be its own function */
+       if (do_rotate) {
+               param_pack_rotate(handle);
+       }
+
        /* we may not use all these boxes */
        boxarray = MEM_mallocN(phandle->ncharts * sizeof(BoxPack), "BoxPack box");
        
@@ -4489,7 +4545,7 @@ void param_pack(ParamHandle *handle, float margin)
                box->index = i; /* warning this index skips PCHART_NOPACK boxes */
                
                if (margin > 0.0f)
-                       area += sqrt(box->w * box->h);
+                       area += (double)sqrtf(box->w * box->h);
        }
        
        if (margin > 0.0f) {
@@ -4515,7 +4571,7 @@ void param_pack(ParamHandle *handle, float margin)
                }
        }
        
-       BLI_box_pack_2D(boxarray, phandle->ncharts - unpacked, &tot_width, &tot_height);
+       BLI_box_pack_2d(boxarray, phandle->ncharts - unpacked, &tot_width, &tot_height);
        
        if (tot_height > tot_width)
                scale = 1.0f / tot_height;
@@ -4591,7 +4647,7 @@ void param_average(ParamHandle *handle)
                        
                        /* Move center to 0,0 */
                        p_chart_uv_translate(chart, trans);
-                       p_chart_uv_scale(chart, sqrt(fac / tot_fac));
+                       p_chart_uv_scale(chart, sqrtf(fac / tot_fac));
                        
                        /* Move to original center */
                        trans[0] = -trans[0];