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