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