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