source/blender/editors/armature/meshlaplacian.c

   1 /**
   2  * $Id$
   3  *
   4  * ***** BEGIN GPL LICENSE BLOCK *****
   5  *
   6  * This program is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU General Public License
   8  * as published by the Free Software Foundation; either version 2
   9  * of the License, or (at your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software Foundation,
  18  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  *
  20  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
  21  * All rights reserved.
  22  *
  23  * The Original Code is: all of this file.
  24  *
  25  * Contributor(s): none yet.
  26  *
  27  * ***** END GPL LICENSE BLOCK *****
  28  * meshlaplacian.c: Algorithms using the mesh laplacian.
  29  */
  30
  31 #include <math.h>
  32 #include <string.h>
  33
  34 #include "MEM_guardedalloc.h"
  35
  36 #include "DNA_object_types.h"
  37 #include "DNA_mesh_types.h"
  38 #include "DNA_meshdata_types.h"
  39 #include "DNA_modifier_types.h"
  40 #include "DNA_scene_types.h"
  41
  42 #include "BLI_math.h"
  43 #include "BLI_edgehash.h"
  44 #include "BLI_memarena.h"
  45
  46 #include "BKE_DerivedMesh.h"
  47 #include "BKE_modifier.h"
  48 #include "BKE_utildefines.h"
  49
  50 #ifdef RIGID_DEFORM
  51 #include "BLI_editVert.h"
  52 #include "BLI_polardecomp.h"
  53 #endif
  54
  55 #include "RE_raytrace.h"
  56
  57 #include "ONL_opennl.h"
  58
  59 #include "BLO_sys_types.h" // for intptr_t support
  60
  61 #include "ED_mesh.h"
  62
  63 #include "meshlaplacian.h"
  64
  65
  66 /* ************* XXX *************** */
  67 static void waitcursor(int val) {}
  68 static void progress_bar() {}
  69 static void start_progress_bar() {}
  70 static void end_progress_bar() {}
  71 static void error(char *str) { printf("error: %s\n", str); }
  72 /* ************* XXX *************** */
  73
  74
  75 /************************** Laplacian System *****************************/
  76
  77 struct LaplacianSystem {
  78         NLContext context;      /* opennl context */
  79
  80         int totvert, totface;
  81
  82         float **verts;                  /* vertex coordinates */
  83         float *varea;                   /* vertex weights for laplacian computation */
  84         char *vpinned;                  /* vertex pinning */
  85         int (*faces)[3];                /* face vertex indices */
  86         float (*fweights)[3];   /* cotangent weights per face */
  87
  88         int areaweights;                /* use area in cotangent weights? */
  89         int storeweights;               /* store cotangent weights in fweights */
  90         int nlbegun;                    /* nlBegin(NL_SYSTEM/NL_MATRIX) done */
  91
  92         EdgeHash *edgehash;             /* edge hash for construction */
  93
  94         struct HeatWeighting {
  95                 MFace *mface;
  96                 int totvert;
  97                 int totface;
  98                 float (*verts)[3];      /* vertex coordinates */
  99                 float (*vnors)[3];      /* vertex normals */
 100
 101                 float (*root)[3];       /* bone root */
 102                 float (*tip)[3];        /* bone tip */
 103                 float (*source)[3]; /* vertex source */
 104                 int numsource;
 105
 106                 float *H;                       /* diagonal H matrix */
 107                 float *p;                       /* values from all p vectors */
 108                 float *mindist;         /* minimum distance to a bone for all vertices */
 109
 110                 RayObject *raytree;     /* ray tracing acceleration structure */
 111                 RayFace   *faces;       /* faces to add to the ray tracing struture */
 112                 MFace     **vface;      /* a face that the vertex belongs to */
 113         } heat;
 114
 115 #ifdef RIGID_DEFORM
 116         struct RigidDeformation {
 117                 EditMesh *mesh;
 118
 119                 float (*R)[3][3];
 120                 float (*rhs)[3];
 121                 float (*origco)[3];
 122                 int thrownerror;
 123         } rigid;
 124 #endif
 125 };
 126
 127 /* Laplacian matrix construction */
 128
 129 /* Computation of these weights for the laplacian is based on:
 130    "Discrete Differential-Geometry Operators for Triangulated 2-Manifolds",
 131    Meyer et al, 2002. Section 3.5, formula (8).
 132
 133    We do it a bit different by going over faces instead of going over each
 134    vertex and adjacent faces, since we don't store this adjacency. Also, the
 135    formulas are tweaked a bit to work for non-manifold meshes. */
 136
 137 static void laplacian_increase_edge_count(EdgeHash *edgehash, int v1, int v2)
 138 {
 139         void **p = BLI_edgehash_lookup_p(edgehash, v1, v2);
 140
 141         if(p)
 142                 *p = (void*)((intptr_t)*p + (intptr_t)1);
 143         else
 144                 BLI_edgehash_insert(edgehash, v1, v2, (void*)(intptr_t)1);
 145 }
 146
 147 static int laplacian_edge_count(EdgeHash *edgehash, int v1, int v2)
 148 {
 149         return (int)(intptr_t)BLI_edgehash_lookup(edgehash, v1, v2);
 150 }
 151
 152 static float cotan_weight(float *v1, float *v2, float *v3)
 153 {
 154         float a[3], b[3], c[3], clen;
 155
 156         sub_v3_v3v3(a, v2, v1);
 157         sub_v3_v3v3(b, v3, v1);
 158         cross_v3_v3v3(c, a, b);
 159
 160         clen = len_v3(c);
 161
 162         if (clen == 0.0f)
 163                 return 0.0f;
 164
 165         return dot_v3v3(a, b)/clen;
 166 }
 167
 168 static void laplacian_triangle_area(LaplacianSystem *sys, int i1, int i2, int i3)
 169 {
 170         float t1, t2, t3, len1, len2, len3, area;
 171         float *varea= sys->varea, *v1, *v2, *v3;
 172         int obtuse = 0;
 173
 174         v1= sys->verts[i1];
 175         v2= sys->verts[i2];
 176         v3= sys->verts[i3];
 177
 178         t1= cotan_weight(v1, v2, v3);
 179         t2= cotan_weight(v2, v3, v1);
 180         t3= cotan_weight(v3, v1, v2);
 181
 182         if(RAD2DEG(angle_v3v3v3(v2, v1, v3)) > 90) obtuse= 1;
 183         else if(RAD2DEG(angle_v3v3v3(v1, v2, v3)) > 90) obtuse= 2;
 184         else if(RAD2DEG(angle_v3v3v3(v1, v3, v2)) > 90) obtuse= 3;
 185
 186         if (obtuse > 0) {
 187                 area= area_tri_v3(v1, v2, v3);
 188
 189                 varea[i1] += (obtuse == 1)? area: area*0.5;
 190                 varea[i2] += (obtuse == 2)? area: area*0.5;
 191                 varea[i3] += (obtuse == 3)? area: area*0.5;
 192         }
 193         else {
 194                 len1= len_v3v3(v2, v3);
 195                 len2= len_v3v3(v1, v3);
 196                 len3= len_v3v3(v1, v2);
 197
 198                 t1 *= len1*len1;
 199                 t2 *= len2*len2;
 200                 t3 *= len3*len3;
 201
 202                 varea[i1] += (t2 + t3)*0.25f;
 203                 varea[i2] += (t1 + t3)*0.25f;
 204                 varea[i3] += (t1 + t2)*0.25f;
 205         }
 206 }
 207
 208 static void laplacian_triangle_weights(LaplacianSystem *sys, int f, int i1, int i2, int i3)
 209 {
 210         float t1, t2, t3;
 211         float *varea= sys->varea, *v1, *v2, *v3;
 212
 213         v1= sys->verts[i1];
 214         v2= sys->verts[i2];
 215         v3= sys->verts[i3];
 216
 217         /* instead of *0.5 we divided by the number of faces of the edge, it still
 218            needs to be verified that this is indeed the correct thing to do! */
 219         t1= cotan_weight(v1, v2, v3)/laplacian_edge_count(sys->edgehash, i2, i3);
 220         t2= cotan_weight(v2, v3, v1)/laplacian_edge_count(sys->edgehash, i3, i1);
 221         t3= cotan_weight(v3, v1, v2)/laplacian_edge_count(sys->edgehash, i1, i2);
 222
 223         nlMatrixAdd(i1, i1, (t2+t3)*varea[i1]);
 224         nlMatrixAdd(i2, i2, (t1+t3)*varea[i2]);
 225         nlMatrixAdd(i3, i3, (t1+t2)*varea[i3]);
 226
 227         nlMatrixAdd(i1, i2, -t3*varea[i1]);
 228         nlMatrixAdd(i2, i1, -t3*varea[i2]);
 229
 230         nlMatrixAdd(i2, i3, -t1*varea[i2]);
 231         nlMatrixAdd(i3, i2, -t1*varea[i3]);
 232
 233         nlMatrixAdd(i3, i1, -t2*varea[i3]);
 234         nlMatrixAdd(i1, i3, -t2*varea[i1]);
 235
 236         if(sys->storeweights) {
 237                 sys->fweights[f][0]= t1*varea[i1];
 238                 sys->fweights[f][1]= t2*varea[i2];
 239                 sys->fweights[f][2]= t3*varea[i3];
 240         }
 241 }
 242
 243 LaplacianSystem *laplacian_system_construct_begin(int totvert, int totface, int lsq)
 244 {
 245         LaplacianSystem *sys;
 246
 247         sys= MEM_callocN(sizeof(LaplacianSystem), "LaplacianSystem");
 248
 249         sys->verts= MEM_callocN(sizeof(float*)*totvert, "LaplacianSystemVerts");
 250         sys->vpinned= MEM_callocN(sizeof(char)*totvert, "LaplacianSystemVpinned");
 251         sys->faces= MEM_callocN(sizeof(int)*3*totface, "LaplacianSystemFaces");
 252
 253         sys->totvert= 0;
 254         sys->totface= 0;
 255
 256         sys->areaweights= 1;
 257         sys->storeweights= 0;
 258
 259         /* create opennl context */
 260         nlNewContext();
 261         nlSolverParameteri(NL_NB_VARIABLES, totvert);
 262         if(lsq)
 263                 nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
 264
 265         sys->context= nlGetCurrent();
 266
 267         return sys;
 268 }
 269
 270 void laplacian_add_vertex(LaplacianSystem *sys, float *co, int pinned)
 271 {
 272         sys->verts[sys->totvert]= co;
 273         sys->vpinned[sys->totvert]= pinned;
 274         sys->totvert++;
 275 }
 276
 277 void laplacian_add_triangle(LaplacianSystem *sys, int v1, int v2, int v3)
 278 {
 279         sys->faces[sys->totface][0]= v1;
 280         sys->faces[sys->totface][1]= v2;
 281         sys->faces[sys->totface][2]= v3;
 282         sys->totface++;
 283 }
 284
 285 void laplacian_system_construct_end(LaplacianSystem *sys)
 286 {
 287         int (*face)[3];
 288         int a, totvert=sys->totvert, totface=sys->totface;
 289
 290         laplacian_begin_solve(sys, 0);
 291
 292         sys->varea= MEM_callocN(sizeof(float)*totvert, "LaplacianSystemVarea");
 293
 294         sys->edgehash= BLI_edgehash_new();
 295         for(a=0, face=sys->faces; a<sys->totface; a++, face++) {
 296                 laplacian_increase_edge_count(sys->edgehash, (*face)[0], (*face)[1]);
 297                 laplacian_increase_edge_count(sys->edgehash, (*face)[1], (*face)[2]);
 298                 laplacian_increase_edge_count(sys->edgehash, (*face)[2], (*face)[0]);
 299         }
 300
 301         if(sys->areaweights)
 302                 for(a=0, face=sys->faces; a<sys->totface; a++, face++)
 303                         laplacian_triangle_area(sys, (*face)[0], (*face)[1], (*face)[2]);
 304
 305         for(a=0; a<totvert; a++) {
 306                 if(sys->areaweights) {
 307                         if(sys->varea[a] != 0.0f)
 308                                 sys->varea[a]= 0.5f/sys->varea[a];
 309                 }
 310                 else
 311                         sys->varea[a]= 1.0f;
 312
 313                 /* for heat weighting */
 314                 if(sys->heat.H)
 315                         nlMatrixAdd(a, a, sys->heat.H[a]);
 316         }
 317
 318         if(sys->storeweights)
 319                 sys->fweights= MEM_callocN(sizeof(float)*3*totface, "LaplacianFWeight");
 320
 321         for(a=0, face=sys->faces; a<totface; a++, face++)
 322                 laplacian_triangle_weights(sys, a, (*face)[0], (*face)[1], (*face)[2]);
 323
 324         MEM_freeN(sys->faces);
 325         sys->faces= NULL;
 326
 327         if(sys->varea) {
 328                 MEM_freeN(sys->varea);
 329                 sys->varea= NULL;
 330         }
 331
 332         BLI_edgehash_free(sys->edgehash, NULL);
 333         sys->edgehash= NULL;
 334 }
 335
 336 void laplacian_system_delete(LaplacianSystem *sys)
 337 {
 338         if(sys->verts) MEM_freeN(sys->verts);
 339         if(sys->varea) MEM_freeN(sys->varea);
 340         if(sys->vpinned) MEM_freeN(sys->vpinned);
 341         if(sys->faces) MEM_freeN(sys->faces);
 342         if(sys->fweights) MEM_freeN(sys->fweights);
 343
 344         nlDeleteContext(sys->context);
 345         MEM_freeN(sys);
 346 }
 347
 348 void laplacian_begin_solve(LaplacianSystem *sys, int index)
 349 {
 350         int a;
 351
 352         if (!sys->nlbegun) {
 353                 nlBegin(NL_SYSTEM);
 354
 355                 if(index >= 0) {
 356                         for(a=0; a<sys->totvert; a++) {
 357                                 if(sys->vpinned[a]) {
 358                                         nlSetVariable(0, a, sys->verts[a][index]);
 359                                         nlLockVariable(a);
 360                                 }
 361                         }
 362                 }
 363
 364                 nlBegin(NL_MATRIX);
 365                 sys->nlbegun = 1;
 366         }
 367 }
 368
 369 void laplacian_add_right_hand_side(LaplacianSystem *sys, int v, float value)
 370 {
 371         nlRightHandSideAdd(0, v, value);
 372 }
 373
 374 int laplacian_system_solve(LaplacianSystem *sys)
 375 {
 376         nlEnd(NL_MATRIX);
 377         nlEnd(NL_SYSTEM);
 378         sys->nlbegun = 0;
 379
 380         //nlPrintMatrix();
 381
 382         return nlSolveAdvanced(NULL, NL_TRUE);
 383 }
 384
 385 float laplacian_system_get_solution(int v)
 386 {
 387         return nlGetVariable(0, v);
 388 }
 389
 390 /************************* Heat Bone Weighting ******************************/
 391 /* From "Automatic Rigging and Animation of 3D Characters"
 392                  Ilya Baran and Jovan Popovic, SIGGRAPH 2007 */
 393
 394 #define C_WEIGHT                        1.0f
 395 #define WEIGHT_LIMIT_START      0.05f
 396 #define WEIGHT_LIMIT_END        0.025f
 397 #define DISTANCE_EPSILON        1e-4f
 398
 399 /* Raytracing for vertex to bone/vertex visibility */
 400 static void heat_ray_tree_create(LaplacianSystem *sys)
 401 {
 402         MFace *mface = sys->heat.mface;
 403         int totface = sys->heat.totface;
 404         int totvert = sys->heat.totvert;
 405         int a;
 406
 407         sys->heat.raytree = RE_rayobject_vbvh_create(totface);
 408         sys->heat.faces = MEM_callocN(sizeof(RayFace)*totface, "Heat RayFaces");
 409         sys->heat.vface = MEM_callocN(sizeof(MFace*)*totvert, "HeatVFaces");
 410
 411         for(a=0; a<totface; a++) {
 412
 413                 MFace *mf = mface+a;
 414                 RayFace *rayface = sys->heat.faces+a;
 415
 416                 RayObject *obj = RE_rayface_from_coords(
 417                                                         rayface, &sys->heat, mf,
 418                                                         sys->heat.verts[mf->v1], sys->heat.verts[mf->v2],
 419                                                         sys->heat.verts[mf->v3], mf->v4 ? sys->heat.verts[mf->v4] : 0
 420                                                 );
 421                 RE_rayobject_add(sys->heat.raytree, obj);
 422
 423                 //Setup inverse pointers to use on isect.orig
 424                 sys->heat.vface[mf->v1]= mf;
 425                 sys->heat.vface[mf->v2]= mf;
 426                 sys->heat.vface[mf->v3]= mf;
 427                 if(mf->v4) sys->heat.vface[mf->v4]= mf;
 428         }
 429         RE_rayobject_done(sys->heat.raytree);
 430 }
 431
 432 static int heat_ray_source_visible(LaplacianSystem *sys, int vertex, int source)
 433 {
 434         Isect isec;
 435         MFace *mface;
 436         float end[3];
 437         int visible;
 438
 439         mface= sys->heat.vface[vertex];
 440         if(!mface)
 441                 return 1;
 442
 443         /* setup isec */
 444         memset(&isec, 0, sizeof(isec));
 445         isec.mode= RE_RAY_SHADOW;
 446         isec.lay= -1;
 447         isec.orig.ob = &sys->heat;
 448         isec.orig.face = mface;
 449         isec.skip = RE_SKIP_CULLFACE;
 450
 451         copy_v3_v3(isec.start, sys->heat.verts[vertex]);
 452
 453         if(sys->heat.root) /* bone */
 454                 closest_to_line_segment_v3(end, isec.start,
 455                         sys->heat.root[source], sys->heat.tip[source]);
 456         else /* vertex */
 457                 copy_v3_v3(end, sys->heat.source[source]);
 458
 459         sub_v3_v3v3(isec.vec, end, isec.start);
 460         isec.labda = 1.0f - 1e-5;
 461         madd_v3_v3v3fl(isec.start, isec.start, isec.vec, 1e-5);
 462
 463         visible= !RE_rayobject_raycast(sys->heat.raytree, &isec);
 464
 465         return visible;
 466 }
 467
 468 static float heat_source_distance(LaplacianSystem *sys, int vertex, int source)
 469 {
 470         float closest[3], d[3], dist, cosine;
 471
 472         /* compute euclidian distance */
 473         if(sys->heat.root) /* bone */
 474                 closest_to_line_segment_v3(closest, sys->heat.verts[vertex],
 475                         sys->heat.root[source], sys->heat.tip[source]);
 476         else /* vertex */
 477                 copy_v3_v3(closest, sys->heat.source[source]);
 478
 479         sub_v3_v3v3(d, sys->heat.verts[vertex], closest);
 480         dist= normalize_v3(d);
 481
 482         /* if the vertex normal does not point along the bone, increase distance */
 483         cosine= INPR(d, sys->heat.vnors[vertex]);
 484
 485         return dist/(0.5f*(cosine + 1.001f));
 486 }
 487
 488 static int heat_source_closest(LaplacianSystem *sys, int vertex, int source)
 489 {
 490         float dist;
 491
 492         dist= heat_source_distance(sys, vertex, source);
 493
 494         if(dist <= sys->heat.mindist[vertex]*(1.0f + DISTANCE_EPSILON))
 495                 if(heat_ray_source_visible(sys, vertex, source))
 496                         return 1;
 497
 498         return 0;
 499 }
 500
 501 static void heat_set_H(LaplacianSystem *sys, int vertex)
 502 {
 503         float dist, mindist, h;
 504         int j, numclosest = 0;
 505
 506         mindist= 1e10;
 507
 508         /* compute minimum distance */
 509         for(j=0; j<sys->heat.numsource; j++) {
 510                 dist= heat_source_distance(sys, vertex, j);
 511
 512                 if(dist < mindist)
 513                         mindist= dist;
 514         }
 515
 516         sys->heat.mindist[vertex]= mindist;
 517
 518         /* count number of sources with approximately this minimum distance */
 519         for(j=0; j<sys->heat.numsource; j++)
 520                 if(heat_source_closest(sys, vertex, j))
 521                         numclosest++;
 522
 523         sys->heat.p[vertex]= (numclosest > 0)? 1.0f/numclosest: 0.0f;
 524
 525         /* compute H entry */
 526         if(numclosest > 0) {
 527                 mindist= maxf(mindist, 1e-4f);
 528                 h= numclosest*C_WEIGHT/(mindist*mindist);
 529         }
 530         else
 531                 h= 0.0f;
 532
 533         sys->heat.H[vertex]= h;
 534 }
 535
 536 void heat_calc_vnormals(LaplacianSystem *sys)
 537 {
 538         float fnor[3];
 539         int a, v1, v2, v3, (*face)[3];
 540
 541         sys->heat.vnors= MEM_callocN(sizeof(float)*3*sys->totvert, "HeatVNors");
 542
 543         for(a=0, face=sys->faces; a<sys->totface; a++, face++) {
 544                 v1= (*face)[0];
 545                 v2= (*face)[1];
 546                 v3= (*face)[2];
 547
 548                 normal_tri_v3( fnor,sys->verts[v1], sys->verts[v2], sys->verts[v3]);
 549
 550                 add_v3_v3(sys->heat.vnors[v1], fnor);
 551                 add_v3_v3(sys->heat.vnors[v2], fnor);
 552                 add_v3_v3(sys->heat.vnors[v3], fnor);
 553         }
 554
 555         for(a=0; a<sys->totvert; a++)
 556                 normalize_v3(sys->heat.vnors[a]);
 557 }
 558
 559 static void heat_laplacian_create(LaplacianSystem *sys)
 560 {
 561         MFace *mface = sys->heat.mface, *mf;
 562         int totface= sys->heat.totface;
 563         int totvert= sys->heat.totvert;
 564         int a;
 565
 566         /* heat specific definitions */
 567         sys->heat.mindist= MEM_callocN(sizeof(float)*totvert, "HeatMinDist");
 568         sys->heat.H= MEM_callocN(sizeof(float)*totvert, "HeatH");
 569         sys->heat.p= MEM_callocN(sizeof(float)*totvert, "HeatP");
 570
 571         /* add verts and faces to laplacian */
 572         for(a=0; a<totvert; a++)
 573                 laplacian_add_vertex(sys, sys->heat.verts[a], 0);
 574
 575         for(a=0, mf=mface; a<totface; a++, mf++) {
 576                 laplacian_add_triangle(sys, mf->v1, mf->v2, mf->v3);
 577                 if(mf->v4)
 578                         laplacian_add_triangle(sys, mf->v1, mf->v3, mf->v4);
 579         }
 580
 581         /* for distance computation in set_H */
 582         heat_calc_vnormals(sys);
 583
 584         for(a=0; a<totvert; a++)
 585                 heat_set_H(sys, a);
 586 }
 587
 588 static void heat_system_free(LaplacianSystem *sys)
 589 {
 590         RE_rayobject_free(sys->heat.raytree);
 591         MEM_freeN(sys->heat.vface);
 592         MEM_freeN(sys->heat.faces);
 593
 594         MEM_freeN(sys->heat.mindist);
 595         MEM_freeN(sys->heat.H);
 596         MEM_freeN(sys->heat.p);
 597         MEM_freeN(sys->heat.vnors);
 598 }
 599
 600 static float heat_limit_weight(float weight)
 601 {
 602         float t;
 603
 604         if(weight < WEIGHT_LIMIT_END) {
 605                 return 0.0f;
 606         }
 607         else if(weight < WEIGHT_LIMIT_START) {
 608                 t= (weight - WEIGHT_LIMIT_END)/(WEIGHT_LIMIT_START - WEIGHT_LIMIT_END);
 609                 return t*WEIGHT_LIMIT_START;
 610         }
 611         else
 612                 return weight;
 613 }
 614
 615 void heat_bone_weighting(Object *ob, Mesh *me, float (*verts)[3], int numsource, bDeformGroup **dgrouplist, bDeformGroup **dgroupflip, float (*root)[3], float (*tip)[3], int *selected)
 616 {
 617         LaplacianSystem *sys;
 618         MFace *mface;
 619         float solution, weight;
 620         int *vertsflipped = NULL, *mask= NULL;
 621         int a, totface, j, bbone, firstsegment, lastsegment, thrownerror = 0;
 622
 623         /* count triangles and create mask */
 624         if(me->editflag & ME_EDIT_PAINT_MASK)
 625                 mask= MEM_callocN(sizeof(int)*me->totvert, "heat_bone_weighting mask");
 626
 627         for(totface=0, a=0, mface=me->mface; a<me->totface; a++, mface++) {
 628                 totface++;
 629                 if(mface->v4) totface++;
 630
 631                 if(mask && (mface->flag & ME_FACE_SEL)) {
 632                         mask[mface->v1]= 1;
 633                         mask[mface->v2]= 1;
 634                         mask[mface->v3]= 1;
 635                         if(mface->v4)
 636                                 mask[mface->v4]= 1;
 637                 }
 638         }
 639
 640         /* create laplacian */
 641         sys = laplacian_system_construct_begin(me->totvert, totface, 1);
 642
 643         sys->heat.mface= me->mface;
 644         sys->heat.totface= me->totface;
 645         sys->heat.totvert= me->totvert;
 646         sys->heat.verts= verts;
 647         sys->heat.root= root;
 648         sys->heat.tip= tip;
 649         sys->heat.numsource= numsource;
 650
 651         heat_ray_tree_create(sys);
 652         heat_laplacian_create(sys);
 653
 654         laplacian_system_construct_end(sys);
 655
 656         if(dgroupflip) {
 657                 vertsflipped = MEM_callocN(sizeof(int)*me->totvert, "vertsflipped");
 658                 for(a=0; a<me->totvert; a++)
 659                         vertsflipped[a] = mesh_get_x_mirror_vert(ob, a);
 660         }
 661
 662         /* compute weights per bone */
 663         for(j=0; j<numsource; j++) {
 664                 if(!selected[j])
 665                         continue;
 666
 667                 firstsegment= (j == 0 || dgrouplist[j-1] != dgrouplist[j]);
 668                 lastsegment= (j == numsource-1 || dgrouplist[j] != dgrouplist[j+1]);
 669                 bbone= !(firstsegment && lastsegment);
 670
 671                 /* clear weights */
 672                 if(bbone && firstsegment) {
 673                         for(a=0; a<me->totvert; a++) {
 674                                 if(mask && !mask[a])
 675                                         continue;
 676
 677                                 ED_vgroup_vert_remove(ob, dgrouplist[j], a);
 678                                 if(vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0)
 679                                         ED_vgroup_vert_remove(ob, dgroupflip[j], vertsflipped[a]);
 680                         }
 681                 }
 682
 683                 /* fill right hand side */
 684                 laplacian_begin_solve(sys, -1);
 685
 686                 for(a=0; a<me->totvert; a++)
 687                         if(heat_source_closest(sys, a, j))
 688                                 laplacian_add_right_hand_side(sys, a,
 689                                         sys->heat.H[a]*sys->heat.p[a]);
 690
 691                 /* solve */
 692                 if(laplacian_system_solve(sys)) {
 693                         /* load solution into vertex groups */
 694                         for(a=0; a<me->totvert; a++) {
 695                                 if(mask && !mask[a])
 696                                         continue;
 697
 698                                 solution= laplacian_system_get_solution(a);
 699
 700                                 if(bbone) {
 701                                         if(solution > 0.0f)
 702                                                 ED_vgroup_vert_add(ob, dgrouplist[j], a, solution,
 703                                                         WEIGHT_ADD);
 704                                 }
 705                                 else {
 706                                         weight= heat_limit_weight(solution);
 707                                         if(weight > 0.0f)
 708                                                 ED_vgroup_vert_add(ob, dgrouplist[j], a, weight,
 709                                                         WEIGHT_REPLACE);
 710                                         else
 711                                                 ED_vgroup_vert_remove(ob, dgrouplist[j], a);
 712                                 }
 713
 714                                 /* do same for mirror */
 715                                 if(vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0) {
 716                                         if(bbone) {
 717                                                 if(solution > 0.0f)
 718                                                         ED_vgroup_vert_add(ob, dgroupflip[j], vertsflipped[a],
 719                                                                 solution, WEIGHT_ADD);
 720                                         }
 721                                         else {
 722                                                 weight= heat_limit_weight(solution);
 723                                                 if(weight > 0.0f)
 724                                                         ED_vgroup_vert_add(ob, dgroupflip[j], vertsflipped[a],
 725                                                                 weight, WEIGHT_REPLACE);
 726                                                 else
 727                                                         ED_vgroup_vert_remove(ob, dgroupflip[j], vertsflipped[a]);
 728                                         }
 729                                 }
 730                         }
 731                 }
 732                 else if(!thrownerror) {
 733                         error("Bone Heat Weighting:"
 734                                 " failed to find solution for one or more bones");
 735                         thrownerror= 1;
 736                         break;
 737                 }
 738
 739                 /* remove too small vertex weights */
 740                 if(bbone && lastsegment) {
 741                         for(a=0; a<me->totvert; a++) {
 742                                 if(mask && !mask[a])
 743                                         continue;
 744
 745                                 weight= ED_vgroup_vert_weight(ob, dgrouplist[j], a);
 746                                 weight= heat_limit_weight(weight);
 747                                 if(weight <= 0.0f)
 748                                         ED_vgroup_vert_remove(ob, dgrouplist[j], a);
 749
 750                                 if(vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0) {
 751                                         weight= ED_vgroup_vert_weight(ob, dgroupflip[j], vertsflipped[a]);
 752                                         weight= heat_limit_weight(weight);
 753                                         if(weight <= 0.0f)
 754                                                 ED_vgroup_vert_remove(ob, dgroupflip[j], vertsflipped[a]);
 755                                 }
 756                         }
 757                 }
 758         }
 759
 760         /* free */
 761         if(vertsflipped) MEM_freeN(vertsflipped);
 762         if(mask) MEM_freeN(mask);
 763
 764         heat_system_free(sys);
 765
 766         laplacian_system_delete(sys);
 767 }
 768
 769 #ifdef RIGID_DEFORM
 770 /********************** As-Rigid-As-Possible Deformation ******************/
 771 /* From "As-Rigid-As-Possible Surface Modeling",
 772                 Olga Sorkine and Marc Alexa, ESGP 2007. */
 773
 774 /* investigate:
 775    - transpose R in orthogonal
 776    - flipped normals and per face adding
 777    - move cancelling to transform, make origco pointer
 778 */
 779
 780 static LaplacianSystem *RigidDeformSystem = NULL;
 781
 782 static void rigid_add_half_edge_to_R(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
 783 {
 784         float e[3], e_[3];
 785         int i;
 786
 787         sub_v3_v3v3(e, sys->rigid.origco[v1->tmp.l], sys->rigid.origco[v2->tmp.l]);
 788         sub_v3_v3v3(e_, v1->co, v2->co);
 789
 790         /* formula (5) */
 791         for (i=0; i<3; i++) {
 792                 sys->rigid.R[v1->tmp.l][i][0] += w*e[0]*e_[i];
 793                 sys->rigid.R[v1->tmp.l][i][1] += w*e[1]*e_[i];
 794                 sys->rigid.R[v1->tmp.l][i][2] += w*e[2]*e_[i];
 795         }
 796 }
 797
 798 static void rigid_add_edge_to_R(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
 799 {
 800         rigid_add_half_edge_to_R(sys, v1, v2, w);
 801         rigid_add_half_edge_to_R(sys, v2, v1, w);
 802 }
 803
 804 static void rigid_orthogonalize_R(float R[][3])
 805 {
 806         HMatrix M, Q, S;
 807
 808         copy_m4_m3(M, R);
 809         polar_decomp(M, Q, S);
 810         copy_m3_m4(R, Q);
 811 }
 812
 813 static void rigid_add_half_edge_to_rhs(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
 814 {
 815         /* formula (8) */
 816         float Rsum[3][3], rhs[3];
 817
 818         if (sys->vpinned[v1->tmp.l])
 819                 return;
 820
 821         add_m3_m3m3(Rsum, sys->rigid.R[v1->tmp.l], sys->rigid.R[v2->tmp.l]);
 822         transpose_m3(Rsum);
 823
 824         sub_v3_v3v3(rhs, sys->rigid.origco[v1->tmp.l], sys->rigid.origco[v2->tmp.l]);
 825         mul_m3_v3(Rsum, rhs);
 826         mul_v3_fl(rhs, 0.5f);
 827         mul_v3_fl(rhs, w);
 828
 829         add_v3_v3(sys->rigid.rhs[v1->tmp.l], rhs);
 830 }
 831
 832 static void rigid_add_edge_to_rhs(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
 833 {
 834         rigid_add_half_edge_to_rhs(sys, v1, v2, w);
 835         rigid_add_half_edge_to_rhs(sys, v2, v1, w);
 836 }
 837
 838 void rigid_deform_iteration()
 839 {
 840         LaplacianSystem *sys= RigidDeformSystem;
 841         EditMesh *em;
 842         EditVert *eve;
 843         EditFace *efa;
 844         int a, i;
 845
 846         if(!sys)
 847                 return;
 848
 849         nlMakeCurrent(sys->context);
 850         em= sys->rigid.mesh;
 851
 852         /* compute R */
 853         memset(sys->rigid.R, 0, sizeof(float)*3*3*sys->totvert);
 854         memset(sys->rigid.rhs, 0, sizeof(float)*3*sys->totvert);
 855
 856         for(a=0, efa=em->faces.first; efa; efa=efa->next, a++) {
 857                 rigid_add_edge_to_R(sys, efa->v1, efa->v2, sys->fweights[a][2]);
 858                 rigid_add_edge_to_R(sys, efa->v2, efa->v3, sys->fweights[a][0]);
 859                 rigid_add_edge_to_R(sys, efa->v3, efa->v1, sys->fweights[a][1]);
 860
 861                 if(efa->v4) {
 862                         a++;
 863                         rigid_add_edge_to_R(sys, efa->v1, efa->v3, sys->fweights[a][2]);
 864                         rigid_add_edge_to_R(sys, efa->v3, efa->v4, sys->fweights[a][0]);
 865                         rigid_add_edge_to_R(sys, efa->v4, efa->v1, sys->fweights[a][1]);
 866                 }
 867         }
 868
 869         for(a=0, eve=em->verts.first; eve; eve=eve->next, a++) {
 870                 rigid_orthogonalize_R(sys->rigid.R[a]);
 871                 eve->tmp.l= a;
 872         }
 873
 874         /* compute right hand sides for solving */
 875         for(a=0, efa=em->faces.first; efa; efa=efa->next, a++) {
 876                 rigid_add_edge_to_rhs(sys, efa->v1, efa->v2, sys->fweights[a][2]);
 877                 rigid_add_edge_to_rhs(sys, efa->v2, efa->v3, sys->fweights[a][0]);
 878                 rigid_add_edge_to_rhs(sys, efa->v3, efa->v1, sys->fweights[a][1]);
 879
 880                 if(efa->v4) {
 881                         a++;
 882                         rigid_add_edge_to_rhs(sys, efa->v1, efa->v3, sys->fweights[a][2]);
 883                         rigid_add_edge_to_rhs(sys, efa->v3, efa->v4, sys->fweights[a][0]);
 884                         rigid_add_edge_to_rhs(sys, efa->v4, efa->v1, sys->fweights[a][1]);
 885                 }
 886         }
 887
 888         /* solve for positions, for X,Y and Z separately */
 889         for(i=0; i<3; i++) {
 890                 laplacian_begin_solve(sys, i);
 891
 892                 for(a=0; a<sys->totvert; a++)
 893                         if(!sys->vpinned[a])
 894                                 laplacian_add_right_hand_side(sys, a, sys->rigid.rhs[a][i]);
 895
 896                 if(laplacian_system_solve(sys)) {
 897                         for(a=0, eve=em->verts.first; eve; eve=eve->next, a++)
 898                                 eve->co[i]= laplacian_system_get_solution(a);
 899                 }
 900                 else {
 901                         if(!sys->rigid.thrownerror) {
 902                                 error("RigidDeform: failed to find solution.");
 903                                 sys->rigid.thrownerror= 1;
 904                         }
 905                         break;
 906                 }
 907         }
 908 }
 909
 910 static void rigid_laplacian_create(LaplacianSystem *sys)
 911 {
 912         EditMesh *em = sys->rigid.mesh;
 913         EditVert *eve;
 914         EditFace *efa;
 915         int a;
 916
 917         /* add verts and faces to laplacian */
 918         for(a=0, eve=em->verts.first; eve; eve=eve->next, a++) {
 919                 laplacian_add_vertex(sys, eve->co, eve->pinned);
 920                 eve->tmp.l= a;
 921         }
 922
 923         for(efa=em->faces.first; efa; efa=efa->next) {
 924                 laplacian_add_triangle(sys,
 925                         efa->v1->tmp.l, efa->v2->tmp.l, efa->v3->tmp.l);
 926                 if(efa->v4)
 927                         laplacian_add_triangle(sys,
 928                                 efa->v1->tmp.l, efa->v3->tmp.l, efa->v4->tmp.l);
 929         }
 930 }
 931
 932 void rigid_deform_begin(EditMesh *em)
 933 {
 934         LaplacianSystem *sys;
 935         EditVert *eve;
 936         EditFace *efa;
 937         int a, totvert, totface;
 938
 939         /* count vertices, triangles */
 940         for(totvert=0, eve=em->verts.first; eve; eve=eve->next)
 941                 totvert++;
 942
 943         for(totface=0, efa=em->faces.first; efa; efa=efa->next) {
 944                 totface++;
 945                 if(efa->v4) totface++;
 946         }
 947
 948         /* create laplacian */
 949         sys = laplacian_system_construct_begin(totvert, totface, 0);
 950
 951         sys->rigid.mesh= em;
 952         sys->rigid.R = MEM_callocN(sizeof(float)*3*3*totvert, "RigidDeformR");
 953         sys->rigid.rhs = MEM_callocN(sizeof(float)*3*totvert, "RigidDeformRHS");
 954         sys->rigid.origco = MEM_callocN(sizeof(float)*3*totvert, "RigidDeformCo");
 955
 956         for(a=0, eve=em->verts.first; eve; eve=eve->next, a++)
 957                 copy_v3_v3(sys->rigid.origco[a], eve->co);
 958
 959         sys->areaweights= 0;
 960         sys->storeweights= 1;
 961
 962         rigid_laplacian_create(sys);
 963
 964         laplacian_system_construct_end(sys);
 965
 966         RigidDeformSystem = sys;
 967 }
 968
 969 void rigid_deform_end(int cancel)
 970 {
 971         LaplacianSystem *sys = RigidDeformSystem;
 972
 973         if(sys) {
 974                 EditMesh *em = sys->rigid.mesh;
 975                 EditVert *eve;
 976                 int a;
 977
 978                 if(cancel)
 979                         for(a=0, eve=em->verts.first; eve; eve=eve->next, a++)
 980                                 if(!eve->pinned)
 981                                         copy_v3_v3(eve->co, sys->rigid.origco[a]);
 982
 983                 if(sys->rigid.R) MEM_freeN(sys->rigid.R);
 984                 if(sys->rigid.rhs) MEM_freeN(sys->rigid.rhs);
 985                 if(sys->rigid.origco) MEM_freeN(sys->rigid.origco);
 986
 987                 /* free */
 988                 laplacian_system_delete(sys);
 989         }
 990
 991         RigidDeformSystem = NULL;
 992 }
 993 #endif
 994
 995 /************************** Harmonic Coordinates ****************************/
 996 /* From "Harmonic Coordinates for Character Articulation",
 997         Pushkar Joshi, Mark Meyer, Tony DeRose, Brian Green and Tom Sanocki,
 998         SIGGRAPH 2007. */
 999
1000 #define EPSILON 0.0001f
1001
1002 #define MESHDEFORM_TAG_UNTYPED  0
1003 #define MESHDEFORM_TAG_BOUNDARY 1
1004 #define MESHDEFORM_TAG_INTERIOR 2
1005 #define MESHDEFORM_TAG_EXTERIOR 3
1006
1007 #define MESHDEFORM_LEN_THRESHOLD 1e-6
1008
1009 #define MESHDEFORM_MIN_INFLUENCE 0.0005
1010
1011 static int MESHDEFORM_OFFSET[7][3] =
1012                 {{0,0,0}, {1,0,0}, {-1,0,0}, {0,1,0}, {0,-1,0}, {0,0,1}, {0,0,-1}};
1013
1014 typedef struct MDefBoundIsect {
1015         float co[3], uvw[4];
1016         int nvert, v[4], facing;
1017         float len;
1018 } MDefBoundIsect;
1019
1020 typedef struct MDefBindInfluence {
1021         struct MDefBindInfluence *next;
1022         float weight;
1023         int vertex;
1024 } MDefBindInfluence;
1025
1026 typedef struct MeshDeformBind {
1027         /* grid dimensions */
1028         float min[3], max[3];
1029         float width[3], halfwidth[3];
1030         int size, size3;
1031
1032         /* meshes */
1033         DerivedMesh *cagedm;
1034         float (*cagecos)[3];
1035         float (*vertexcos)[3];
1036         int totvert, totcagevert;
1037
1038         /* grids */
1039         MemArena *memarena;
1040         MDefBoundIsect *(*boundisect)[6];
1041         int *semibound;
1042         int *tag;
1043         float *phi, *totalphi;
1044
1045         /* mesh stuff */
1046         int *inside;
1047         float *weights;
1048         MDefBindInfluence **dyngrid;
1049         float cagemat[4][4];
1050
1051         /* direct solver */
1052         int *varidx;
1053
1054         /* raytrace */
1055         RayObject *raytree;
1056 } MeshDeformBind;
1057
1058 /* ray intersection */
1059
1060 /* our own triangle intersection, so we can fully control the epsilons and
1061  * prevent corner case from going wrong*/
1062 static int meshdeform_tri_intersect(float orig[3], float end[3], float vert0[3],
1063         float vert1[3], float vert2[3], float *isectco, float *uvw)
1064 {
1065         float edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
1066         float det,inv_det, u, v, dir[3], isectdir[3];
1067
1068         sub_v3_v3v3(dir, end, orig);
1069
1070         /* find vectors for two edges sharing vert0 */
1071         sub_v3_v3v3(edge1, vert1, vert0);
1072         sub_v3_v3v3(edge2, vert2, vert0);
1073
1074         /* begin calculating determinant - also used to calculate U parameter */
1075         cross_v3_v3v3(pvec, dir, edge2);
1076
1077         /* if determinant is near zero, ray lies in plane of triangle */
1078         det = INPR(edge1, pvec);
1079
1080         if (det == 0.0f)
1081           return 0;
1082         inv_det = 1.0f / det;
1083
1084         /* calculate distance from vert0 to ray origin */
1085         sub_v3_v3v3(tvec, orig, vert0);
1086
1087         /* calculate U parameter and test bounds */
1088         u = INPR(tvec, pvec) * inv_det;
1089         if (u < -EPSILON || u > 1.0f+EPSILON)
1090           return 0;
1091
1092         /* prepare to test V parameter */
1093         cross_v3_v3v3(qvec, tvec, edge1);
1094
1095         /* calculate V parameter and test bounds */
1096         v = INPR(dir, qvec) * inv_det;
1097         if (v < -EPSILON || u + v > 1.0f+EPSILON)
1098           return 0;
1099
1100         isectco[0]= (1.0f - u - v)*vert0[0] + u*vert1[0] + v*vert2[0];
1101         isectco[1]= (1.0f - u - v)*vert0[1] + u*vert1[1] + v*vert2[1];
1102         isectco[2]= (1.0f - u - v)*vert0[2] + u*vert1[2] + v*vert2[2];
1103
1104         uvw[0]= 1.0 - u - v;
1105         uvw[1]= u;
1106         uvw[2]= v;
1107
1108         /* check if it is within the length of the line segment */
1109         sub_v3_v3v3(isectdir, isectco, orig);
1110
1111         if(INPR(dir, isectdir) < -EPSILON)
1112                 return 0;
1113
1114         if(INPR(dir, dir) + EPSILON < INPR(isectdir, isectdir))
1115                 return 0;
1116
1117         return 1;
1118 }
1119
1120 /* blender's raytracer is not use now, even though it is much faster. it can
1121  * give problems with rays falling through, so we use our own intersection
1122  * function above with tweaked epsilons */
1123
1124 #if 0
1125 static MeshDeformBind *MESHDEFORM_BIND = NULL;
1126
1127 static void meshdeform_ray_coords_func(RayFace *face, float **v1, float **v2, float **v3, float **v4)
1128 {
1129         MFace *mface= (MFace*)face;
1130         float (*cagecos)[3]= MESHDEFORM_BIND->cagecos;
1131
1132         *v1= cagecos[mface->v1];
1133         *v2= cagecos[mface->v2];
1134         *v3= cagecos[mface->v3];
1135         *v4= (mface->v4)? cagecos[mface->v4]: NULL;
1136 }
1137
1138 static int meshdeform_ray_check_func(Isect *is, RayFace *face)
1139 {
1140         return 1;
1141 }
1142
1143 static void meshdeform_ray_tree_create(MeshDeformBind *mdb)
1144 {
1145         MFace *mface;
1146         float min[3], max[3];
1147         int a, totface;
1148
1149         /* create a raytrace tree from the mesh */
1150         INIT_MINMAX(min, max);
1151
1152         for(a=0; a<mdb->totcagevert; a++)
1153                 DO_MINMAX(mdb->cagecos[a], min, max)
1154
1155         MESHDEFORM_BIND= mdb;
1156
1157         mface= mdb->cagedm->getFaceArray(mdb->cagedm);
1158         totface= mdb->cagedm->getNumFaces(mdb->cagedm);
1159
1160         mdb->raytree= RE_ray_tree_create(64, totface, min, max,
1161                 meshdeform_ray_coords_func, meshdeform_ray_check_func);
1162
1163         for(a=0; a<totface; a++, mface++)
1164                 RE_ray_tree_add_face(mdb->raytree, mface);
1165
1166         RE_ray_tree_done(mdb->raytree);
1167 }
1168
1169 static void meshdeform_ray_tree_free(MeshDeformBind *mdb)
1170 {
1171         MESHDEFORM_BIND= NULL;
1172         RE_ray_tree_free(mdb->raytree);
1173 }
1174 #endif
1175
1176 static int meshdeform_intersect(MeshDeformBind *mdb, Isect *isec)
1177 {
1178         MFace *mface;
1179         float face[4][3], co[3], uvw[3], len, nor[3], end[3];
1180         int f, hit, is= 0, totface;
1181
1182         isec->labda= 1e10;
1183
1184         mface= mdb->cagedm->getFaceArray(mdb->cagedm);
1185         totface= mdb->cagedm->getNumFaces(mdb->cagedm);
1186
1187         add_v3_v3v3(end, isec->start, isec->vec);
1188
1189         for(f=0; f<totface; f++, mface++) {
1190                 copy_v3_v3(face[0], mdb->cagecos[mface->v1]);
1191                 copy_v3_v3(face[1], mdb->cagecos[mface->v2]);
1192                 copy_v3_v3(face[2], mdb->cagecos[mface->v3]);
1193
1194                 if(mface->v4) {
1195                         copy_v3_v3(face[3], mdb->cagecos[mface->v4]);
1196                         hit = meshdeform_tri_intersect(isec->start, end, face[0], face[1], face[2], co, uvw);
1197
1198                         if(hit) {
1199                                 normal_tri_v3( nor,face[0], face[1], face[2]);
1200                         }
1201                         else {
1202                                 hit= meshdeform_tri_intersect(isec->start, end, face[0], face[2], face[3], co, uvw);
1203                                 normal_tri_v3( nor,face[0], face[2], face[3]);
1204                         }
1205                 }
1206                 else {
1207                         hit= meshdeform_tri_intersect(isec->start, end, face[0], face[1], face[2], co, uvw);
1208                         normal_tri_v3( nor,face[0], face[1], face[2]);
1209                 }
1210
1211                 if(hit) {
1212                         len= len_v3v3(isec->start, co)/len_v3v3(isec->start, end);
1213                         if(len < isec->labda) {
1214                                 isec->labda= len;
1215                                 isec->hit.face = mface;
1216                                 isec->isect= (INPR(isec->vec, nor) <= 0.0f);
1217                                 is= 1;
1218                         }
1219                 }
1220         }
1221
1222         return is;
1223 }
1224
1225 static MDefBoundIsect *meshdeform_ray_tree_intersect(MeshDeformBind *mdb, float *co1, float *co2)
1226 {
1227         MDefBoundIsect *isect;
1228         Isect isec;
1229         float (*cagecos)[3];
1230         MFace *mface;
1231         float vert[4][3], len, end[3];
1232         static float epsilon[3]= {0, 0, 0}; //1e-4, 1e-4, 1e-4};
1233
1234         /* setup isec */
1235         memset(&isec, 0, sizeof(isec));
1236         isec.mode= RE_RAY_MIRROR; /* we want the closest intersection */
1237         isec.lay= -1;
1238         isec.labda= 1e10f;
1239
1240         VECADD(isec.start, co1, epsilon);
1241         VECADD(end, co2, epsilon);
1242         sub_v3_v3v3(isec.vec, end, isec.start);
1243
1244 #if 0
1245         /*if(RE_ray_tree_intersect(mdb->raytree, &isec)) {*/
1246 #endif
1247
1248         if(meshdeform_intersect(mdb, &isec)) {
1249                 len= isec.labda;
1250                 mface=(MFace*)isec.hit.face;
1251
1252                 /* create MDefBoundIsect */
1253                 isect= BLI_memarena_alloc(mdb->memarena, sizeof(*isect));
1254
1255                 /* compute intersection coordinate */
1256                 isect->co[0]= co1[0] + isec.vec[0]*len;
1257                 isect->co[1]= co1[1] + isec.vec[1]*len;
1258                 isect->co[2]= co1[2] + isec.vec[2]*len;
1259
1260                 isect->len= len_v3v3(co1, isect->co);
1261                 if(isect->len < MESHDEFORM_LEN_THRESHOLD)
1262                         isect->len= MESHDEFORM_LEN_THRESHOLD;
1263
1264                 isect->v[0]= mface->v1;
1265                 isect->v[1]= mface->v2;
1266                 isect->v[2]= mface->v3;
1267                 isect->v[3]= mface->v4;
1268                 isect->nvert= (mface->v4)? 4: 3;
1269
1270                 isect->facing= isec.isect;
1271
1272                 /* compute mean value coordinates for interpolation */
1273                 cagecos= mdb->cagecos;
1274                 copy_v3_v3(vert[0], cagecos[mface->v1]);
1275                 copy_v3_v3(vert[1], cagecos[mface->v2]);
1276                 copy_v3_v3(vert[2], cagecos[mface->v3]);
1277                 if(mface->v4) copy_v3_v3(vert[3], cagecos[mface->v4]);
1278                 interp_weights_poly_v3( isect->uvw,vert, isect->nvert, isect->co);
1279
1280                 return isect;
1281         }
1282
1283         return NULL;
1284 }
1285
1286 static int meshdeform_inside_cage(MeshDeformBind *mdb, float *co)
1287 {
1288         MDefBoundIsect *isect;
1289         float outside[3], start[3], dir[3];
1290         int i, counter;
1291
1292         for(i=1; i<=6; i++) {
1293                 counter = 0;
1294
1295                 outside[0] = co[0] + (mdb->max[0] - mdb->min[0] + 1.0f)*MESHDEFORM_OFFSET[i][0];
1296                 outside[1] = co[1] + (mdb->max[1] - mdb->min[1] + 1.0f)*MESHDEFORM_OFFSET[i][1];
1297                 outside[2] = co[2] + (mdb->max[2] - mdb->min[2] + 1.0f)*MESHDEFORM_OFFSET[i][2];
1298
1299                 copy_v3_v3(start, co);
1300                 sub_v3_v3v3(dir, outside, start);
1301                 normalize_v3(dir);
1302
1303                 isect = meshdeform_ray_tree_intersect(mdb, start, outside);
1304                 if(isect && !isect->facing)
1305                         return 1;
1306         }
1307
1308         return 0;
1309 }
1310
1311 /* solving */
1312
1313 static int meshdeform_index(MeshDeformBind *mdb, int x, int y, int z, int n)
1314 {
1315         int size= mdb->size;
1316
1317         x += MESHDEFORM_OFFSET[n][0];
1318         y += MESHDEFORM_OFFSET[n][1];
1319         z += MESHDEFORM_OFFSET[n][2];
1320
1321         if(x < 0 || x >= mdb->size)
1322                 return -1;
1323         if(y < 0 || y >= mdb->size)
1324                 return -1;
1325         if(z < 0 || z >= mdb->size)
1326                 return -1;
1327
1328         return x + y*size + z*size*size;
1329 }
1330
1331 static void meshdeform_cell_center(MeshDeformBind *mdb, int x, int y, int z, int n, float *center)
1332 {
1333         x += MESHDEFORM_OFFSET[n][0];
1334         y += MESHDEFORM_OFFSET[n][1];
1335         z += MESHDEFORM_OFFSET[n][2];
1336
1337         center[0]= mdb->min[0] + x*mdb->width[0] + mdb->halfwidth[0];
1338         center[1]= mdb->min[1] + y*mdb->width[1] + mdb->halfwidth[1];
1339         center[2]= mdb->min[2] + z*mdb->width[2] + mdb->halfwidth[2];
1340 }
1341
1342 static void meshdeform_add_intersections(MeshDeformBind *mdb, int x, int y, int z)
1343 {
1344         MDefBoundIsect *isect;
1345         float center[3], ncenter[3];
1346         int i, a;
1347
1348         a= meshdeform_index(mdb, x, y, z, 0);
1349         meshdeform_cell_center(mdb, x, y, z, 0, center);
1350
1351         /* check each outgoing edge for intersection */
1352         for(i=1; i<=6; i++) {
1353                 if(meshdeform_index(mdb, x, y, z, i) == -1)
1354                         continue;
1355
1356                 meshdeform_cell_center(mdb, x, y, z, i, ncenter);
1357
1358                 isect= meshdeform_ray_tree_intersect(mdb, center, ncenter);
1359                 if(isect) {
1360                         mdb->boundisect[a][i-1]= isect;
1361                         mdb->tag[a]= MESHDEFORM_TAG_BOUNDARY;
1362                 }
1363         }
1364 }
1365
1366 static void meshdeform_bind_floodfill(MeshDeformBind *mdb)
1367 {
1368         int *stack, *tag= mdb->tag;
1369         int a, b, i, xyz[3], stacksize, size= mdb->size;
1370
1371         stack= MEM_callocN(sizeof(int)*mdb->size3, "MeshDeformBindStack");
1372
1373         /* we know lower left corner is EXTERIOR because of padding */
1374         tag[0]= MESHDEFORM_TAG_EXTERIOR;
1375         stack[0]= 0;
1376         stacksize= 1;
1377
1378         /* floodfill exterior tag */
1379         while(stacksize > 0) {
1380                 a= stack[--stacksize];
1381
1382                 xyz[2]= a/(size*size);
1383                 xyz[1]= (a - xyz[2]*size*size)/size;
1384                 xyz[0]= a - xyz[1]*size - xyz[2]*size*size;
1385
1386                 for(i=1; i<=6; i++) {
1387                         b= meshdeform_index(mdb, xyz[0], xyz[1], xyz[2], i);
1388
1389                         if(b != -1) {
1390                                 if(tag[b] == MESHDEFORM_TAG_UNTYPED ||
1391                                    (tag[b] == MESHDEFORM_TAG_BOUNDARY && !mdb->boundisect[a][i-1])) {
1392                                         tag[b]= MESHDEFORM_TAG_EXTERIOR;
1393                                         stack[stacksize++]= b;
1394                                 }
1395                         }
1396                 }
1397         }
1398
1399         /* other cells are interior */
1400         for(a=0; a<size*size*size; a++)
1401                 if(tag[a]==MESHDEFORM_TAG_UNTYPED)
1402                         tag[a]= MESHDEFORM_TAG_INTERIOR;
1403
1404 #if 0
1405         {
1406                 int tb, ti, te, ts;
1407                 tb= ti= te= ts= 0;
1408                 for(a=0; a<size*size*size; a++)
1409                         if(tag[a]==MESHDEFORM_TAG_BOUNDARY)
1410                                 tb++;
1411                         else if(tag[a]==MESHDEFORM_TAG_INTERIOR)
1412                                 ti++;
1413                         else if(tag[a]==MESHDEFORM_TAG_EXTERIOR) {
1414                                 te++;
1415
1416                                 if(mdb->semibound[a])
1417                                         ts++;
1418                         }
1419
1420                 printf("interior %d exterior %d boundary %d semi-boundary %d\n", ti, te, tb, ts);
1421         }
1422 #endif
1423
1424         MEM_freeN(stack);
1425 }
1426
1427 static float meshdeform_boundary_phi(MeshDeformBind *mdb, MDefBoundIsect *isect, int cagevert)
1428 {
1429         int a;
1430
1431         for(a=0; a<isect->nvert; a++)
1432                 if(isect->v[a] == cagevert)
1433                         return isect->uvw[a];
1434
1435         return 0.0f;
1436 }
1437
1438 static float meshdeform_interp_w(MeshDeformBind *mdb, float *gridvec, float *vec, int cagevert)
1439 {
1440         float dvec[3], ivec[3], wx, wy, wz, result=0.0f;
1441         float weight, totweight= 0.0f;
1442         int i, a, x, y, z;
1443
1444         for(i=0; i<3; i++) {
1445                 ivec[i]= (int)gridvec[i];
1446                 dvec[i]= gridvec[i] - ivec[i];
1447         }
1448
1449         for(i=0; i<8; i++) {
1450                 if(i & 1) { x= ivec[0]+1; wx= dvec[0]; }
1451                 else { x= ivec[0]; wx= 1.0f-dvec[0]; }
1452
1453                 if(i & 2) { y= ivec[1]+1; wy= dvec[1]; }
1454                 else { y= ivec[1]; wy= 1.0f-dvec[1]; }
1455
1456                 if(i & 4) { z= ivec[2]+1; wz= dvec[2]; }
1457                 else { z= ivec[2]; wz= 1.0f-dvec[2]; }
1458
1459                 CLAMP(x, 0, mdb->size-1);
1460                 CLAMP(y, 0, mdb->size-1);
1461                 CLAMP(z, 0, mdb->size-1);
1462
1463                 a= meshdeform_index(mdb, x, y, z, 0);
1464                 weight= wx*wy*wz;
1465                 result += weight*mdb->phi[a];
1466                 totweight += weight;
1467         }
1468
1469         if(totweight > 0.0f)
1470                 result /= totweight;
1471
1472         return result;
1473 }
1474
1475 static void meshdeform_check_semibound(MeshDeformBind *mdb, int x, int y, int z)
1476 {
1477         int i, a;
1478
1479         a= meshdeform_index(mdb, x, y, z, 0);
1480         if(mdb->tag[a] != MESHDEFORM_TAG_EXTERIOR)
1481                 return;
1482
1483         for(i=1; i<=6; i++)
1484                 if(mdb->boundisect[a][i-1])
1485                         mdb->semibound[a]= 1;
1486 }
1487
1488 static float meshdeform_boundary_total_weight(MeshDeformBind *mdb, int x, int y, int z)
1489 {
1490         float weight, totweight= 0.0f;
1491         int i, a;
1492
1493         a= meshdeform_index(mdb, x, y, z, 0);
1494
1495         /* count weight for neighbour cells */
1496         for(i=1; i<=6; i++) {
1497                 if(meshdeform_index(mdb, x, y, z, i) == -1)
1498                         continue;
1499
1500                 if(mdb->boundisect[a][i-1])
1501                         weight= 1.0f/mdb->boundisect[a][i-1]->len;
1502                 else if(!mdb->semibound[a])
1503                         weight= 1.0f/mdb->width[0];
1504                 else
1505                         weight= 0.0f;
1506
1507                 totweight += weight;
1508         }
1509
1510         return totweight;
1511 }
1512
1513 static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, int x, int y, int z)
1514 {
1515         MDefBoundIsect *isect;
1516         float weight, totweight;
1517         int i, a, acenter;
1518
1519         acenter= meshdeform_index(mdb, x, y, z, 0);
1520         if(mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR)
1521                 return;
1522
1523         nlMatrixAdd(mdb->varidx[acenter], mdb->varidx[acenter], 1.0f);
1524
1525         totweight= meshdeform_boundary_total_weight(mdb, x, y, z);
1526         for(i=1; i<=6; i++) {
1527                 a= meshdeform_index(mdb, x, y, z, i);
1528                 if(a == -1 || mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR)
1529                         continue;
1530
1531                 isect= mdb->boundisect[acenter][i-1];
1532                 if (!isect) {
1533                         weight= (1.0f/mdb->width[0])/totweight;
1534                         nlMatrixAdd(mdb->varidx[acenter], mdb->varidx[a], -weight);
1535                 }
1536         }
1537 }
1538
1539 static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, int x, int y, int z, int cagevert)
1540 {
1541         MDefBoundIsect *isect;
1542         float rhs, weight, totweight;
1543         int i, a, acenter;
1544
1545         acenter= meshdeform_index(mdb, x, y, z, 0);
1546         if(mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR)
1547                 return;
1548
1549         totweight= meshdeform_boundary_total_weight(mdb, x, y, z);
1550         for(i=1; i<=6; i++) {
1551                 a= meshdeform_index(mdb, x, y, z, i);
1552                 if(a == -1)
1553                         continue;
1554
1555                 isect= mdb->boundisect[acenter][i-1];
1556
1557                 if (isect) {
1558                         weight= (1.0f/isect->len)/totweight;
1559                         rhs= weight*meshdeform_boundary_phi(mdb, isect, cagevert);
1560                         nlRightHandSideAdd(0, mdb->varidx[acenter], rhs);
1561                 }
1562         }
1563 }
1564
1565 static void meshdeform_matrix_add_semibound_phi(MeshDeformBind *mdb, int x, int y, int z, int cagevert)
1566 {
1567         MDefBoundIsect *isect;
1568         float rhs, weight, totweight;
1569         int i, a;
1570
1571         a= meshdeform_index(mdb, x, y, z, 0);
1572         if(!mdb->semibound[a])
1573                 return;
1574
1575         mdb->phi[a]= 0.0f;
1576
1577         totweight= meshdeform_boundary_total_weight(mdb, x, y, z);
1578         for(i=1; i<=6; i++) {
1579                 isect= mdb->boundisect[a][i-1];
1580
1581                 if (isect) {
1582                         weight= (1.0f/isect->len)/totweight;
1583                         rhs= weight*meshdeform_boundary_phi(mdb, isect, cagevert);
1584                         mdb->phi[a] += rhs;
1585                 }
1586         }
1587 }
1588
1589 static void meshdeform_matrix_add_exterior_phi(MeshDeformBind *mdb, int x, int y, int z, int cagevert)
1590 {
1591         float phi, totweight;
1592         int i, a, acenter;
1593
1594         acenter= meshdeform_index(mdb, x, y, z, 0);
1595         if(mdb->tag[acenter] != MESHDEFORM_TAG_EXTERIOR || mdb->semibound[acenter])
1596                 return;
1597
1598         phi= 0.0f;
1599         totweight= 0.0f;
1600         for(i=1; i<=6; i++) {
1601                 a= meshdeform_index(mdb, x, y, z, i);
1602
1603                 if(a != -1 && mdb->semibound[a]) {
1604                         phi += mdb->phi[a];
1605                         totweight += 1.0f;
1606                 }
1607         }
1608
1609         if(totweight != 0.0f)
1610                 mdb->phi[acenter]= phi/totweight;
1611 }
1612
1613 static void meshdeform_matrix_solve(MeshDeformBind *mdb)
1614 {
1615         NLContext *context;
1616         float vec[3], gridvec[3];
1617         int a, b, x, y, z, totvar;
1618         char message[1024];
1619
1620         /* setup variable indices */
1621         mdb->varidx= MEM_callocN(sizeof(int)*mdb->size3, "MeshDeformDSvaridx");
1622         for(a=0, totvar=0; a<mdb->size3; a++)
1623                 mdb->varidx[a]= (mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR)? -1: totvar++;
1624
1625         if(totvar == 0) {
1626                 MEM_freeN(mdb->varidx);
1627                 return;
1628         }
1629
1630         progress_bar(0, "Starting mesh deform solve");
1631
1632         /* setup opennl solver */
1633         nlNewContext();
1634         context= nlGetCurrent();
1635
1636         nlSolverParameteri(NL_NB_VARIABLES, totvar);
1637         nlSolverParameteri(NL_NB_ROWS, totvar);
1638         nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 1);
1639
1640         nlBegin(NL_SYSTEM);
1641         nlBegin(NL_MATRIX);
1642
1643         /* build matrix */
1644         for(z=0; z<mdb->size; z++)
1645                 for(y=0; y<mdb->size; y++)
1646                         for(x=0; x<mdb->size; x++)
1647                                 meshdeform_matrix_add_cell(mdb, x, y, z);
1648
1649         /* solve for each cage vert */
1650         for(a=0; a<mdb->totcagevert; a++) {
1651                 if(a != 0) {
1652                         nlBegin(NL_SYSTEM);
1653                         nlBegin(NL_MATRIX);
1654                 }
1655
1656                 /* fill in right hand side and solve */
1657                 for(z=0; z<mdb->size; z++)
1658                         for(y=0; y<mdb->size; y++)
1659                                 for(x=0; x<mdb->size; x++)
1660                                         meshdeform_matrix_add_rhs(mdb, x, y, z, a);
1661
1662                 nlEnd(NL_MATRIX);
1663                 nlEnd(NL_SYSTEM);
1664
1665 #if 0
1666                 nlPrintMatrix();
1667 #endif
1668
1669                 if(nlSolveAdvanced(NULL, NL_TRUE)) {
1670                         for(z=0; z<mdb->size; z++)
1671                                 for(y=0; y<mdb->size; y++)
1672                                         for(x=0; x<mdb->size; x++)
1673                                                 meshdeform_matrix_add_semibound_phi(mdb, x, y, z, a);
1674
1675                         for(z=0; z<mdb->size; z++)
1676                                 for(y=0; y<mdb->size; y++)
1677                                         for(x=0; x<mdb->size; x++)
1678                                                 meshdeform_matrix_add_exterior_phi(mdb, x, y, z, a);
1679
1680                         for(b=0; b<mdb->size3; b++) {
1681                                 if(mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR)
1682                                         mdb->phi[b]= nlGetVariable(0, mdb->varidx[b]);
1683                                 mdb->totalphi[b] += mdb->phi[b];
1684                         }
1685
1686                         if(mdb->weights) {
1687                                 /* static bind : compute weights for each vertex */
1688                                 for(b=0; b<mdb->totvert; b++) {
1689                                         if(mdb->inside[b]) {
1690                                                 copy_v3_v3(vec, mdb->vertexcos[b]);
1691                                                 gridvec[0]= (vec[0] - mdb->min[0] - mdb->halfwidth[0])/mdb->width[0];
1692                                                 gridvec[1]= (vec[1] - mdb->min[1] - mdb->halfwidth[1])/mdb->width[1];
1693                                                 gridvec[2]= (vec[2] - mdb->min[2] - mdb->halfwidth[2])/mdb->width[2];
1694
1695                                                 mdb->weights[b*mdb->totcagevert + a]= meshdeform_interp_w(mdb, gridvec, vec, a);
1696                                         }
1697                                 }
1698                         }
1699                         else {
1700                                 MDefBindInfluence *inf;
1701
1702                                 /* dynamic bind */
1703                                 for(b=0; b<mdb->size3; b++) {
1704                                         if(mdb->phi[b] >= MESHDEFORM_MIN_INFLUENCE) {
1705                                                 inf= BLI_memarena_alloc(mdb->memarena, sizeof(*inf));
1706                                                 inf->vertex= a;
1707                                                 inf->weight= mdb->phi[b];
1708                                                 inf->next= mdb->dyngrid[b];
1709                                                 mdb->dyngrid[b]= inf;
1710                                         }
1711                                 }
1712                         }
1713                 }
1714                 else {
1715                         error("Mesh Deform: failed to find solution.");
1716                         break;
1717                 }
1718
1719                 sprintf(message, "Mesh deform solve %d / %d       |||", a+1, mdb->totcagevert);
1720                 progress_bar((float)(a+1)/(float)(mdb->totcagevert), message);
1721         }
1722
1723 #if 0
1724         /* sanity check */
1725         for(b=0; b<mdb->size3; b++)
1726                 if(mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR)
1727                         if(fabs(mdb->totalphi[b] - 1.0f) > 1e-4)
1728                                 printf("totalphi deficiency [%s|%d] %d: %.10f\n",
1729                                         (mdb->tag[b] == MESHDEFORM_TAG_INTERIOR)? "interior": "boundary", mdb->semibound[b], mdb->varidx[b], mdb->totalphi[b]);
1730 #endif
1731
1732         /* free */
1733         MEM_freeN(mdb->varidx);
1734
1735         nlDeleteContext(context);
1736 }
1737
1738 static void harmonic_coordinates_bind(Scene *scene, MeshDeformModifierData *mmd, MeshDeformBind *mdb)
1739 {
1740         MDefBindInfluence *inf;
1741         MDefInfluence *mdinf;
1742         MDefCell *cell;
1743         float center[3], vec[3], maxwidth, totweight;
1744         int a, b, x, y, z, totinside, offset;
1745
1746         /* compute bounding box of the cage mesh */
1747         INIT_MINMAX(mdb->min, mdb->max);
1748
1749         for(a=0; a<mdb->totcagevert; a++)
1750                 DO_MINMAX(mdb->cagecos[a], mdb->min, mdb->max);
1751
1752         /* allocate memory */
1753         mdb->size= (2<<(mmd->gridsize-1)) + 2;
1754         mdb->size3= mdb->size*mdb->size*mdb->size;
1755         mdb->tag= MEM_callocN(sizeof(int)*mdb->size3, "MeshDeformBindTag");
1756         mdb->phi= MEM_callocN(sizeof(float)*mdb->size3, "MeshDeformBindPhi");
1757         mdb->totalphi= MEM_callocN(sizeof(float)*mdb->size3, "MeshDeformBindTotalPhi");
1758         mdb->boundisect= MEM_callocN(sizeof(*mdb->boundisect)*mdb->size3, "MDefBoundIsect");
1759         mdb->semibound= MEM_callocN(sizeof(int)*mdb->size3, "MDefSemiBound");
1760
1761         mdb->inside= MEM_callocN(sizeof(int)*mdb->totvert, "MDefInside");
1762
1763         if(mmd->flag & MOD_MDEF_DYNAMIC_BIND)
1764                 mdb->dyngrid= MEM_callocN(sizeof(MDefBindInfluence*)*mdb->size3, "MDefDynGrid");
1765         else
1766                 mdb->weights= MEM_callocN(sizeof(float)*mdb->totvert*mdb->totcagevert, "MDefWeights");
1767
1768         mdb->memarena= BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "harmonic coords arena");
1769         BLI_memarena_use_calloc(mdb->memarena);
1770
1771         /* make bounding box equal size in all directions, add padding, and compute
1772          * width of the cells */
1773         maxwidth = -1.0f;
1774         for(a=0; a<3; a++)
1775                 if(mdb->max[a]-mdb->min[a] > maxwidth)
1776                         maxwidth= mdb->max[a]-mdb->min[a];
1777
1778         for(a=0; a<3; a++) {
1779                 center[a]= (mdb->min[a]+mdb->max[a])*0.5f;
1780                 mdb->min[a]= center[a] - maxwidth*0.5f;
1781                 mdb->max[a]= center[a] + maxwidth*0.5f;
1782
1783                 mdb->width[a]= (mdb->max[a]-mdb->min[a])/(mdb->size-4);
1784                 mdb->min[a] -= 2.1f*mdb->width[a];
1785                 mdb->max[a] += 2.1f*mdb->width[a];
1786
1787                 mdb->width[a]= (mdb->max[a]-mdb->min[a])/mdb->size;
1788                 mdb->halfwidth[a]= mdb->width[a]*0.5f;
1789         }
1790
1791         progress_bar(0, "Setting up mesh deform system");
1792
1793 #if 0
1794         /* create ray tree */
1795         meshdeform_ray_tree_create(mdb);
1796 #endif
1797
1798         totinside= 0;
1799         for(a=0; a<mdb->totvert; a++) {
1800                 copy_v3_v3(vec, mdb->vertexcos[a]);
1801                 mdb->inside[a]= meshdeform_inside_cage(mdb, vec);
1802                 if(mdb->inside[a])
1803                         totinside++;
1804         }
1805
1806         /* free temporary MDefBoundIsects */
1807         BLI_memarena_free(mdb->memarena);
1808         mdb->memarena= BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "harmonic coords arena");
1809
1810         /* start with all cells untyped */
1811         for(a=0; a<mdb->size3; a++)
1812                 mdb->tag[a]= MESHDEFORM_TAG_UNTYPED;
1813
1814         /* detect intersections and tag boundary cells */
1815         for(z=0; z<mdb->size; z++)
1816                 for(y=0; y<mdb->size; y++)
1817                         for(x=0; x<mdb->size; x++)
1818                                 meshdeform_add_intersections(mdb, x, y, z);
1819
1820 #if 0
1821         /* free ray tree */
1822         meshdeform_ray_tree_free(mdb);
1823 #endif
1824
1825         /* compute exterior and interior tags */
1826         meshdeform_bind_floodfill(mdb);
1827
1828         for(z=0; z<mdb->size; z++)
1829                 for(y=0; y<mdb->size; y++)
1830                         for(x=0; x<mdb->size; x++)
1831                                 meshdeform_check_semibound(mdb, x, y, z);
1832
1833         /* solve */
1834         meshdeform_matrix_solve(mdb);
1835
1836         /* assign results */
1837         if(mmd->flag & MOD_MDEF_DYNAMIC_BIND) {
1838                 mmd->totinfluence= 0;
1839                 for(a=0; a<mdb->size3; a++)
1840                         for(inf=mdb->dyngrid[a]; inf; inf=inf->next)
1841                                 mmd->totinfluence++;
1842
1843                 /* convert MDefBindInfluences to smaller MDefInfluences */
1844                 mmd->dyngrid= MEM_callocN(sizeof(MDefCell)*mdb->size3, "MDefDynGrid");
1845                 mmd->dyninfluences= MEM_callocN(sizeof(MDefInfluence)*mmd->totinfluence, "MDefInfluence");
1846                 offset= 0;
1847                 for(a=0; a<mdb->size3; a++) {
1848                         cell= &mmd->dyngrid[a];
1849                         cell->offset= offset;
1850
1851                         totweight= 0.0f;
1852                         mdinf= mmd->dyninfluences + cell->offset;
1853                         for(inf=mdb->dyngrid[a]; inf; inf=inf->next, mdinf++) {
1854                                 mdinf->weight= inf->weight;
1855                                 mdinf->vertex= inf->vertex;
1856                                 totweight += mdinf->weight;
1857                                 cell->totinfluence++;
1858                         }
1859
1860                         if(totweight > 0.0f) {
1861                                 mdinf= mmd->dyninfluences + cell->offset;
1862                                 for(b=0; b<cell->totinfluence; b++, mdinf++)
1863                                         mdinf->weight /= totweight;
1864                         }
1865
1866                         offset += cell->totinfluence;
1867                 }
1868
1869                 mmd->dynverts= mdb->inside;
1870                 mmd->dyngridsize= mdb->size;
1871                 copy_v3_v3(mmd->dyncellmin, mdb->min);
1872                 mmd->dyncellwidth= mdb->width[0];
1873                 MEM_freeN(mdb->dyngrid);
1874         }
1875         else {
1876                 mmd->bindweights= mdb->weights;
1877                 MEM_freeN(mdb->inside);
1878         }
1879
1880         MEM_freeN(mdb->tag);
1881         MEM_freeN(mdb->phi);
1882         MEM_freeN(mdb->totalphi);
1883         MEM_freeN(mdb->boundisect);
1884         MEM_freeN(mdb->semibound);
1885         BLI_memarena_free(mdb->memarena);
1886 }
1887
1888 #if 0
1889 static void heat_weighting_bind(Scene *scene, DerivedMesh *dm, MeshDeformModifierData *mmd, MeshDeformBind *mdb)
1890 {
1891         LaplacianSystem *sys;
1892         MFace *mface= dm->getFaceArray(dm), *mf;
1893         int totvert= dm->getNumVerts(dm);
1894         int totface= dm->getNumFaces(dm);
1895         float solution, weight;
1896         int a, tottri, j, thrownerror = 0;
1897
1898         mdb->weights= MEM_callocN(sizeof(float)*mdb->totvert*mdb->totcagevert, "MDefWeights");
1899
1900         /* count triangles */
1901         for(tottri=0, a=0, mf=mface; a<totface; a++, mf++) {
1902                 tottri++;
1903                 if(mf->v4) tottri++;
1904         }
1905
1906         /* create laplacian */
1907         sys = laplacian_system_construct_begin(totvert, tottri, 1);
1908
1909         sys->heat.mface= mface;
1910         sys->heat.totface= totface;
1911         sys->heat.totvert= totvert;
1912         sys->heat.verts= mdb->vertexcos;
1913         sys->heat.source = mdb->cagecos;
1914         sys->heat.numsource= mdb->totcagevert;
1915
1916         heat_ray_tree_create(sys);
1917         heat_laplacian_create(sys);
1918
1919         laplacian_system_construct_end(sys);
1920
1921         /* compute weights per bone */
1922         for(j=0; j<mdb->totcagevert; j++) {
1923                 /* fill right hand side */
1924                 laplacian_begin_solve(sys, -1);
1925
1926                 for(a=0; a<totvert; a++)
1927                         if(heat_source_closest(sys, a, j))
1928                                 laplacian_add_right_hand_side(sys, a,
1929                                         sys->heat.H[a]*sys->heat.p[a]);
1930
1931                 /* solve */
1932                 if(laplacian_system_solve(sys)) {
1933                         /* load solution into vertex groups */
1934                         for(a=0; a<totvert; a++) {
1935                                 solution= laplacian_system_get_solution(a);
1936
1937                                 weight= heat_limit_weight(solution);
1938                                 if(weight > 0.0f)
1939                                         mdb->weights[a*mdb->totcagevert + j] = weight;
1940                         }
1941                 }
1942                 else if(!thrownerror) {
1943                         error("Mesh Deform Heat Weighting:"
1944                                 " failed to find solution for one or more vertices");
1945                         thrownerror= 1;
1946                         break;
1947                 }
1948         }
1949
1950         /* free */
1951         heat_system_free(sys);
1952         laplacian_system_delete(sys);
1953
1954         mmd->bindweights= mdb->weights;
1955 }
1956 #endif
1957
1958 void mesh_deform_bind(Scene *scene, DerivedMesh *dm, MeshDeformModifierData *mmd, float *vertexcos, int totvert, float cagemat[][4])
1959 {
1960         MeshDeformBind mdb;
1961         MVert *mvert;
1962         int a;
1963
1964         waitcursor(1);
1965         start_progress_bar();
1966
1967         memset(&mdb, 0, sizeof(MeshDeformBind));
1968
1969         /* get mesh and cage mesh */
1970         mdb.vertexcos= MEM_callocN(sizeof(float)*3*totvert, "MeshDeformCos");
1971         mdb.totvert= totvert;
1972
1973         mdb.cagedm= mesh_create_derived_no_deform(scene, mmd->object, NULL, CD_MASK_BAREMESH);
1974         mdb.totcagevert= mdb.cagedm->getNumVerts(mdb.cagedm);
1975         mdb.cagecos= MEM_callocN(sizeof(*mdb.cagecos)*mdb.totcagevert, "MeshDeformBindCos");
1976         copy_m4_m4(mdb.cagemat, cagemat);
1977
1978         mvert= mdb.cagedm->getVertArray(mdb.cagedm);
1979         for(a=0; a<mdb.totcagevert; a++)
1980                 copy_v3_v3(mdb.cagecos[a], mvert[a].co);
1981         for(a=0; a<mdb.totvert; a++)
1982                 mul_v3_m4v3(mdb.vertexcos[a], mdb.cagemat, vertexcos + a*3);
1983
1984         /* solve */
1985 #if 0
1986         if(mmd->mode == MOD_MDEF_VOLUME)
1987                 harmonic_coordinates_bind(scene, mmd, &mdb);
1988         else
1989                 heat_weighting_bind(scene, dm, mmd, &mdb);
1990 #else
1991         harmonic_coordinates_bind(scene, mmd, &mdb);
1992 #endif
1993
1994         /* assign bind variables */
1995         mmd->bindcagecos= (float*)mdb.cagecos;
1996         mmd->totvert= mdb.totvert;
1997         mmd->totcagevert= mdb.totcagevert;
1998         copy_m4_m4(mmd->bindmat, mmd->object->obmat);
1999
2000         /* transform bindcagecos to world space */
2001         for(a=0; a<mdb.totcagevert; a++)
2002                 mul_m4_v3(mmd->object->obmat, mmd->bindcagecos+a*3);
2003
2004         /* free */
2005         mdb.cagedm->release(mdb.cagedm);
2006         MEM_freeN(mdb.vertexcos);
2007
2008         /* compact weights */
2009         modifier_mdef_compact_influences((ModifierData*)mmd);
2010
2011         end_progress_bar();
2012         waitcursor(0);
2013 }
2014