Cleanup/refactor binding code for MeshDeform modifier.
[blender.git] / source / blender / editors / armature / meshlaplacian.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * Contributor(s): none yet.
19  *
20  * ***** END GPL LICENSE BLOCK *****
21  * meshlaplacian.c: Algorithms using the mesh laplacian.
22  */
23
24 /** \file blender/editors/armature/meshlaplacian.c
25  *  \ingroup edarmature
26  */
27
28 #include "MEM_guardedalloc.h"
29
30 #include "DNA_object_types.h"
31 #include "DNA_mesh_types.h"
32 #include "DNA_meshdata_types.h"
33 #include "DNA_scene_types.h"
34
35 #include "BLI_math.h"
36 #include "BLI_edgehash.h"
37 #include "BLI_memarena.h"
38 #include "BLI_string.h"
39 #include "BLI_alloca.h"
40
41 #include "BLT_translation.h"
42
43 #include "BKE_bvhutils.h"
44 #include "BKE_mesh.h"
45 #include "BKE_mesh_runtime.h"
46 #include "BKE_modifier.h"
47
48 #include "ED_mesh.h"
49 #include "ED_armature.h"
50
51 #include "DEG_depsgraph.h"
52
53 #include "eigen_capi.h"
54
55 #include "meshlaplacian.h"
56
57 /* ************* XXX *************** */
58 static void waitcursor(int UNUSED(val)) {}
59 static void progress_bar(int UNUSED(dummy_val), const char *UNUSED(dummy)) {}
60 static void start_progress_bar(void) {}
61 static void end_progress_bar(void) {}
62 static void error(const char *str) { printf("error: %s\n", str); }
63 /* ************* XXX *************** */
64
65
66 /************************** Laplacian System *****************************/
67
68 struct LaplacianSystem {
69         LinearSolver *context;  /* linear solver */
70
71         int totvert, totface;
72
73         float **verts;          /* vertex coordinates */
74         float *varea;           /* vertex weights for laplacian computation */
75         char *vpinned;          /* vertex pinning */
76         int (*faces)[3];        /* face vertex indices */
77         float (*fweights)[3];   /* cotangent weights per face */
78
79         int areaweights;        /* use area in cotangent weights? */
80         int storeweights;       /* store cotangent weights in fweights */
81         bool variablesdone;     /* variables set in linear system */
82
83         EdgeHash *edgehash;     /* edge hash for construction */
84
85         struct HeatWeighting {
86                 const MLoopTri *mlooptri;
87                 const MLoop *mloop;     /* needed to find vertices by index */
88                 int totvert;
89                 int tottri;
90                 float (*verts)[3];  /* vertex coordinates */
91                 float (*vnors)[3];  /* vertex normals */
92
93                 float (*root)[3];   /* bone root */
94                 float (*tip)[3];    /* bone tip */
95                 float (*source)[3]; /* vertex source */
96                 int numsource;
97
98                 float *H;           /* diagonal H matrix */
99                 float *p;           /* values from all p vectors */
100                 float *mindist;     /* minimum distance to a bone for all vertices */
101
102                 BVHTree   *bvhtree; /* ray tracing acceleration structure */
103                 const MLoopTri **vltree;  /* a looptri that the vertex belongs to */
104         } heat;
105 };
106
107 /* Laplacian matrix construction */
108
109 /* Computation of these weights for the laplacian is based on:
110  * "Discrete Differential-Geometry Operators for Triangulated 2-Manifolds",
111  * Meyer et al, 2002. Section 3.5, formula (8).
112  *
113  * We do it a bit different by going over faces instead of going over each
114  * vertex and adjacent faces, since we don't store this adjacency. Also, the
115  * formulas are tweaked a bit to work for non-manifold meshes. */
116
117 static void laplacian_increase_edge_count(EdgeHash *edgehash, int v1, int v2)
118 {
119         void **p;
120
121         if (BLI_edgehash_ensure_p(edgehash, v1, v2, &p))
122                 *p = (void *)((intptr_t)*p + (intptr_t)1);
123         else
124                 *p = (void *)((intptr_t)1);
125 }
126
127 static int laplacian_edge_count(EdgeHash *edgehash, int v1, int v2)
128 {
129         return (int)(intptr_t)BLI_edgehash_lookup(edgehash, v1, v2);
130 }
131
132 static void laplacian_triangle_area(LaplacianSystem *sys, int i1, int i2, int i3)
133 {
134         float t1, t2, t3, len1, len2, len3, area;
135         float *varea = sys->varea, *v1, *v2, *v3;
136         int obtuse = 0;
137
138         v1 = sys->verts[i1];
139         v2 = sys->verts[i2];
140         v3 = sys->verts[i3];
141
142         t1 = cotangent_tri_weight_v3(v1, v2, v3);
143         t2 = cotangent_tri_weight_v3(v2, v3, v1);
144         t3 = cotangent_tri_weight_v3(v3, v1, v2);
145
146         if      (angle_v3v3v3(v2, v1, v3) > DEG2RADF(90.0f)) obtuse = 1;
147         else if (angle_v3v3v3(v1, v2, v3) > DEG2RADF(90.0f)) obtuse = 2;
148         else if (angle_v3v3v3(v1, v3, v2) > DEG2RADF(90.0f)) obtuse = 3;
149
150         if (obtuse > 0) {
151                 area = area_tri_v3(v1, v2, v3);
152
153                 varea[i1] += (obtuse == 1) ? area : area * 0.5f;
154                 varea[i2] += (obtuse == 2) ? area : area * 0.5f;
155                 varea[i3] += (obtuse == 3) ? area : area * 0.5f;
156         }
157         else {
158                 len1 = len_v3v3(v2, v3);
159                 len2 = len_v3v3(v1, v3);
160                 len3 = len_v3v3(v1, v2);
161
162                 t1 *= len1 * len1;
163                 t2 *= len2 * len2;
164                 t3 *= len3 * len3;
165
166                 varea[i1] += (t2 + t3) * 0.25f;
167                 varea[i2] += (t1 + t3) * 0.25f;
168                 varea[i3] += (t1 + t2) * 0.25f;
169         }
170 }
171
172 static void laplacian_triangle_weights(LaplacianSystem *sys, int f, int i1, int i2, int i3)
173 {
174         float t1, t2, t3;
175         float *varea = sys->varea, *v1, *v2, *v3;
176
177         v1 = sys->verts[i1];
178         v2 = sys->verts[i2];
179         v3 = sys->verts[i3];
180
181         /* instead of *0.5 we divided by the number of faces of the edge, it still
182          * needs to be verified that this is indeed the correct thing to do! */
183         t1 = cotangent_tri_weight_v3(v1, v2, v3) / laplacian_edge_count(sys->edgehash, i2, i3);
184         t2 = cotangent_tri_weight_v3(v2, v3, v1) / laplacian_edge_count(sys->edgehash, i3, i1);
185         t3 = cotangent_tri_weight_v3(v3, v1, v2) / laplacian_edge_count(sys->edgehash, i1, i2);
186
187         EIG_linear_solver_matrix_add(sys->context, i1, i1, (t2 + t3) * varea[i1]);
188         EIG_linear_solver_matrix_add(sys->context, i2, i2, (t1 + t3) * varea[i2]);
189         EIG_linear_solver_matrix_add(sys->context, i3, i3, (t1 + t2) * varea[i3]);
190
191         EIG_linear_solver_matrix_add(sys->context, i1, i2, -t3 * varea[i1]);
192         EIG_linear_solver_matrix_add(sys->context, i2, i1, -t3 * varea[i2]);
193
194         EIG_linear_solver_matrix_add(sys->context, i2, i3, -t1 * varea[i2]);
195         EIG_linear_solver_matrix_add(sys->context, i3, i2, -t1 * varea[i3]);
196
197         EIG_linear_solver_matrix_add(sys->context, i3, i1, -t2 * varea[i3]);
198         EIG_linear_solver_matrix_add(sys->context, i1, i3, -t2 * varea[i1]);
199
200         if (sys->storeweights) {
201                 sys->fweights[f][0] = t1 * varea[i1];
202                 sys->fweights[f][1] = t2 * varea[i2];
203                 sys->fweights[f][2] = t3 * varea[i3];
204         }
205 }
206
207 static LaplacianSystem *laplacian_system_construct_begin(int totvert, int totface, int lsq)
208 {
209         LaplacianSystem *sys;
210
211         sys = MEM_callocN(sizeof(LaplacianSystem), "LaplacianSystem");
212
213         sys->verts = MEM_callocN(sizeof(float *) * totvert, "LaplacianSystemVerts");
214         sys->vpinned = MEM_callocN(sizeof(char) * totvert, "LaplacianSystemVpinned");
215         sys->faces = MEM_callocN(sizeof(int) * 3 * totface, "LaplacianSystemFaces");
216
217         sys->totvert = 0;
218         sys->totface = 0;
219
220         sys->areaweights = 1;
221         sys->storeweights = 0;
222
223         /* create linear solver */
224         if (lsq)
225                 sys->context = EIG_linear_least_squares_solver_new(0, totvert, 1);
226         else
227                 sys->context = EIG_linear_solver_new(0, totvert, 1);
228
229         return sys;
230 }
231
232 void laplacian_add_vertex(LaplacianSystem *sys, float *co, int pinned)
233 {
234         sys->verts[sys->totvert] = co;
235         sys->vpinned[sys->totvert] = pinned;
236         sys->totvert++;
237 }
238
239 void laplacian_add_triangle(LaplacianSystem *sys, int v1, int v2, int v3)
240 {
241         sys->faces[sys->totface][0] = v1;
242         sys->faces[sys->totface][1] = v2;
243         sys->faces[sys->totface][2] = v3;
244         sys->totface++;
245 }
246
247 static void laplacian_system_construct_end(LaplacianSystem *sys)
248 {
249         int (*face)[3];
250         int a, totvert = sys->totvert, totface = sys->totface;
251
252         laplacian_begin_solve(sys, 0);
253
254         sys->varea = MEM_callocN(sizeof(float) * totvert, "LaplacianSystemVarea");
255
256         sys->edgehash = BLI_edgehash_new_ex(__func__, BLI_EDGEHASH_SIZE_GUESS_FROM_POLYS(sys->totface));
257         for (a = 0, face = sys->faces; a < sys->totface; a++, face++) {
258                 laplacian_increase_edge_count(sys->edgehash, (*face)[0], (*face)[1]);
259                 laplacian_increase_edge_count(sys->edgehash, (*face)[1], (*face)[2]);
260                 laplacian_increase_edge_count(sys->edgehash, (*face)[2], (*face)[0]);
261         }
262
263         if (sys->areaweights)
264                 for (a = 0, face = sys->faces; a < sys->totface; a++, face++)
265                         laplacian_triangle_area(sys, (*face)[0], (*face)[1], (*face)[2]);
266
267         for (a = 0; a < totvert; a++) {
268                 if (sys->areaweights) {
269                         if (sys->varea[a] != 0.0f)
270                                 sys->varea[a] = 0.5f / sys->varea[a];
271                 }
272                 else
273                         sys->varea[a] = 1.0f;
274
275                 /* for heat weighting */
276                 if (sys->heat.H)
277                         EIG_linear_solver_matrix_add(sys->context, a, a, sys->heat.H[a]);
278         }
279
280         if (sys->storeweights)
281                 sys->fweights = MEM_callocN(sizeof(float) * 3 * totface, "LaplacianFWeight");
282
283         for (a = 0, face = sys->faces; a < totface; a++, face++)
284                 laplacian_triangle_weights(sys, a, (*face)[0], (*face)[1], (*face)[2]);
285
286         MEM_freeN(sys->faces);
287         sys->faces = NULL;
288
289         if (sys->varea) {
290                 MEM_freeN(sys->varea);
291                 sys->varea = NULL;
292         }
293
294         BLI_edgehash_free(sys->edgehash, NULL);
295         sys->edgehash = NULL;
296 }
297
298 static void laplacian_system_delete(LaplacianSystem *sys)
299 {
300         if (sys->verts) MEM_freeN(sys->verts);
301         if (sys->varea) MEM_freeN(sys->varea);
302         if (sys->vpinned) MEM_freeN(sys->vpinned);
303         if (sys->faces) MEM_freeN(sys->faces);
304         if (sys->fweights) MEM_freeN(sys->fweights);
305
306         EIG_linear_solver_delete(sys->context);
307         MEM_freeN(sys);
308 }
309
310 void laplacian_begin_solve(LaplacianSystem *sys, int index)
311 {
312         int a;
313
314         if (!sys->variablesdone) {
315                 if (index >= 0) {
316                         for (a = 0; a < sys->totvert; a++) {
317                                 if (sys->vpinned[a]) {
318                                         EIG_linear_solver_variable_set(sys->context, 0, a, sys->verts[a][index]);
319                                         EIG_linear_solver_variable_lock(sys->context, a);
320                                 }
321                         }
322                 }
323
324                 sys->variablesdone = true;
325         }
326 }
327
328 void laplacian_add_right_hand_side(LaplacianSystem *sys, int v, float value)
329 {
330         EIG_linear_solver_right_hand_side_add(sys->context, 0, v, value);
331 }
332
333 int laplacian_system_solve(LaplacianSystem *sys)
334 {
335         sys->variablesdone = false;
336
337         //EIG_linear_solver_print_matrix(sys->context, );
338
339         return EIG_linear_solver_solve(sys->context);
340 }
341
342 float laplacian_system_get_solution(LaplacianSystem *sys, int v)
343 {
344         return EIG_linear_solver_variable_get(sys->context, 0, v);
345 }
346
347 /************************* Heat Bone Weighting ******************************/
348 /* From "Automatic Rigging and Animation of 3D Characters"
349  * Ilya Baran and Jovan Popovic, SIGGRAPH 2007 */
350
351 #define C_WEIGHT            1.0f
352 #define WEIGHT_LIMIT_START  0.05f
353 #define WEIGHT_LIMIT_END    0.025f
354 #define DISTANCE_EPSILON    1e-4f
355
356 typedef struct BVHCallbackUserData {
357         float start[3];
358         float vec[3];
359         LaplacianSystem *sys;
360 } BVHCallbackUserData;
361
362 static void bvh_callback(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
363 {
364         BVHCallbackUserData *data = (struct BVHCallbackUserData *)userdata;
365         const MLoopTri *lt = &data->sys->heat.mlooptri[index];
366         const MLoop *mloop = data->sys->heat.mloop;
367         float (*verts)[3] = data->sys->heat.verts;
368         const float *vtri_co[3];
369         float dist_test;
370
371         vtri_co[0] = verts[mloop[lt->tri[0]].v];
372         vtri_co[1] = verts[mloop[lt->tri[1]].v];
373         vtri_co[2] = verts[mloop[lt->tri[2]].v];
374
375 #ifdef USE_KDOPBVH_WATERTIGHT
376         if (isect_ray_tri_watertight_v3(data->start, ray->isect_precalc, UNPACK3(vtri_co), &dist_test, NULL))
377 #else
378         UNUSED_VARS(ray);
379         if (isect_ray_tri_v3(data->start, data->vec, UNPACK3(vtri_co), &dist_test, NULL))
380 #endif
381         {
382                 if (dist_test < hit->dist) {
383                         float n[3];
384                         normal_tri_v3(n, UNPACK3(vtri_co));
385                         if (dot_v3v3(n, data->vec) < -1e-5f) {
386                                 hit->index = index;
387                                 hit->dist = dist_test;
388                         }
389                 }
390         }
391 }
392
393 /* Raytracing for vertex to bone/vertex visibility */
394 static void heat_ray_tree_create(LaplacianSystem *sys)
395 {
396         const MLoopTri *looptri = sys->heat.mlooptri;
397         const MLoop *mloop = sys->heat.mloop;
398         float (*verts)[3] = sys->heat.verts;
399         int tottri = sys->heat.tottri;
400         int totvert = sys->heat.totvert;
401         int a;
402
403         sys->heat.bvhtree = BLI_bvhtree_new(tottri, 0.0f, 4, 6);
404         sys->heat.vltree = MEM_callocN(sizeof(MLoopTri *) * totvert, "HeatVFaces");
405
406         for (a = 0; a < tottri; a++) {
407                 const MLoopTri *lt = &looptri[a];
408                 float bb[6];
409                 int vtri[3];
410
411                 vtri[0] = mloop[lt->tri[0]].v;
412                 vtri[1] = mloop[lt->tri[1]].v;
413                 vtri[2] = mloop[lt->tri[2]].v;
414
415                 INIT_MINMAX(bb, bb + 3);
416                 minmax_v3v3_v3(bb, bb + 3, verts[vtri[0]]);
417                 minmax_v3v3_v3(bb, bb + 3, verts[vtri[1]]);
418                 minmax_v3v3_v3(bb, bb + 3, verts[vtri[2]]);
419
420                 BLI_bvhtree_insert(sys->heat.bvhtree, a, bb, 2);
421
422                 //Setup inverse pointers to use on isect.orig
423                 sys->heat.vltree[vtri[0]] = lt;
424                 sys->heat.vltree[vtri[1]] = lt;
425                 sys->heat.vltree[vtri[2]] = lt;
426         }
427
428         BLI_bvhtree_balance(sys->heat.bvhtree);
429 }
430
431 static int heat_ray_source_visible(LaplacianSystem *sys, int vertex, int source)
432 {
433         BVHTreeRayHit hit;
434         BVHCallbackUserData data;
435         const MLoopTri *lt;
436         float end[3];
437         int visible;
438
439         lt = sys->heat.vltree[vertex];
440         if (lt == NULL)
441                 return 1;
442
443         data.sys = sys;
444         copy_v3_v3(data.start, sys->heat.verts[vertex]);
445
446         closest_to_line_segment_v3(end, data.start, sys->heat.root[source], sys->heat.tip[source]);
447
448         sub_v3_v3v3(data.vec, end, data.start);
449         madd_v3_v3v3fl(data.start, data.start, data.vec, 1e-5);
450         mul_v3_fl(data.vec, 1.0f - 2e-5f);
451
452         /* pass normalized vec + distance to bvh */
453         hit.index = -1;
454         hit.dist = normalize_v3(data.vec);
455
456         visible = BLI_bvhtree_ray_cast(sys->heat.bvhtree, data.start, data.vec, 0.0f, &hit, bvh_callback, (void *)&data) == -1;
457
458         return visible;
459 }
460
461 static float heat_source_distance(LaplacianSystem *sys, int vertex, int source)
462 {
463         float closest[3], d[3], dist, cosine;
464
465         /* compute euclidian distance */
466         closest_to_line_segment_v3(closest, sys->heat.verts[vertex], sys->heat.root[source], sys->heat.tip[source]);
467
468         sub_v3_v3v3(d, sys->heat.verts[vertex], closest);
469         dist = normalize_v3(d);
470
471         /* if the vertex normal does not point along the bone, increase distance */
472         cosine = dot_v3v3(d, sys->heat.vnors[vertex]);
473
474         return dist / (0.5f * (cosine + 1.001f));
475 }
476
477 static int heat_source_closest(LaplacianSystem *sys, int vertex, int source)
478 {
479         float dist;
480
481         dist = heat_source_distance(sys, vertex, source);
482
483         if (dist <= sys->heat.mindist[vertex] * (1.0f + DISTANCE_EPSILON))
484                 if (heat_ray_source_visible(sys, vertex, source))
485                         return 1;
486
487         return 0;
488 }
489
490 static void heat_set_H(LaplacianSystem *sys, int vertex)
491 {
492         float dist, mindist, h;
493         int j, numclosest = 0;
494
495         mindist = 1e10;
496
497         /* compute minimum distance */
498         for (j = 0; j < sys->heat.numsource; j++) {
499                 dist = heat_source_distance(sys, vertex, j);
500
501                 if (dist < mindist)
502                         mindist = dist;
503         }
504
505         sys->heat.mindist[vertex] = mindist;
506
507         /* count number of sources with approximately this minimum distance */
508         for (j = 0; j < sys->heat.numsource; j++)
509                 if (heat_source_closest(sys, vertex, j))
510                         numclosest++;
511
512         sys->heat.p[vertex] = (numclosest > 0) ? 1.0f / numclosest : 0.0f;
513
514         /* compute H entry */
515         if (numclosest > 0) {
516                 mindist = max_ff(mindist, 1e-4f);
517                 h = numclosest * C_WEIGHT / (mindist * mindist);
518         }
519         else
520                 h = 0.0f;
521
522         sys->heat.H[vertex] = h;
523 }
524
525 static void heat_calc_vnormals(LaplacianSystem *sys)
526 {
527         float fnor[3];
528         int a, v1, v2, v3, (*face)[3];
529
530         sys->heat.vnors = MEM_callocN(sizeof(float) * 3 * sys->totvert, "HeatVNors");
531
532         for (a = 0, face = sys->faces; a < sys->totface; a++, face++) {
533                 v1 = (*face)[0];
534                 v2 = (*face)[1];
535                 v3 = (*face)[2];
536
537                 normal_tri_v3(fnor, sys->verts[v1], sys->verts[v2], sys->verts[v3]);
538
539                 add_v3_v3(sys->heat.vnors[v1], fnor);
540                 add_v3_v3(sys->heat.vnors[v2], fnor);
541                 add_v3_v3(sys->heat.vnors[v3], fnor);
542         }
543
544         for (a = 0; a < sys->totvert; a++)
545                 normalize_v3(sys->heat.vnors[a]);
546 }
547
548 static void heat_laplacian_create(LaplacianSystem *sys)
549 {
550         const MLoopTri *mlooptri = sys->heat.mlooptri, *lt;
551         const MLoop *mloop = sys->heat.mloop;
552         int tottri = sys->heat.tottri;
553         int totvert = sys->heat.totvert;
554         int a;
555
556         /* heat specific definitions */
557         sys->heat.mindist = MEM_callocN(sizeof(float) * totvert, "HeatMinDist");
558         sys->heat.H = MEM_callocN(sizeof(float) * totvert, "HeatH");
559         sys->heat.p = MEM_callocN(sizeof(float) * totvert, "HeatP");
560
561         /* add verts and faces to laplacian */
562         for (a = 0; a < totvert; a++)
563                 laplacian_add_vertex(sys, sys->heat.verts[a], 0);
564
565         for (a = 0, lt = mlooptri; a < tottri; a++, lt++) {
566                 int vtri[3];
567                 vtri[0] = mloop[lt->tri[0]].v;
568                 vtri[1] = mloop[lt->tri[1]].v;
569                 vtri[2] = mloop[lt->tri[2]].v;
570                 laplacian_add_triangle(sys, UNPACK3(vtri));
571         }
572
573         /* for distance computation in set_H */
574         heat_calc_vnormals(sys);
575
576         for (a = 0; a < totvert; a++)
577                 heat_set_H(sys, a);
578 }
579
580 static void heat_system_free(LaplacianSystem *sys)
581 {
582         BLI_bvhtree_free(sys->heat.bvhtree);
583         MEM_freeN((void *)sys->heat.vltree);
584         MEM_freeN((void *)sys->heat.mlooptri);
585
586         MEM_freeN(sys->heat.mindist);
587         MEM_freeN(sys->heat.H);
588         MEM_freeN(sys->heat.p);
589         MEM_freeN(sys->heat.vnors);
590 }
591
592 static float heat_limit_weight(float weight)
593 {
594         float t;
595
596         if (weight < WEIGHT_LIMIT_END) {
597                 return 0.0f;
598         }
599         else if (weight < WEIGHT_LIMIT_START) {
600                 t = (weight - WEIGHT_LIMIT_END) / (WEIGHT_LIMIT_START - WEIGHT_LIMIT_END);
601                 return t * WEIGHT_LIMIT_START;
602         }
603         else
604                 return weight;
605 }
606
607 void heat_bone_weighting(
608         Object *ob, Mesh *me, float (*verts)[3], int numsource,
609         bDeformGroup **dgrouplist, bDeformGroup **dgroupflip,
610         float (*root)[3], float (*tip)[3], int *selected, const char **err_str)
611 {
612         LaplacianSystem *sys;
613         MLoopTri *mlooptri;
614         MPoly *mp;
615         MLoop *ml;
616         float solution, weight;
617         int *vertsflipped = NULL, *mask = NULL;
618         int a, tottri, j, bbone, firstsegment, lastsegment;
619         bool use_topology = (me->editflag & ME_EDIT_MIRROR_TOPO) != 0;
620
621         MVert *mvert = me->mvert;
622         bool use_vert_sel = (me->editflag & ME_EDIT_PAINT_VERT_SEL) != 0;
623         bool use_face_sel = (me->editflag & ME_EDIT_PAINT_FACE_SEL) != 0;
624
625         *err_str = NULL;
626
627         /* bone heat needs triangulated faces */
628         tottri = poly_to_tri_count(me->totpoly, me->totloop);
629
630         /* count triangles and create mask */
631         if (ob->mode & OB_MODE_WEIGHT_PAINT &&
632             (use_face_sel || use_vert_sel))
633         {
634                 mask = MEM_callocN(sizeof(int) * me->totvert, "heat_bone_weighting mask");
635
636                 /*  (added selectedVerts content for vertex mask, they used to just equal 1) */
637                 if (use_vert_sel) {
638                         for (a = 0, mp = me->mpoly; a < me->totpoly; mp++, a++) {
639                                 for (j = 0, ml = me->mloop + mp->loopstart; j < mp->totloop; j++, ml++) {
640                                         mask[ml->v] = (mvert[ml->v].flag & SELECT) != 0;
641                                 }
642                         }
643                 }
644                 else if (use_face_sel) {
645                         for (a = 0, mp = me->mpoly; a < me->totpoly; mp++, a++) {
646                                 if (mp->flag & ME_FACE_SEL) {
647                                         for (j = 0, ml = me->mloop + mp->loopstart; j < mp->totloop; j++, ml++) {
648                                                 mask[ml->v] = 1;
649                                         }
650                                 }
651                         }
652                 }
653         }
654
655         /* create laplacian */
656         sys = laplacian_system_construct_begin(me->totvert, tottri, 1);
657
658         sys->heat.tottri = poly_to_tri_count(me->totpoly, me->totloop);
659         mlooptri = MEM_mallocN(sizeof(*sys->heat.mlooptri) * sys->heat.tottri, __func__);
660
661         BKE_mesh_recalc_looptri(
662                 me->mloop, me->mpoly,
663                 me->mvert,
664                 me->totloop, me->totpoly,
665                 mlooptri);
666
667         sys->heat.mlooptri = mlooptri;
668         sys->heat.mloop = me->mloop;
669         sys->heat.totvert = me->totvert;
670         sys->heat.verts = verts;
671         sys->heat.root = root;
672         sys->heat.tip = tip;
673         sys->heat.numsource = numsource;
674
675         heat_ray_tree_create(sys);
676         heat_laplacian_create(sys);
677
678         laplacian_system_construct_end(sys);
679
680         if (dgroupflip) {
681                 vertsflipped = MEM_callocN(sizeof(int) * me->totvert, "vertsflipped");
682                 for (a = 0; a < me->totvert; a++)
683                         vertsflipped[a] = mesh_get_x_mirror_vert(ob, NULL, a, use_topology);
684         }
685
686         /* compute weights per bone */
687         for (j = 0; j < numsource; j++) {
688                 if (!selected[j])
689                         continue;
690
691                 firstsegment = (j == 0 || dgrouplist[j - 1] != dgrouplist[j]);
692                 lastsegment = (j == numsource - 1 || dgrouplist[j] != dgrouplist[j + 1]);
693                 bbone = !(firstsegment && lastsegment);
694
695                 /* clear weights */
696                 if (bbone && firstsegment) {
697                         for (a = 0; a < me->totvert; a++) {
698                                 if (mask && !mask[a])
699                                         continue;
700
701                                 ED_vgroup_vert_remove(ob, dgrouplist[j], a);
702                                 if (vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0)
703                                         ED_vgroup_vert_remove(ob, dgroupflip[j], vertsflipped[a]);
704                         }
705                 }
706
707                 /* fill right hand side */
708                 laplacian_begin_solve(sys, -1);
709
710                 for (a = 0; a < me->totvert; a++)
711                         if (heat_source_closest(sys, a, j))
712                                 laplacian_add_right_hand_side(sys, a,
713                                                               sys->heat.H[a] * sys->heat.p[a]);
714
715                 /* solve */
716                 if (laplacian_system_solve(sys)) {
717                         /* load solution into vertex groups */
718                         for (a = 0; a < me->totvert; a++) {
719                                 if (mask && !mask[a])
720                                         continue;
721
722                                 solution = laplacian_system_get_solution(sys, a);
723
724                                 if (bbone) {
725                                         if (solution > 0.0f)
726                                                 ED_vgroup_vert_add(ob, dgrouplist[j], a, solution,
727                                                                    WEIGHT_ADD);
728                                 }
729                                 else {
730                                         weight = heat_limit_weight(solution);
731                                         if (weight > 0.0f)
732                                                 ED_vgroup_vert_add(ob, dgrouplist[j], a, weight,
733                                                                    WEIGHT_REPLACE);
734                                         else
735                                                 ED_vgroup_vert_remove(ob, dgrouplist[j], a);
736                                 }
737
738                                 /* do same for mirror */
739                                 if (vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0) {
740                                         if (bbone) {
741                                                 if (solution > 0.0f)
742                                                         ED_vgroup_vert_add(ob, dgroupflip[j], vertsflipped[a],
743                                                                            solution, WEIGHT_ADD);
744                                         }
745                                         else {
746                                                 weight = heat_limit_weight(solution);
747                                                 if (weight > 0.0f)
748                                                         ED_vgroup_vert_add(ob, dgroupflip[j], vertsflipped[a],
749                                                                            weight, WEIGHT_REPLACE);
750                                                 else
751                                                         ED_vgroup_vert_remove(ob, dgroupflip[j], vertsflipped[a]);
752                                         }
753                                 }
754                         }
755                 }
756                 else if (*err_str == NULL) {
757                         *err_str = N_("Bone Heat Weighting: failed to find solution for one or more bones");
758                         break;
759                 }
760
761                 /* remove too small vertex weights */
762                 if (bbone && lastsegment) {
763                         for (a = 0; a < me->totvert; a++) {
764                                 if (mask && !mask[a])
765                                         continue;
766
767                                 weight = ED_vgroup_vert_weight(ob, dgrouplist[j], a);
768                                 weight = heat_limit_weight(weight);
769                                 if (weight <= 0.0f)
770                                         ED_vgroup_vert_remove(ob, dgrouplist[j], a);
771
772                                 if (vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0) {
773                                         weight = ED_vgroup_vert_weight(ob, dgroupflip[j], vertsflipped[a]);
774                                         weight = heat_limit_weight(weight);
775                                         if (weight <= 0.0f)
776                                                 ED_vgroup_vert_remove(ob, dgroupflip[j], vertsflipped[a]);
777                                 }
778                         }
779                 }
780         }
781
782         /* free */
783         if (vertsflipped) MEM_freeN(vertsflipped);
784         if (mask) MEM_freeN(mask);
785
786         heat_system_free(sys);
787
788         laplacian_system_delete(sys);
789 }
790
791 /************************** Harmonic Coordinates ****************************/
792 /* From "Harmonic Coordinates for Character Articulation",
793  * Pushkar Joshi, Mark Meyer, Tony DeRose, Brian Green and Tom Sanocki,
794  * SIGGRAPH 2007. */
795
796 #define EPSILON 0.0001f
797
798 #define MESHDEFORM_TAG_UNTYPED  0
799 #define MESHDEFORM_TAG_BOUNDARY 1
800 #define MESHDEFORM_TAG_INTERIOR 2
801 #define MESHDEFORM_TAG_EXTERIOR 3
802
803 /** minimum length for #MDefBoundIsect.len */
804 #define MESHDEFORM_LEN_THRESHOLD 1e-6f
805
806 #define MESHDEFORM_MIN_INFLUENCE 0.0005f
807
808 static const int MESHDEFORM_OFFSET[7][3] = {
809         {0, 0, 0}, {1, 0, 0}, {-1, 0, 0}, {0, 1, 0}, {0, -1, 0}, {0, 0, 1}, {0, 0, -1}
810 };
811
812 typedef struct MDefBoundIsect {
813         /* intersection on the cage 'cagecos' */
814         float co[3];
815         /* non-facing intersections are considered interior */
816         bool facing;
817         /* ray-cast index aligned with MPoly (ray-hit-triangle isn't needed) */
818         int poly_index;
819         /* distance from 'co' to the ray-cast start (clamped to avoid zero division) */
820         float len;
821         /* weights aligned with the MPoly's loop indices */
822         float poly_weights[0];
823 } MDefBoundIsect;
824
825 typedef struct MDefBindInfluence {
826         struct MDefBindInfluence *next;
827         float weight;
828         int vertex;
829 } MDefBindInfluence;
830
831 typedef struct MeshDeformBind {
832         /* grid dimensions */
833         float min[3], max[3];
834         float width[3], halfwidth[3];
835         int size, size3;
836
837         /* meshes */
838         Mesh *cagemesh;
839         float (*cagecos)[3];
840         float (*vertexcos)[3];
841         int totvert, totcagevert;
842
843         /* grids */
844         MemArena *memarena;
845         MDefBoundIsect *(*boundisect)[6];
846         int *semibound;
847         int *tag;
848         float *phi, *totalphi;
849
850         /* mesh stuff */
851         int *inside;
852         float *weights;
853         MDefBindInfluence **dyngrid;
854         float cagemat[4][4];
855
856         /* direct solver */
857         int *varidx;
858
859         BVHTree *bvhtree;
860         BVHTreeFromMesh bvhdata;
861
862         /* avoid DM function calls during intersections */
863         struct {
864                 const MPoly *mpoly;
865                 const MLoop *mloop;
866                 const MLoopTri *looptri;
867                 const float (*poly_nors)[3];
868         } cagemesh_cache;
869 } MeshDeformBind;
870
871 typedef struct MeshDeformIsect {
872         float start[3];
873         float vec[3];
874         float vec_length;
875         float lambda;
876
877         bool isect;
878         float u, v;
879
880 } MeshDeformIsect;
881
882 /* ray intersection */
883
884 struct MeshRayCallbackData {
885         MeshDeformBind *mdb;
886         MeshDeformIsect *isec;
887 };
888
889 static void harmonic_ray_callback(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
890 {
891         struct MeshRayCallbackData *data = userdata;
892         MeshDeformBind *mdb = data->mdb;
893         const MLoop *mloop = mdb->cagemesh_cache.mloop;
894         const MLoopTri *looptri = mdb->cagemesh_cache.looptri, *lt;
895         const float (*poly_nors)[3] = mdb->cagemesh_cache.poly_nors;
896         MeshDeformIsect *isec = data->isec;
897         float no[3], co[3], dist;
898         float *face[3];
899
900         lt = &looptri[index];
901
902         face[0] = mdb->cagecos[mloop[lt->tri[0]].v];
903         face[1] = mdb->cagecos[mloop[lt->tri[1]].v];
904         face[2] = mdb->cagecos[mloop[lt->tri[2]].v];
905
906         if (!isect_ray_tri_watertight_v3(
907                 ray->origin, ray->isect_precalc, UNPACK3(face), &dist, NULL))
908         {
909                 return;
910         }
911
912         if (poly_nors) {
913                 copy_v3_v3(no, poly_nors[lt->poly]);
914         }
915         else {
916                 normal_tri_v3(no, UNPACK3(face));
917         }
918
919         madd_v3_v3v3fl(co, ray->origin, ray->direction, dist);
920         dist /= isec->vec_length;
921         if (dist < hit->dist) {
922                 hit->index = index;
923                 hit->dist = dist;
924                 copy_v3_v3(hit->co, co);
925
926                 isec->isect = (dot_v3v3(no, ray->direction) <= 0.0f);
927                 isec->lambda = dist;
928         }
929 }
930
931 static MDefBoundIsect *meshdeform_ray_tree_intersect(MeshDeformBind *mdb, const float co1[3], const float co2[3])
932 {
933         BVHTreeRayHit hit;
934         MeshDeformIsect isect_mdef;
935         struct MeshRayCallbackData data = {
936                 mdb,
937                 &isect_mdef,
938         };
939         float end[3], vec_normal[3];
940
941         /* happens binding when a cage has no faces */
942         if (UNLIKELY(mdb->bvhtree == NULL))
943                 return NULL;
944
945         /* setup isec */
946         memset(&isect_mdef, 0, sizeof(isect_mdef));
947         isect_mdef.lambda = 1e10f;
948
949         copy_v3_v3(isect_mdef.start, co1);
950         copy_v3_v3(end, co2);
951         sub_v3_v3v3(isect_mdef.vec, end, isect_mdef.start);
952         isect_mdef.vec_length = normalize_v3_v3(vec_normal, isect_mdef.vec);
953
954         hit.index = -1;
955         hit.dist = BVH_RAYCAST_DIST_MAX;
956         if (BLI_bvhtree_ray_cast_ex(mdb->bvhtree, isect_mdef.start, vec_normal,
957                                     0.0, &hit, harmonic_ray_callback, &data, BVH_RAYCAST_WATERTIGHT) != -1)
958         {
959                 const MLoop *mloop = mdb->cagemesh_cache.mloop;
960                 const MLoopTri *lt = &mdb->cagemesh_cache.looptri[hit.index];
961                 const MPoly *mp = &mdb->cagemesh_cache.mpoly[lt->poly];
962                 const float (*cagecos)[3] = mdb->cagecos;
963                 const float len = isect_mdef.lambda;
964                 MDefBoundIsect *isect;
965
966                 float (*mp_cagecos)[3] = BLI_array_alloca(mp_cagecos, mp->totloop);
967                 int i;
968
969                 /* create MDefBoundIsect, and extra for 'poly_weights[]' */
970                 isect = BLI_memarena_alloc(mdb->memarena, sizeof(*isect) + (sizeof(float) * mp->totloop));
971
972                 /* compute intersection coordinate */
973                 madd_v3_v3v3fl(isect->co, co1, isect_mdef.vec, len);
974
975                 isect->facing = isect_mdef.isect;
976
977                 isect->poly_index = lt->poly;
978
979                 isect->len = max_ff(len_v3v3(co1, isect->co), MESHDEFORM_LEN_THRESHOLD);
980
981                 /* compute mean value coordinates for interpolation */
982                 for (i = 0; i < mp->totloop; i++) {
983                         copy_v3_v3(mp_cagecos[i], cagecos[mloop[mp->loopstart + i].v]);
984                 }
985
986                 interp_weights_poly_v3(isect->poly_weights, mp_cagecos, mp->totloop, isect->co);
987
988                 return isect;
989         }
990
991         return NULL;
992 }
993
994 static int meshdeform_inside_cage(MeshDeformBind *mdb, float *co)
995 {
996         MDefBoundIsect *isect;
997         float outside[3], start[3], dir[3];
998         int i;
999
1000         for (i = 1; i <= 6; i++) {
1001                 outside[0] = co[0] + (mdb->max[0] - mdb->min[0] + 1.0f) * MESHDEFORM_OFFSET[i][0];
1002                 outside[1] = co[1] + (mdb->max[1] - mdb->min[1] + 1.0f) * MESHDEFORM_OFFSET[i][1];
1003                 outside[2] = co[2] + (mdb->max[2] - mdb->min[2] + 1.0f) * MESHDEFORM_OFFSET[i][2];
1004
1005                 copy_v3_v3(start, co);
1006                 sub_v3_v3v3(dir, outside, start);
1007                 normalize_v3(dir);
1008
1009                 isect = meshdeform_ray_tree_intersect(mdb, start, outside);
1010                 if (isect && !isect->facing)
1011                         return 1;
1012         }
1013
1014         return 0;
1015 }
1016
1017 /* solving */
1018
1019 BLI_INLINE int meshdeform_index(MeshDeformBind *mdb, int x, int y, int z, int n)
1020 {
1021         int size = mdb->size;
1022
1023         x += MESHDEFORM_OFFSET[n][0];
1024         y += MESHDEFORM_OFFSET[n][1];
1025         z += MESHDEFORM_OFFSET[n][2];
1026
1027         if (x < 0 || x >= mdb->size)
1028                 return -1;
1029         if (y < 0 || y >= mdb->size)
1030                 return -1;
1031         if (z < 0 || z >= mdb->size)
1032                 return -1;
1033
1034         return x + y * size + z * size * size;
1035 }
1036
1037 BLI_INLINE void meshdeform_cell_center(MeshDeformBind *mdb, int x, int y, int z, int n, float *center)
1038 {
1039         x += MESHDEFORM_OFFSET[n][0];
1040         y += MESHDEFORM_OFFSET[n][1];
1041         z += MESHDEFORM_OFFSET[n][2];
1042
1043         center[0] = mdb->min[0] + x * mdb->width[0] + mdb->halfwidth[0];
1044         center[1] = mdb->min[1] + y * mdb->width[1] + mdb->halfwidth[1];
1045         center[2] = mdb->min[2] + z * mdb->width[2] + mdb->halfwidth[2];
1046 }
1047
1048 static void meshdeform_add_intersections(MeshDeformBind *mdb, int x, int y, int z)
1049 {
1050         MDefBoundIsect *isect;
1051         float center[3], ncenter[3];
1052         int i, a;
1053
1054         a = meshdeform_index(mdb, x, y, z, 0);
1055         meshdeform_cell_center(mdb, x, y, z, 0, center);
1056
1057         /* check each outgoing edge for intersection */
1058         for (i = 1; i <= 6; i++) {
1059                 if (meshdeform_index(mdb, x, y, z, i) == -1)
1060                         continue;
1061
1062                 meshdeform_cell_center(mdb, x, y, z, i, ncenter);
1063
1064                 isect = meshdeform_ray_tree_intersect(mdb, center, ncenter);
1065                 if (isect) {
1066                         mdb->boundisect[a][i - 1] = isect;
1067                         mdb->tag[a] = MESHDEFORM_TAG_BOUNDARY;
1068                 }
1069         }
1070 }
1071
1072 static void meshdeform_bind_floodfill(MeshDeformBind *mdb)
1073 {
1074         int *stack, *tag = mdb->tag;
1075         int a, b, i, xyz[3], stacksize, size = mdb->size;
1076
1077         stack = MEM_callocN(sizeof(int) * mdb->size3, "MeshDeformBindStack");
1078
1079         /* we know lower left corner is EXTERIOR because of padding */
1080         tag[0] = MESHDEFORM_TAG_EXTERIOR;
1081         stack[0] = 0;
1082         stacksize = 1;
1083
1084         /* floodfill exterior tag */
1085         while (stacksize > 0) {
1086                 a = stack[--stacksize];
1087
1088                 xyz[2] = a / (size * size);
1089                 xyz[1] = (a - xyz[2] * size * size) / size;
1090                 xyz[0] = a - xyz[1] * size - xyz[2] * size * size;
1091
1092                 for (i = 1; i <= 6; i++) {
1093                         b = meshdeform_index(mdb, xyz[0], xyz[1], xyz[2], i);
1094
1095                         if (b != -1) {
1096                                 if (tag[b] == MESHDEFORM_TAG_UNTYPED ||
1097                                     (tag[b] == MESHDEFORM_TAG_BOUNDARY && !mdb->boundisect[a][i - 1]))
1098                                 {
1099                                         tag[b] = MESHDEFORM_TAG_EXTERIOR;
1100                                         stack[stacksize++] = b;
1101                                 }
1102                         }
1103                 }
1104         }
1105
1106         /* other cells are interior */
1107         for (a = 0; a < size * size * size; a++)
1108                 if (tag[a] == MESHDEFORM_TAG_UNTYPED)
1109                         tag[a] = MESHDEFORM_TAG_INTERIOR;
1110
1111 #if 0
1112         {
1113                 int tb, ti, te, ts;
1114                 tb = ti = te = ts = 0;
1115                 for (a = 0; a < size * size * size; a++)
1116                         if (tag[a] == MESHDEFORM_TAG_BOUNDARY)
1117                                 tb++;
1118                         else if (tag[a] == MESHDEFORM_TAG_INTERIOR)
1119                                 ti++;
1120                         else if (tag[a] == MESHDEFORM_TAG_EXTERIOR) {
1121                                 te++;
1122
1123                                 if (mdb->semibound[a])
1124                                         ts++;
1125                         }
1126
1127                 printf("interior %d exterior %d boundary %d semi-boundary %d\n", ti, te, tb, ts);
1128         }
1129 #endif
1130
1131         MEM_freeN(stack);
1132 }
1133
1134 static float meshdeform_boundary_phi(const MeshDeformBind *mdb, const MDefBoundIsect *isect, int cagevert)
1135 {
1136         const MLoop *mloop = mdb->cagemesh_cache.mloop;
1137         const MPoly *mp = &mdb->cagemesh_cache.mpoly[isect->poly_index];
1138         int i;
1139
1140         for (i = 0; i < mp->totloop; i++) {
1141                 if (mloop[mp->loopstart + i].v == cagevert) {
1142                         return isect->poly_weights[i];
1143                 }
1144         }
1145
1146         return 0.0f;
1147 }
1148
1149 static float meshdeform_interp_w(MeshDeformBind *mdb, float *gridvec, float *UNUSED(vec), int UNUSED(cagevert))
1150 {
1151         float dvec[3], ivec[3], wx, wy, wz, result = 0.0f;
1152         float weight, totweight = 0.0f;
1153         int i, a, x, y, z;
1154
1155         for (i = 0; i < 3; i++) {
1156                 ivec[i] = (int)gridvec[i];
1157                 dvec[i] = gridvec[i] - ivec[i];
1158         }
1159
1160         for (i = 0; i < 8; i++) {
1161                 if (i & 1) { x = ivec[0] + 1; wx = dvec[0]; }
1162                 else       { x = ivec[0];     wx = 1.0f - dvec[0]; }
1163
1164                 if (i & 2) { y = ivec[1] + 1; wy = dvec[1]; }
1165                 else       { y = ivec[1];     wy = 1.0f - dvec[1]; }
1166
1167                 if (i & 4) { z = ivec[2] + 1; wz = dvec[2]; }
1168                 else       { z = ivec[2];     wz = 1.0f - dvec[2]; }
1169
1170                 CLAMP(x, 0, mdb->size - 1);
1171                 CLAMP(y, 0, mdb->size - 1);
1172                 CLAMP(z, 0, mdb->size - 1);
1173
1174                 a = meshdeform_index(mdb, x, y, z, 0);
1175                 weight = wx * wy * wz;
1176                 result += weight * mdb->phi[a];
1177                 totweight += weight;
1178         }
1179
1180         if (totweight > 0.0f)
1181                 result /= totweight;
1182
1183         return result;
1184 }
1185
1186 static void meshdeform_check_semibound(MeshDeformBind *mdb, int x, int y, int z)
1187 {
1188         int i, a;
1189
1190         a = meshdeform_index(mdb, x, y, z, 0);
1191         if (mdb->tag[a] != MESHDEFORM_TAG_EXTERIOR)
1192                 return;
1193
1194         for (i = 1; i <= 6; i++)
1195                 if (mdb->boundisect[a][i - 1])
1196                         mdb->semibound[a] = 1;
1197 }
1198
1199 static float meshdeform_boundary_total_weight(MeshDeformBind *mdb, int x, int y, int z)
1200 {
1201         float weight, totweight = 0.0f;
1202         int i, a;
1203
1204         a = meshdeform_index(mdb, x, y, z, 0);
1205
1206         /* count weight for neighbor cells */
1207         for (i = 1; i <= 6; i++) {
1208                 if (meshdeform_index(mdb, x, y, z, i) == -1)
1209                         continue;
1210
1211                 if (mdb->boundisect[a][i - 1])
1212                         weight = 1.0f / mdb->boundisect[a][i - 1]->len;
1213                 else if (!mdb->semibound[a])
1214                         weight = 1.0f / mdb->width[0];
1215                 else
1216                         weight = 0.0f;
1217
1218                 totweight += weight;
1219         }
1220
1221         return totweight;
1222 }
1223
1224 static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, LinearSolver *context, int x, int y, int z)
1225 {
1226         MDefBoundIsect *isect;
1227         float weight, totweight;
1228         int i, a, acenter;
1229
1230         acenter = meshdeform_index(mdb, x, y, z, 0);
1231         if (mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR)
1232                 return;
1233
1234         EIG_linear_solver_matrix_add(context, mdb->varidx[acenter], mdb->varidx[acenter], 1.0f);
1235
1236         totweight = meshdeform_boundary_total_weight(mdb, x, y, z);
1237         for (i = 1; i <= 6; i++) {
1238                 a = meshdeform_index(mdb, x, y, z, i);
1239                 if (a == -1 || mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR)
1240                         continue;
1241
1242                 isect = mdb->boundisect[acenter][i - 1];
1243                 if (!isect) {
1244                         weight = (1.0f / mdb->width[0]) / totweight;
1245                         EIG_linear_solver_matrix_add(context, mdb->varidx[acenter], mdb->varidx[a], -weight);
1246                 }
1247         }
1248 }
1249
1250 static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, LinearSolver *context, int x, int y, int z, int cagevert)
1251 {
1252         MDefBoundIsect *isect;
1253         float rhs, weight, totweight;
1254         int i, a, acenter;
1255
1256         acenter = meshdeform_index(mdb, x, y, z, 0);
1257         if (mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR)
1258                 return;
1259
1260         totweight = meshdeform_boundary_total_weight(mdb, x, y, z);
1261         for (i = 1; i <= 6; i++) {
1262                 a = meshdeform_index(mdb, x, y, z, i);
1263                 if (a == -1)
1264                         continue;
1265
1266                 isect = mdb->boundisect[acenter][i - 1];
1267
1268                 if (isect) {
1269                         weight = (1.0f / isect->len) / totweight;
1270                         rhs = weight * meshdeform_boundary_phi(mdb, isect, cagevert);
1271                         EIG_linear_solver_right_hand_side_add(context, 0, mdb->varidx[acenter], rhs);
1272                 }
1273         }
1274 }
1275
1276 static void meshdeform_matrix_add_semibound_phi(MeshDeformBind *mdb, int x, int y, int z, int cagevert)
1277 {
1278         MDefBoundIsect *isect;
1279         float rhs, weight, totweight;
1280         int i, a;
1281
1282         a = meshdeform_index(mdb, x, y, z, 0);
1283         if (!mdb->semibound[a])
1284                 return;
1285
1286         mdb->phi[a] = 0.0f;
1287
1288         totweight = meshdeform_boundary_total_weight(mdb, x, y, z);
1289         for (i = 1; i <= 6; i++) {
1290                 isect = mdb->boundisect[a][i - 1];
1291
1292                 if (isect) {
1293                         weight = (1.0f / isect->len) / totweight;
1294                         rhs = weight * meshdeform_boundary_phi(mdb, isect, cagevert);
1295                         mdb->phi[a] += rhs;
1296                 }
1297         }
1298 }
1299
1300 static void meshdeform_matrix_add_exterior_phi(MeshDeformBind *mdb, int x, int y, int z, int UNUSED(cagevert))
1301 {
1302         float phi, totweight;
1303         int i, a, acenter;
1304
1305         acenter = meshdeform_index(mdb, x, y, z, 0);
1306         if (mdb->tag[acenter] != MESHDEFORM_TAG_EXTERIOR || mdb->semibound[acenter])
1307                 return;
1308
1309         phi = 0.0f;
1310         totweight = 0.0f;
1311         for (i = 1; i <= 6; i++) {
1312                 a = meshdeform_index(mdb, x, y, z, i);
1313
1314                 if (a != -1 && mdb->semibound[a]) {
1315                         phi += mdb->phi[a];
1316                         totweight += 1.0f;
1317                 }
1318         }
1319
1320         if (totweight != 0.0f)
1321                 mdb->phi[acenter] = phi / totweight;
1322 }
1323
1324 static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind *mdb)
1325 {
1326         LinearSolver *context;
1327         float vec[3], gridvec[3];
1328         int a, b, x, y, z, totvar;
1329         char message[256];
1330
1331         /* setup variable indices */
1332         mdb->varidx = MEM_callocN(sizeof(int) * mdb->size3, "MeshDeformDSvaridx");
1333         for (a = 0, totvar = 0; a < mdb->size3; a++)
1334                 mdb->varidx[a] = (mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR) ? -1 : totvar++;
1335
1336         if (totvar == 0) {
1337                 MEM_freeN(mdb->varidx);
1338                 return;
1339         }
1340
1341         progress_bar(0, "Starting mesh deform solve");
1342
1343         /* setup linear solver */
1344         context = EIG_linear_solver_new(totvar, totvar, 1);
1345
1346         /* build matrix */
1347         for (z = 0; z < mdb->size; z++)
1348                 for (y = 0; y < mdb->size; y++)
1349                         for (x = 0; x < mdb->size; x++)
1350                                 meshdeform_matrix_add_cell(mdb, context, x, y, z);
1351
1352         /* solve for each cage vert */
1353         for (a = 0; a < mdb->totcagevert; a++) {
1354                 /* fill in right hand side and solve */
1355                 for (z = 0; z < mdb->size; z++)
1356                         for (y = 0; y < mdb->size; y++)
1357                                 for (x = 0; x < mdb->size; x++)
1358                                         meshdeform_matrix_add_rhs(mdb, context, x, y, z, a);
1359
1360                 if (EIG_linear_solver_solve(context)) {
1361                         for (z = 0; z < mdb->size; z++)
1362                                 for (y = 0; y < mdb->size; y++)
1363                                         for (x = 0; x < mdb->size; x++)
1364                                                 meshdeform_matrix_add_semibound_phi(mdb, x, y, z, a);
1365
1366                         for (z = 0; z < mdb->size; z++)
1367                                 for (y = 0; y < mdb->size; y++)
1368                                         for (x = 0; x < mdb->size; x++)
1369                                                 meshdeform_matrix_add_exterior_phi(mdb, x, y, z, a);
1370
1371                         for (b = 0; b < mdb->size3; b++) {
1372                                 if (mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR)
1373                                         mdb->phi[b] = EIG_linear_solver_variable_get(context, 0, mdb->varidx[b]);
1374                                 mdb->totalphi[b] += mdb->phi[b];
1375                         }
1376
1377                         if (mdb->weights) {
1378                                 /* static bind : compute weights for each vertex */
1379                                 for (b = 0; b < mdb->totvert; b++) {
1380                                         if (mdb->inside[b]) {
1381                                                 copy_v3_v3(vec, mdb->vertexcos[b]);
1382                                                 gridvec[0] = (vec[0] - mdb->min[0] - mdb->halfwidth[0]) / mdb->width[0];
1383                                                 gridvec[1] = (vec[1] - mdb->min[1] - mdb->halfwidth[1]) / mdb->width[1];
1384                                                 gridvec[2] = (vec[2] - mdb->min[2] - mdb->halfwidth[2]) / mdb->width[2];
1385
1386                                                 mdb->weights[b * mdb->totcagevert + a] = meshdeform_interp_w(mdb, gridvec, vec, a);
1387                                         }
1388                                 }
1389                         }
1390                         else {
1391                                 MDefBindInfluence *inf;
1392
1393                                 /* dynamic bind */
1394                                 for (b = 0; b < mdb->size3; b++) {
1395                                         if (mdb->phi[b] >= MESHDEFORM_MIN_INFLUENCE) {
1396                                                 inf = BLI_memarena_alloc(mdb->memarena, sizeof(*inf));
1397                                                 inf->vertex = a;
1398                                                 inf->weight = mdb->phi[b];
1399                                                 inf->next = mdb->dyngrid[b];
1400                                                 mdb->dyngrid[b] = inf;
1401                                         }
1402                                 }
1403                         }
1404                 }
1405                 else {
1406                         modifier_setError(&mmd->modifier, "Failed to find bind solution (increase precision?)");
1407                         error("Mesh Deform: failed to find bind solution.");
1408                         break;
1409                 }
1410
1411                 BLI_snprintf(message, sizeof(message), "Mesh deform solve %d / %d       |||", a + 1, mdb->totcagevert);
1412                 progress_bar((float)(a + 1) / (float)(mdb->totcagevert), message);
1413         }
1414
1415 #if 0
1416         /* sanity check */
1417         for (b = 0; b < mdb->size3; b++)
1418                 if (mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR)
1419                         if (fabsf(mdb->totalphi[b] - 1.0f) > 1e-4f)
1420                                 printf("totalphi deficiency [%s|%d] %d: %.10f\n",
1421                                        (mdb->tag[b] == MESHDEFORM_TAG_INTERIOR) ? "interior" : "boundary", mdb->semibound[b], mdb->varidx[b], mdb->totalphi[b]);
1422 #endif
1423
1424         /* free */
1425         MEM_freeN(mdb->varidx);
1426
1427         EIG_linear_solver_delete(context);
1428 }
1429
1430 static void harmonic_coordinates_bind(MeshDeformModifierData *mmd, MeshDeformBind *mdb)
1431 {
1432         MDefBindInfluence *inf;
1433         MDefInfluence *mdinf;
1434         MDefCell *cell;
1435         float center[3], vec[3], maxwidth, totweight;
1436         int a, b, x, y, z, totinside, offset;
1437
1438         /* compute bounding box of the cage mesh */
1439         INIT_MINMAX(mdb->min, mdb->max);
1440
1441         for (a = 0; a < mdb->totcagevert; a++)
1442                 minmax_v3v3_v3(mdb->min, mdb->max, mdb->cagecos[a]);
1443
1444         /* allocate memory */
1445         mdb->size = (2 << (mmd->gridsize - 1)) + 2;
1446         mdb->size3 = mdb->size * mdb->size * mdb->size;
1447         mdb->tag = MEM_callocN(sizeof(int) * mdb->size3, "MeshDeformBindTag");
1448         mdb->phi = MEM_callocN(sizeof(float) * mdb->size3, "MeshDeformBindPhi");
1449         mdb->totalphi = MEM_callocN(sizeof(float) * mdb->size3, "MeshDeformBindTotalPhi");
1450         mdb->boundisect = MEM_callocN(sizeof(*mdb->boundisect) * mdb->size3, "MDefBoundIsect");
1451         mdb->semibound = MEM_callocN(sizeof(int) * mdb->size3, "MDefSemiBound");
1452         mdb->bvhtree = BKE_bvhtree_from_mesh_get(&mdb->bvhdata, mdb->cagemesh, BVHTREE_FROM_LOOPTRI, 4);
1453         mdb->inside = MEM_callocN(sizeof(int) * mdb->totvert, "MDefInside");
1454
1455         if (mmd->flag & MOD_MDEF_DYNAMIC_BIND)
1456                 mdb->dyngrid = MEM_callocN(sizeof(MDefBindInfluence *) * mdb->size3, "MDefDynGrid");
1457         else
1458                 mdb->weights = MEM_callocN(sizeof(float) * mdb->totvert * mdb->totcagevert, "MDefWeights");
1459
1460         mdb->memarena = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "harmonic coords arena");
1461         BLI_memarena_use_calloc(mdb->memarena);
1462
1463         /* initialize data from 'cagedm' for reuse */
1464         {
1465                 Mesh *me = mdb->cagemesh;
1466                 mdb->cagemesh_cache.mpoly = me->mpoly;
1467                 mdb->cagemesh_cache.mloop = me->mloop;
1468                 mdb->cagemesh_cache.looptri = BKE_mesh_runtime_looptri_ensure(me);
1469                 mdb->cagemesh_cache.poly_nors = CustomData_get_layer(&me->pdata, CD_NORMAL);  /* can be NULL */
1470         }
1471
1472         /* make bounding box equal size in all directions, add padding, and compute
1473          * width of the cells */
1474         maxwidth = -1.0f;
1475         for (a = 0; a < 3; a++)
1476                 if (mdb->max[a] - mdb->min[a] > maxwidth)
1477                         maxwidth = mdb->max[a] - mdb->min[a];
1478
1479         for (a = 0; a < 3; a++) {
1480                 center[a] = (mdb->min[a] + mdb->max[a]) * 0.5f;
1481                 mdb->min[a] = center[a] - maxwidth * 0.5f;
1482                 mdb->max[a] = center[a] + maxwidth * 0.5f;
1483
1484                 mdb->width[a] = (mdb->max[a] - mdb->min[a]) / (mdb->size - 4);
1485                 mdb->min[a] -= 2.1f * mdb->width[a];
1486                 mdb->max[a] += 2.1f * mdb->width[a];
1487
1488                 mdb->width[a] = (mdb->max[a] - mdb->min[a]) / mdb->size;
1489                 mdb->halfwidth[a] = mdb->width[a] * 0.5f;
1490         }
1491
1492         progress_bar(0, "Setting up mesh deform system");
1493
1494         totinside = 0;
1495         for (a = 0; a < mdb->totvert; a++) {
1496                 copy_v3_v3(vec, mdb->vertexcos[a]);
1497                 mdb->inside[a] = meshdeform_inside_cage(mdb, vec);
1498                 if (mdb->inside[a])
1499                         totinside++;
1500         }
1501
1502         /* free temporary MDefBoundIsects */
1503         BLI_memarena_free(mdb->memarena);
1504         mdb->memarena = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "harmonic coords arena");
1505
1506         /* start with all cells untyped */
1507         for (a = 0; a < mdb->size3; a++)
1508                 mdb->tag[a] = MESHDEFORM_TAG_UNTYPED;
1509
1510         /* detect intersections and tag boundary cells */
1511         for (z = 0; z < mdb->size; z++)
1512                 for (y = 0; y < mdb->size; y++)
1513                         for (x = 0; x < mdb->size; x++)
1514                                 meshdeform_add_intersections(mdb, x, y, z);
1515
1516         /* compute exterior and interior tags */
1517         meshdeform_bind_floodfill(mdb);
1518
1519         for (z = 0; z < mdb->size; z++)
1520                 for (y = 0; y < mdb->size; y++)
1521                         for (x = 0; x < mdb->size; x++)
1522                                 meshdeform_check_semibound(mdb, x, y, z);
1523
1524         /* solve */
1525         meshdeform_matrix_solve(mmd, mdb);
1526
1527         /* assign results */
1528         if (mmd->flag & MOD_MDEF_DYNAMIC_BIND) {
1529                 mmd->totinfluence = 0;
1530                 for (a = 0; a < mdb->size3; a++)
1531                         for (inf = mdb->dyngrid[a]; inf; inf = inf->next)
1532                                 mmd->totinfluence++;
1533
1534                 /* convert MDefBindInfluences to smaller MDefInfluences */
1535                 mmd->dyngrid = MEM_callocN(sizeof(MDefCell) * mdb->size3, "MDefDynGrid");
1536                 mmd->dyninfluences = MEM_callocN(sizeof(MDefInfluence) * mmd->totinfluence, "MDefInfluence");
1537                 offset = 0;
1538                 for (a = 0; a < mdb->size3; a++) {
1539                         cell = &mmd->dyngrid[a];
1540                         cell->offset = offset;
1541
1542                         totweight = 0.0f;
1543                         mdinf = mmd->dyninfluences + cell->offset;
1544                         for (inf = mdb->dyngrid[a]; inf; inf = inf->next, mdinf++) {
1545                                 mdinf->weight = inf->weight;
1546                                 mdinf->vertex = inf->vertex;
1547                                 totweight += mdinf->weight;
1548                                 cell->totinfluence++;
1549                         }
1550
1551                         if (totweight > 0.0f) {
1552                                 mdinf = mmd->dyninfluences + cell->offset;
1553                                 for (b = 0; b < cell->totinfluence; b++, mdinf++)
1554                                         mdinf->weight /= totweight;
1555                         }
1556
1557                         offset += cell->totinfluence;
1558                 }
1559
1560                 mmd->dynverts = mdb->inside;
1561                 mmd->dyngridsize = mdb->size;
1562                 copy_v3_v3(mmd->dyncellmin, mdb->min);
1563                 mmd->dyncellwidth = mdb->width[0];
1564                 MEM_freeN(mdb->dyngrid);
1565         }
1566         else {
1567                 mmd->bindweights = mdb->weights;
1568                 MEM_freeN(mdb->inside);
1569         }
1570
1571         MEM_freeN(mdb->tag);
1572         MEM_freeN(mdb->phi);
1573         MEM_freeN(mdb->totalphi);
1574         MEM_freeN(mdb->boundisect);
1575         MEM_freeN(mdb->semibound);
1576         BLI_memarena_free(mdb->memarena);
1577         free_bvhtree_from_mesh(&mdb->bvhdata);
1578 }
1579
1580 void ED_mesh_deform_bind_callback(
1581         MeshDeformModifierData *mmd, Mesh *cagemesh,
1582         float *vertexcos, int totvert, float cagemat[4][4])
1583 {
1584         MeshDeformBind mdb;
1585         MVert *mvert;
1586         int a;
1587
1588         waitcursor(1);
1589         start_progress_bar();
1590
1591         memset(&mdb, 0, sizeof(MeshDeformBind));
1592
1593         /* get mesh and cage mesh */
1594         mdb.vertexcos = MEM_callocN(sizeof(float) * 3 * totvert, "MeshDeformCos");
1595         mdb.totvert = totvert;
1596
1597         mdb.cagemesh = cagemesh;
1598         mdb.totcagevert = mdb.cagemesh->totvert;
1599         mdb.cagecos = MEM_callocN(sizeof(*mdb.cagecos) * mdb.totcagevert, "MeshDeformBindCos");
1600         copy_m4_m4(mdb.cagemat, cagemat);
1601
1602         mvert = mdb.cagemesh->mvert;
1603         for (a = 0; a < mdb.totcagevert; a++)
1604                 copy_v3_v3(mdb.cagecos[a], mvert[a].co);
1605         for (a = 0; a < mdb.totvert; a++)
1606                 mul_v3_m4v3(mdb.vertexcos[a], mdb.cagemat, vertexcos + a * 3);
1607
1608         /* solve */
1609         harmonic_coordinates_bind(mmd, &mdb);
1610
1611         /* assign bind variables */
1612         mmd->bindcagecos = (float *)mdb.cagecos;
1613         mmd->totvert = mdb.totvert;
1614         mmd->totcagevert = mdb.totcagevert;
1615         copy_m4_m4(mmd->bindmat, mmd->object->obmat);
1616
1617         /* transform bindcagecos to world space */
1618         for (a = 0; a < mdb.totcagevert; a++)
1619                 mul_m4_v3(mmd->object->obmat, mmd->bindcagecos + a * 3);
1620
1621         /* free */
1622         MEM_freeN(mdb.vertexcos);
1623
1624         /* compact weights */
1625         modifier_mdef_compact_influences((ModifierData *)mmd);
1626
1627         end_progress_bar();
1628         waitcursor(0);
1629 }