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