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