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