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