Fix T38941: Laplacian Deform crashes on OSX
[blender.git] / source / blender / modifiers / intern / MOD_laplaciandeform.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  * The Original Code is Copyright (C) 2013 by the Blender Foundation.
19  * All rights reserved.
20  *
21  * Contributor(s): Alexander Pinzon Fernandez
22  *
23  * ***** END GPL LICENSE BLOCK *****
24  *
25  */
26
27 /** \file blender/modifiers/intern/MOD_laplaciandeform.c
28  *  \ingroup modifiers
29  */
30
31 #include "BLI_math.h"
32 #include "BLI_utildefines.h"
33 #include "BLI_string.h"
34
35 #include "MEM_guardedalloc.h"
36
37 #include "BKE_mesh.h"
38 #include "BKE_mesh_mapping.h"
39 #include "BKE_cdderivedmesh.h"
40 #include "BKE_particle.h"
41 #include "BKE_deform.h"
42
43 #include "MOD_util.h"
44
45
46 enum {
47         LAPDEFORM_SYSTEM_NOT_CHANGE = 0,
48         LAPDEFORM_SYSTEM_IS_DIFFERENT,
49         LAPDEFORM_SYSTEM_ONLY_CHANGE_ANCHORS,
50         LAPDEFORM_SYSTEM_ONLY_CHANGE_GROUP,
51         LAPDEFORM_SYSTEM_ONLY_CHANGE_MESH,
52         LAPDEFORM_SYSTEM_CHANGE_VERTEXES,
53         LAPDEFORM_SYSTEM_CHANGE_EDGES,
54         LAPDEFORM_SYSTEM_CHANGE_NOT_VALID_GROUP,
55 };
56
57 #ifdef WITH_OPENNL
58
59 #include "ONL_opennl.h"
60
61 typedef struct LaplacianSystem {
62         bool is_matrix_computed;
63         bool has_solution;
64         int total_verts;
65         int total_edges;
66         int total_faces;
67         int total_anchors;
68         int repeat;
69         char anchor_grp_name[64];       /* Vertex Group name */
70         float (*co)[3];                         /* Original vertex coordinates */
71         float (*no)[3];                         /* Original vertex normal */
72         float (*delta)[3];                      /* Differential Coordinates */
73         unsigned int (*faces)[4];       /* Copy of MFace (tessface) v1-v4 */
74         int *index_anchors;                     /* Static vertex index list */
75         int *unit_verts;                        /* Unit vectors of projected edges onto the plane orthogonal to n */
76         int *ringf_indices;                     /* Indices of faces per vertex */
77         int *ringv_indices;                     /* Indices of neighbors(vertex) per vertex */
78         NLContext *context;                     /* System for solve general implicit rotations */
79         MeshElemMap *ringf_map;         /* Map of faces per vertex */
80         MeshElemMap *ringv_map;         /* Map of vertex per vertex */
81 } LaplacianSystem;
82
83 static LaplacianSystem *newLaplacianSystem(void)
84 {
85         LaplacianSystem *sys;
86         sys = MEM_callocN(sizeof(LaplacianSystem), "DeformCache");
87
88         sys->is_matrix_computed = false;
89         sys->has_solution = false;
90         sys->total_verts = 0;
91         sys->total_edges = 0;
92         sys->total_anchors = 0;
93         sys->total_faces = 0;
94         sys->repeat = 1;
95         sys->anchor_grp_name[0] = '\0';
96
97         return sys;
98 }
99
100 static LaplacianSystem *initLaplacianSystem(int totalVerts, int totalEdges, int totalFaces, int totalAnchors,
101                                             const char defgrpName[64], int iterations)
102 {
103         LaplacianSystem *sys = newLaplacianSystem();
104
105         sys->is_matrix_computed = false;
106         sys->has_solution = false;
107         sys->total_verts = totalVerts;
108         sys->total_edges = totalEdges;
109         sys->total_faces = totalFaces;
110         sys->total_anchors = totalAnchors;
111         sys->repeat = iterations;
112         BLI_strncpy(sys->anchor_grp_name, defgrpName, sizeof(sys->anchor_grp_name));
113         sys->co = MEM_mallocN(sizeof(float[3]) * totalVerts, "DeformCoordinates");
114         sys->no = MEM_callocN(sizeof(float[3]) * totalVerts, "DeformNormals");
115         sys->delta = MEM_callocN(sizeof(float[3]) * totalVerts, "DeformDeltas");
116         sys->faces = MEM_mallocN(sizeof(int[4]) * totalFaces, "DeformFaces");
117         sys->index_anchors = MEM_mallocN(sizeof(int) * (totalAnchors), "DeformAnchors");
118         sys->unit_verts = MEM_callocN(sizeof(int) * totalVerts, "DeformUnitVerts");
119         return sys;
120 }
121
122 static void deleteLaplacianSystem(LaplacianSystem *sys)
123 {
124         MEM_SAFE_FREE(sys->co);
125         MEM_SAFE_FREE(sys->no);
126         MEM_SAFE_FREE(sys->delta);
127         MEM_SAFE_FREE(sys->faces);
128         MEM_SAFE_FREE(sys->index_anchors);
129         MEM_SAFE_FREE(sys->unit_verts);
130         MEM_SAFE_FREE(sys->ringf_indices);
131         MEM_SAFE_FREE(sys->ringv_indices);
132         MEM_SAFE_FREE(sys->ringf_map);
133         MEM_SAFE_FREE(sys->ringv_map);
134
135         if (sys->context) {
136                 nlDeleteContext(sys->context);
137         }
138         MEM_SAFE_FREE(sys);
139 }
140
141 static float cotan_weight(const float v1[3], const float v2[3], const float v3[3])
142 {
143         float a[3], b[3], c[3], clen;
144
145         sub_v3_v3v3(a, v2, v1);
146         sub_v3_v3v3(b, v3, v1);
147         cross_v3_v3v3(c, a, b);
148
149         clen = len_v3(c);
150
151         if (clen > FLT_EPSILON) {
152                 return dot_v3v3(a, b) / clen;
153         }
154         else {
155                 return 0.0f;
156         }
157 }
158
159 static void createFaceRingMap(
160         const int mvert_tot, const MFace *mface, const int mface_tot,
161         MeshElemMap **r_map, int **r_indices)
162 {
163         int i, j, totalr = 0;
164         int *indices, *index_iter;
165         MeshElemMap *map = MEM_callocN(sizeof(MeshElemMap) * mvert_tot, "DeformRingMap");
166         const MFace *mf;
167
168         for (i = 0, mf = mface; i < mface_tot; i++, mf++) {
169                 bool has_4_vert;
170
171                 has_4_vert = mf->v4 ? 1 : 0;
172
173                 for (j = 0; j < (has_4_vert ? 4 : 3); j++) {
174                         const unsigned int v_index = (*(&mf->v1 + j));
175                         map[v_index].count++;
176                         totalr++;
177                 }
178         }
179         indices = MEM_callocN(sizeof(int) * totalr, "DeformRingIndex");
180         index_iter = indices;
181         for (i = 0; i < mvert_tot; i++) {
182                 map[i].indices = index_iter;
183                 index_iter += map[i].count;
184                 map[i].count = 0;
185         }
186         for (i = 0, mf = mface; i < mface_tot; i++, mf++) {
187                 bool has_4_vert;
188
189                 has_4_vert = mf->v4 ? 1 : 0;
190
191                 for (j = 0; j < (has_4_vert ? 4 : 3); j++) {
192                         const unsigned int v_index = (*(&mf->v1 + j));
193                         map[v_index].indices[map[v_index].count] = i;
194                         map[v_index].count++;
195                 }
196         }
197         *r_map = map;
198         *r_indices = indices;
199 }
200
201 static void createVertRingMap(
202         const int mvert_tot, const MEdge *medge, const int medge_tot,
203         MeshElemMap **r_map, int **r_indices)
204 {
205         MeshElemMap *map = MEM_callocN(sizeof(MeshElemMap) * mvert_tot, "DeformNeighborsMap");
206         int i, vid[2], totalr = 0;
207         int *indices, *index_iter;
208         const MEdge *me;
209
210         for (i = 0, me = medge; i < medge_tot; i++, me++) {
211                 vid[0] = me->v1;
212                 vid[1] = me->v2;
213                 map[vid[0]].count++;
214                 map[vid[1]].count++;
215                 totalr += 2;
216         }
217         indices = MEM_callocN(sizeof(int) * totalr, "DeformNeighborsIndex");
218         index_iter = indices;
219         for (i = 0; i < mvert_tot; i++) {
220                 map[i].indices = index_iter;
221                 index_iter += map[i].count;
222                 map[i].count = 0;
223         }
224         for (i = 0, me = medge; i < medge_tot; i++, me++) {
225                 vid[0] = me->v1;
226                 vid[1] = me->v2;
227                 map[vid[0]].indices[map[vid[0]].count] = vid[1];
228                 map[vid[0]].count++;
229                 map[vid[1]].indices[map[vid[1]].count] = vid[0];
230                 map[vid[1]].count++;
231         }
232         *r_map = map;
233         *r_indices = indices;
234 }
235
236 /**
237  * This method computes the Laplacian Matrix and Differential Coordinates for all vertex in the mesh.
238  * The Linear system is LV = d
239  * Where L is Laplacian Matrix, V as the vertexes in Mesh, d is the differential coordinates
240  * The Laplacian Matrix is computes as a
241  * Lij = sum(Wij) (if i == j)
242  * Lij = Wij (if i != j)
243  * Wij is weight between vertex Vi and vertex Vj, we use cotangent weight
244  *
245  * The Differential Coordinate is computes as a
246  * di = Vi * sum(Wij) - sum(Wij * Vj)
247  * Where :
248  * di is the Differential Coordinate i
249  * sum (Wij) is the sum of all weights between vertex Vi and its vertexes neighbors (Vj)
250  * sum (Wij * Vj) is the sum of the product between vertex neighbor Vj and weight Wij for all neighborhood.
251  *
252  * This Laplacian Matrix is described in the paper:
253  * Desbrun M. et.al, Implicit fairing of irregular meshes using diffusion and curvature flow, SIGGRAPH '99, pag 317-324,
254  * New York, USA
255  *
256  * The computation of Laplace Beltrami operator on Hybrid Triangle/Quad Meshes is described in the paper:
257  * Pinzon A., Romero E., Shape Inflation With an Adapted Laplacian Operator For Hybrid Quad/Triangle Meshes,
258  * Conference on Graphics Patterns and Images, SIBGRAPI, 2013
259  *
260  * The computation of Differential Coordinates is described in the paper:
261  * Sorkine, O. Laplacian Surface Editing. Proceedings of the EUROGRAPHICS/ACM SIGGRAPH Symposium on Geometry Processing,
262  * 2004. p. 179-188.
263  */
264 static void initLaplacianMatrix(LaplacianSystem *sys)
265 {
266         float v1[3], v2[3], v3[3], v4[3], no[3];
267         float w2, w3, w4;
268         int i, j, fi;
269         bool has_4_vert;
270         unsigned int idv1, idv2, idv3, idv4;
271
272         for (fi = 0; fi < sys->total_faces; fi++) {
273                 const unsigned int *vidf = sys->faces[fi];
274
275                 idv1 = vidf[0];
276                 idv2 = vidf[1];
277                 idv3 = vidf[2];
278                 idv4 = vidf[3];
279
280                 has_4_vert = vidf[3] ? 1 : 0;
281                 if (has_4_vert) {
282                         normal_quad_v3(no, sys->co[idv1], sys->co[idv2], sys->co[idv3], sys->co[idv4]);
283                         add_v3_v3(sys->no[idv4], no);
284                         i = 4;
285                 }
286                 else {
287                         normal_tri_v3(no, sys->co[idv1], sys->co[idv2], sys->co[idv3]);
288                         i = 3;
289                 }
290                 add_v3_v3(sys->no[idv1], no);
291                 add_v3_v3(sys->no[idv2], no);
292                 add_v3_v3(sys->no[idv3], no);
293
294                 for (j = 0; j < i; j++) {
295                         idv1 = vidf[j];
296                         idv2 = vidf[(j + 1) % i];
297                         idv3 = vidf[(j + 2) % i];
298                         idv4 = has_4_vert ? vidf[(j + 3) % i] : 0;
299
300                         copy_v3_v3(v1, sys->co[idv1]);
301                         copy_v3_v3(v2, sys->co[idv2]);
302                         copy_v3_v3(v3, sys->co[idv3]);
303                         if (has_4_vert) {
304                                 copy_v3_v3(v4, sys->co[idv4]);
305                         }
306
307                         if (has_4_vert) {
308
309                                 w2 = (cotan_weight(v4, v1, v2) + cotan_weight(v3, v1, v2)) / 2.0f;
310                                 w3 = (cotan_weight(v2, v3, v1) + cotan_weight(v4, v1, v3)) / 2.0f;
311                                 w4 = (cotan_weight(v2, v4, v1) + cotan_weight(v3, v4, v1)) / 2.0f;
312
313                                 sys->delta[idv1][0] -= v4[0] * w4;
314                                 sys->delta[idv1][1] -= v4[1] * w4;
315                                 sys->delta[idv1][2] -= v4[2] * w4;
316
317                                 nlRightHandSideAdd(0, idv1, -v4[0] * w4);
318                                 nlRightHandSideAdd(1, idv1, -v4[1] * w4);
319                                 nlRightHandSideAdd(2, idv1, -v4[2] * w4);
320
321                                 nlMatrixAdd(idv1, idv4, -w4);
322                         }
323                         else {
324                                 w2 = cotan_weight(v3, v1, v2);
325                                 w3 = cotan_weight(v2, v3, v1);
326                                 w4 = 0.0f;
327                         }
328
329                         sys->delta[idv1][0] += v1[0] * (w2 + w3 + w4);
330                         sys->delta[idv1][1] += v1[1] * (w2 + w3 + w4);
331                         sys->delta[idv1][2] += v1[2] * (w2 + w3 + w4);
332
333                         sys->delta[idv1][0] -= v2[0] * w2;
334                         sys->delta[idv1][1] -= v2[1] * w2;
335                         sys->delta[idv1][2] -= v2[2] * w2;
336
337                         sys->delta[idv1][0] -= v3[0] * w3;
338                         sys->delta[idv1][1] -= v3[1] * w3;
339                         sys->delta[idv1][2] -= v3[2] * w3;
340
341                         nlMatrixAdd(idv1, idv2, -w2);
342                         nlMatrixAdd(idv1, idv3, -w3);
343                         nlMatrixAdd(idv1, idv1, w2 + w3 + w4);
344
345                 }
346         }
347 }
348
349 static void computeImplictRotations(LaplacianSystem *sys)
350 {
351         int vid, *vidn = NULL;
352         float minj, mjt, qj[3], vj[3];
353         int i, j, ln;
354
355         for (i = 0; i < sys->total_verts; i++) {
356                 normalize_v3(sys->no[i]);
357                 vidn = sys->ringv_map[i].indices;
358                 ln = sys->ringv_map[i].count;
359                 minj = 1000000.0f;
360                 for (j = 0; j < ln; j++) {
361                         vid = vidn[j];
362                         copy_v3_v3(qj, sys->co[vid]);
363                         sub_v3_v3v3(vj, qj, sys->co[i]);
364                         normalize_v3(vj);
365                         mjt = fabsf(dot_v3v3(vj, sys->no[i]));
366                         if (mjt < minj) {
367                                 minj = mjt;
368                                 sys->unit_verts[i] = vidn[j];
369                         }
370                 }
371         }
372 }
373
374 static void rotateDifferentialCoordinates(LaplacianSystem *sys)
375 {
376         float alpha, beta, gamma;
377         float pj[3], ni[3], di[3];
378         float uij[3], dun[3], e2[3], pi[3], fni[3], vn[4][3];
379         int i, j, lvin, num_fni, k, fi;
380         int *fidn;
381
382         for (i = 0; i < sys->total_verts; i++) {
383                 copy_v3_v3(pi, sys->co[i]);
384                 copy_v3_v3(ni, sys->no[i]);
385                 k = sys->unit_verts[i];
386                 copy_v3_v3(pj, sys->co[k]);
387                 sub_v3_v3v3(uij, pj, pi);
388                 mul_v3_v3fl(dun, ni, dot_v3v3(uij, ni));
389                 sub_v3_v3(uij, dun);
390                 normalize_v3(uij);
391                 cross_v3_v3v3(e2, ni, uij);
392                 copy_v3_v3(di, sys->delta[i]);
393                 alpha = dot_v3v3(ni, di);
394                 beta = dot_v3v3(uij, di);
395                 gamma = dot_v3v3(e2, di);
396
397                 pi[0] = nlGetVariable(0, i);
398                 pi[1] = nlGetVariable(1, i);
399                 pi[2] = nlGetVariable(2, i);
400                 zero_v3(ni);
401                 num_fni = 0;
402                 num_fni = sys->ringf_map[i].count;
403                 for (fi = 0; fi < num_fni; fi++) {
404                         const unsigned int *vin;
405                         fidn = sys->ringf_map[i].indices;
406                         vin = sys->faces[fidn[fi]];
407                         lvin = vin[3] ? 4 : 3;
408                         for (j = 0; j < lvin; j++) {
409                                 vn[j][0] = nlGetVariable(0, vin[j]);
410                                 vn[j][1] = nlGetVariable(1, vin[j]);
411                                 vn[j][2] = nlGetVariable(2, vin[j]);
412                                 if (vin[j] == sys->unit_verts[i]) {
413                                         copy_v3_v3(pj, vn[j]);
414                                 }
415                         }
416
417                         if (lvin == 3) {
418                                 normal_tri_v3(fni, vn[0], vn[1], vn[2]);
419                         }
420                         else if (lvin == 4) {
421                                 normal_quad_v3(fni, vn[0], vn[1], vn[2], vn[3]);
422                         }
423                         add_v3_v3(ni, fni);
424                 }
425
426                 normalize_v3(ni);
427                 sub_v3_v3v3(uij, pj, pi);
428                 mul_v3_v3fl(dun, ni, dot_v3v3(uij, ni));
429                 sub_v3_v3(uij, dun);
430                 normalize_v3(uij);
431                 cross_v3_v3v3(e2, ni, uij);
432                 fni[0] = alpha * ni[0] + beta * uij[0] + gamma * e2[0];
433                 fni[1] = alpha * ni[1] + beta * uij[1] + gamma * e2[1];
434                 fni[2] = alpha * ni[2] + beta * uij[2] + gamma * e2[2];
435
436                 if (len_squared_v3(fni) > FLT_EPSILON) {
437                         nlRightHandSideSet(0, i, fni[0]);
438                         nlRightHandSideSet(1, i, fni[1]);
439                         nlRightHandSideSet(2, i, fni[2]);
440                 }
441                 else {
442                         nlRightHandSideSet(0, i, sys->delta[i][0]);
443                         nlRightHandSideSet(1, i, sys->delta[i][1]);
444                         nlRightHandSideSet(2, i, sys->delta[i][2]);
445                 }
446         }
447 }
448
449 static void laplacianDeformPreview(LaplacianSystem *sys, float (*vertexCos)[3])
450 {
451         int vid, i, j, n, na;
452         n = sys->total_verts;
453         na = sys->total_anchors;
454
455 #ifdef OPENNL_THREADING_HACK
456         modifier_opennl_lock();
457 #endif
458
459         if (!sys->is_matrix_computed) {
460                 nlNewContext();
461                 sys->context = nlGetCurrent();
462
463                 nlSolverParameteri(NL_NB_VARIABLES, n);
464                 nlSolverParameteri(NL_SYMMETRIC, NL_FALSE);
465                 nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
466                 nlSolverParameteri(NL_NB_ROWS, n + na);
467                 nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 3);
468                 nlBegin(NL_SYSTEM);
469                 for (i = 0; i < n; i++) {
470                         nlSetVariable(0, i, sys->co[i][0]);
471                         nlSetVariable(1, i, sys->co[i][1]);
472                         nlSetVariable(2, i, sys->co[i][2]);
473                 }
474                 for (i = 0; i < na; i++) {
475                         vid = sys->index_anchors[i];
476                         nlSetVariable(0, vid, vertexCos[vid][0]);
477                         nlSetVariable(1, vid, vertexCos[vid][1]);
478                         nlSetVariable(2, vid, vertexCos[vid][2]);
479                 }
480                 nlBegin(NL_MATRIX);
481
482                 initLaplacianMatrix(sys);
483                 computeImplictRotations(sys);
484
485                 for (i = 0; i < n; i++) {
486                         nlRightHandSideSet(0, i, sys->delta[i][0]);
487                         nlRightHandSideSet(1, i, sys->delta[i][1]);
488                         nlRightHandSideSet(2, i, sys->delta[i][2]);
489                 }
490                 for (i = 0; i < na; i++) {
491                         vid = sys->index_anchors[i];
492                         nlRightHandSideSet(0, n + i, vertexCos[vid][0]);
493                         nlRightHandSideSet(1, n + i, vertexCos[vid][1]);
494                         nlRightHandSideSet(2, n + i, vertexCos[vid][2]);
495                         nlMatrixAdd(n + i, vid, 1.0f);
496                 }
497                 nlEnd(NL_MATRIX);
498                 nlEnd(NL_SYSTEM);
499                 if (nlSolveAdvanced(NULL, NL_TRUE)) {
500                         sys->has_solution = true;
501
502                         for (j = 1; j <= sys->repeat; j++) {
503                                 nlBegin(NL_SYSTEM);
504                                 nlBegin(NL_MATRIX);
505                                 rotateDifferentialCoordinates(sys);
506
507                                 for (i = 0; i < na; i++) {
508                                         vid = sys->index_anchors[i];
509                                         nlRightHandSideSet(0, n + i, vertexCos[vid][0]);
510                                         nlRightHandSideSet(1, n + i, vertexCos[vid][1]);
511                                         nlRightHandSideSet(2, n + i, vertexCos[vid][2]);
512                                 }
513
514                                 nlEnd(NL_MATRIX);
515                                 nlEnd(NL_SYSTEM);
516                                 if (!nlSolveAdvanced(NULL, NL_FALSE)) {
517                                         sys->has_solution = false;
518                                         break;
519                                 }
520                         }
521                         if (sys->has_solution) {
522                                 for (vid = 0; vid < sys->total_verts; vid++) {
523                                         vertexCos[vid][0] = nlGetVariable(0, vid);
524                                         vertexCos[vid][1] = nlGetVariable(1, vid);
525                                         vertexCos[vid][2] = nlGetVariable(2, vid);
526                                 }
527                         }
528                         else {
529                                 sys->has_solution = false;
530                         }
531
532                 }
533                 else {
534                         sys->has_solution = false;
535                 }
536                 sys->is_matrix_computed = true;
537
538         }
539         else if (sys->has_solution) {
540                 nlMakeCurrent(sys->context);
541
542                 nlBegin(NL_SYSTEM);
543                 nlBegin(NL_MATRIX);
544
545                 for (i = 0; i < n; i++) {
546                         nlRightHandSideSet(0, i, sys->delta[i][0]);
547                         nlRightHandSideSet(1, i, sys->delta[i][1]);
548                         nlRightHandSideSet(2, i, sys->delta[i][2]);
549                 }
550                 for (i = 0; i < na; i++) {
551                         vid = sys->index_anchors[i];
552                         nlRightHandSideSet(0, n + i, vertexCos[vid][0]);
553                         nlRightHandSideSet(1, n + i, vertexCos[vid][1]);
554                         nlRightHandSideSet(2, n + i, vertexCos[vid][2]);
555                         nlMatrixAdd(n + i, vid, 1.0f);
556                 }
557
558                 nlEnd(NL_MATRIX);
559                 nlEnd(NL_SYSTEM);
560                 if (nlSolveAdvanced(NULL, NL_FALSE)) {
561                         sys->has_solution = true;
562                         for (j = 1; j <= sys->repeat; j++) {
563                                 nlBegin(NL_SYSTEM);
564                                 nlBegin(NL_MATRIX);
565                                 rotateDifferentialCoordinates(sys);
566
567                                 for (i = 0; i < na; i++) {
568                                         vid = sys->index_anchors[i];
569                                         nlRightHandSideSet(0, n + i, vertexCos[vid][0]);
570                                         nlRightHandSideSet(1, n + i, vertexCos[vid][1]);
571                                         nlRightHandSideSet(2, n + i, vertexCos[vid][2]);
572                                 }
573                                 nlEnd(NL_MATRIX);
574                                 nlEnd(NL_SYSTEM);
575                                 if (!nlSolveAdvanced(NULL, NL_FALSE)) {
576                                         sys->has_solution = false;
577                                         break;
578                                 }
579                         }
580                         if (sys->has_solution) {
581                                 for (vid = 0; vid < sys->total_verts; vid++) {
582                                         vertexCos[vid][0] = nlGetVariable(0, vid);
583                                         vertexCos[vid][1] = nlGetVariable(1, vid);
584                                         vertexCos[vid][2] = nlGetVariable(2, vid);
585                                 }
586                         }
587                         else {
588                                 sys->has_solution = false;
589                         }
590                 }
591                 else {
592                         sys->has_solution = false;
593                 }
594         }
595
596 #ifdef OPENNL_THREADING_HACK
597         modifier_opennl_unlock();
598 #endif
599 }
600
601 static bool isValidVertexGroup(LaplacianDeformModifierData *lmd, Object *ob, DerivedMesh *dm)
602 {
603         int defgrp_index;
604         MDeformVert *dvert = NULL;
605
606         modifier_get_vgroup(ob, dm, lmd->anchor_grp_name, &dvert, &defgrp_index);
607
608         return  (dvert != NULL);
609 }
610
611 static void initSystem(LaplacianDeformModifierData *lmd, Object *ob, DerivedMesh *dm,
612                        float (*vertexCos)[3], int numVerts)
613 {
614         int i;
615         int defgrp_index;
616         int total_anchors;
617         float wpaint;
618         MDeformVert *dvert = NULL;
619         MDeformVert *dv = NULL;
620         LaplacianSystem *sys;
621         if (isValidVertexGroup(lmd, ob, dm)) {
622                 int *index_anchors = MEM_mallocN(sizeof(int) * numVerts, __func__);  /* over-alloc */
623                 MFace *tessface;
624                 STACK_DECLARE(index_anchors);
625                 STACK_INIT(index_anchors);
626
627                 modifier_get_vgroup(ob, dm, lmd->anchor_grp_name, &dvert, &defgrp_index);
628                 BLI_assert(dvert != NULL);
629                 dv = dvert;
630                 for (i = 0; i < numVerts; i++) {
631                         wpaint = defvert_find_weight(dv, defgrp_index);
632                         dv++;
633                         if (wpaint > 0.0f) {
634                                 STACK_PUSH(index_anchors, i);
635                         }
636                 }
637                 DM_ensure_tessface(dm);
638                 total_anchors = STACK_SIZE(index_anchors);
639                 lmd->cache_system = initLaplacianSystem(numVerts, dm->getNumEdges(dm), dm->getNumTessFaces(dm),
640                                                        total_anchors, lmd->anchor_grp_name, lmd->repeat);
641                 sys = (LaplacianSystem *)lmd->cache_system;
642                 memcpy(sys->index_anchors, index_anchors, sizeof(int) * total_anchors);
643                 memcpy(sys->co, vertexCos, sizeof(float[3]) * numVerts);
644                 MEM_freeN(index_anchors);
645                 STACK_FREE(index_anchors);
646                 lmd->vertexco = MEM_mallocN(sizeof(float[3]) * numVerts, "ModDeformCoordinates");
647                 memcpy(lmd->vertexco, vertexCos, sizeof(float[3]) * numVerts);
648                 lmd->total_verts = numVerts;
649
650                 createFaceRingMap(
651                             dm->getNumVerts(dm), dm->getTessFaceArray(dm), dm->getNumTessFaces(dm),
652                             &sys->ringf_map, &sys->ringf_indices);
653                 createVertRingMap(
654                             dm->getNumVerts(dm), dm->getEdgeArray(dm), dm->getNumEdges(dm),
655                             &sys->ringv_map, &sys->ringv_indices);
656
657
658                 tessface = dm->getTessFaceArray(dm);
659
660                 for (i = 0; i < sys->total_faces; i++) {
661                         memcpy(&sys->faces[i], &tessface[i].v1, sizeof(*sys->faces));
662                 }
663         }
664 }
665
666 static int isSystemDifferent(LaplacianDeformModifierData *lmd, Object *ob, DerivedMesh *dm, int numVerts)
667 {
668         int i;
669         int defgrp_index;
670         int total_anchors = 0;
671         float wpaint;
672         MDeformVert *dvert = NULL;
673         MDeformVert *dv = NULL;
674         LaplacianSystem *sys = (LaplacianSystem *)lmd->cache_system;
675
676         if (sys->total_verts != numVerts) {
677                 return LAPDEFORM_SYSTEM_CHANGE_VERTEXES;
678         }
679         if (sys->total_edges != dm->getNumEdges(dm)) {
680                 return LAPDEFORM_SYSTEM_CHANGE_EDGES;
681         }
682         if (!STREQ(lmd->anchor_grp_name, sys->anchor_grp_name)) {
683                 return LAPDEFORM_SYSTEM_ONLY_CHANGE_GROUP;
684         }
685         modifier_get_vgroup(ob, dm, lmd->anchor_grp_name, &dvert, &defgrp_index);
686         if (!dvert) {
687                 return LAPDEFORM_SYSTEM_CHANGE_NOT_VALID_GROUP;
688         }
689         dv = dvert;
690         for (i = 0; i < numVerts; i++) {
691                 wpaint = defvert_find_weight(dv, defgrp_index);
692                 dv++;
693                 if (wpaint > 0.0f) {
694                         total_anchors++;
695                 }
696         }
697         if (sys->total_anchors != total_anchors) {
698                 return LAPDEFORM_SYSTEM_ONLY_CHANGE_ANCHORS;
699         }
700
701         return LAPDEFORM_SYSTEM_NOT_CHANGE;
702 }
703
704 static void LaplacianDeformModifier_do(
705         LaplacianDeformModifierData *lmd, Object *ob, DerivedMesh *dm,
706         float (*vertexCos)[3], int numVerts)
707 {
708         float (*filevertexCos)[3];
709         int sysdif;
710         LaplacianSystem *sys = NULL;
711         filevertexCos = NULL;
712         if (!(lmd->flag & MOD_LAPLACIANDEFORM_BIND)) {
713                 if (lmd->cache_system) {
714                         sys = lmd->cache_system;
715                         deleteLaplacianSystem(sys);
716                         lmd->cache_system = NULL;
717                 }
718                 lmd->total_verts = 0;
719                 MEM_SAFE_FREE(lmd->vertexco);
720                 return;
721         }
722         if (lmd->cache_system) {
723                 sysdif = isSystemDifferent(lmd, ob, dm, numVerts);
724                 sys = lmd->cache_system;
725                 if (sysdif) {
726                         if (sysdif == LAPDEFORM_SYSTEM_ONLY_CHANGE_ANCHORS || sysdif == LAPDEFORM_SYSTEM_ONLY_CHANGE_GROUP) {
727                                 filevertexCos = MEM_mallocN(sizeof(float[3]) * numVerts, "TempModDeformCoordinates");
728                                 memcpy(filevertexCos, lmd->vertexco, sizeof(float[3]) * numVerts);
729                                 MEM_SAFE_FREE(lmd->vertexco);
730                                 lmd->total_verts = 0;
731                                 deleteLaplacianSystem(sys);
732                                 lmd->cache_system = NULL;
733                                 initSystem(lmd, ob, dm, filevertexCos, numVerts);
734                                 sys = lmd->cache_system; /* may have been reallocated */
735                                 MEM_SAFE_FREE(filevertexCos);
736                                 if (sys) {
737                                         laplacianDeformPreview(sys, vertexCos);
738                                 }
739                         }
740                         else {
741                                 if (sysdif == LAPDEFORM_SYSTEM_CHANGE_VERTEXES) {
742                                         modifier_setError(&lmd->modifier, "Vertices changed from %d to %d", lmd->total_verts, numVerts);
743                                 }
744                                 else if (sysdif == LAPDEFORM_SYSTEM_CHANGE_EDGES) {
745                                         modifier_setError(&lmd->modifier, "Edges changed from %d to %d", sys->total_edges, dm->getNumEdges(dm));
746                                 }
747                                 else if (sysdif == LAPDEFORM_SYSTEM_CHANGE_NOT_VALID_GROUP) {
748                                         modifier_setError(&lmd->modifier, "Vertex group '%s' is not valid", sys->anchor_grp_name);
749                                 }
750                         }
751                 }
752                 else {
753                         sys->repeat = lmd->repeat;
754                         laplacianDeformPreview(sys, vertexCos);
755                 }
756         }
757         else {
758                 if (!isValidVertexGroup(lmd, ob, dm)) {
759                         modifier_setError(&lmd->modifier, "Vertex group '%s' is not valid", lmd->anchor_grp_name);
760                         lmd->flag &= ~MOD_LAPLACIANDEFORM_BIND;
761                 }
762                 else if (lmd->total_verts > 0 && lmd->total_verts == numVerts) {
763                         filevertexCos = MEM_mallocN(sizeof(float[3]) * numVerts, "TempDeformCoordinates");
764                         memcpy(filevertexCos, lmd->vertexco, sizeof(float[3]) * numVerts);
765                         MEM_SAFE_FREE(lmd->vertexco);
766                         lmd->total_verts = 0;
767                         initSystem(lmd, ob, dm, filevertexCos, numVerts);
768                         sys = lmd->cache_system;
769                         MEM_SAFE_FREE(filevertexCos);
770                         laplacianDeformPreview(sys, vertexCos);
771                 }
772                 else {
773                         initSystem(lmd, ob, dm, vertexCos, numVerts);
774                         sys = lmd->cache_system;
775                         laplacianDeformPreview(sys, vertexCos);
776                 }
777         }
778         if (sys && sys->is_matrix_computed && !sys->has_solution) {
779                 modifier_setError(&lmd->modifier, "The system did not find a solution");
780         }
781 }
782
783 #else  /* WITH_OPENNL */
784 static void LaplacianDeformModifier_do(
785         LaplacianDeformModifierData *lmd, Object *ob, DerivedMesh *dm,
786         float (*vertexCos)[3], int numVerts)
787 {
788         (void)lmd, (void)ob, (void)dm, (void)vertexCos, (void)numVerts;
789 }
790 #endif  /* WITH_OPENNL */
791
792 static void initData(ModifierData *md)
793 {
794         LaplacianDeformModifierData *lmd = (LaplacianDeformModifierData *)md;
795         lmd->anchor_grp_name[0] = '\0';
796         lmd->total_verts = 0;
797         lmd->repeat = 1;
798         lmd->vertexco = NULL;
799         lmd->cache_system = NULL;
800         lmd->flag = 0;
801 }
802
803 static void copyData(ModifierData *md, ModifierData *target)
804 {
805         LaplacianDeformModifierData *lmd = (LaplacianDeformModifierData *)md;
806         LaplacianDeformModifierData *tlmd = (LaplacianDeformModifierData *)target;
807
808         modifier_copyData_generic(md, target);
809
810         tlmd->vertexco = MEM_dupallocN(lmd->vertexco);
811         tlmd->cache_system = MEM_dupallocN(lmd->cache_system);
812 }
813
814 static bool isDisabled(ModifierData *md, int UNUSED(useRenderParams))
815 {
816         LaplacianDeformModifierData *lmd = (LaplacianDeformModifierData *)md;
817         if (lmd->anchor_grp_name[0]) return 0;
818         return 1;
819 }
820
821 static CustomDataMask requiredDataMask(Object *UNUSED(ob), ModifierData *md)
822 {
823         LaplacianDeformModifierData *lmd = (LaplacianDeformModifierData *)md;
824         CustomDataMask dataMask = 0;
825         if (lmd->anchor_grp_name[0]) dataMask |= CD_MASK_MDEFORMVERT;
826         return dataMask;
827 }
828
829 static void deformVerts(ModifierData *md, Object *ob, DerivedMesh *derivedData,
830                         float (*vertexCos)[3], int numVerts, ModifierApplyFlag UNUSED(flag))
831 {
832         DerivedMesh *dm = get_dm(ob, NULL, derivedData, NULL, false, false);
833
834         LaplacianDeformModifier_do((LaplacianDeformModifierData *)md, ob, dm, vertexCos, numVerts);
835         if (dm != derivedData) {
836                 dm->release(dm);
837         }
838 }
839
840 static void deformVertsEM(
841         ModifierData *md, Object *ob, struct BMEditMesh *editData,
842         DerivedMesh *derivedData, float (*vertexCos)[3], int numVerts)
843 {
844         DerivedMesh *dm = get_dm(ob, editData, derivedData, NULL, false, false);
845         LaplacianDeformModifier_do((LaplacianDeformModifierData *)md, ob, dm,
846                                    vertexCos, numVerts);
847         if (dm != derivedData) {
848                 dm->release(dm);
849         }
850 }
851
852 static void freeData(ModifierData *md)
853 {
854         LaplacianDeformModifierData *lmd = (LaplacianDeformModifierData *)md;
855 #ifdef WITH_OPENNL
856         LaplacianSystem *sys = (LaplacianSystem *)lmd->cache_system;
857         if (sys) {
858                 deleteLaplacianSystem(sys);
859         }
860 #endif
861         MEM_SAFE_FREE(lmd->vertexco);
862         lmd->total_verts = 0;
863 }
864
865 ModifierTypeInfo modifierType_LaplacianDeform = {
866         /* name */              "LaplacianDeform",
867         /* structName */        "LaplacianDeformModifierData",
868         /* structSize */        sizeof(LaplacianDeformModifierData),
869         /* type */              eModifierTypeType_OnlyDeform,
870         /* flags */             eModifierTypeFlag_AcceptsMesh | eModifierTypeFlag_SupportsEditmode,
871         /* copyData */          copyData,
872         /* deformVerts */       deformVerts,
873         /* deformMatrices */    NULL,
874         /* deformVertsEM */     deformVertsEM,
875         /* deformMatricesEM */  NULL,
876         /* applyModifier */     NULL,
877         /* applyModifierEM */   NULL,
878         /* initData */          initData,
879         /* requiredDataMask */  requiredDataMask,
880         /* freeData */          freeData,
881         /* isDisabled */        isDisabled,
882         /* updateDepgraph */    NULL,
883         /* dependsOnTime */     NULL,
884         /* dependsOnNormals */  NULL,
885         /* foreachObjectLink */ NULL,
886         /* foreachIDLink */     NULL,
887         /* foreachTexLink */    NULL,
888 };