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