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