Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat (limited to 'source/blender/bmesh/operators/bmo_smooth_laplacian.c')
-rw-r--r--source/blender/bmesh/operators/bmo_smooth_laplacian.c619
1 files changed, 619 insertions, 0 deletions
diff --git a/source/blender/bmesh/operators/bmo_smooth_laplacian.c b/source/blender/bmesh/operators/bmo_smooth_laplacian.c
new file mode 100644
index 00000000000..362123e412d
--- /dev/null
+++ b/source/blender/bmesh/operators/bmo_smooth_laplacian.c
@@ -0,0 +1,619 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * Contributor(s): Alexander Pinzon
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+/** \file blender/bmesh/operators/bmo_smooth_laplacian.c
+ * \ingroup bmesh
+ */
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_meshdata_types.h"
+
+#include "BLI_array.h"
+#include "BLI_heap.h"
+#include "BLI_math.h"
+#include "BLI_math_geom.h"
+#include "BLI_smallhash.h"
+
+#include "BKE_customdata.h"
+#include "BKE_mesh.h"
+
+#include "bmesh.h"
+
+#include "ONL_opennl.h"
+
+#include "intern/bmesh_operators_private.h" /* own include */
+
+// #define SMOOTH_LAPLACIAN_AREA_FACTOR 4.0f /* UNUSED */
+// #define SMOOTH_LAPLACIAN_EDGE_FACTOR 2.0f /* UNUSED */
+#define SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE 1.8f
+#define SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE 0.15f
+
+struct BLaplacianSystem {
+ float *eweights; /* Length weights per Edge */
+ float (*fweights)[3]; /* Cotangent weights per face */
+ float *ring_areas; /* Total area per ring*/
+ float *vlengths; /* Total sum of lengths(edges) per vertice*/
+ float *vweights; /* Total sum of weights per vertice*/
+ int numEdges; /* Number of edges*/
+ int numFaces; /* Number of faces*/
+ int numVerts; /* Number of verts*/
+ short *zerola; /* Is zero area or length*/
+
+ /* Pointers to data*/
+ BMesh *bm;
+ BMOperator *op;
+ NLContext *context;
+
+ /*Data*/
+ float min_area;
+};
+typedef struct BLaplacianSystem LaplacianSystem;
+
+static float compute_volume(BMesh *bm);
+static float cotan_weight(float *v1, float *v2, float *v3);
+static int vert_is_boundary(BMVert *v);
+static LaplacianSystem *init_laplacian_system(int a_numEdges, int a_numFaces, int a_numVerts);
+static void init_laplacian_matrix(LaplacianSystem *sys);
+static void delete_laplacian_system(LaplacianSystem *sys);
+static void delete_void_pointer(void *data);
+static void fill_laplacian_matrix(LaplacianSystem *sys);
+static void memset_laplacian_system(LaplacianSystem *sys, int val);
+static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez, int volumepreservation);
+static void volume_preservation(BMesh *bm, BMOperator *op, float vini, float vend, int usex, int usey, int usez);
+
+static void delete_void_pointer(void *data)
+{
+ if (data) {
+ MEM_freeN(data);
+ data = NULL;
+ }
+}
+
+static void delete_laplacian_system(LaplacianSystem *sys)
+{
+ delete_void_pointer(sys->eweights);
+ delete_void_pointer(sys->fweights);
+ delete_void_pointer(sys->ring_areas);
+ delete_void_pointer(sys->vlengths);
+ delete_void_pointer(sys->vweights);
+ delete_void_pointer(sys->zerola);
+ if (sys->context) {
+ nlDeleteContext(sys->context);
+ }
+ sys->bm = NULL;
+ sys->op = NULL;
+ MEM_freeN(sys);
+}
+
+static void memset_laplacian_system(LaplacianSystem *sys, int val)
+{
+ memset(sys->eweights, val, sizeof(float) * sys->numEdges);
+ memset(sys->fweights, val, sizeof(float) * sys->numFaces * 3);
+ memset(sys->ring_areas, val, sizeof(float) * sys->numVerts);
+ memset(sys->vlengths, val, sizeof(float) * sys->numVerts);
+ memset(sys->vweights, val, sizeof(float) * sys->numVerts);
+ memset(sys->zerola, val, sizeof(short) * sys->numVerts);
+}
+
+static LaplacianSystem *init_laplacian_system(int a_numEdges, int a_numFaces, int a_numVerts)
+{
+ LaplacianSystem *sys;
+ sys = MEM_callocN(sizeof(LaplacianSystem), "ModLaplSmoothSystem");
+ sys->numEdges = a_numEdges;
+ sys->numFaces = a_numFaces;
+ sys->numVerts = a_numVerts;
+
+ sys->eweights = MEM_callocN(sizeof(float) * sys->numEdges, "ModLaplSmoothEWeight");
+ if (!sys->eweights) {
+ delete_laplacian_system(sys);
+ return NULL;
+ }
+
+ sys->fweights = MEM_callocN(sizeof(float) * 3 * sys->numFaces, "ModLaplSmoothFWeight");
+ if (!sys->fweights) {
+ delete_laplacian_system(sys);
+ return NULL;
+ }
+
+ sys->ring_areas = MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothRingAreas");
+ if (!sys->ring_areas) {
+ delete_laplacian_system(sys);
+ return NULL;
+ }
+
+ sys->vlengths = MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVlengths");
+ if (!sys->vlengths) {
+ delete_laplacian_system(sys);
+ return NULL;
+ }
+
+ sys->vweights = MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVweights");
+ if (!sys->vweights) {
+ delete_laplacian_system(sys);
+ return NULL;
+ }
+
+ sys->zerola = MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothZeloa");
+ if (!sys->zerola) {
+ delete_laplacian_system(sys);
+ return NULL;
+ }
+
+ return sys;
+}
+
+/* Compute weigth between vertice v_i and all your neighbors
+ * weight between v_i and v_neighbor
+ * Wij = cot(alpha) + cot(beta) / (4.0 * total area of all faces * sum all weight)
+ * v_i *
+ * / | \
+ * / | \
+ * v_beta* | * v_alpha
+ * \ | /
+ * \ | /
+ * * v_neighbor
+ */
+
+static void init_laplacian_matrix(LaplacianSystem *sys)
+{
+ float areaf;
+ float *v1, *v2, *v3, *v4;
+ float w1, w2, w3, w4;
+ int i, j;
+ int has_4_vert;
+ unsigned int idv1, idv2, idv3, idv4, idv[4];
+ BMEdge *e;
+ BMFace *f;
+ BMIter eiter;
+ BMIter fiter;
+ BMIter vi;
+ BMVert *vn;
+ BMVert *vf[4];
+
+ BM_ITER_MESH_INDEX (e, &eiter, sys->bm, BM_EDGES_OF_MESH, j) {
+ if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e)) {
+ v1 = e->v1->co;
+ v2 = e->v2->co;
+ idv1 = BM_elem_index_get(e->v1);
+ idv2 = BM_elem_index_get(e->v2);
+
+ w1 = len_v3v3(v1, v2);
+ if (w1 > sys->min_area) {
+ w1 = 1.0f / w1;
+ i = BM_elem_index_get(e);
+ sys->eweights[i] = w1;
+ sys->vlengths[idv1] += w1;
+ sys->vlengths[idv2] += w1;
+ }
+ else {
+ sys->zerola[idv1] = 1;
+ sys->zerola[idv2] = 1;
+ }
+ }
+ }
+
+ BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) {
+ if (BM_elem_flag_test(f, BM_ELEM_SELECT)) {
+
+ BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) {
+ vf[i] = vn;
+ }
+ has_4_vert = (i == 4) ? 1 : 0;
+ idv1 = BM_elem_index_get(vf[0]);
+ idv2 = BM_elem_index_get(vf[1]);
+ idv3 = BM_elem_index_get(vf[2]);
+ idv4 = has_4_vert ? BM_elem_index_get(vf[3]) : 0;
+
+ v1 = vf[0]->co;
+ v2 = vf[1]->co;
+ v3 = vf[2]->co;
+ v4 = has_4_vert ? vf[3]->co : 0;
+
+ if (has_4_vert) {
+ areaf = area_quad_v3(v1, v2, v3, v4);
+ }
+ else {
+ areaf = area_tri_v3(v1, v2, v3);
+ }
+
+ if (fabsf(areaf) < sys->min_area) {
+ sys->zerola[idv1] = 1;
+ sys->zerola[idv2] = 1;
+ sys->zerola[idv3] = 1;
+ if (has_4_vert) sys->zerola[idv4] = 1;
+ }
+
+ sys->ring_areas[idv1] += areaf;
+ sys->ring_areas[idv2] += areaf;
+ sys->ring_areas[idv3] += areaf;
+ if (has_4_vert) sys->ring_areas[idv4] += areaf;
+
+ if (has_4_vert) {
+
+ idv[0] = idv1;
+ idv[1] = idv2;
+ idv[2] = idv3;
+ idv[3] = idv4;
+
+ for (j = 0; j < 4; j++) {
+ idv1 = idv[j];
+ idv2 = idv[(j + 1) % 4];
+ idv3 = idv[(j + 2) % 4];
+ idv4 = idv[(j + 3) % 4];
+
+ v1 = vf[j]->co;
+ v2 = vf[(j + 1) % 4]->co;
+ v3 = vf[(j + 2) % 4]->co;
+ v4 = vf[(j + 3) % 4]->co;
+
+ w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
+ w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
+ w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
+
+ sys->vweights[idv1] += (w2 + w3 + w4) / 4.0f;
+ }
+ }
+ else {
+ i = BM_elem_index_get(f);
+
+ w1 = cotan_weight(v1, v2, v3);
+ w2 = cotan_weight(v2, v3, v1);
+ w3 = cotan_weight(v3, v1, v2);
+
+ sys->fweights[i][0] += w1;
+ sys->fweights[i][1] += w2;
+ sys->fweights[i][2] += w3;
+
+ sys->vweights[idv1] += w2 + w3;
+ sys->vweights[idv2] += w1 + w3;
+ sys->vweights[idv3] += w1 + w2;
+ }
+ }
+ }
+}
+
+static void fill_laplacian_matrix(LaplacianSystem *sys)
+{
+ float *v1, *v2, *v3, *v4;
+ float w2, w3, w4;
+ int i, j;
+ int has_4_vert;
+ unsigned int idv1, idv2, idv3, idv4, idv[4];
+
+ BMEdge *e;
+ BMFace *f;
+ BMIter eiter;
+ BMIter fiter;
+ BMIter vi;
+ BMVert *vn;
+ BMVert *vf[4];
+
+ BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) {
+ if (BM_elem_flag_test(f, BM_ELEM_SELECT)) {
+ BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) {
+ vf[i] = vn;
+ }
+ has_4_vert = (i == 4) ? 1 : 0;
+ if (has_4_vert) {
+ idv[0] = BM_elem_index_get(vf[0]);
+ idv[1] = BM_elem_index_get(vf[1]);
+ idv[2] = BM_elem_index_get(vf[2]);
+ idv[3] = BM_elem_index_get(vf[3]);
+ for (j = 0; j < 4; j++) {
+ idv1 = idv[j];
+ idv2 = idv[(j + 1) % 4];
+ idv3 = idv[(j + 2) % 4];
+ idv4 = idv[(j + 3) % 4];
+
+ v1 = vf[j]->co;
+ v2 = vf[(j + 1) % 4]->co;
+ v3 = vf[(j + 2) % 4]->co;
+ v4 = vf[(j + 3) % 4]->co;
+
+ w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
+ w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
+ w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
+
+ w2 = w2 / 4.0f;
+ w3 = w3 / 4.0f;
+ w4 = w4 / 4.0f;
+
+ if (!vert_is_boundary(vf[j]) && sys->zerola[idv1] == 0) {
+ nlMatrixAdd(idv1, idv2, w2 * sys->vweights[idv1]);
+ nlMatrixAdd(idv1, idv3, w3 * sys->vweights[idv1]);
+ nlMatrixAdd(idv1, idv4, w4 * sys->vweights[idv1]);
+ }
+ }
+ }
+ else {
+ idv1 = BM_elem_index_get(vf[0]);
+ idv2 = BM_elem_index_get(vf[1]);
+ idv3 = BM_elem_index_get(vf[2]);
+ /* Is ring if number of faces == number of edges around vertice*/
+ i = BM_elem_index_get(f);
+ if (!vert_is_boundary(vf[0]) && sys->zerola[idv1] == 0) {
+ nlMatrixAdd(idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
+ nlMatrixAdd(idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
+ }
+ if (!vert_is_boundary(vf[1]) && sys->zerola[idv2] == 0) {
+ nlMatrixAdd(idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
+ nlMatrixAdd(idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
+ }
+ if (!vert_is_boundary(vf[2]) && sys->zerola[idv3] == 0) {
+ nlMatrixAdd(idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
+ nlMatrixAdd(idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
+ }
+ }
+ }
+ }
+ BM_ITER_MESH (e, &eiter, sys->bm, BM_EDGES_OF_MESH) {
+ if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e) ) {
+ v1 = e->v1->co;
+ v2 = e->v2->co;
+ idv1 = BM_elem_index_get(e->v1);
+ idv2 = BM_elem_index_get(e->v2);
+ if (sys->zerola[idv1] == 0 && sys->zerola[idv2] == 0) {
+ i = BM_elem_index_get(e);
+ nlMatrixAdd(idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
+ nlMatrixAdd(idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
+ }
+ }
+ }
+}
+
+static float cotan_weight(float *v1, float *v2, float *v3)
+{
+ 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 int vert_is_boundary(BMVert *v)
+{
+ BMEdge *ed;
+ BMFace *f;
+ BMIter ei;
+ BMIter fi;
+ BM_ITER_ELEM (ed, &ei, v, BM_EDGES_OF_VERT) {
+ if (BM_edge_is_boundary(ed)) {
+ return 1;
+ }
+ }
+ BM_ITER_ELEM (f, &fi, v, BM_FACES_OF_VERT) {
+ if (!BM_elem_flag_test(f, BM_ELEM_SELECT)) {
+ return 1;
+ }
+ }
+ return 0;
+}
+
+static float compute_volume(BMesh *bm)
+{
+ float vol = 0.0f;
+ float x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
+ int i;
+ BMFace *f;
+ BMIter fiter;
+ BMIter vi;
+ BMVert *vn;
+ BMVert *vf[4];
+
+ BM_ITER_MESH (f, &fiter, bm, BM_FACES_OF_MESH) {
+ BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) {
+ vf[i] = vn;
+ }
+ x1 = vf[0]->co[0];
+ y1 = vf[0]->co[1];
+ z1 = vf[0]->co[2];
+
+ x2 = vf[1]->co[0];
+ y2 = vf[1]->co[1];
+ z2 = vf[1]->co[2];
+
+ x3 = vf[2]->co[0];
+ y3 = vf[2]->co[1];
+ z3 = vf[2]->co[2];
+
+ vol += (1.0f / 6.0f) * (0.0f - x3 * y2 * z1 + x2 * y3 * z1 + x3 * y1 * z2 - x1 * y3 * z2 - x2 * y1 * z3 + x1 * y2 * z3);
+
+ if (i == 4) {
+ x4 = vf[3]->co[0];
+ y4 = vf[3]->co[1];
+ z4 = vf[3]->co[2];
+ vol += (1.0f / 6.0f) * (x1 * y3 * z4 - x1 * y4 * z3 - x3 * y1 * z4 + x3 * z1 * y4 + y1 * x4 * z3 - x4 * y3 * z1);
+ }
+ }
+ return fabs(vol);
+}
+
+static void volume_preservation(BMesh *bm, BMOperator *op, float vini, float vend, int usex, int usey, int usez)
+{
+ float beta;
+ BMOIter siter;
+ BMVert *v;
+
+ if (vend != 0.0f) {
+ beta = pow(vini / vend, 1.0f / 3.0f);
+ BMO_ITER (v, &siter, bm, op, "verts", BM_VERT) {
+ if (usex) {
+ v->co[0] *= beta;
+ }
+ if (usey) {
+ v->co[1] *= beta;
+ }
+ if (usez) {
+ v->co[2] *= beta;
+ }
+
+ }
+ }
+}
+
+static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez, int volumepreservation)
+{
+ int m_vertex_id;
+ float leni, lene;
+ float vini, vend;
+ float *vi1, *vi2, ve1[3], ve2[3];
+ unsigned int idv1, idv2;
+ BMOIter siter;
+ BMVert *v;
+ BMEdge *e;
+ BMIter eiter;
+
+ BM_ITER_MESH (e, &eiter, sys->bm, BM_EDGES_OF_MESH) {
+ idv1 = BM_elem_index_get(e->v1);
+ idv2 = BM_elem_index_get(e->v2);
+ vi1 = e->v1->co;
+ vi2 = e->v2->co;
+ ve1[0] = nlGetVariable(0, idv1);
+ ve1[1] = nlGetVariable(1, idv1);
+ ve1[2] = nlGetVariable(2, idv1);
+ ve2[0] = nlGetVariable(0, idv2);
+ ve2[1] = nlGetVariable(1, idv2);
+ ve2[2] = nlGetVariable(2, idv2);
+ leni = len_v3v3(vi1, vi2);
+ lene = len_v3v3(ve1, ve2);
+ if (lene > leni * SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE || lene < leni * SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE) {
+ sys->zerola[idv1] = 1;
+ sys->zerola[idv2] = 1;
+ }
+ }
+
+ if (volumepreservation) {
+ vini = compute_volume(sys->bm);
+ }
+ BMO_ITER (v, &siter, sys->bm, sys->op, "verts", BM_VERT) {
+ m_vertex_id = BM_elem_index_get(v);
+ if (sys->zerola[m_vertex_id] == 0) {
+ if (usex) {
+ v->co[0] = nlGetVariable(0, m_vertex_id);
+ }
+ if (usey) {
+ v->co[1] = nlGetVariable(1, m_vertex_id);
+ }
+ if (usez) {
+ v->co[2] = nlGetVariable(2, m_vertex_id);
+ }
+ }
+ }
+ if (volumepreservation) {
+ vend = compute_volume(sys->bm);
+ volume_preservation(sys->bm, sys->op, vini, vend, usex, usey, usez);
+ }
+
+}
+
+void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op)
+{
+ int i;
+ int m_vertex_id;
+ int usex, usey, usez, volumepreservation;
+ float lambda, lambda_border;
+ float w;
+ BMOIter siter;
+ BMVert *v;
+ LaplacianSystem *sys;
+
+ sys = init_laplacian_system(bm->totedge, bm->totface, bm->totvert);
+ if (!sys) return;
+ sys->bm = bm;
+ sys->op = op;
+
+ memset_laplacian_system(sys, 0);
+
+ BM_mesh_elem_index_ensure(bm, BM_VERT);
+ lambda = BMO_slot_float_get(op, "lambda");
+ lambda_border = BMO_slot_float_get(op, "lambda_border");
+ sys->min_area = 0.00001f;
+ usex = BMO_slot_bool_get(op, "use_x");
+ usey = BMO_slot_bool_get(op, "use_y");
+ usez = BMO_slot_bool_get(op, "use_z");
+ volumepreservation = BMO_slot_bool_get(op, "volume_preservation");
+
+
+ nlNewContext();
+ sys->context = nlGetCurrent();
+
+ nlSolverParameteri(NL_NB_VARIABLES, bm->totvert);
+ nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
+ nlSolverParameteri(NL_NB_ROWS, bm->totvert);
+ nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 3);
+
+ nlBegin(NL_SYSTEM);
+ for (i = 0; i < bm->totvert; i++) {
+ nlLockVariable(i);
+ }
+ BMO_ITER (v, &siter, bm, op, "verts", BM_VERT) {
+ m_vertex_id = BM_elem_index_get(v);
+ nlUnlockVariable(m_vertex_id);
+ nlSetVariable(0, m_vertex_id, v->co[0]);
+ nlSetVariable(1, m_vertex_id, v->co[1]);
+ nlSetVariable(2, m_vertex_id, v->co[2]);
+ }
+
+ nlBegin(NL_MATRIX);
+ init_laplacian_matrix(sys);
+ BMO_ITER (v, &siter, bm, op, "verts", BM_VERT) {
+ m_vertex_id = BM_elem_index_get(v);
+ nlRightHandSideAdd(0, m_vertex_id, v->co[0]);
+ nlRightHandSideAdd(1, m_vertex_id, v->co[1]);
+ nlRightHandSideAdd(2, m_vertex_id, v->co[2]);
+ i = m_vertex_id;
+ if (sys->zerola[i] == 0) {
+ w = sys->vweights[i] * sys->ring_areas[i];
+ sys->vweights[i] = (w == 0.0f) ? 0.0f : -lambda / (4.0f * w);
+ w = sys->vlengths[i];
+ sys->vlengths[i] = (w == 0.0f) ? 0.0f : -lambda_border * 2.0f / w;
+
+ if (!vert_is_boundary(v)) {
+ nlMatrixAdd(i, i, 1.0f + lambda / (4.0f * sys->ring_areas[i]));
+ }
+ else {
+ nlMatrixAdd(i, i, 1.0f + lambda_border * 2.0f);
+ }
+ }
+ else {
+ nlMatrixAdd(i, i, 1.0f);
+ }
+ }
+ fill_laplacian_matrix(sys);
+
+ nlEnd(NL_MATRIX);
+ nlEnd(NL_SYSTEM);
+
+ if (nlSolveAdvanced(NULL, NL_TRUE) ) {
+ validate_solution(sys, usex, usey, usez, volumepreservation);
+ }
+
+ delete_laplacian_system(sys);
+}