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