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