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