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