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