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