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