4 * ***** BEGIN GPL LICENSE BLOCK *****
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.
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.
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.
20 * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
21 * All rights reserved.
23 * The Original Code is: all of this file.
25 * Contributor(s): none yet.
27 * ***** END GPL LICENSE BLOCK *****
28 * meshlaplacian.c: Algorithms using the mesh laplacian.
34 #include "MEM_guardedalloc.h"
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"
43 #include "BLI_arithb.h"
44 #include "BLI_edgehash.h"
45 #include "BLI_memarena.h"
47 #include "BKE_DerivedMesh.h"
48 #include "BKE_utildefines.h"
51 #include "BLI_editVert.h"
52 #include "BLI_polardecomp.h"
55 #include "RE_raytrace.h"
57 #include "ONL_opennl.h"
59 #include "BLO_sys_types.h" // for intptr_t support
61 #include "ED_armature.h"
64 #include "meshlaplacian.h"
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 *************** */
76 /************************** Laplacian System *****************************/
78 struct LaplacianSystem {
79 NLContext context; /* opennl context */
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 */
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 */
93 EdgeHash *edgehash; /* edge hash for construction */
95 struct HeatWeighting {
97 float (*verts)[3]; /* vertex coordinates */
98 float (*vnors)[3]; /* vertex normals */
100 float (*root)[3]; /* bone root */
101 float (*tip)[3]; /* bone tip */
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 */
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 */
114 struct RigidDeformation {
125 /* Laplacian matrix construction */
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).
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. */
135 static void laplacian_increase_edge_count(EdgeHash *edgehash, int v1, int v2)
137 void **p = BLI_edgehash_lookup_p(edgehash, v1, v2);
140 *p = (void*)((intptr_t)*p + (intptr_t)1);
142 BLI_edgehash_insert(edgehash, v1, v2, (void*)(intptr_t)1);
145 static int laplacian_edge_count(EdgeHash *edgehash, int v1, int v2)
147 return (int)(intptr_t)BLI_edgehash_lookup(edgehash, v1, v2);
150 static float cotan_weight(float *v1, float *v2, float *v3)
152 float a[3], b[3], c[3], clen;
163 return Inpf(a, b)/clen;
166 static void laplacian_triangle_area(LaplacianSystem *sys, int i1, int i2, int i3)
168 float t1, t2, t3, len1, len2, len3, area;
169 float *varea= sys->varea, *v1, *v2, *v3;
176 t1= cotan_weight(v1, v2, v3);
177 t2= cotan_weight(v2, v3, v1);
178 t3= cotan_weight(v3, v1, v2);
180 if(RAD2DEG(VecAngle3(v2, v1, v3)) > 90) obtuse= 1;
181 else if(RAD2DEG(VecAngle3(v1, v2, v3)) > 90) obtuse= 2;
182 else if(RAD2DEG(VecAngle3(v1, v3, v2)) > 90) obtuse= 3;
185 area= AreaT3Dfl(v1, v2, v3);
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;
192 len1= VecLenf(v2, v3);
193 len2= VecLenf(v1, v3);
194 len3= VecLenf(v1, v2);
200 varea[i1] += (t2 + t3)*0.25f;
201 varea[i2] += (t1 + t3)*0.25f;
202 varea[i3] += (t1 + t2)*0.25f;
206 static void laplacian_triangle_weights(LaplacianSystem *sys, int f, int i1, int i2, int i3)
209 float *varea= sys->varea, *v1, *v2, *v3;
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);
221 nlMatrixAdd(i1, i1, (t2+t3)*varea[i1]);
222 nlMatrixAdd(i2, i2, (t1+t3)*varea[i2]);
223 nlMatrixAdd(i3, i3, (t1+t2)*varea[i3]);
225 nlMatrixAdd(i1, i2, -t3*varea[i1]);
226 nlMatrixAdd(i2, i1, -t3*varea[i2]);
228 nlMatrixAdd(i2, i3, -t1*varea[i2]);
229 nlMatrixAdd(i3, i2, -t1*varea[i3]);
231 nlMatrixAdd(i3, i1, -t2*varea[i3]);
232 nlMatrixAdd(i1, i3, -t2*varea[i1]);
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];
241 LaplacianSystem *laplacian_system_construct_begin(int totvert, int totface, int lsq)
243 LaplacianSystem *sys;
245 sys= MEM_callocN(sizeof(LaplacianSystem), "LaplacianSystem");
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");
255 sys->storeweights= 0;
257 /* create opennl context */
259 nlSolverParameteri(NL_NB_VARIABLES, totvert);
261 nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
263 sys->context= nlGetCurrent();
268 void laplacian_add_vertex(LaplacianSystem *sys, float *co, int pinned)
270 sys->verts[sys->totvert]= co;
271 sys->vpinned[sys->totvert]= pinned;
275 void laplacian_add_triangle(LaplacianSystem *sys, int v1, int v2, int v3)
277 sys->faces[sys->totface][0]= v1;
278 sys->faces[sys->totface][1]= v2;
279 sys->faces[sys->totface][2]= v3;
283 void laplacian_system_construct_end(LaplacianSystem *sys)
286 int a, totvert=sys->totvert, totface=sys->totface;
288 laplacian_begin_solve(sys, 0);
290 sys->varea= MEM_callocN(sizeof(float)*totvert, "LaplacianSystemVarea");
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]);
300 for(a=0, face=sys->faces; a<sys->totface; a++, face++)
301 laplacian_triangle_area(sys, (*face)[0], (*face)[1], (*face)[2]);
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];
311 /* for heat weighting */
313 nlMatrixAdd(a, a, sys->heat.H[a]);
316 if(sys->storeweights)
317 sys->fweights= MEM_callocN(sizeof(float)*3*totface, "LaplacianFWeight");
319 for(a=0, face=sys->faces; a<totface; a++, face++)
320 laplacian_triangle_weights(sys, a, (*face)[0], (*face)[1], (*face)[2]);
322 MEM_freeN(sys->faces);
326 MEM_freeN(sys->varea);
330 BLI_edgehash_free(sys->edgehash, NULL);
334 void laplacian_system_delete(LaplacianSystem *sys)
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);
342 nlDeleteContext(sys->context);
346 void laplacian_begin_solve(LaplacianSystem *sys, int index)
354 for(a=0; a<sys->totvert; a++) {
355 if(sys->vpinned[a]) {
356 nlSetVariable(0, a, sys->verts[a][index]);
367 void laplacian_add_right_hand_side(LaplacianSystem *sys, int v, float value)
369 nlRightHandSideAdd(0, v, value);
372 int laplacian_system_solve(LaplacianSystem *sys)
380 return nlSolveAdvanced(NULL, NL_TRUE);
383 float laplacian_system_get_solution(int v)
385 return nlGetVariable(0, v);
388 /************************* Heat Bone Weighting ******************************/
389 /* From "Automatic Rigging and Animation of 3D Characters"
390 Ilya Baran and Jovan Popovic, SIGGRAPH 2007 */
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
397 /* Raytracing for vertex to bone visibility */
398 static void heat_ray_tree_create(LaplacianSystem *sys)
400 Mesh *me = sys->heat.mesh;
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");
407 for(a=0; a<me->totface; a++) {
409 MFace *mface = me->mface+a;
410 RayFace *rayface = sys->heat.faces+a;
412 RayObject *obj = RE_rayface_from_coords(
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
417 RE_rayobject_add(sys->heat.raytree, obj);
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;
425 RE_rayobject_done(sys->heat.raytree);
428 static int heat_ray_bone_visible(LaplacianSystem *sys, int vertex, int bone)
435 mface= sys->heat.vface[vertex];
440 memset(&isec, 0, sizeof(isec));
441 isec.mode= RE_RAY_SHADOW;
443 isec.orig.ob = sys->heat.mesh;
444 isec.orig.face = mface;
445 isec.skip = RE_SKIP_CULLFACE;
448 VECCOPY(isec.start, sys->heat.verts[vertex]);
449 PclosestVL3Dfl(end, isec.start, sys->heat.root[bone], sys->heat.tip[bone]);
451 VECSUB(isec.vec, end, isec.start);
452 isec.labda = 1.0f - 1e-5;
453 VECADDFAC( isec.start, isec.start, isec.vec, 1e-5);
455 visible= !RE_rayobject_raycast(sys->heat.raytree, &isec);
460 static float heat_bone_distance(LaplacianSystem *sys, int vertex, int bone)
462 float closest[3], d[3], dist, cosine;
464 /* compute euclidian distance */
465 PclosestVL3Dfl(closest, sys->heat.verts[vertex],
466 sys->heat.root[bone], sys->heat.tip[bone]);
468 VecSubf(d, sys->heat.verts[vertex], closest);
471 /* if the vertex normal does not point along the bone, increase distance */
472 cosine= INPR(d, sys->heat.vnors[vertex]);
474 return dist/(0.5f*(cosine + 1.001f));
477 static int heat_bone_closest(LaplacianSystem *sys, int vertex, int bone)
481 dist= heat_bone_distance(sys, vertex, bone);
483 if(dist <= sys->heat.mindist[vertex]*(1.0f + DISTANCE_EPSILON))
484 if(heat_ray_bone_visible(sys, vertex, bone))
490 static void heat_set_H(LaplacianSystem *sys, int vertex)
492 float dist, mindist, h;
493 int j, numclosest = 0;
497 /* compute minimum distance */
498 for(j=0; j<sys->heat.numbones; j++) {
499 dist= heat_bone_distance(sys, vertex, j);
505 sys->heat.mindist[vertex]= mindist;
507 /* count number of bones with approximately this minimum distance */
508 for(j=0; j<sys->heat.numbones; j++)
509 if(heat_bone_closest(sys, vertex, j))
512 sys->heat.p[vertex]= (numclosest > 0)? 1.0f/numclosest: 0.0f;
514 /* compute H entry */
517 h= numclosest*C_WEIGHT/(mindist*mindist);
524 sys->heat.H[vertex]= h;
527 void heat_calc_vnormals(LaplacianSystem *sys)
530 int a, v1, v2, v3, (*face)[3];
532 sys->heat.vnors= MEM_callocN(sizeof(float)*3*sys->totvert, "HeatVNors");
534 for(a=0, face=sys->faces; a<sys->totface; a++, face++) {
539 CalcNormFloat(sys->verts[v1], sys->verts[v2], sys->verts[v3], fnor);
541 VecAddf(sys->heat.vnors[v1], sys->heat.vnors[v1], fnor);
542 VecAddf(sys->heat.vnors[v2], sys->heat.vnors[v2], fnor);
543 VecAddf(sys->heat.vnors[v3], sys->heat.vnors[v3], fnor);
546 for(a=0; a<sys->totvert; a++)
547 Normalize(sys->heat.vnors[a]);
550 static void heat_laplacian_create(LaplacianSystem *sys)
552 Mesh *me = sys->heat.mesh;
556 /* heat specific definitions */
557 sys->heat.mindist= MEM_callocN(sizeof(float)*me->totvert, "HeatMinDist");
558 sys->heat.H= MEM_callocN(sizeof(float)*me->totvert, "HeatH");
559 sys->heat.p= MEM_callocN(sizeof(float)*me->totvert, "HeatP");
561 /* add verts and faces to laplacian */
562 for(a=0; a<me->totvert; a++)
563 laplacian_add_vertex(sys, sys->heat.verts[a], 0);
565 for(a=0, mface=me->mface; a<me->totface; a++, mface++) {
566 laplacian_add_triangle(sys, mface->v1, mface->v2, mface->v3);
568 laplacian_add_triangle(sys, mface->v1, mface->v3, mface->v4);
571 /* for distance computation in set_H */
572 heat_calc_vnormals(sys);
574 for(a=0; a<me->totvert; a++)
578 static float heat_limit_weight(float weight)
582 if(weight < WEIGHT_LIMIT_END) {
585 else if(weight < WEIGHT_LIMIT_START) {
586 t= (weight - WEIGHT_LIMIT_END)/(WEIGHT_LIMIT_START - WEIGHT_LIMIT_END);
587 return t*WEIGHT_LIMIT_START;
593 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)
595 LaplacianSystem *sys;
597 float solution, weight;
598 int *vertsflipped = NULL;
599 int a, totface, j, bbone, firstsegment, lastsegment, thrownerror = 0;
601 /* count triangles */
602 for(totface=0, a=0, mface=me->mface; a<me->totface; a++, mface++) {
604 if(mface->v4) totface++;
607 /* create laplacian */
608 sys = laplacian_system_construct_begin(me->totvert, totface, 1);
611 sys->heat.verts= verts;
612 sys->heat.root= root;
614 sys->heat.numbones= numbones;
616 heat_ray_tree_create(sys);
617 heat_laplacian_create(sys);
619 laplacian_system_construct_end(sys);
622 vertsflipped = MEM_callocN(sizeof(int)*me->totvert, "vertsflipped");
623 for(a=0; a<me->totvert; a++)
624 vertsflipped[a] = mesh_get_x_mirror_vert(ob, a);
627 /* compute weights per bone */
628 for(j=0; j<numbones; j++) {
632 firstsegment= (j == 0 || dgrouplist[j-1] != dgrouplist[j]);
633 lastsegment= (j == numbones-1 || dgrouplist[j] != dgrouplist[j+1]);
634 bbone= !(firstsegment && lastsegment);
637 if(bbone && firstsegment) {
638 for(a=0; a<me->totvert; a++) {
639 ED_vgroup_vert_remove(ob, dgrouplist[j], a);
640 if(vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0)
641 ED_vgroup_vert_remove(ob, dgroupflip[j], vertsflipped[a]);
645 /* fill right hand side */
646 laplacian_begin_solve(sys, -1);
648 for(a=0; a<me->totvert; a++)
649 if(heat_bone_closest(sys, a, j))
650 laplacian_add_right_hand_side(sys, a,
651 sys->heat.H[a]*sys->heat.p[a]);
654 if(laplacian_system_solve(sys)) {
655 /* load solution into vertex groups */
656 for(a=0; a<me->totvert; a++) {
657 solution= laplacian_system_get_solution(a);
661 ED_vgroup_vert_add(ob, dgrouplist[j], a, solution,
665 weight= heat_limit_weight(solution);
667 ED_vgroup_vert_add(ob, dgrouplist[j], a, weight,
670 ED_vgroup_vert_remove(ob, dgrouplist[j], a);
673 /* do same for mirror */
674 if(vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0) {
677 ED_vgroup_vert_add(ob, dgroupflip[j], vertsflipped[a],
678 solution, WEIGHT_ADD);
681 weight= heat_limit_weight(solution);
683 ED_vgroup_vert_add(ob, dgroupflip[j], vertsflipped[a],
684 weight, WEIGHT_REPLACE);
686 ED_vgroup_vert_remove(ob, dgroupflip[j], vertsflipped[a]);
691 else if(!thrownerror) {
692 error("Bone Heat Weighting:"
693 " failed to find solution for one or more bones");
698 /* remove too small vertex weights */
699 if(bbone && lastsegment) {
700 for(a=0; a<me->totvert; a++) {
701 weight= ED_vgroup_vert_weight(ob, dgrouplist[j], a);
702 weight= heat_limit_weight(weight);
704 ED_vgroup_vert_remove(ob, dgrouplist[j], a);
706 if(vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0) {
707 weight= ED_vgroup_vert_weight(ob, dgroupflip[j], vertsflipped[a]);
708 weight= heat_limit_weight(weight);
710 ED_vgroup_vert_remove(ob, dgroupflip[j], vertsflipped[a]);
717 if(vertsflipped) MEM_freeN(vertsflipped);
719 RE_rayobject_free(sys->heat.raytree);
720 MEM_freeN(sys->heat.vface);
721 MEM_freeN(sys->heat.faces);
723 MEM_freeN(sys->heat.mindist);
724 MEM_freeN(sys->heat.H);
725 MEM_freeN(sys->heat.p);
726 MEM_freeN(sys->heat.vnors);
728 laplacian_system_delete(sys);
732 /********************** As-Rigid-As-Possible Deformation ******************/
733 /* From "As-Rigid-As-Possible Surface Modeling",
734 Olga Sorkine and Marc Alexa, ESGP 2007. */
737 - transpose R in orthogonal
738 - flipped normals and per face adding
739 - move cancelling to transform, make origco pointer
742 static LaplacianSystem *RigidDeformSystem = NULL;
744 static void rigid_add_half_edge_to_R(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
749 VecSubf(e, sys->rigid.origco[v1->tmp.l], sys->rigid.origco[v2->tmp.l]);
750 VecSubf(e_, v1->co, v2->co);
753 for (i=0; i<3; i++) {
754 sys->rigid.R[v1->tmp.l][i][0] += w*e[0]*e_[i];
755 sys->rigid.R[v1->tmp.l][i][1] += w*e[1]*e_[i];
756 sys->rigid.R[v1->tmp.l][i][2] += w*e[2]*e_[i];
760 static void rigid_add_edge_to_R(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
762 rigid_add_half_edge_to_R(sys, v1, v2, w);
763 rigid_add_half_edge_to_R(sys, v2, v1, w);
766 static void rigid_orthogonalize_R(float R[][3])
771 polar_decomp(M, Q, S);
775 static void rigid_add_half_edge_to_rhs(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
778 float Rsum[3][3], rhs[3];
780 if (sys->vpinned[v1->tmp.l])
783 Mat3AddMat3(Rsum, sys->rigid.R[v1->tmp.l], sys->rigid.R[v2->tmp.l]);
786 VecSubf(rhs, sys->rigid.origco[v1->tmp.l], sys->rigid.origco[v2->tmp.l]);
787 Mat3MulVecfl(Rsum, rhs);
791 VecAddf(sys->rigid.rhs[v1->tmp.l], sys->rigid.rhs[v1->tmp.l], rhs);
794 static void rigid_add_edge_to_rhs(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w)
796 rigid_add_half_edge_to_rhs(sys, v1, v2, w);
797 rigid_add_half_edge_to_rhs(sys, v2, v1, w);
800 void rigid_deform_iteration()
802 LaplacianSystem *sys= RigidDeformSystem;
811 nlMakeCurrent(sys->context);
815 memset(sys->rigid.R, 0, sizeof(float)*3*3*sys->totvert);
816 memset(sys->rigid.rhs, 0, sizeof(float)*3*sys->totvert);
818 for(a=0, efa=em->faces.first; efa; efa=efa->next, a++) {
819 rigid_add_edge_to_R(sys, efa->v1, efa->v2, sys->fweights[a][2]);
820 rigid_add_edge_to_R(sys, efa->v2, efa->v3, sys->fweights[a][0]);
821 rigid_add_edge_to_R(sys, efa->v3, efa->v1, sys->fweights[a][1]);
825 rigid_add_edge_to_R(sys, efa->v1, efa->v3, sys->fweights[a][2]);
826 rigid_add_edge_to_R(sys, efa->v3, efa->v4, sys->fweights[a][0]);
827 rigid_add_edge_to_R(sys, efa->v4, efa->v1, sys->fweights[a][1]);
831 for(a=0, eve=em->verts.first; eve; eve=eve->next, a++) {
832 rigid_orthogonalize_R(sys->rigid.R[a]);
836 /* compute right hand sides for solving */
837 for(a=0, efa=em->faces.first; efa; efa=efa->next, a++) {
838 rigid_add_edge_to_rhs(sys, efa->v1, efa->v2, sys->fweights[a][2]);
839 rigid_add_edge_to_rhs(sys, efa->v2, efa->v3, sys->fweights[a][0]);
840 rigid_add_edge_to_rhs(sys, efa->v3, efa->v1, sys->fweights[a][1]);
844 rigid_add_edge_to_rhs(sys, efa->v1, efa->v3, sys->fweights[a][2]);
845 rigid_add_edge_to_rhs(sys, efa->v3, efa->v4, sys->fweights[a][0]);
846 rigid_add_edge_to_rhs(sys, efa->v4, efa->v1, sys->fweights[a][1]);
850 /* solve for positions, for X,Y and Z separately */
852 laplacian_begin_solve(sys, i);
854 for(a=0; a<sys->totvert; a++)
856 laplacian_add_right_hand_side(sys, a, sys->rigid.rhs[a][i]);
858 if(laplacian_system_solve(sys)) {
859 for(a=0, eve=em->verts.first; eve; eve=eve->next, a++)
860 eve->co[i]= laplacian_system_get_solution(a);
863 if(!sys->rigid.thrownerror) {
864 error("RigidDeform: failed to find solution.");
865 sys->rigid.thrownerror= 1;
872 static void rigid_laplacian_create(LaplacianSystem *sys)
874 EditMesh *em = sys->rigid.mesh;
879 /* add verts and faces to laplacian */
880 for(a=0, eve=em->verts.first; eve; eve=eve->next, a++) {
881 laplacian_add_vertex(sys, eve->co, eve->pinned);
885 for(efa=em->faces.first; efa; efa=efa->next) {
886 laplacian_add_triangle(sys,
887 efa->v1->tmp.l, efa->v2->tmp.l, efa->v3->tmp.l);
889 laplacian_add_triangle(sys,
890 efa->v1->tmp.l, efa->v3->tmp.l, efa->v4->tmp.l);
894 void rigid_deform_begin(EditMesh *em)
896 LaplacianSystem *sys;
899 int a, totvert, totface;
901 /* count vertices, triangles */
902 for(totvert=0, eve=em->verts.first; eve; eve=eve->next)
905 for(totface=0, efa=em->faces.first; efa; efa=efa->next) {
907 if(efa->v4) totface++;
910 /* create laplacian */
911 sys = laplacian_system_construct_begin(totvert, totface, 0);
914 sys->rigid.R = MEM_callocN(sizeof(float)*3*3*totvert, "RigidDeformR");
915 sys->rigid.rhs = MEM_callocN(sizeof(float)*3*totvert, "RigidDeformRHS");
916 sys->rigid.origco = MEM_callocN(sizeof(float)*3*totvert, "RigidDeformCo");
918 for(a=0, eve=em->verts.first; eve; eve=eve->next, a++)
919 VecCopyf(sys->rigid.origco[a], eve->co);
922 sys->storeweights= 1;
924 rigid_laplacian_create(sys);
926 laplacian_system_construct_end(sys);
928 RigidDeformSystem = sys;
931 void rigid_deform_end(int cancel)
933 LaplacianSystem *sys = RigidDeformSystem;
936 EditMesh *em = sys->rigid.mesh;
941 for(a=0, eve=em->verts.first; eve; eve=eve->next, a++)
943 VecCopyf(eve->co, sys->rigid.origco[a]);
945 if(sys->rigid.R) MEM_freeN(sys->rigid.R);
946 if(sys->rigid.rhs) MEM_freeN(sys->rigid.rhs);
947 if(sys->rigid.origco) MEM_freeN(sys->rigid.origco);
950 laplacian_system_delete(sys);
953 RigidDeformSystem = NULL;
957 /************************** Harmonic Coordinates ****************************/
958 /* From "Harmonic Coordinates for Character Articulation",
959 Pushkar Joshi, Mark Meyer, Tony DeRose, Brian Green and Tom Sanocki,
962 #define EPSILON 0.0001f
964 #define MESHDEFORM_TAG_UNTYPED 0
965 #define MESHDEFORM_TAG_BOUNDARY 1
966 #define MESHDEFORM_TAG_INTERIOR 2
967 #define MESHDEFORM_TAG_EXTERIOR 3
969 #define MESHDEFORM_LEN_THRESHOLD 1e-6
971 #define MESHDEFORM_MIN_INFLUENCE 0.0005
973 static int MESHDEFORM_OFFSET[7][3] =
974 {{0,0,0}, {1,0,0}, {-1,0,0}, {0,1,0}, {0,-1,0}, {0,0,1}, {0,0,-1}};
976 typedef struct MDefBoundIsect {
978 int nvert, v[4], facing;
982 typedef struct MDefBindInfluence {
983 struct MDefBindInfluence *next;
988 typedef struct MeshDeformBind {
989 /* grid dimensions */
990 float min[3], max[3];
991 float width[3], halfwidth[3];
997 float (*vertexcos)[3];
998 int totvert, totcagevert;
1002 MDefBoundIsect *(*boundisect)[6];
1005 float *phi, *totalphi;
1010 MDefBindInfluence **dyngrid;
1011 float cagemat[4][4];
1020 /* ray intersection */
1022 /* our own triangle intersection, so we can fully control the epsilons and
1023 * prevent corner case from going wrong*/
1024 static int meshdeform_tri_intersect(float orig[3], float end[3], float vert0[3],
1025 float vert1[3], float vert2[3], float *isectco, float *uvw)
1027 float edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
1028 float det,inv_det, u, v, dir[3], isectdir[3];
1030 VECSUB(dir, end, orig);
1032 /* find vectors for two edges sharing vert0 */
1033 VECSUB(edge1, vert1, vert0);
1034 VECSUB(edge2, vert2, vert0);
1036 /* begin calculating determinant - also used to calculate U parameter */
1037 Crossf(pvec, dir, edge2);
1039 /* if determinant is near zero, ray lies in plane of triangle */
1040 det = INPR(edge1, pvec);
1044 inv_det = 1.0f / det;
1046 /* calculate distance from vert0 to ray origin */
1047 VECSUB(tvec, orig, vert0);
1049 /* calculate U parameter and test bounds */
1050 u = INPR(tvec, pvec) * inv_det;
1051 if (u < -EPSILON || u > 1.0f+EPSILON)
1054 /* prepare to test V parameter */
1055 Crossf(qvec, tvec, edge1);
1057 /* calculate V parameter and test bounds */
1058 v = INPR(dir, qvec) * inv_det;
1059 if (v < -EPSILON || u + v > 1.0f+EPSILON)
1062 isectco[0]= (1.0f - u - v)*vert0[0] + u*vert1[0] + v*vert2[0];
1063 isectco[1]= (1.0f - u - v)*vert0[1] + u*vert1[1] + v*vert2[1];
1064 isectco[2]= (1.0f - u - v)*vert0[2] + u*vert1[2] + v*vert2[2];
1066 uvw[0]= 1.0 - u - v;
1070 /* check if it is within the length of the line segment */
1071 VECSUB(isectdir, isectco, orig);
1073 if(INPR(dir, isectdir) < -EPSILON)
1076 if(INPR(dir, dir) + EPSILON < INPR(isectdir, isectdir))
1082 /* blender's raytracer is not use now, even though it is much faster. it can
1083 * give problems with rays falling through, so we use our own intersection
1084 * function above with tweaked epsilons */
1087 static MeshDeformBind *MESHDEFORM_BIND = NULL;
1089 static void meshdeform_ray_coords_func(RayFace *face, float **v1, float **v2, float **v3, float **v4)
1091 MFace *mface= (MFace*)face;
1092 float (*cagecos)[3]= MESHDEFORM_BIND->cagecos;
1094 *v1= cagecos[mface->v1];
1095 *v2= cagecos[mface->v2];
1096 *v3= cagecos[mface->v3];
1097 *v4= (mface->v4)? cagecos[mface->v4]: NULL;
1100 static int meshdeform_ray_check_func(Isect *is, RayFace *face)
1105 static void meshdeform_ray_tree_create(MeshDeformBind *mdb)
1108 float min[3], max[3];
1111 /* create a raytrace tree from the mesh */
1112 INIT_MINMAX(min, max);
1114 for(a=0; a<mdb->totcagevert; a++)
1115 DO_MINMAX(mdb->cagecos[a], min, max)
1117 MESHDEFORM_BIND= mdb;
1119 mface= mdb->cagedm->getFaceArray(mdb->cagedm);
1120 totface= mdb->cagedm->getNumFaces(mdb->cagedm);
1122 mdb->raytree= RE_ray_tree_create(64, totface, min, max,
1123 meshdeform_ray_coords_func, meshdeform_ray_check_func);
1125 for(a=0; a<totface; a++, mface++)
1126 RE_ray_tree_add_face(mdb->raytree, mface);
1128 RE_ray_tree_done(mdb->raytree);
1131 static void meshdeform_ray_tree_free(MeshDeformBind *mdb)
1133 MESHDEFORM_BIND= NULL;
1134 RE_ray_tree_free(mdb->raytree);
1138 static int meshdeform_intersect(MeshDeformBind *mdb, Isect *isec)
1141 float face[4][3], co[3], uvw[3], len, nor[3], end[3];
1142 int f, hit, is= 0, totface;
1146 mface= mdb->cagedm->getFaceArray(mdb->cagedm);
1147 totface= mdb->cagedm->getNumFaces(mdb->cagedm);
1149 VECADDFAC( end, isec->start, isec->vec, isec->labda );
1151 for(f=0; f<totface; f++, mface++) {
1152 VECCOPY(face[0], mdb->cagecos[mface->v1]);
1153 VECCOPY(face[1], mdb->cagecos[mface->v2]);
1154 VECCOPY(face[2], mdb->cagecos[mface->v3]);
1157 VECCOPY(face[3], mdb->cagecos[mface->v4]);
1158 hit = meshdeform_tri_intersect(isec->start, end, face[0], face[1], face[2], co, uvw);
1161 CalcNormFloat(face[0], face[1], face[2], nor);
1164 hit= meshdeform_tri_intersect(isec->start, end, face[0], face[2], face[3], co, uvw);
1165 CalcNormFloat(face[0], face[2], face[3], nor);
1169 hit= meshdeform_tri_intersect(isec->start, end, face[0], face[1], face[2], co, uvw);
1170 CalcNormFloat(face[0], face[1], face[2], nor);
1174 len= VecLenf(isec->start, co)/VecLenf(isec->start, end);
1175 if(len < isec->labda) {
1177 isec->hit.face = mface;
1178 isec->isect= (INPR(isec->vec, nor) <= 0.0f);
1187 static MDefBoundIsect *meshdeform_ray_tree_intersect(MeshDeformBind *mdb, float *co1, float *co2)
1189 MDefBoundIsect *isect;
1191 float (*cagecos)[3];
1193 float vert[4][3], len, end[3];
1194 static float epsilon[3]= {0, 0, 0}; //1e-4, 1e-4, 1e-4};
1197 memset(&isec, 0, sizeof(isec));
1198 isec.mode= RE_RAY_MIRROR; /* we want the closest intersection */
1202 VECADD(isec.start, co1, epsilon);
1203 VECADD(end, co2, epsilon);
1204 VECSUB(isec.vec, end, isec.start);
1207 /*if(RE_ray_tree_intersect(mdb->raytree, &isec)) {*/
1210 if(meshdeform_intersect(mdb, &isec)) {
1212 mface=(MFace*)isec.hit.face;
1214 /* create MDefBoundIsect */
1215 isect= BLI_memarena_alloc(mdb->memarena, sizeof(*isect));
1217 /* compute intersection coordinate */
1218 isect->co[0]= co1[0] + isec.vec[0]*len;
1219 isect->co[1]= co1[1] + isec.vec[1]*len;
1220 isect->co[2]= co1[2] + isec.vec[2]*len;
1222 isect->len= VecLenf(co1, isect->co);
1223 if(isect->len < MESHDEFORM_LEN_THRESHOLD)
1224 isect->len= MESHDEFORM_LEN_THRESHOLD;
1226 isect->v[0]= mface->v1;
1227 isect->v[1]= mface->v2;
1228 isect->v[2]= mface->v3;
1229 isect->v[3]= mface->v4;
1230 isect->nvert= (mface->v4)? 4: 3;
1232 isect->facing= isec.isect;
1234 /* compute mean value coordinates for interpolation */
1235 cagecos= mdb->cagecos;
1236 VECCOPY(vert[0], cagecos[mface->v1]);
1237 VECCOPY(vert[1], cagecos[mface->v2]);
1238 VECCOPY(vert[2], cagecos[mface->v3]);
1239 if(mface->v4) VECCOPY(vert[3], cagecos[mface->v4]);
1240 MeanValueWeights(vert, isect->nvert, isect->co, isect->uvw);
1248 static int meshdeform_inside_cage(MeshDeformBind *mdb, float *co)
1250 MDefBoundIsect *isect;
1251 float outside[3], start[3], dir[3];
1254 for(i=1; i<=6; i++) {
1257 outside[0] = co[0] + (mdb->max[0] - mdb->min[0] + 1.0f)*MESHDEFORM_OFFSET[i][0];
1258 outside[1] = co[1] + (mdb->max[1] - mdb->min[1] + 1.0f)*MESHDEFORM_OFFSET[i][1];
1259 outside[2] = co[2] + (mdb->max[2] - mdb->min[2] + 1.0f)*MESHDEFORM_OFFSET[i][2];
1262 VECSUB(dir, outside, start);
1265 isect = meshdeform_ray_tree_intersect(mdb, start, outside);
1266 if(isect && !isect->facing)
1275 static int meshdeform_index(MeshDeformBind *mdb, int x, int y, int z, int n)
1277 int size= mdb->size;
1279 x += MESHDEFORM_OFFSET[n][0];
1280 y += MESHDEFORM_OFFSET[n][1];
1281 z += MESHDEFORM_OFFSET[n][2];
1283 if(x < 0 || x >= mdb->size)
1285 if(y < 0 || y >= mdb->size)
1287 if(z < 0 || z >= mdb->size)
1290 return x + y*size + z*size*size;
1293 static void meshdeform_cell_center(MeshDeformBind *mdb, int x, int y, int z, int n, float *center)
1295 x += MESHDEFORM_OFFSET[n][0];
1296 y += MESHDEFORM_OFFSET[n][1];
1297 z += MESHDEFORM_OFFSET[n][2];
1299 center[0]= mdb->min[0] + x*mdb->width[0] + mdb->halfwidth[0];
1300 center[1]= mdb->min[1] + y*mdb->width[1] + mdb->halfwidth[1];
1301 center[2]= mdb->min[2] + z*mdb->width[2] + mdb->halfwidth[2];
1304 static void meshdeform_add_intersections(MeshDeformBind *mdb, int x, int y, int z)
1306 MDefBoundIsect *isect;
1307 float center[3], ncenter[3];
1310 a= meshdeform_index(mdb, x, y, z, 0);
1311 meshdeform_cell_center(mdb, x, y, z, 0, center);
1313 /* check each outgoing edge for intersection */
1314 for(i=1; i<=6; i++) {
1315 if(meshdeform_index(mdb, x, y, z, i) == -1)
1318 meshdeform_cell_center(mdb, x, y, z, i, ncenter);
1320 isect= meshdeform_ray_tree_intersect(mdb, center, ncenter);
1322 mdb->boundisect[a][i-1]= isect;
1323 mdb->tag[a]= MESHDEFORM_TAG_BOUNDARY;
1328 static void meshdeform_bind_floodfill(MeshDeformBind *mdb)
1330 int *stack, *tag= mdb->tag;
1331 int a, b, i, xyz[3], stacksize, size= mdb->size;
1333 stack= MEM_callocN(sizeof(int)*mdb->size3, "MeshDeformBindStack");
1335 /* we know lower left corner is EXTERIOR because of padding */
1336 tag[0]= MESHDEFORM_TAG_EXTERIOR;
1340 /* floodfill exterior tag */
1341 while(stacksize > 0) {
1342 a= stack[--stacksize];
1344 xyz[2]= a/(size*size);
1345 xyz[1]= (a - xyz[2]*size*size)/size;
1346 xyz[0]= a - xyz[1]*size - xyz[2]*size*size;
1348 for(i=1; i<=6; i++) {
1349 b= meshdeform_index(mdb, xyz[0], xyz[1], xyz[2], i);
1352 if(tag[b] == MESHDEFORM_TAG_UNTYPED ||
1353 (tag[b] == MESHDEFORM_TAG_BOUNDARY && !mdb->boundisect[a][i-1])) {
1354 tag[b]= MESHDEFORM_TAG_EXTERIOR;
1355 stack[stacksize++]= b;
1361 /* other cells are interior */
1362 for(a=0; a<size*size*size; a++)
1363 if(tag[a]==MESHDEFORM_TAG_UNTYPED)
1364 tag[a]= MESHDEFORM_TAG_INTERIOR;
1370 for(a=0; a<size*size*size; a++)
1371 if(tag[a]==MESHDEFORM_TAG_BOUNDARY)
1373 else if(tag[a]==MESHDEFORM_TAG_INTERIOR)
1375 else if(tag[a]==MESHDEFORM_TAG_EXTERIOR) {
1378 if(mdb->semibound[a])
1382 printf("interior %d exterior %d boundary %d semi-boundary %d\n", ti, te, tb, ts);
1389 static float meshdeform_boundary_phi(MeshDeformBind *mdb, MDefBoundIsect *isect, int cagevert)
1393 for(a=0; a<isect->nvert; a++)
1394 if(isect->v[a] == cagevert)
1395 return isect->uvw[a];
1400 static float meshdeform_interp_w(MeshDeformBind *mdb, float *gridvec, float *vec, int cagevert)
1402 float dvec[3], ivec[3], wx, wy, wz, result=0.0f;
1403 float weight, totweight= 0.0f;
1406 for(i=0; i<3; i++) {
1407 ivec[i]= (int)gridvec[i];
1408 dvec[i]= gridvec[i] - ivec[i];
1411 for(i=0; i<8; i++) {
1412 if(i & 1) { x= ivec[0]+1; wx= dvec[0]; }
1413 else { x= ivec[0]; wx= 1.0f-dvec[0]; }
1415 if(i & 2) { y= ivec[1]+1; wy= dvec[1]; }
1416 else { y= ivec[1]; wy= 1.0f-dvec[1]; }
1418 if(i & 4) { z= ivec[2]+1; wz= dvec[2]; }
1419 else { z= ivec[2]; wz= 1.0f-dvec[2]; }
1421 CLAMP(x, 0, mdb->size-1);
1422 CLAMP(y, 0, mdb->size-1);
1423 CLAMP(z, 0, mdb->size-1);
1425 a= meshdeform_index(mdb, x, y, z, 0);
1427 result += weight*mdb->phi[a];
1428 totweight += weight;
1431 if(totweight > 0.0f)
1432 result /= totweight;
1437 static void meshdeform_check_semibound(MeshDeformBind *mdb, int x, int y, int z)
1441 a= meshdeform_index(mdb, x, y, z, 0);
1442 if(mdb->tag[a] != MESHDEFORM_TAG_EXTERIOR)
1446 if(mdb->boundisect[a][i-1])
1447 mdb->semibound[a]= 1;
1450 static float meshdeform_boundary_total_weight(MeshDeformBind *mdb, int x, int y, int z)
1452 float weight, totweight= 0.0f;
1455 a= meshdeform_index(mdb, x, y, z, 0);
1457 /* count weight for neighbour cells */
1458 for(i=1; i<=6; i++) {
1459 if(meshdeform_index(mdb, x, y, z, i) == -1)
1462 if(mdb->boundisect[a][i-1])
1463 weight= 1.0f/mdb->boundisect[a][i-1]->len;
1464 else if(!mdb->semibound[a])
1465 weight= 1.0f/mdb->width[0];
1469 totweight += weight;
1475 static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, int x, int y, int z)
1477 MDefBoundIsect *isect;
1478 float weight, totweight;
1481 acenter= meshdeform_index(mdb, x, y, z, 0);
1482 if(mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR)
1485 nlMatrixAdd(mdb->varidx[acenter], mdb->varidx[acenter], 1.0f);
1487 totweight= meshdeform_boundary_total_weight(mdb, x, y, z);
1488 for(i=1; i<=6; i++) {
1489 a= meshdeform_index(mdb, x, y, z, i);
1490 if(a == -1 || mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR)
1493 isect= mdb->boundisect[acenter][i-1];
1495 weight= (1.0f/mdb->width[0])/totweight;
1496 nlMatrixAdd(mdb->varidx[acenter], mdb->varidx[a], -weight);
1501 static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, int x, int y, int z, int cagevert)
1503 MDefBoundIsect *isect;
1504 float rhs, weight, totweight;
1507 acenter= meshdeform_index(mdb, x, y, z, 0);
1508 if(mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR)
1511 totweight= meshdeform_boundary_total_weight(mdb, x, y, z);
1512 for(i=1; i<=6; i++) {
1513 a= meshdeform_index(mdb, x, y, z, i);
1517 isect= mdb->boundisect[acenter][i-1];
1520 weight= (1.0f/isect->len)/totweight;
1521 rhs= weight*meshdeform_boundary_phi(mdb, isect, cagevert);
1522 nlRightHandSideAdd(0, mdb->varidx[acenter], rhs);
1527 static void meshdeform_matrix_add_semibound_phi(MeshDeformBind *mdb, int x, int y, int z, int cagevert)
1529 MDefBoundIsect *isect;
1530 float rhs, weight, totweight;
1533 a= meshdeform_index(mdb, x, y, z, 0);
1534 if(!mdb->semibound[a])
1539 totweight= meshdeform_boundary_total_weight(mdb, x, y, z);
1540 for(i=1; i<=6; i++) {
1541 isect= mdb->boundisect[a][i-1];
1544 weight= (1.0f/isect->len)/totweight;
1545 rhs= weight*meshdeform_boundary_phi(mdb, isect, cagevert);
1551 static void meshdeform_matrix_add_exterior_phi(MeshDeformBind *mdb, int x, int y, int z, int cagevert)
1553 float phi, totweight;
1556 acenter= meshdeform_index(mdb, x, y, z, 0);
1557 if(mdb->tag[acenter] != MESHDEFORM_TAG_EXTERIOR || mdb->semibound[acenter])
1562 for(i=1; i<=6; i++) {
1563 a= meshdeform_index(mdb, x, y, z, i);
1565 if(a != -1 && mdb->semibound[a]) {
1571 if(totweight != 0.0f)
1572 mdb->phi[acenter]= phi/totweight;
1575 static void meshdeform_matrix_solve(MeshDeformBind *mdb)
1578 float vec[3], gridvec[3];
1579 int a, b, x, y, z, totvar;
1582 /* setup variable indices */
1583 mdb->varidx= MEM_callocN(sizeof(int)*mdb->size3, "MeshDeformDSvaridx");
1584 for(a=0, totvar=0; a<mdb->size3; a++)
1585 mdb->varidx[a]= (mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR)? -1: totvar++;
1588 MEM_freeN(mdb->varidx);
1592 progress_bar(0, "Starting mesh deform solve");
1594 /* setup opennl solver */
1596 context= nlGetCurrent();
1598 nlSolverParameteri(NL_NB_VARIABLES, totvar);
1599 nlSolverParameteri(NL_NB_ROWS, totvar);
1600 nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 1);
1606 for(z=0; z<mdb->size; z++)
1607 for(y=0; y<mdb->size; y++)
1608 for(x=0; x<mdb->size; x++)
1609 meshdeform_matrix_add_cell(mdb, x, y, z);
1611 /* solve for each cage vert */
1612 for(a=0; a<mdb->totcagevert; a++) {
1618 /* fill in right hand side and solve */
1619 for(z=0; z<mdb->size; z++)
1620 for(y=0; y<mdb->size; y++)
1621 for(x=0; x<mdb->size; x++)
1622 meshdeform_matrix_add_rhs(mdb, x, y, z, a);
1631 if(nlSolveAdvanced(NULL, NL_TRUE)) {
1632 for(z=0; z<mdb->size; z++)
1633 for(y=0; y<mdb->size; y++)
1634 for(x=0; x<mdb->size; x++)
1635 meshdeform_matrix_add_semibound_phi(mdb, x, y, z, a);
1637 for(z=0; z<mdb->size; z++)
1638 for(y=0; y<mdb->size; y++)
1639 for(x=0; x<mdb->size; x++)
1640 meshdeform_matrix_add_exterior_phi(mdb, x, y, z, a);
1642 for(b=0; b<mdb->size3; b++) {
1643 if(mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR)
1644 mdb->phi[b]= nlGetVariable(0, mdb->varidx[b]);
1645 mdb->totalphi[b] += mdb->phi[b];
1649 /* static bind : compute weights for each vertex */
1650 for(b=0; b<mdb->totvert; b++) {
1651 if(mdb->inside[b]) {
1652 VECCOPY(vec, mdb->vertexcos[b]);
1653 Mat4MulVecfl(mdb->cagemat, vec);
1654 gridvec[0]= (vec[0] - mdb->min[0] - mdb->halfwidth[0])/mdb->width[0];
1655 gridvec[1]= (vec[1] - mdb->min[1] - mdb->halfwidth[1])/mdb->width[1];
1656 gridvec[2]= (vec[2] - mdb->min[2] - mdb->halfwidth[2])/mdb->width[2];
1658 mdb->weights[b*mdb->totcagevert + a]= meshdeform_interp_w(mdb, gridvec, vec, a);
1663 MDefBindInfluence *inf;
1666 for(b=0; b<mdb->size3; b++) {
1667 if(mdb->phi[b] >= MESHDEFORM_MIN_INFLUENCE) {
1668 inf= BLI_memarena_alloc(mdb->memarena, sizeof(*inf));
1670 inf->weight= mdb->phi[b];
1671 inf->next= mdb->dyngrid[b];
1672 mdb->dyngrid[b]= inf;
1678 error("Mesh Deform: failed to find solution.");
1682 sprintf(message, "Mesh deform solve %d / %d |||", a+1, mdb->totcagevert);
1683 progress_bar((float)(a+1)/(float)(mdb->totcagevert), message);
1688 for(b=0; b<mdb->size3; b++)
1689 if(mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR)
1690 if(fabs(mdb->totalphi[b] - 1.0f) > 1e-4)
1691 printf("totalphi deficiency [%s|%d] %d: %.10f\n",
1692 (mdb->tag[b] == MESHDEFORM_TAG_INTERIOR)? "interior": "boundary", mdb->semibound[b], mdb->varidx[b], mdb->totalphi[b]);
1696 MEM_freeN(mdb->varidx);
1698 nlDeleteContext(context);
1701 void harmonic_coordinates_bind(Scene *scene, MeshDeformModifierData *mmd, float (*vertexcos)[3], int totvert, float cagemat[][4])
1704 MDefBindInfluence *inf;
1705 MDefInfluence *mdinf;
1708 float center[3], vec[3], maxwidth, totweight;
1709 int a, b, x, y, z, totinside, offset;
1712 start_progress_bar();
1714 memset(&mdb, 0, sizeof(MeshDeformBind));
1716 /* get mesh and cage mesh */
1717 mdb.vertexcos= vertexcos;
1718 mdb.totvert= totvert;
1720 mdb.cagedm= mesh_create_derived_no_deform(scene, mmd->object, NULL, CD_MASK_BAREMESH);
1721 mdb.totcagevert= mdb.cagedm->getNumVerts(mdb.cagedm);
1722 mdb.cagecos= MEM_callocN(sizeof(*mdb.cagecos)*mdb.totcagevert, "MeshDeformBindCos");
1723 Mat4CpyMat4(mdb.cagemat, cagemat);
1725 mvert= mdb.cagedm->getVertArray(mdb.cagedm);
1726 for(a=0; a<mdb.totcagevert; a++)
1727 VECCOPY(mdb.cagecos[a], mvert[a].co)
1729 /* compute bounding box of the cage mesh */
1730 INIT_MINMAX(mdb.min, mdb.max);
1732 for(a=0; a<mdb.totcagevert; a++)
1733 DO_MINMAX(mdb.cagecos[a], mdb.min, mdb.max);
1735 /* allocate memory */
1736 mdb.size= (2<<(mmd->gridsize-1)) + 2;
1737 mdb.size3= mdb.size*mdb.size*mdb.size;
1738 mdb.tag= MEM_callocN(sizeof(int)*mdb.size3, "MeshDeformBindTag");
1739 mdb.phi= MEM_callocN(sizeof(float)*mdb.size3, "MeshDeformBindPhi");
1740 mdb.totalphi= MEM_callocN(sizeof(float)*mdb.size3, "MeshDeformBindTotalPhi");
1741 mdb.boundisect= MEM_callocN(sizeof(*mdb.boundisect)*mdb.size3, "MDefBoundIsect");
1742 mdb.semibound= MEM_callocN(sizeof(int)*mdb.size3, "MDefSemiBound");
1744 mdb.inside= MEM_callocN(sizeof(int)*mdb.totvert, "MDefInside");
1746 if(mmd->flag & MOD_MDEF_DYNAMIC_BIND)
1747 mdb.dyngrid= MEM_callocN(sizeof(MDefBindInfluence*)*mdb.size3, "MDefDynGrid");
1749 mdb.weights= MEM_callocN(sizeof(float)*mdb.totvert*mdb.totcagevert, "MDefWeights");
1751 mdb.memarena= BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE);
1752 BLI_memarena_use_calloc(mdb.memarena);
1754 /* make bounding box equal size in all directions, add padding, and compute
1755 * width of the cells */
1758 if(mdb.max[a]-mdb.min[a] > maxwidth)
1759 maxwidth= mdb.max[a]-mdb.min[a];
1761 for(a=0; a<3; a++) {
1762 center[a]= (mdb.min[a]+mdb.max[a])*0.5f;
1763 mdb.min[a]= center[a] - maxwidth*0.5f;
1764 mdb.max[a]= center[a] + maxwidth*0.5f;
1766 mdb.width[a]= (mdb.max[a]-mdb.min[a])/(mdb.size-4);
1767 mdb.min[a] -= 2.1f*mdb.width[a];
1768 mdb.max[a] += 2.1f*mdb.width[a];
1770 mdb.width[a]= (mdb.max[a]-mdb.min[a])/mdb.size;
1771 mdb.halfwidth[a]= mdb.width[a]*0.5f;
1774 progress_bar(0, "Setting up mesh deform system");
1777 /* create ray tree */
1778 meshdeform_ray_tree_create(&mdb);
1782 for(a=0; a<mdb.totvert; a++) {
1783 VECCOPY(vec, mdb.vertexcos[a]);
1784 Mat4MulVecfl(mdb.cagemat, vec);
1785 mdb.inside[a]= meshdeform_inside_cage(&mdb, vec);
1790 /* free temporary MDefBoundIsects */
1791 BLI_memarena_free(mdb.memarena);
1792 mdb.memarena= BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE);
1794 /* start with all cells untyped */
1795 for(a=0; a<mdb.size3; a++)
1796 mdb.tag[a]= MESHDEFORM_TAG_UNTYPED;
1798 /* detect intersections and tag boundary cells */
1799 for(z=0; z<mdb.size; z++)
1800 for(y=0; y<mdb.size; y++)
1801 for(x=0; x<mdb.size; x++)
1802 meshdeform_add_intersections(&mdb, x, y, z);
1806 meshdeform_ray_tree_free(&mdb);
1809 /* compute exterior and interior tags */
1810 meshdeform_bind_floodfill(&mdb);
1812 for(z=0; z<mdb.size; z++)
1813 for(y=0; y<mdb.size; y++)
1814 for(x=0; x<mdb.size; x++)
1815 meshdeform_check_semibound(&mdb, x, y, z);
1818 meshdeform_matrix_solve(&mdb);
1820 /* assign results */
1821 mmd->bindcos= (float*)mdb.cagecos;
1822 mmd->totvert= mdb.totvert;
1823 mmd->totcagevert= mdb.totcagevert;
1824 Mat4CpyMat4(mmd->bindmat, mmd->object->obmat);
1826 if(mmd->flag & MOD_MDEF_DYNAMIC_BIND) {
1827 mmd->totinfluence= 0;
1828 for(a=0; a<mdb.size3; a++)
1829 for(inf=mdb.dyngrid[a]; inf; inf=inf->next)
1830 mmd->totinfluence++;
1832 /* convert MDefBindInfluences to smaller MDefInfluences */
1833 mmd->dyngrid= MEM_callocN(sizeof(MDefCell)*mdb.size3, "MDefDynGrid");
1834 mmd->dyninfluences= MEM_callocN(sizeof(MDefInfluence)*mmd->totinfluence, "MDefInfluence");
1836 for(a=0; a<mdb.size3; a++) {
1837 cell= &mmd->dyngrid[a];
1838 cell->offset= offset;
1841 mdinf= mmd->dyninfluences + cell->offset;
1842 for(inf=mdb.dyngrid[a]; inf; inf=inf->next, mdinf++) {
1843 mdinf->weight= inf->weight;
1844 mdinf->vertex= inf->vertex;
1845 totweight += mdinf->weight;
1846 cell->totinfluence++;
1849 if(totweight > 0.0f) {
1850 mdinf= mmd->dyninfluences + cell->offset;
1851 for(b=0; b<cell->totinfluence; b++, mdinf++)
1852 mdinf->weight /= totweight;
1855 offset += cell->totinfluence;
1858 mmd->dynverts= mdb.inside;
1859 mmd->dyngridsize= mdb.size;
1860 VECCOPY(mmd->dyncellmin, mdb.min);
1861 mmd->dyncellwidth= mdb.width[0];
1862 MEM_freeN(mdb.dyngrid);
1865 mmd->bindweights= mdb.weights;
1866 MEM_freeN(mdb.inside);
1869 /* transform bindcos to world space */
1870 for(a=0; a<mdb.totcagevert; a++)
1871 Mat4MulVecfl(mmd->object->obmat, mmd->bindcos+a*3);
1874 mdb.cagedm->release(mdb.cagedm);
1877 MEM_freeN(mdb.totalphi);
1878 MEM_freeN(mdb.boundisect);
1879 MEM_freeN(mdb.semibound);
1880 BLI_memarena_free(mdb.memarena);