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