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_mesh_types.h"
  25 #include "DNA_meshdata_types.h"
  26 #include "DNA_object_types.h"
  27 #include "DNA_scene_types.h"
  28
  29 #include "BLI_alloca.h"
  30 #include "BLI_edgehash.h"
  31 #include "BLI_math.h"
  32 #include "BLI_memarena.h"
  33 #include "BLI_string.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_armature.h"
  43 #include "ED_mesh.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       }
1205       else if (tag[a] == MESHDEFORM_TAG_INTERIOR) {
1206         ti++;
1207       }
1208       else if (tag[a] == MESHDEFORM_TAG_EXTERIOR) {
1209         te++;
1210
1211         if (mdb->semibound[a]) {
1212           ts++;
1213         }
1214       }
1215     }
1216
1217     printf("interior %d exterior %d boundary %d semi-boundary %d\n", ti, te, tb, ts);
1218   }
1219 #endif
1220
1221   MEM_freeN(stack);
1222 }
1223
1224 static float meshdeform_boundary_phi(const MeshDeformBind *mdb,
1225                                      const MDefBoundIsect *isect,
1226                                      int cagevert)
1227 {
1228   const MLoop *mloop = mdb->cagemesh_cache.mloop;
1229   const MPoly *mp = &mdb->cagemesh_cache.mpoly[isect->poly_index];
1230   int i;
1231
1232   for (i = 0; i < mp->totloop; i++) {
1233     if (mloop[mp->loopstart + i].v == cagevert) {
1234       return isect->poly_weights[i];
1235     }
1236   }
1237
1238   return 0.0f;
1239 }
1240
1241 static float meshdeform_interp_w(MeshDeformBind *mdb,
1242                                  float *gridvec,
1243                                  float *UNUSED(vec),
1244                                  int UNUSED(cagevert))
1245 {
1246   float dvec[3], ivec[3], wx, wy, wz, result = 0.0f;
1247   float weight, totweight = 0.0f;
1248   int i, a, x, y, z;
1249
1250   for (i = 0; i < 3; i++) {
1251     ivec[i] = (int)gridvec[i];
1252     dvec[i] = gridvec[i] - ivec[i];
1253   }
1254
1255   for (i = 0; i < 8; i++) {
1256     if (i & 1) {
1257       x = ivec[0] + 1;
1258       wx = dvec[0];
1259     }
1260     else {
1261       x = ivec[0];
1262       wx = 1.0f - dvec[0];
1263     }
1264
1265     if (i & 2) {
1266       y = ivec[1] + 1;
1267       wy = dvec[1];
1268     }
1269     else {
1270       y = ivec[1];
1271       wy = 1.0f - dvec[1];
1272     }
1273
1274     if (i & 4) {
1275       z = ivec[2] + 1;
1276       wz = dvec[2];
1277     }
1278     else {
1279       z = ivec[2];
1280       wz = 1.0f - dvec[2];
1281     }
1282
1283     CLAMP(x, 0, mdb->size - 1);
1284     CLAMP(y, 0, mdb->size - 1);
1285     CLAMP(z, 0, mdb->size - 1);
1286
1287     a = meshdeform_index(mdb, x, y, z, 0);
1288     weight = wx * wy * wz;
1289     result += weight * mdb->phi[a];
1290     totweight += weight;
1291   }
1292
1293   if (totweight > 0.0f) {
1294     result /= totweight;
1295   }
1296
1297   return result;
1298 }
1299
1300 static void meshdeform_check_semibound(MeshDeformBind *mdb, int x, int y, int z)
1301 {
1302   int i, a;
1303
1304   a = meshdeform_index(mdb, x, y, z, 0);
1305   if (mdb->tag[a] != MESHDEFORM_TAG_EXTERIOR) {
1306     return;
1307   }
1308
1309   for (i = 1; i <= 6; i++) {
1310     if (mdb->boundisect[a][i - 1]) {
1311       mdb->semibound[a] = 1;
1312     }
1313   }
1314 }
1315
1316 static float meshdeform_boundary_total_weight(MeshDeformBind *mdb, int x, int y, int z)
1317 {
1318   float weight, totweight = 0.0f;
1319   int i, a;
1320
1321   a = meshdeform_index(mdb, x, y, z, 0);
1322
1323   /* count weight for neighbor cells */
1324   for (i = 1; i <= 6; i++) {
1325     if (meshdeform_index(mdb, x, y, z, i) == -1) {
1326       continue;
1327     }
1328
1329     if (mdb->boundisect[a][i - 1]) {
1330       weight = 1.0f / mdb->boundisect[a][i - 1]->len;
1331     }
1332     else if (!mdb->semibound[a]) {
1333       weight = 1.0f / mdb->width[0];
1334     }
1335     else {
1336       weight = 0.0f;
1337     }
1338
1339     totweight += weight;
1340   }
1341
1342   return totweight;
1343 }
1344
1345 static void meshdeform_matrix_add_cell(
1346     MeshDeformBind *mdb, LinearSolver *context, int x, int y, int z)
1347 {
1348   MDefBoundIsect *isect;
1349   float weight, totweight;
1350   int i, a, acenter;
1351
1352   acenter = meshdeform_index(mdb, x, y, z, 0);
1353   if (mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR) {
1354     return;
1355   }
1356
1357   EIG_linear_solver_matrix_add(context, mdb->varidx[acenter], mdb->varidx[acenter], 1.0f);
1358
1359   totweight = meshdeform_boundary_total_weight(mdb, x, y, z);
1360   for (i = 1; i <= 6; i++) {
1361     a = meshdeform_index(mdb, x, y, z, i);
1362     if (a == -1 || mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR) {
1363       continue;
1364     }
1365
1366     isect = mdb->boundisect[acenter][i - 1];
1367     if (!isect) {
1368       weight = (1.0f / mdb->width[0]) / totweight;
1369       EIG_linear_solver_matrix_add(context, mdb->varidx[acenter], mdb->varidx[a], -weight);
1370     }
1371   }
1372 }
1373
1374 static void meshdeform_matrix_add_rhs(
1375     MeshDeformBind *mdb, LinearSolver *context, int x, int y, int z, int cagevert)
1376 {
1377   MDefBoundIsect *isect;
1378   float rhs, weight, totweight;
1379   int i, a, acenter;
1380
1381   acenter = meshdeform_index(mdb, x, y, z, 0);
1382   if (mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR) {
1383     return;
1384   }
1385
1386   totweight = meshdeform_boundary_total_weight(mdb, x, y, z);
1387   for (i = 1; i <= 6; i++) {
1388     a = meshdeform_index(mdb, x, y, z, i);
1389     if (a == -1) {
1390       continue;
1391     }
1392
1393     isect = mdb->boundisect[acenter][i - 1];
1394
1395     if (isect) {
1396       weight = (1.0f / isect->len) / totweight;
1397       rhs = weight * meshdeform_boundary_phi(mdb, isect, cagevert);
1398       EIG_linear_solver_right_hand_side_add(context, 0, mdb->varidx[acenter], rhs);
1399     }
1400   }
1401 }
1402
1403 static void meshdeform_matrix_add_semibound_phi(
1404     MeshDeformBind *mdb, int x, int y, int z, int cagevert)
1405 {
1406   MDefBoundIsect *isect;
1407   float rhs, weight, totweight;
1408   int i, a;
1409
1410   a = meshdeform_index(mdb, x, y, z, 0);
1411   if (!mdb->semibound[a]) {
1412     return;
1413   }
1414
1415   mdb->phi[a] = 0.0f;
1416
1417   totweight = meshdeform_boundary_total_weight(mdb, x, y, z);
1418   for (i = 1; i <= 6; i++) {
1419     isect = mdb->boundisect[a][i - 1];
1420
1421     if (isect) {
1422       weight = (1.0f / isect->len) / totweight;
1423       rhs = weight * meshdeform_boundary_phi(mdb, isect, cagevert);
1424       mdb->phi[a] += rhs;
1425     }
1426   }
1427 }
1428
1429 static void meshdeform_matrix_add_exterior_phi(
1430     MeshDeformBind *mdb, int x, int y, int z, int UNUSED(cagevert))
1431 {
1432   float phi, totweight;
1433   int i, a, acenter;
1434
1435   acenter = meshdeform_index(mdb, x, y, z, 0);
1436   if (mdb->tag[acenter] != MESHDEFORM_TAG_EXTERIOR || mdb->semibound[acenter]) {
1437     return;
1438   }
1439
1440   phi = 0.0f;
1441   totweight = 0.0f;
1442   for (i = 1; i <= 6; i++) {
1443     a = meshdeform_index(mdb, x, y, z, i);
1444
1445     if (a != -1 && mdb->semibound[a]) {
1446       phi += mdb->phi[a];
1447       totweight += 1.0f;
1448     }
1449   }
1450
1451   if (totweight != 0.0f) {
1452     mdb->phi[acenter] = phi / totweight;
1453   }
1454 }
1455
1456 static void meshdeform_matrix_solve(MeshDeformModifierData *mmd, MeshDeformBind *mdb)
1457 {
1458   LinearSolver *context;
1459   float vec[3], gridvec[3];
1460   int a, b, x, y, z, totvar;
1461   char message[256];
1462
1463   /* setup variable indices */
1464   mdb->varidx = MEM_callocN(sizeof(int) * mdb->size3, "MeshDeformDSvaridx");
1465   for (a = 0, totvar = 0; a < mdb->size3; a++) {
1466     mdb->varidx[a] = (mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR) ? -1 : totvar++;
1467   }
1468
1469   if (totvar == 0) {
1470     MEM_freeN(mdb->varidx);
1471     return;
1472   }
1473
1474   progress_bar(0, "Starting mesh deform solve");
1475
1476   /* setup linear solver */
1477   context = EIG_linear_solver_new(totvar, totvar, 1);
1478
1479   /* build matrix */
1480   for (z = 0; z < mdb->size; z++) {
1481     for (y = 0; y < mdb->size; y++) {
1482       for (x = 0; x < mdb->size; x++) {
1483         meshdeform_matrix_add_cell(mdb, context, x, y, z);
1484       }
1485     }
1486   }
1487
1488   /* solve for each cage vert */
1489   for (a = 0; a < mdb->totcagevert; a++) {
1490     /* fill in right hand side and solve */
1491     for (z = 0; z < mdb->size; z++) {
1492       for (y = 0; y < mdb->size; y++) {
1493         for (x = 0; x < mdb->size; x++) {
1494           meshdeform_matrix_add_rhs(mdb, context, x, y, z, a);
1495         }
1496       }
1497     }
1498
1499     if (EIG_linear_solver_solve(context)) {
1500       for (z = 0; z < mdb->size; z++) {
1501         for (y = 0; y < mdb->size; y++) {
1502           for (x = 0; x < mdb->size; x++) {
1503             meshdeform_matrix_add_semibound_phi(mdb, x, y, z, a);
1504           }
1505         }
1506       }
1507
1508       for (z = 0; z < mdb->size; z++) {
1509         for (y = 0; y < mdb->size; y++) {
1510           for (x = 0; x < mdb->size; x++) {
1511             meshdeform_matrix_add_exterior_phi(mdb, x, y, z, a);
1512           }
1513         }
1514       }
1515
1516       for (b = 0; b < mdb->size3; b++) {
1517         if (mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR) {
1518           mdb->phi[b] = EIG_linear_solver_variable_get(context, 0, mdb->varidx[b]);
1519         }
1520         mdb->totalphi[b] += mdb->phi[b];
1521       }
1522
1523       if (mdb->weights) {
1524         /* static bind : compute weights for each vertex */
1525         for (b = 0; b < mdb->totvert; b++) {
1526           if (mdb->inside[b]) {
1527             copy_v3_v3(vec, mdb->vertexcos[b]);
1528             gridvec[0] = (vec[0] - mdb->min[0] - mdb->halfwidth[0]) / mdb->width[0];
1529             gridvec[1] = (vec[1] - mdb->min[1] - mdb->halfwidth[1]) / mdb->width[1];
1530             gridvec[2] = (vec[2] - mdb->min[2] - mdb->halfwidth[2]) / mdb->width[2];
1531
1532             mdb->weights[b * mdb->totcagevert + a] = meshdeform_interp_w(mdb, gridvec, vec, a);
1533           }
1534         }
1535       }
1536       else {
1537         MDefBindInfluence *inf;
1538
1539         /* dynamic bind */
1540         for (b = 0; b < mdb->size3; b++) {
1541           if (mdb->phi[b] >= MESHDEFORM_MIN_INFLUENCE) {
1542             inf = BLI_memarena_alloc(mdb->memarena, sizeof(*inf));
1543             inf->vertex = a;
1544             inf->weight = mdb->phi[b];
1545             inf->next = mdb->dyngrid[b];
1546             mdb->dyngrid[b] = inf;
1547           }
1548         }
1549       }
1550     }
1551     else {
1552       modifier_setError(&mmd->modifier, "Failed to find bind solution (increase precision?)");
1553       error("Mesh Deform: failed to find bind solution.");
1554       break;
1555     }
1556
1557     BLI_snprintf(
1558         message, sizeof(message), "Mesh deform solve %d / %d       |||", a + 1, mdb->totcagevert);
1559     progress_bar((float)(a + 1) / (float)(mdb->totcagevert), message);
1560   }
1561
1562 #if 0
1563   /* sanity check */
1564   for (b = 0; b < mdb->size3; b++) {
1565     if (mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR) {
1566       if (fabsf(mdb->totalphi[b] - 1.0f) > 1e-4f) {
1567         printf("totalphi deficiency [%s|%d] %d: %.10f\n",
1568                (mdb->tag[b] == MESHDEFORM_TAG_INTERIOR) ? "interior" : "boundary",
1569                mdb->semibound[b],
1570                mdb->varidx[b],
1571                mdb->totalphi[b]);
1572       }
1573     }
1574   }
1575 #endif
1576
1577   /* free */
1578   MEM_freeN(mdb->varidx);
1579
1580   EIG_linear_solver_delete(context);
1581 }
1582
1583 static void harmonic_coordinates_bind(MeshDeformModifierData *mmd, MeshDeformBind *mdb)
1584 {
1585   MDefBindInfluence *inf;
1586   MDefInfluence *mdinf;
1587   MDefCell *cell;
1588   float center[3], vec[3], maxwidth, totweight;
1589   int a, b, x, y, z, totinside, offset;
1590
1591   /* compute bounding box of the cage mesh */
1592   INIT_MINMAX(mdb->min, mdb->max);
1593
1594   for (a = 0; a < mdb->totcagevert; a++) {
1595     minmax_v3v3_v3(mdb->min, mdb->max, mdb->cagecos[a]);
1596   }
1597
1598   /* allocate memory */
1599   mdb->size = (2 << (mmd->gridsize - 1)) + 2;
1600   mdb->size3 = mdb->size * mdb->size * mdb->size;
1601   mdb->tag = MEM_callocN(sizeof(int) * mdb->size3, "MeshDeformBindTag");
1602   mdb->phi = MEM_callocN(sizeof(float) * mdb->size3, "MeshDeformBindPhi");
1603   mdb->totalphi = MEM_callocN(sizeof(float) * mdb->size3, "MeshDeformBindTotalPhi");
1604   mdb->boundisect = MEM_callocN(sizeof(*mdb->boundisect) * mdb->size3, "MDefBoundIsect");
1605   mdb->semibound = MEM_callocN(sizeof(int) * mdb->size3, "MDefSemiBound");
1606   mdb->bvhtree = BKE_bvhtree_from_mesh_get(&mdb->bvhdata, mdb->cagemesh, BVHTREE_FROM_LOOPTRI, 4);
1607   mdb->inside = MEM_callocN(sizeof(int) * mdb->totvert, "MDefInside");
1608
1609   if (mmd->flag & MOD_MDEF_DYNAMIC_BIND) {
1610     mdb->dyngrid = MEM_callocN(sizeof(MDefBindInfluence *) * mdb->size3, "MDefDynGrid");
1611   }
1612   else {
1613     mdb->weights = MEM_callocN(sizeof(float) * mdb->totvert * mdb->totcagevert, "MDefWeights");
1614   }
1615
1616   mdb->memarena = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "harmonic coords arena");
1617   BLI_memarena_use_calloc(mdb->memarena);
1618
1619   /* initialize data from 'cagedm' for reuse */
1620   {
1621     Mesh *me = mdb->cagemesh;
1622     mdb->cagemesh_cache.mpoly = me->mpoly;
1623     mdb->cagemesh_cache.mloop = me->mloop;
1624     mdb->cagemesh_cache.looptri = BKE_mesh_runtime_looptri_ensure(me);
1625     /* can be NULL */
1626     mdb->cagemesh_cache.poly_nors = CustomData_get_layer(&me->pdata, CD_NORMAL);
1627   }
1628
1629   /* make bounding box equal size in all directions, add padding, and compute
1630    * width of the cells */
1631   maxwidth = -1.0f;
1632   for (a = 0; a < 3; a++) {
1633     if (mdb->max[a] - mdb->min[a] > maxwidth) {
1634       maxwidth = mdb->max[a] - mdb->min[a];
1635     }
1636   }
1637
1638   for (a = 0; a < 3; a++) {
1639     center[a] = (mdb->min[a] + mdb->max[a]) * 0.5f;
1640     mdb->min[a] = center[a] - maxwidth * 0.5f;
1641     mdb->max[a] = center[a] + maxwidth * 0.5f;
1642
1643     mdb->width[a] = (mdb->max[a] - mdb->min[a]) / (mdb->size - 4);
1644     mdb->min[a] -= 2.1f * mdb->width[a];
1645     mdb->max[a] += 2.1f * mdb->width[a];
1646
1647     mdb->width[a] = (mdb->max[a] - mdb->min[a]) / mdb->size;
1648     mdb->halfwidth[a] = mdb->width[a] * 0.5f;
1649   }
1650
1651   progress_bar(0, "Setting up mesh deform system");
1652
1653   totinside = 0;
1654   for (a = 0; a < mdb->totvert; a++) {
1655     copy_v3_v3(vec, mdb->vertexcos[a]);
1656     mdb->inside[a] = meshdeform_inside_cage(mdb, vec);
1657     if (mdb->inside[a]) {
1658       totinside++;
1659     }
1660   }
1661
1662   /* free temporary MDefBoundIsects */
1663   BLI_memarena_free(mdb->memarena);
1664   mdb->memarena = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "harmonic coords arena");
1665
1666   /* start with all cells untyped */
1667   for (a = 0; a < mdb->size3; a++) {
1668     mdb->tag[a] = MESHDEFORM_TAG_UNTYPED;
1669   }
1670
1671   /* detect intersections and tag boundary cells */
1672   for (z = 0; z < mdb->size; z++) {
1673     for (y = 0; y < mdb->size; y++) {
1674       for (x = 0; x < mdb->size; x++) {
1675         meshdeform_add_intersections(mdb, x, y, z);
1676       }
1677     }
1678   }
1679
1680   /* compute exterior and interior tags */
1681   meshdeform_bind_floodfill(mdb);
1682
1683   for (z = 0; z < mdb->size; z++) {
1684     for (y = 0; y < mdb->size; y++) {
1685       for (x = 0; x < mdb->size; x++) {
1686         meshdeform_check_semibound(mdb, x, y, z);
1687       }
1688     }
1689   }
1690
1691   /* solve */
1692   meshdeform_matrix_solve(mmd, mdb);
1693
1694   /* assign results */
1695   if (mmd->flag & MOD_MDEF_DYNAMIC_BIND) {
1696     mmd->totinfluence = 0;
1697     for (a = 0; a < mdb->size3; a++) {
1698       for (inf = mdb->dyngrid[a]; inf; inf = inf->next) {
1699         mmd->totinfluence++;
1700       }
1701     }
1702
1703     /* convert MDefBindInfluences to smaller MDefInfluences */
1704     mmd->dyngrid = MEM_callocN(sizeof(MDefCell) * mdb->size3, "MDefDynGrid");
1705     mmd->dyninfluences = MEM_callocN(sizeof(MDefInfluence) * mmd->totinfluence, "MDefInfluence");
1706     offset = 0;
1707     for (a = 0; a < mdb->size3; a++) {
1708       cell = &mmd->dyngrid[a];
1709       cell->offset = offset;
1710
1711       totweight = 0.0f;
1712       mdinf = mmd->dyninfluences + cell->offset;
1713       for (inf = mdb->dyngrid[a]; inf; inf = inf->next, mdinf++) {
1714         mdinf->weight = inf->weight;
1715         mdinf->vertex = inf->vertex;
1716         totweight += mdinf->weight;
1717         cell->totinfluence++;
1718       }
1719
1720       if (totweight > 0.0f) {
1721         mdinf = mmd->dyninfluences + cell->offset;
1722         for (b = 0; b < cell->totinfluence; b++, mdinf++) {
1723           mdinf->weight /= totweight;
1724         }
1725       }
1726
1727       offset += cell->totinfluence;
1728     }
1729
1730     mmd->dynverts = mdb->inside;
1731     mmd->dyngridsize = mdb->size;
1732     copy_v3_v3(mmd->dyncellmin, mdb->min);
1733     mmd->dyncellwidth = mdb->width[0];
1734     MEM_freeN(mdb->dyngrid);
1735   }
1736   else {
1737     mmd->bindweights = mdb->weights;
1738     MEM_freeN(mdb->inside);
1739   }
1740
1741   MEM_freeN(mdb->tag);
1742   MEM_freeN(mdb->phi);
1743   MEM_freeN(mdb->totalphi);
1744   MEM_freeN(mdb->boundisect);
1745   MEM_freeN(mdb->semibound);
1746   BLI_memarena_free(mdb->memarena);
1747   free_bvhtree_from_mesh(&mdb->bvhdata);
1748 }
1749
1750 void ED_mesh_deform_bind_callback(MeshDeformModifierData *mmd,
1751                                   Mesh *cagemesh,
1752                                   float *vertexcos,
1753                                   int totvert,
1754                                   float cagemat[4][4])
1755 {
1756   MeshDeformModifierData *mmd_orig = (MeshDeformModifierData *)modifier_get_original(
1757       &mmd->modifier);
1758   MeshDeformBind mdb;
1759   MVert *mvert;
1760   int a;
1761
1762   waitcursor(1);
1763   start_progress_bar();
1764
1765   memset(&mdb, 0, sizeof(MeshDeformBind));
1766
1767   /* get mesh and cage mesh */
1768   mdb.vertexcos = MEM_callocN(sizeof(float) * 3 * totvert, "MeshDeformCos");
1769   mdb.totvert = totvert;
1770
1771   mdb.cagemesh = cagemesh;
1772   mdb.totcagevert = mdb.cagemesh->totvert;
1773   mdb.cagecos = MEM_callocN(sizeof(*mdb.cagecos) * mdb.totcagevert, "MeshDeformBindCos");
1774   copy_m4_m4(mdb.cagemat, cagemat);
1775
1776   mvert = mdb.cagemesh->mvert;
1777   for (a = 0; a < mdb.totcagevert; a++) {
1778     copy_v3_v3(mdb.cagecos[a], mvert[a].co);
1779   }
1780   for (a = 0; a < mdb.totvert; a++) {
1781     mul_v3_m4v3(mdb.vertexcos[a], mdb.cagemat, vertexcos + a * 3);
1782   }
1783
1784   /* solve */
1785   harmonic_coordinates_bind(mmd_orig, &mdb);
1786
1787   /* assign bind variables */
1788   mmd_orig->bindcagecos = (float *)mdb.cagecos;
1789   mmd_orig->totvert = mdb.totvert;
1790   mmd_orig->totcagevert = mdb.totcagevert;
1791   copy_m4_m4(mmd_orig->bindmat, mmd_orig->object->obmat);
1792
1793   /* transform bindcagecos to world space */
1794   for (a = 0; a < mdb.totcagevert; a++) {
1795     mul_m4_v3(mmd_orig->object->obmat, mmd_orig->bindcagecos + a * 3);
1796   }
1797
1798   /* free */
1799   MEM_freeN(mdb.vertexcos);
1800
1801   /* compact weights */
1802   modifier_mdef_compact_influences((ModifierData *)mmd_orig);
1803
1804   end_progress_bar();
1805   waitcursor(0);
1806 }