c36c7467b89308b750876cce8ea476e6d35ef4d2
[blender.git] / source / blender / bmesh / operators / bmo_smooth_laplacian.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * Contributor(s): Alexander Pinzon
19  *
20  * ***** END GPL LICENSE BLOCK *****
21  */
22
23 /** \file blender/bmesh/operators/bmo_smooth_laplacian.c
24  *  \ingroup bmesh
25  *
26  * Advanced smoothing.
27  */
28
29 #include "MEM_guardedalloc.h"
30
31
32 #include "BLI_math.h"
33
34
35 #include "bmesh.h"
36
37 #include "intern/bmesh_operators_private.h" /* own include */
38
39 #ifdef WITH_OPENNL
40
41 #include "ONL_opennl.h"
42
43 // #define SMOOTH_LAPLACIAN_AREA_FACTOR 4.0f  /* UNUSED */
44 // #define SMOOTH_LAPLACIAN_EDGE_FACTOR 2.0f  /* UNUSED */
45 #define SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE 1.8f
46 #define SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE 0.15f
47
48 struct BLaplacianSystem {
49         float *eweights;        /* Length weights per Edge */
50         float (*fweights)[3];   /* Cotangent weights per face */
51         float *ring_areas;      /* Total area per ring*/
52         float *vlengths;        /* Total sum of lengths(edges) per vertice*/
53         float *vweights;        /* Total sum of weights per vertice*/
54         int numEdges;           /* Number of edges*/
55         int numFaces;           /* Number of faces*/
56         int numVerts;           /* Number of verts*/
57         short *zerola;          /* Is zero area or length*/
58
59         /* Pointers to data*/
60         BMesh *bm;
61         BMOperator *op;
62         NLContext *context;
63
64         /*Data*/
65         float min_area;
66 };
67 typedef struct BLaplacianSystem LaplacianSystem;
68
69 static float cotan_weight(float *v1, float *v2, float *v3);
70 static int vert_is_boundary(BMVert *v);
71 static LaplacianSystem *init_laplacian_system(int a_numEdges, int a_numFaces, int a_numVerts);
72 static void init_laplacian_matrix(LaplacianSystem *sys);
73 static void delete_laplacian_system(LaplacianSystem *sys);
74 static void delete_void_pointer(void *data);
75 static void fill_laplacian_matrix(LaplacianSystem *sys);
76 static void memset_laplacian_system(LaplacianSystem *sys, int val);
77 static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez, int preserve_volume);
78 static void volume_preservation(BMOperator *op, float vini, float vend, int usex, int usey, int usez);
79
80 static void delete_void_pointer(void *data)
81 {
82         if (data) {
83                 MEM_freeN(data);
84         }
85 }
86
87 static void delete_laplacian_system(LaplacianSystem *sys)
88 {
89         delete_void_pointer(sys->eweights);
90         delete_void_pointer(sys->fweights);
91         delete_void_pointer(sys->ring_areas);
92         delete_void_pointer(sys->vlengths);
93         delete_void_pointer(sys->vweights);
94         delete_void_pointer(sys->zerola);
95         if (sys->context) {
96                 nlDeleteContext(sys->context);
97         }
98         sys->bm = NULL;
99         sys->op = NULL;
100         MEM_freeN(sys);
101 }
102
103 static void memset_laplacian_system(LaplacianSystem *sys, int val)
104 {
105         memset(sys->eweights,     val, sizeof(float) * sys->numEdges);
106         memset(sys->fweights,     val, sizeof(float) * sys->numFaces * 3);
107         memset(sys->ring_areas,   val, sizeof(float) * sys->numVerts);
108         memset(sys->vlengths,     val, sizeof(float) * sys->numVerts);
109         memset(sys->vweights,     val, sizeof(float) * sys->numVerts);
110         memset(sys->zerola,       val, sizeof(short) * sys->numVerts);
111 }
112
113 static LaplacianSystem *init_laplacian_system(int a_numEdges, int a_numFaces, int a_numVerts)
114 {
115         LaplacianSystem *sys;
116         sys = MEM_callocN(sizeof(LaplacianSystem), "ModLaplSmoothSystem");
117         sys->numEdges = a_numEdges;
118         sys->numFaces = a_numFaces;
119         sys->numVerts = a_numVerts;
120
121         sys->eweights =  MEM_callocN(sizeof(float) * sys->numEdges, "ModLaplSmoothEWeight");
122         if (!sys->eweights) {
123                 delete_laplacian_system(sys);
124                 return NULL;
125         }
126
127         sys->fweights =  MEM_callocN(sizeof(float) * 3 * sys->numFaces, "ModLaplSmoothFWeight");
128         if (!sys->fweights) {
129                 delete_laplacian_system(sys);
130                 return NULL;
131         }
132
133         sys->ring_areas =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothRingAreas");
134         if (!sys->ring_areas) {
135                 delete_laplacian_system(sys);
136                 return NULL;
137         }
138
139         sys->vlengths =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVlengths");
140         if (!sys->vlengths) {
141                 delete_laplacian_system(sys);
142                 return NULL;
143         }
144
145         sys->vweights =  MEM_callocN(sizeof(float) * sys->numVerts, "ModLaplSmoothVweights");
146         if (!sys->vweights) {
147                 delete_laplacian_system(sys);
148                 return NULL;
149         }
150
151         sys->zerola =  MEM_callocN(sizeof(short) * sys->numVerts, "ModLaplSmoothZeloa");
152         if (!sys->zerola) {
153                 delete_laplacian_system(sys);
154                 return NULL;
155         }
156
157         return sys;
158 }
159
160 /* Compute weigth between vertice v_i and all your neighbors
161  * weight between v_i and v_neighbor
162  * Wij = cot(alpha) + cot(beta) / (4.0 * total area of all faces  * sum all weight)
163  *        v_i *
164  *          / | \
165  *         /  |  \
166  *  v_beta*   |   * v_alpha
167  *         \  |  /
168  *          \ | /
169  *            * v_neighbor
170  */
171
172 static void init_laplacian_matrix(LaplacianSystem *sys)
173 {
174         float areaf;
175         float *v1, *v2, *v3, *v4;
176         float w1, w2, w3, w4;
177         int i, j;
178         bool has_4_vert;
179         unsigned int idv1, idv2, idv3, idv4, idv[4];
180         BMEdge *e;
181         BMFace *f;
182         BMIter eiter;
183         BMIter fiter;
184         BMIter vi;
185         BMVert *vn;
186         BMVert *vf[4];
187
188         BM_ITER_MESH_INDEX (e, &eiter, sys->bm, BM_EDGES_OF_MESH, j) {
189                 if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e)) {
190                         v1 = e->v1->co;
191                         v2 =  e->v2->co;
192                         idv1 = BM_elem_index_get(e->v1);
193                         idv2 = BM_elem_index_get(e->v2);
194
195                         w1 = len_v3v3(v1, v2);
196                         if (w1 > sys->min_area) {
197                                 w1 = 1.0f / w1;
198                                 i = BM_elem_index_get(e);
199                                 sys->eweights[i] = w1;
200                                 sys->vlengths[idv1] += w1;
201                                 sys->vlengths[idv2] += w1;
202                         }
203                         else {
204                                 sys->zerola[idv1] = 1;
205                                 sys->zerola[idv2] = 1;
206                         }
207                 }
208         }
209
210         BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) {
211                 if (BM_elem_flag_test(f, BM_ELEM_SELECT)) {
212
213                         BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) {
214                                 vf[i] = vn;
215                         }
216                         has_4_vert = (i == 4) ? 1 : 0;
217                         idv1 = BM_elem_index_get(vf[0]);
218                         idv2 = BM_elem_index_get(vf[1]);
219                         idv3 = BM_elem_index_get(vf[2]);
220                         idv4 = has_4_vert ? BM_elem_index_get(vf[3]) : 0;
221
222                         v1 = vf[0]->co;
223                         v2 = vf[1]->co;
224                         v3 = vf[2]->co;
225                         v4 = has_4_vert ? vf[3]->co : NULL;
226
227                         if (has_4_vert) {
228                                 areaf = area_quad_v3(v1, v2, v3, v4);
229                         }
230                         else {
231                                 areaf = area_tri_v3(v1, v2, v3);
232                         }
233
234                         if (fabsf(areaf) < sys->min_area) {
235                                 sys->zerola[idv1] = 1;
236                                 sys->zerola[idv2] = 1;
237                                 sys->zerola[idv3] = 1;
238                                 if (has_4_vert) sys->zerola[idv4] = 1;
239                         }
240
241                         sys->ring_areas[idv1] += areaf;
242                         sys->ring_areas[idv2] += areaf;
243                         sys->ring_areas[idv3] += areaf;
244                         if (has_4_vert) sys->ring_areas[idv4] += areaf;
245
246                         if (has_4_vert) {
247
248                                 idv[0] = idv1;
249                                 idv[1] = idv2;
250                                 idv[2] = idv3;
251                                 idv[3] = idv4;
252
253                                 for (j = 0; j < 4; j++) {
254                                         idv1 = idv[j];
255                                         idv2 = idv[(j + 1) % 4];
256                                         idv3 = idv[(j + 2) % 4];
257                                         idv4 = idv[(j + 3) % 4];
258
259                                         v1 = vf[j]->co;
260                                         v2 = vf[(j + 1) % 4]->co;
261                                         v3 = vf[(j + 2) % 4]->co;
262                                         v4 = vf[(j + 3) % 4]->co;
263
264                                         w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
265                                         w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
266                                         w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
267
268                                         sys->vweights[idv1] += (w2 + w3 + w4) / 4.0f;
269                                 }
270                         }
271                         else {
272                                 i = BM_elem_index_get(f);
273
274                                 w1 = cotan_weight(v1, v2, v3);
275                                 w2 = cotan_weight(v2, v3, v1);
276                                 w3 = cotan_weight(v3, v1, v2);
277
278                                 sys->fweights[i][0] += w1;
279                                 sys->fweights[i][1] += w2;
280                                 sys->fweights[i][2] += w3;
281
282                                 sys->vweights[idv1] += w2 + w3;
283                                 sys->vweights[idv2] += w1 + w3;
284                                 sys->vweights[idv3] += w1 + w2;
285                         }
286                 }
287         }
288 }
289
290 static void fill_laplacian_matrix(LaplacianSystem *sys)
291 {
292         float *v1, *v2, *v3, *v4;
293         float w2, w3, w4;
294         int i, j;
295         bool has_4_vert;
296         unsigned int idv1, idv2, idv3, idv4, idv[4];
297
298         BMEdge *e;
299         BMFace *f;
300         BMIter eiter;
301         BMIter fiter;
302         BMIter vi;
303         BMVert *vn;
304         BMVert *vf[4];
305
306         BM_ITER_MESH (f, &fiter, sys->bm, BM_FACES_OF_MESH) {
307                 if (BM_elem_flag_test(f, BM_ELEM_SELECT)) {
308                         BM_ITER_ELEM_INDEX (vn, &vi, f, BM_VERTS_OF_FACE, i) {
309                                 vf[i] = vn;
310                         }
311                         has_4_vert = (i == 4) ? 1 : 0;
312                         if (has_4_vert) {
313                                 idv[0] = BM_elem_index_get(vf[0]);
314                                 idv[1] = BM_elem_index_get(vf[1]);
315                                 idv[2] = BM_elem_index_get(vf[2]);
316                                 idv[3] = BM_elem_index_get(vf[3]);
317                                 for (j = 0; j < 4; j++) {
318                                         idv1 = idv[j];
319                                         idv2 = idv[(j + 1) % 4];
320                                         idv3 = idv[(j + 2) % 4];
321                                         idv4 = idv[(j + 3) % 4];
322
323                                         v1 = vf[j]->co;
324                                         v2 = vf[(j + 1) % 4]->co;
325                                         v3 = vf[(j + 2) % 4]->co;
326                                         v4 = vf[(j + 3) % 4]->co;
327
328                                         w2 = cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2);
329                                         w3 = cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3);
330                                         w4 = cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1);
331
332                                         w2 = w2 / 4.0f;
333                                         w3 = w3 / 4.0f;
334                                         w4 = w4 / 4.0f;
335
336                                         if (!vert_is_boundary(vf[j]) && sys->zerola[idv1] == 0) {
337                                                 nlMatrixAdd(idv1, idv2, w2 * sys->vweights[idv1]);
338                                                 nlMatrixAdd(idv1, idv3, w3 * sys->vweights[idv1]);
339                                                 nlMatrixAdd(idv1, idv4, w4 * sys->vweights[idv1]);
340                                         }
341                                 }
342                         }
343                         else {
344                                 idv1 = BM_elem_index_get(vf[0]);
345                                 idv2 = BM_elem_index_get(vf[1]);
346                                 idv3 = BM_elem_index_get(vf[2]);
347                                 /* Is ring if number of faces == number of edges around vertice*/
348                                 i = BM_elem_index_get(f);
349                                 if (!vert_is_boundary(vf[0]) && sys->zerola[idv1] == 0) {
350                                         nlMatrixAdd(idv1, idv2, sys->fweights[i][2] * sys->vweights[idv1]);
351                                         nlMatrixAdd(idv1, idv3, sys->fweights[i][1] * sys->vweights[idv1]);
352                                 }
353                                 if (!vert_is_boundary(vf[1]) && sys->zerola[idv2] == 0) {
354                                         nlMatrixAdd(idv2, idv1, sys->fweights[i][2] * sys->vweights[idv2]);
355                                         nlMatrixAdd(idv2, idv3, sys->fweights[i][0] * sys->vweights[idv2]);
356                                 }
357                                 if (!vert_is_boundary(vf[2]) && sys->zerola[idv3] == 0) {
358                                         nlMatrixAdd(idv3, idv1, sys->fweights[i][1] * sys->vweights[idv3]);
359                                         nlMatrixAdd(idv3, idv2, sys->fweights[i][0] * sys->vweights[idv3]);
360                                 }
361                         }
362                 }
363         }
364         BM_ITER_MESH (e, &eiter, sys->bm, BM_EDGES_OF_MESH) {
365                 if (!BM_elem_flag_test(e, BM_ELEM_SELECT) && BM_edge_is_boundary(e)) {
366                         v1 = e->v1->co;
367                         v2 =  e->v2->co;
368                         idv1 = BM_elem_index_get(e->v1);
369                         idv2 = BM_elem_index_get(e->v2);
370                         if (sys->zerola[idv1] == 0 && sys->zerola[idv2] == 0) {
371                                 i = BM_elem_index_get(e);
372                                 nlMatrixAdd(idv1, idv2, sys->eweights[i] * sys->vlengths[idv1]);
373                                 nlMatrixAdd(idv2, idv1, sys->eweights[i] * sys->vlengths[idv2]);
374                         }
375                 }
376         }
377 }
378
379 static float cotan_weight(float *v1, float *v2, float *v3)
380 {
381         float a[3], b[3], c[3], clen;
382
383         sub_v3_v3v3(a, v2, v1);
384         sub_v3_v3v3(b, v3, v1);
385         cross_v3_v3v3(c, a, b);
386
387         clen = len_v3(c);
388
389         if (clen == 0.0f)
390                 return 0.0f;
391
392         return dot_v3v3(a, b) / clen;
393 }
394
395 static int vert_is_boundary(BMVert *v)
396 {
397         BMEdge *ed;
398         BMFace *f;
399         BMIter ei;
400         BMIter fi;
401         BM_ITER_ELEM (ed, &ei, v, BM_EDGES_OF_VERT) {
402                 if (BM_edge_is_boundary(ed)) {
403                         return 1;
404                 }
405         }
406         BM_ITER_ELEM (f, &fi, v, BM_FACES_OF_VERT) {
407                 if (!BM_elem_flag_test(f, BM_ELEM_SELECT)) {
408                         return 1;
409                 }
410         }
411         return 0;
412 }
413
414 static void volume_preservation(BMOperator *op, float vini, float vend, int usex, int usey, int usez)
415 {
416         float beta;
417         BMOIter siter;
418         BMVert *v;
419
420         if (vend != 0.0f) {
421                 beta  = pow(vini / vend, 1.0f / 3.0f);
422                 BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
423                         if (usex) {
424                                 v->co[0] *= beta;
425                         }
426                         if (usey) {
427                                 v->co[1] *= beta;
428                         }
429                         if (usez) {
430                                 v->co[2] *= beta;
431                         }
432
433                 }
434         }
435 }
436
437 static void validate_solution(LaplacianSystem *sys, int usex, int usey, int usez, int preserve_volume)
438 {
439         int m_vertex_id;
440         float leni, lene;
441         float vini, vend;
442         float *vi1, *vi2, ve1[3], ve2[3];
443         unsigned int idv1, idv2;
444         BMOIter siter;
445         BMVert *v;
446         BMEdge *e;
447         BMIter eiter;
448
449         BM_ITER_MESH  (e, &eiter, sys->bm, BM_EDGES_OF_MESH) {
450                 idv1 = BM_elem_index_get(e->v1);
451                 idv2 = BM_elem_index_get(e->v2);
452                 vi1 = e->v1->co;
453                 vi2 =  e->v2->co;
454                 ve1[0] = nlGetVariable(0, idv1);
455                 ve1[1] = nlGetVariable(1, idv1);
456                 ve1[2] = nlGetVariable(2, idv1);
457                 ve2[0] = nlGetVariable(0, idv2);
458                 ve2[1] = nlGetVariable(1, idv2);
459                 ve2[2] = nlGetVariable(2, idv2);
460                 leni = len_v3v3(vi1, vi2);
461                 lene = len_v3v3(ve1, ve2);
462                 if (lene > leni * SMOOTH_LAPLACIAN_MAX_EDGE_PERCENTAGE || lene < leni * SMOOTH_LAPLACIAN_MIN_EDGE_PERCENTAGE) {
463                         sys->zerola[idv1] = 1;
464                         sys->zerola[idv2] = 1;
465                 }
466         }
467
468         if (preserve_volume) {
469                 vini = BM_mesh_calc_volume(sys->bm, false);
470         }
471         BMO_ITER (v, &siter, sys->op->slots_in, "verts", BM_VERT) {
472                 m_vertex_id = BM_elem_index_get(v);
473                 if (sys->zerola[m_vertex_id] == 0) {
474                         if (usex) {
475                                 v->co[0] =  nlGetVariable(0, m_vertex_id);
476                         }
477                         if (usey) {
478                                 v->co[1] =  nlGetVariable(1, m_vertex_id);
479                         }
480                         if (usez) {
481                                 v->co[2] =  nlGetVariable(2, m_vertex_id);
482                         }
483                 }
484         }
485         if (preserve_volume) {
486                 vend = BM_mesh_calc_volume(sys->bm, false);
487                 volume_preservation(sys->op, vini, vend, usex, usey, usez);
488         }
489
490 }
491
492 void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op)
493 {
494         int i;
495         int m_vertex_id;
496         bool usex, usey, usez, preserve_volume;
497         float lambda_factor, lambda_border;
498         float w;
499         BMOIter siter;
500         BMVert *v;
501         LaplacianSystem *sys;
502
503         if (bm->totface == 0) return;
504         sys = init_laplacian_system(bm->totedge, bm->totface, bm->totvert);
505         if (!sys) return;
506         sys->bm = bm;
507         sys->op = op;
508
509         memset_laplacian_system(sys, 0);
510
511         BM_mesh_elem_index_ensure(bm, BM_VERT);
512         lambda_factor = BMO_slot_float_get(op->slots_in, "lambda_factor");
513         lambda_border = BMO_slot_float_get(op->slots_in, "lambda_border");
514         sys->min_area = 0.00001f;
515         usex = BMO_slot_bool_get(op->slots_in, "use_x");
516         usey = BMO_slot_bool_get(op->slots_in, "use_y");
517         usez = BMO_slot_bool_get(op->slots_in, "use_z");
518         preserve_volume = BMO_slot_bool_get(op->slots_in, "preserve_volume");
519
520
521         nlNewContext();
522         sys->context = nlGetCurrent();
523
524         nlSolverParameteri(NL_NB_VARIABLES, bm->totvert);
525         nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
526         nlSolverParameteri(NL_NB_ROWS, bm->totvert);
527         nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 3);
528
529         nlBegin(NL_SYSTEM);
530         for (i = 0; i < bm->totvert; i++) {
531                 nlLockVariable(i);
532         }
533         BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
534                 m_vertex_id = BM_elem_index_get(v);
535                 nlUnlockVariable(m_vertex_id);
536                 nlSetVariable(0, m_vertex_id, v->co[0]);
537                 nlSetVariable(1, m_vertex_id, v->co[1]);
538                 nlSetVariable(2, m_vertex_id, v->co[2]);
539         }
540
541         nlBegin(NL_MATRIX);
542         init_laplacian_matrix(sys);
543         BMO_ITER (v, &siter, op->slots_in, "verts", BM_VERT) {
544                 m_vertex_id = BM_elem_index_get(v);
545                 nlRightHandSideAdd(0, m_vertex_id, v->co[0]);
546                 nlRightHandSideAdd(1, m_vertex_id, v->co[1]);
547                 nlRightHandSideAdd(2, m_vertex_id, v->co[2]);
548                 i = m_vertex_id;
549                 if (sys->zerola[i] == 0) {
550                         w = sys->vweights[i] * sys->ring_areas[i];
551                         sys->vweights[i] = (w == 0.0f) ? 0.0f : -lambda_factor  / (4.0f * w);
552                         w = sys->vlengths[i];
553                         sys->vlengths[i] = (w == 0.0f) ? 0.0f : -lambda_border  * 2.0f / w;
554
555                         if (!vert_is_boundary(v)) {
556                                 nlMatrixAdd(i, i,  1.0f + lambda_factor / (4.0f * sys->ring_areas[i]));
557                         }
558                         else {
559                                 nlMatrixAdd(i, i,  1.0f + lambda_border * 2.0f);
560                         }
561                 }
562                 else {
563                         nlMatrixAdd(i, i, 1.0f);
564                 }
565         }
566         fill_laplacian_matrix(sys);
567
568         nlEnd(NL_MATRIX);
569         nlEnd(NL_SYSTEM);
570
571         if (nlSolveAdvanced(NULL, NL_TRUE) ) {
572                 validate_solution(sys, usex, usey, usez, preserve_volume);
573         }
574
575         delete_laplacian_system(sys);
576 }
577
578 #else  /* WITH_OPENNL */
579
580 #ifdef __GNUC__
581 #  pragma GCC diagnostic ignored "-Wunused-parameter"
582 #endif
583
584 void bmo_smooth_laplacian_vert_exec(BMesh *bm, BMOperator *op) {}
585
586 #endif  /* WITH_OPENNL */