diff options
Diffstat (limited to 'source/blender/editors/armature/meshlaplacian.c')
-rw-r--r-- | source/blender/editors/armature/meshlaplacian.c | 1920 |
1 files changed, 1920 insertions, 0 deletions
diff --git a/source/blender/editors/armature/meshlaplacian.c b/source/blender/editors/armature/meshlaplacian.c new file mode 100644 index 00000000000..8807b21e653 --- /dev/null +++ b/source/blender/editors/armature/meshlaplacian.c @@ -0,0 +1,1920 @@ +/** + * $Id$ + * + * ***** BEGIN GPL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + * + * The Original Code is: all of this file. + * + * Contributor(s): none yet. + * + * ***** END GPL LICENSE BLOCK ***** + * meshlaplacian.c: Algorithms using the mesh laplacian. + */ + +#include <math.h> +#include <string.h> + +#include "MEM_guardedalloc.h" + +#include "DNA_listBase.h" +#include "DNA_object_types.h" +#include "DNA_mesh_types.h" +#include "DNA_meshdata_types.h" +#include "DNA_modifier_types.h" +#include "DNA_scene_types.h" + +#include "BLI_arithb.h" +#include "BLI_edgehash.h" +#include "BLI_memarena.h" + +#include "BKE_DerivedMesh.h" +#include "BKE_utildefines.h" + +#ifdef RIGID_DEFORM +#include "BLI_editVert.h" +#include "BLI_polardecomp.h" +#endif + +#include "RE_raytrace.h" + +#include "ONL_opennl.h" + +#include "BLO_sys_types.h" // for intptr_t support + +#include "ED_armature.h" +#include "ED_mesh.h" + +#include "meshlaplacian.h" + + +/* ************* XXX *************** */ +static void waitcursor() {} +static void progress_bar() {} +static void start_progress_bar() {} +static void end_progress_bar() {} +static void error() {} +/* ************* XXX *************** */ + + +/************************** Laplacian System *****************************/ + +struct LaplacianSystem { + NLContext context; /* opennl context */ + + int totvert, totface; + + float **verts; /* vertex coordinates */ + float *varea; /* vertex weights for laplacian computation */ + char *vpinned; /* vertex pinning */ + int (*faces)[3]; /* face vertex indices */ + float (*fweights)[3]; /* cotangent weights per face */ + + int areaweights; /* use area in cotangent weights? */ + int storeweights; /* store cotangent weights in fweights */ + int nlbegun; /* nlBegin(NL_SYSTEM/NL_MATRIX) done */ + + EdgeHash *edgehash; /* edge hash for construction */ + + struct HeatWeighting { + Mesh *mesh; + float (*verts)[3]; /* vertex coordinates */ + float (*vnors)[3]; /* vertex normals */ + + float (*root)[3]; /* bone root */ + float (*tip)[3]; /* bone tip */ + int numbones; + + float *H; /* diagonal H matrix */ + float *p; /* values from all p vectors */ + float *mindist; /* minimum distance to a bone for all vertices */ + + RayTree *raytree; /* ray tracing acceleration structure */ + MFace **vface; /* a face that the vertex belongs to */ + } heat; + +#ifdef RIGID_DEFORM + struct RigidDeformation { + EditMesh *mesh; + + float (*R)[3][3]; + float (*rhs)[3]; + float (*origco)[3]; + int thrownerror; + } rigid; +#endif +}; + +/* Laplacian matrix construction */ + +/* Computation of these weights for the laplacian is based on: + "Discrete Differential-Geometry Operators for Triangulated 2-Manifolds", + Meyer et al, 2002. Section 3.5, formula (8). + + We do it a bit different by going over faces instead of going over each + vertex and adjacent faces, since we don't store this adjacency. Also, the + formulas are tweaked a bit to work for non-manifold meshes. */ + +static void laplacian_increase_edge_count(EdgeHash *edgehash, int v1, int v2) +{ + void **p = BLI_edgehash_lookup_p(edgehash, v1, v2); + + if(p) + *p = (void*)((intptr_t)*p + (intptr_t)1); + else + BLI_edgehash_insert(edgehash, v1, v2, (void*)(intptr_t)1); +} + +static int laplacian_edge_count(EdgeHash *edgehash, int v1, int v2) +{ + return (int)(intptr_t)BLI_edgehash_lookup(edgehash, v1, v2); +} + +static float cotan_weight(float *v1, float *v2, float *v3) +{ + float a[3], b[3], c[3], clen; + + VecSubf(a, v2, v1); + VecSubf(b, v3, v1); + Crossf(c, a, b); + + clen = VecLength(c); + + if (clen == 0.0f) + return 0.0f; + + return Inpf(a, b)/clen; +} + +static void laplacian_triangle_area(LaplacianSystem *sys, int i1, int i2, int i3) +{ + float t1, t2, t3, len1, len2, len3, area; + float *varea= sys->varea, *v1, *v2, *v3; + int obtuse = 0; + + v1= sys->verts[i1]; + v2= sys->verts[i2]; + v3= sys->verts[i3]; + + t1= cotan_weight(v1, v2, v3); + t2= cotan_weight(v2, v3, v1); + t3= cotan_weight(v3, v1, v2); + + if(VecAngle3(v2, v1, v3) > 90) obtuse= 1; + else if(VecAngle3(v1, v2, v3) > 90) obtuse= 2; + else if(VecAngle3(v1, v3, v2) > 90) obtuse= 3; + + if (obtuse > 0) { + area= AreaT3Dfl(v1, v2, v3); + + varea[i1] += (obtuse == 1)? area: area*0.5; + varea[i2] += (obtuse == 2)? area: area*0.5; + varea[i3] += (obtuse == 3)? area: area*0.5; + } + else { + len1= VecLenf(v2, v3); + len2= VecLenf(v1, v3); + len3= VecLenf(v1, v2); + + t1 *= len1*len1; + t2 *= len2*len2; + t3 *= len3*len3; + + varea[i1] += (t2 + t3)*0.25f; + varea[i2] += (t1 + t3)*0.25f; + varea[i3] += (t1 + t2)*0.25f; + } +} + +static void laplacian_triangle_weights(LaplacianSystem *sys, int f, int i1, int i2, int i3) +{ + float t1, t2, t3; + float *varea= sys->varea, *v1, *v2, *v3; + + v1= sys->verts[i1]; + v2= sys->verts[i2]; + v3= sys->verts[i3]; + + /* instead of *0.5 we divided by the number of faces of the edge, it still + needs to be verified that this is indeed the correct thing to do! */ + t1= cotan_weight(v1, v2, v3)/laplacian_edge_count(sys->edgehash, i2, i3); + t2= cotan_weight(v2, v3, v1)/laplacian_edge_count(sys->edgehash, i3, i1); + t3= cotan_weight(v3, v1, v2)/laplacian_edge_count(sys->edgehash, i1, i2); + + nlMatrixAdd(i1, i1, (t2+t3)*varea[i1]); + nlMatrixAdd(i2, i2, (t1+t3)*varea[i2]); + nlMatrixAdd(i3, i3, (t1+t2)*varea[i3]); + + nlMatrixAdd(i1, i2, -t3*varea[i1]); + nlMatrixAdd(i2, i1, -t3*varea[i2]); + + nlMatrixAdd(i2, i3, -t1*varea[i2]); + nlMatrixAdd(i3, i2, -t1*varea[i3]); + + nlMatrixAdd(i3, i1, -t2*varea[i3]); + nlMatrixAdd(i1, i3, -t2*varea[i1]); + + if(sys->storeweights) { + sys->fweights[f][0]= t1*varea[i1]; + sys->fweights[f][1]= t2*varea[i2]; + sys->fweights[f][2]= t3*varea[i3]; + } +} + +LaplacianSystem *laplacian_system_construct_begin(int totvert, int totface, int lsq) +{ + LaplacianSystem *sys; + + sys= MEM_callocN(sizeof(LaplacianSystem), "LaplacianSystem"); + + sys->verts= MEM_callocN(sizeof(float*)*totvert, "LaplacianSystemVerts"); + sys->vpinned= MEM_callocN(sizeof(char)*totvert, "LaplacianSystemVpinned"); + sys->faces= MEM_callocN(sizeof(int)*3*totface, "LaplacianSystemFaces"); + + sys->totvert= 0; + sys->totface= 0; + + sys->areaweights= 1; + sys->storeweights= 0; + + /* create opennl context */ + nlNewContext(); + nlSolverParameteri(NL_NB_VARIABLES, totvert); + if(lsq) + nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE); + + sys->context= nlGetCurrent(); + + return sys; +} + +void laplacian_add_vertex(LaplacianSystem *sys, float *co, int pinned) +{ + sys->verts[sys->totvert]= co; + sys->vpinned[sys->totvert]= pinned; + sys->totvert++; +} + +void laplacian_add_triangle(LaplacianSystem *sys, int v1, int v2, int v3) +{ + sys->faces[sys->totface][0]= v1; + sys->faces[sys->totface][1]= v2; + sys->faces[sys->totface][2]= v3; + sys->totface++; +} + +void laplacian_system_construct_end(LaplacianSystem *sys) +{ + int (*face)[3]; + int a, totvert=sys->totvert, totface=sys->totface; + + laplacian_begin_solve(sys, 0); + + sys->varea= MEM_callocN(sizeof(float)*totvert, "LaplacianSystemVarea"); + + sys->edgehash= BLI_edgehash_new(); + for(a=0, face=sys->faces; a<sys->totface; a++, face++) { + laplacian_increase_edge_count(sys->edgehash, (*face)[0], (*face)[1]); + laplacian_increase_edge_count(sys->edgehash, (*face)[1], (*face)[2]); + laplacian_increase_edge_count(sys->edgehash, (*face)[2], (*face)[0]); + } + + if(sys->areaweights) + for(a=0, face=sys->faces; a<sys->totface; a++, face++) + laplacian_triangle_area(sys, (*face)[0], (*face)[1], (*face)[2]); + + for(a=0; a<totvert; a++) { + if(sys->areaweights) { + if(sys->varea[a] != 0.0f) + sys->varea[a]= 0.5f/sys->varea[a]; + } + else + sys->varea[a]= 1.0f; + + /* for heat weighting */ + if(sys->heat.H) + nlMatrixAdd(a, a, sys->heat.H[a]); + } + + if(sys->storeweights) + sys->fweights= MEM_callocN(sizeof(float)*3*totface, "LaplacianFWeight"); + + for(a=0, face=sys->faces; a<totface; a++, face++) + laplacian_triangle_weights(sys, a, (*face)[0], (*face)[1], (*face)[2]); + + MEM_freeN(sys->faces); + sys->faces= NULL; + + if(sys->varea) { + MEM_freeN(sys->varea); + sys->varea= NULL; + } + + BLI_edgehash_free(sys->edgehash, NULL); + sys->edgehash= NULL; +} + +void laplacian_system_delete(LaplacianSystem *sys) +{ + if(sys->verts) MEM_freeN(sys->verts); + if(sys->varea) MEM_freeN(sys->varea); + if(sys->vpinned) MEM_freeN(sys->vpinned); + if(sys->faces) MEM_freeN(sys->faces); + if(sys->fweights) MEM_freeN(sys->fweights); + + nlDeleteContext(sys->context); + MEM_freeN(sys); +} + +void laplacian_begin_solve(LaplacianSystem *sys, int index) +{ + int a; + + if (!sys->nlbegun) { + nlBegin(NL_SYSTEM); + + if(index >= 0) { + for(a=0; a<sys->totvert; a++) { + if(sys->vpinned[a]) { + nlSetVariable(0, a, sys->verts[a][index]); + nlLockVariable(a); + } + } + } + + nlBegin(NL_MATRIX); + sys->nlbegun = 1; + } +} + +void laplacian_add_right_hand_side(LaplacianSystem *sys, int v, float value) +{ + nlRightHandSideAdd(0, v, value); +} + +int laplacian_system_solve(LaplacianSystem *sys) +{ + nlEnd(NL_MATRIX); + nlEnd(NL_SYSTEM); + sys->nlbegun = 0; + + //nlPrintMatrix(); + + return nlSolveAdvanced(NULL, NL_TRUE); +} + +float laplacian_system_get_solution(int v) +{ + return nlGetVariable(0, v); +} + +/************************* Heat Bone Weighting ******************************/ +/* From "Automatic Rigging and Animation of 3D Characters" + Ilya Baran and Jovan Popovic, SIGGRAPH 2007 */ + +#define C_WEIGHT 1.0f +#define WEIGHT_LIMIT_START 0.05f +#define WEIGHT_LIMIT_END 0.025f +#define DISTANCE_EPSILON 1e-4f + +/* Raytracing for vertex to bone visibility */ + +static LaplacianSystem *HeatSys = NULL; + +static void heat_ray_coords_func(RayFace *face, float **v1, float **v2, float **v3, float **v4) +{ + MFace *mface= (MFace*)face; + float (*verts)[3]= HeatSys->heat.verts; + + *v1= verts[mface->v1]; + *v2= verts[mface->v2]; + *v3= verts[mface->v3]; + *v4= (mface->v4)? verts[mface->v4]: NULL; +} + +static int heat_ray_check_func(Isect *is, int ob, RayFace *face) +{ + float *v1, *v2, *v3, *v4, nor[3]; + + /* don't intersect if the ray faces along the face normal */ + heat_ray_coords_func(face, &v1, &v2, &v3, &v4); + + if(v4) CalcNormFloat4(v1, v2, v3, v4, nor); + else CalcNormFloat(v1, v2, v3, nor); + + return (INPR(nor, is->vec) < 0); +} + +static void heat_ray_tree_create(LaplacianSystem *sys) +{ + Mesh *me = sys->heat.mesh; + RayTree *tree; + MFace *mface; + float min[3], max[3]; + int a; + + /* create a raytrace tree from the mesh */ + INIT_MINMAX(min, max); + + for(a=0; a<me->totvert; a++) + DO_MINMAX(sys->heat.verts[a], min, max); + + tree= RE_ray_tree_create(64, me->totface, min, max, + heat_ray_coords_func, heat_ray_check_func, NULL, NULL); + + sys->heat.vface= MEM_callocN(sizeof(MFace*)*me->totvert, "HeatVFaces"); + + HeatSys= sys; + + for(a=0, mface=me->mface; a<me->totface; a++, mface++) { + RE_ray_tree_add_face(tree, 0, mface); + + sys->heat.vface[mface->v1]= mface; + sys->heat.vface[mface->v2]= mface; + sys->heat.vface[mface->v3]= mface; + if(mface->v4) sys->heat.vface[mface->v4]= mface; + } + + HeatSys= NULL; + + RE_ray_tree_done(tree); + + sys->heat.raytree= tree; +} + +static int heat_ray_bone_visible(LaplacianSystem *sys, int vertex, int bone) +{ + Isect isec; + MFace *mface; + float dir[3]; + int visible; + + mface= sys->heat.vface[vertex]; + if(!mface) + return 1; + + /* setup isec */ + memset(&isec, 0, sizeof(isec)); + isec.mode= RE_RAY_SHADOW; + isec.lay= -1; + isec.face_last= NULL; + isec.faceorig= mface; + + VECCOPY(isec.start, sys->heat.verts[vertex]); + PclosestVL3Dfl(isec.end, isec.start, + sys->heat.root[bone], sys->heat.tip[bone]); + + /* add an extra offset to the start position to avoid self intersection */ + VECSUB(dir, isec.end, isec.start); + Normalize(dir); + VecMulf(dir, 1e-5); + VecAddf(isec.start, isec.start, dir); + + HeatSys= sys; + visible= !RE_ray_tree_intersect(sys->heat.raytree, &isec); + HeatSys= NULL; + + return visible; +} + +static float heat_bone_distance(LaplacianSystem *sys, int vertex, int bone) +{ + float closest[3], d[3], dist, cosine; + + /* compute euclidian distance */ + PclosestVL3Dfl(closest, sys->heat.verts[vertex], + sys->heat.root[bone], sys->heat.tip[bone]); + + VecSubf(d, sys->heat.verts[vertex], closest); + dist= Normalize(d); + + /* if the vertex normal does not point along the bone, increase distance */ + cosine= INPR(d, sys->heat.vnors[vertex]); + + return dist/(0.5f*(cosine + 1.001f)); +} + +static int heat_bone_closest(LaplacianSystem *sys, int vertex, int bone) +{ + float dist; + + dist= heat_bone_distance(sys, vertex, bone); + + if(dist <= sys->heat.mindist[vertex]*(1.0f + DISTANCE_EPSILON)) + if(heat_ray_bone_visible(sys, vertex, bone)) + return 1; + + return 0; +} + +static void heat_set_H(LaplacianSystem *sys, int vertex) +{ + float dist, mindist, h; + int j, numclosest = 0; + + mindist= 1e10; + + /* compute minimum distance */ + for(j=0; j<sys->heat.numbones; j++) { + dist= heat_bone_distance(sys, vertex, j); + + if(dist < mindist) + mindist= dist; + } + + sys->heat.mindist[vertex]= mindist; + + /* count number of bones with approximately this minimum distance */ + for(j=0; j<sys->heat.numbones; j++) + if(heat_bone_closest(sys, vertex, j)) + numclosest++; + + sys->heat.p[vertex]= (numclosest > 0)? 1.0f/numclosest: 0.0f; + + /* compute H entry */ + if(numclosest > 0) { + if(mindist > 1e-5) + h= numclosest*C_WEIGHT/(mindist*mindist); + else + h= 1e10f; + } + else + h= 0.0f; + + sys->heat.H[vertex]= h; +} + +void heat_calc_vnormals(LaplacianSystem *sys) +{ + float fnor[3]; + int a, v1, v2, v3, (*face)[3]; + + sys->heat.vnors= MEM_callocN(sizeof(float)*3*sys->totvert, "HeatVNors"); + + for(a=0, face=sys->faces; a<sys->totface; a++, face++) { + v1= (*face)[0]; + v2= (*face)[1]; + v3= (*face)[2]; + + CalcNormFloat(sys->verts[v1], sys->verts[v2], sys->verts[v3], fnor); + + VecAddf(sys->heat.vnors[v1], sys->heat.vnors[v1], fnor); + VecAddf(sys->heat.vnors[v2], sys->heat.vnors[v2], fnor); + VecAddf(sys->heat.vnors[v3], sys->heat.vnors[v3], fnor); + } + + for(a=0; a<sys->totvert; a++) + Normalize(sys->heat.vnors[a]); +} + +static void heat_laplacian_create(LaplacianSystem *sys) +{ + Mesh *me = sys->heat.mesh; + MFace *mface; + int a; + + /* heat specific definitions */ + sys->heat.mindist= MEM_callocN(sizeof(float)*me->totvert, "HeatMinDist"); + sys->heat.H= MEM_callocN(sizeof(float)*me->totvert, "HeatH"); + sys->heat.p= MEM_callocN(sizeof(float)*me->totvert, "HeatP"); + + /* add verts and faces to laplacian */ + for(a=0; a<me->totvert; a++) + laplacian_add_vertex(sys, sys->heat.verts[a], 0); + + for(a=0, mface=me->mface; a<me->totface; a++, mface++) { + laplacian_add_triangle(sys, mface->v1, mface->v2, mface->v3); + if(mface->v4) + laplacian_add_triangle(sys, mface->v1, mface->v3, mface->v4); + } + + /* for distance computation in set_H */ + heat_calc_vnormals(sys); + + for(a=0; a<me->totvert; a++) + heat_set_H(sys, a); +} + +static float heat_limit_weight(float weight) +{ + float t; + + if(weight < WEIGHT_LIMIT_END) { + return 0.0f; + } + else if(weight < WEIGHT_LIMIT_START) { + t= (weight - WEIGHT_LIMIT_END)/(WEIGHT_LIMIT_START - WEIGHT_LIMIT_END); + return t*WEIGHT_LIMIT_START; + } + else + return weight; +} + +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) +{ + LaplacianSystem *sys; + MFace *mface; + float solution, weight; + int *vertsflipped = NULL; + int a, totface, j, bbone, firstsegment, lastsegment, thrownerror = 0; + + /* count triangles */ + for(totface=0, a=0, mface=me->mface; a<me->totface; a++, mface++) { + totface++; + if(mface->v4) totface++; + } + + /* create laplacian */ + sys = laplacian_system_construct_begin(me->totvert, totface, 1); + + sys->heat.mesh= me; + sys->heat.verts= verts; + sys->heat.root= root; + sys->heat.tip= tip; + sys->heat.numbones= numbones; + + heat_ray_tree_create(sys); + heat_laplacian_create(sys); + + laplacian_system_construct_end(sys); + + if(dgroupflip) { + vertsflipped = MEM_callocN(sizeof(int)*me->totvert, "vertsflipped"); + for(a=0; a<me->totvert; a++) + vertsflipped[a] = mesh_get_x_mirror_vert(ob, a); + } + + /* compute weights per bone */ + for(j=0; j<numbones; j++) { + if(!selected[j]) + continue; + + firstsegment= (j == 0 || dgrouplist[j-1] != dgrouplist[j]); + lastsegment= (j == numbones-1 || dgrouplist[j] != dgrouplist[j+1]); + bbone= !(firstsegment && lastsegment); + + /* clear weights */ + if(bbone && firstsegment) { + for(a=0; a<me->totvert; a++) { + remove_vert_defgroup(ob, dgrouplist[j], a); + if(vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0) + remove_vert_defgroup(ob, dgroupflip[j], vertsflipped[a]); + } + } + + /* fill right hand side */ + laplacian_begin_solve(sys, -1); + + for(a=0; a<me->totvert; a++) + if(heat_bone_closest(sys, a, j)) + laplacian_add_right_hand_side(sys, a, + sys->heat.H[a]*sys->heat.p[a]); + + /* solve */ + if(laplacian_system_solve(sys)) { + /* load solution into vertex groups */ + for(a=0; a<me->totvert; a++) { + solution= laplacian_system_get_solution(a); + + if(bbone) { + if(solution > 0.0f) + add_vert_to_defgroup(ob, dgrouplist[j], a, solution, + WEIGHT_ADD); + } + else { + weight= heat_limit_weight(solution); + if(weight > 0.0f) + add_vert_to_defgroup(ob, dgrouplist[j], a, weight, + WEIGHT_REPLACE); + else + remove_vert_defgroup(ob, dgrouplist[j], a); + } + + /* do same for mirror */ + if(vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0) { + if(bbone) { + if(solution > 0.0f) + add_vert_to_defgroup(ob, dgroupflip[j], vertsflipped[a], + solution, WEIGHT_ADD); + } + else { + weight= heat_limit_weight(solution); + if(weight > 0.0f) + add_vert_to_defgroup(ob, dgroupflip[j], vertsflipped[a], + weight, WEIGHT_REPLACE); + else + remove_vert_defgroup(ob, dgroupflip[j], vertsflipped[a]); + } + } + } + } + else if(!thrownerror) { + error("Bone Heat Weighting:" + " failed to find solution for one or more bones"); + thrownerror= 1; + break; + } + + /* remove too small vertex weights */ + if(bbone && lastsegment) { + for(a=0; a<me->totvert; a++) { + weight= get_vert_defgroup(ob, dgrouplist[j], a); + weight= heat_limit_weight(weight); + if(weight <= 0.0f) + remove_vert_defgroup(ob, dgrouplist[j], a); + + if(vertsflipped && dgroupflip[j] && vertsflipped[a] >= 0) { + weight= get_vert_defgroup(ob, dgroupflip[j], vertsflipped[a]); + weight= heat_limit_weight(weight); + if(weight <= 0.0f) + remove_vert_defgroup(ob, dgroupflip[j], vertsflipped[a]); + } + } + } + } + + /* free */ + if(vertsflipped) MEM_freeN(vertsflipped); + + RE_ray_tree_free(sys->heat.raytree); + MEM_freeN(sys->heat.vface); + + MEM_freeN(sys->heat.mindist); + MEM_freeN(sys->heat.H); + MEM_freeN(sys->heat.p); + MEM_freeN(sys->heat.vnors); + + laplacian_system_delete(sys); +} + +#ifdef RIGID_DEFORM +/********************** As-Rigid-As-Possible Deformation ******************/ +/* From "As-Rigid-As-Possible Surface Modeling", + Olga Sorkine and Marc Alexa, ESGP 2007. */ + +/* investigate: + - transpose R in orthogonal + - flipped normals and per face adding + - move cancelling to transform, make origco pointer +*/ + +static LaplacianSystem *RigidDeformSystem = NULL; + +static void rigid_add_half_edge_to_R(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w) +{ + float e[3], e_[3]; + int i; + + VecSubf(e, sys->rigid.origco[v1->tmp.l], sys->rigid.origco[v2->tmp.l]); + VecSubf(e_, v1->co, v2->co); + + /* formula (5) */ + for (i=0; i<3; i++) { + sys->rigid.R[v1->tmp.l][i][0] += w*e[0]*e_[i]; + sys->rigid.R[v1->tmp.l][i][1] += w*e[1]*e_[i]; + sys->rigid.R[v1->tmp.l][i][2] += w*e[2]*e_[i]; + } +} + +static void rigid_add_edge_to_R(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w) +{ + rigid_add_half_edge_to_R(sys, v1, v2, w); + rigid_add_half_edge_to_R(sys, v2, v1, w); +} + +static void rigid_orthogonalize_R(float R[][3]) +{ + HMatrix M, Q, S; + + Mat4CpyMat3(M, R); + polar_decomp(M, Q, S); + Mat3CpyMat4(R, Q); +} + +static void rigid_add_half_edge_to_rhs(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w) +{ + /* formula (8) */ + float Rsum[3][3], rhs[3]; + + if (sys->vpinned[v1->tmp.l]) + return; + + Mat3AddMat3(Rsum, sys->rigid.R[v1->tmp.l], sys->rigid.R[v2->tmp.l]); + Mat3Transp(Rsum); + + VecSubf(rhs, sys->rigid.origco[v1->tmp.l], sys->rigid.origco[v2->tmp.l]); + Mat3MulVecfl(Rsum, rhs); + VecMulf(rhs, 0.5f); + VecMulf(rhs, w); + + VecAddf(sys->rigid.rhs[v1->tmp.l], sys->rigid.rhs[v1->tmp.l], rhs); +} + +static void rigid_add_edge_to_rhs(LaplacianSystem *sys, EditVert *v1, EditVert *v2, float w) +{ + rigid_add_half_edge_to_rhs(sys, v1, v2, w); + rigid_add_half_edge_to_rhs(sys, v2, v1, w); +} + +void rigid_deform_iteration() +{ + LaplacianSystem *sys= RigidDeformSystem; + EditMesh *em; + EditVert *eve; + EditFace *efa; + int a, i; + + if(!sys) + return; + + nlMakeCurrent(sys->context); + em= sys->rigid.mesh; + + /* compute R */ + memset(sys->rigid.R, 0, sizeof(float)*3*3*sys->totvert); + memset(sys->rigid.rhs, 0, sizeof(float)*3*sys->totvert); + + for(a=0, efa=em->faces.first; efa; efa=efa->next, a++) { + rigid_add_edge_to_R(sys, efa->v1, efa->v2, sys->fweights[a][2]); + rigid_add_edge_to_R(sys, efa->v2, efa->v3, sys->fweights[a][0]); + rigid_add_edge_to_R(sys, efa->v3, efa->v1, sys->fweights[a][1]); + + if(efa->v4) { + a++; + rigid_add_edge_to_R(sys, efa->v1, efa->v3, sys->fweights[a][2]); + rigid_add_edge_to_R(sys, efa->v3, efa->v4, sys->fweights[a][0]); + rigid_add_edge_to_R(sys, efa->v4, efa->v1, sys->fweights[a][1]); + } + } + + for(a=0, eve=em->verts.first; eve; eve=eve->next, a++) { + rigid_orthogonalize_R(sys->rigid.R[a]); + eve->tmp.l= a; + } + + /* compute right hand sides for solving */ + for(a=0, efa=em->faces.first; efa; efa=efa->next, a++) { + rigid_add_edge_to_rhs(sys, efa->v1, efa->v2, sys->fweights[a][2]); + rigid_add_edge_to_rhs(sys, efa->v2, efa->v3, sys->fweights[a][0]); + rigid_add_edge_to_rhs(sys, efa->v3, efa->v1, sys->fweights[a][1]); + + if(efa->v4) { + a++; + rigid_add_edge_to_rhs(sys, efa->v1, efa->v3, sys->fweights[a][2]); + rigid_add_edge_to_rhs(sys, efa->v3, efa->v4, sys->fweights[a][0]); + rigid_add_edge_to_rhs(sys, efa->v4, efa->v1, sys->fweights[a][1]); + } + } + + /* solve for positions, for X,Y and Z separately */ + for(i=0; i<3; i++) { + laplacian_begin_solve(sys, i); + + for(a=0; a<sys->totvert; a++) + if(!sys->vpinned[a]) + laplacian_add_right_hand_side(sys, a, sys->rigid.rhs[a][i]); + + if(laplacian_system_solve(sys)) { + for(a=0, eve=em->verts.first; eve; eve=eve->next, a++) + eve->co[i]= laplacian_system_get_solution(a); + } + else { + if(!sys->rigid.thrownerror) { + error("RigidDeform: failed to find solution."); + sys->rigid.thrownerror= 1; + } + break; + } + } +} + +static void rigid_laplacian_create(LaplacianSystem *sys) +{ + EditMesh *em = sys->rigid.mesh; + EditVert *eve; + EditFace *efa; + int a; + + /* add verts and faces to laplacian */ + for(a=0, eve=em->verts.first; eve; eve=eve->next, a++) { + laplacian_add_vertex(sys, eve->co, eve->pinned); + eve->tmp.l= a; + } + + for(efa=em->faces.first; efa; efa=efa->next) { + laplacian_add_triangle(sys, + efa->v1->tmp.l, efa->v2->tmp.l, efa->v3->tmp.l); + if(efa->v4) + laplacian_add_triangle(sys, + efa->v1->tmp.l, efa->v3->tmp.l, efa->v4->tmp.l); + } +} + +void rigid_deform_begin(EditMesh *em) +{ + LaplacianSystem *sys; + EditVert *eve; + EditFace *efa; + int a, totvert, totface; + + /* count vertices, triangles */ + for(totvert=0, eve=em->verts.first; eve; eve=eve->next) + totvert++; + + for(totface=0, efa=em->faces.first; efa; efa=efa->next) { + totface++; + if(efa->v4) totface++; + } + + /* create laplacian */ + sys = laplacian_system_construct_begin(totvert, totface, 0); + + sys->rigid.mesh= em; + sys->rigid.R = MEM_callocN(sizeof(float)*3*3*totvert, "RigidDeformR"); + sys->rigid.rhs = MEM_callocN(sizeof(float)*3*totvert, "RigidDeformRHS"); + sys->rigid.origco = MEM_callocN(sizeof(float)*3*totvert, "RigidDeformCo"); + + for(a=0, eve=em->verts.first; eve; eve=eve->next, a++) + VecCopyf(sys->rigid.origco[a], eve->co); + + sys->areaweights= 0; + sys->storeweights= 1; + + rigid_laplacian_create(sys); + + laplacian_system_construct_end(sys); + + RigidDeformSystem = sys; +} + +void rigid_deform_end(int cancel) +{ + LaplacianSystem *sys = RigidDeformSystem; + + if(sys) { + EditMesh *em = sys->rigid.mesh; + EditVert *eve; + int a; + + if(cancel) + for(a=0, eve=em->verts.first; eve; eve=eve->next, a++) + if(!eve->pinned) + VecCopyf(eve->co, sys->rigid.origco[a]); + + if(sys->rigid.R) MEM_freeN(sys->rigid.R); + if(sys->rigid.rhs) MEM_freeN(sys->rigid.rhs); + if(sys->rigid.origco) MEM_freeN(sys->rigid.origco); + + /* free */ + laplacian_system_delete(sys); + } + + RigidDeformSystem = NULL; +} +#endif + +/************************** Harmonic Coordinates ****************************/ +/* From "Harmonic Coordinates for Character Articulation", + Pushkar Joshi, Mark Meyer, Tony DeRose, Brian Green and Tom Sanocki, + SIGGRAPH 2007. */ + +#define EPSILON 0.0001f + +#define MESHDEFORM_TAG_UNTYPED 0 +#define MESHDEFORM_TAG_BOUNDARY 1 +#define MESHDEFORM_TAG_INTERIOR 2 +#define MESHDEFORM_TAG_EXTERIOR 3 + +#define MESHDEFORM_LEN_THRESHOLD 1e-6 + +#define MESHDEFORM_MIN_INFLUENCE 0.0005 + +static int MESHDEFORM_OFFSET[7][3] = + {{0,0,0}, {1,0,0}, {-1,0,0}, {0,1,0}, {0,-1,0}, {0,0,1}, {0,0,-1}}; + +typedef struct MDefBoundIsect { + float co[3], uvw[4]; + int nvert, v[4], facing; + float len; +} MDefBoundIsect; + +typedef struct MDefBindInfluence { + struct MDefBindInfluence *next; + float weight; + int vertex; +} MDefBindInfluence; + +typedef struct MeshDeformBind { + /* grid dimensions */ + float min[3], max[3]; + float width[3], halfwidth[3]; + int size, size3; + + /* meshes */ + DerivedMesh *cagedm; + float (*cagecos)[3]; + float (*vertexcos)[3]; + int totvert, totcagevert; + + /* grids */ + MemArena *memarena; + MDefBoundIsect *(*boundisect)[6]; + int *semibound; + int *tag; + float *phi, *totalphi; + + /* mesh stuff */ + int *inside; + float *weights; + MDefBindInfluence **dyngrid; + float cagemat[4][4]; + + /* direct solver */ + int *varidx; + + /* raytrace */ + RayTree *raytree; +} MeshDeformBind; + +/* ray intersection */ + +/* our own triangle intersection, so we can fully control the epsilons and + * prevent corner case from going wrong*/ +static int meshdeform_tri_intersect(float orig[3], float end[3], float vert0[3], + float vert1[3], float vert2[3], float *isectco, float *uvw) +{ + float edge1[3], edge2[3], tvec[3], pvec[3], qvec[3]; + float det,inv_det, u, v, dir[3], isectdir[3]; + + VECSUB(dir, end, orig); + + /* find vectors for two edges sharing vert0 */ + VECSUB(edge1, vert1, vert0); + VECSUB(edge2, vert2, vert0); + + /* begin calculating determinant - also used to calculate U parameter */ + Crossf(pvec, dir, edge2); + + /* if determinant is near zero, ray lies in plane of triangle */ + det = INPR(edge1, pvec); + + if (det == 0.0f) + return 0; + inv_det = 1.0f / det; + + /* calculate distance from vert0 to ray origin */ + VECSUB(tvec, orig, vert0); + + /* calculate U parameter and test bounds */ + u = INPR(tvec, pvec) * inv_det; + if (u < -EPSILON || u > 1.0f+EPSILON) + return 0; + + /* prepare to test V parameter */ + Crossf(qvec, tvec, edge1); + + /* calculate V parameter and test bounds */ + v = INPR(dir, qvec) * inv_det; + if (v < -EPSILON || u + v > 1.0f+EPSILON) + return 0; + + isectco[0]= (1.0f - u - v)*vert0[0] + u*vert1[0] + v*vert2[0]; + isectco[1]= (1.0f - u - v)*vert0[1] + u*vert1[1] + v*vert2[1]; + isectco[2]= (1.0f - u - v)*vert0[2] + u*vert1[2] + v*vert2[2]; + + uvw[0]= 1.0 - u - v; + uvw[1]= u; + uvw[2]= v; + + /* check if it is within the length of the line segment */ + VECSUB(isectdir, isectco, orig); + + if(INPR(dir, isectdir) < -EPSILON) + return 0; + + if(INPR(dir, dir) + EPSILON < INPR(isectdir, isectdir)) + return 0; + + return 1; +} + +/* blender's raytracer is not use now, even though it is much faster. it can + * give problems with rays falling through, so we use our own intersection + * function above with tweaked epsilons */ + +#if 0 +static MeshDeformBind *MESHDEFORM_BIND = NULL; + +static void meshdeform_ray_coords_func(RayFace *face, float **v1, float **v2, float **v3, float **v4) +{ + MFace *mface= (MFace*)face; + float (*cagecos)[3]= MESHDEFORM_BIND->cagecos; + + *v1= cagecos[mface->v1]; + *v2= cagecos[mface->v2]; + *v3= cagecos[mface->v3]; + *v4= (mface->v4)? cagecos[mface->v4]: NULL; +} + +static int meshdeform_ray_check_func(Isect *is, RayFace *face) +{ + return 1; +} + +static void meshdeform_ray_tree_create(MeshDeformBind *mdb) +{ + MFace *mface; + float min[3], max[3]; + int a, totface; + + /* create a raytrace tree from the mesh */ + INIT_MINMAX(min, max); + + for(a=0; a<mdb->totcagevert; a++) + DO_MINMAX(mdb->cagecos[a], min, max) + + MESHDEFORM_BIND= mdb; + + mface= mdb->cagedm->getFaceArray(mdb->cagedm); + totface= mdb->cagedm->getNumFaces(mdb->cagedm); + + mdb->raytree= RE_ray_tree_create(64, totface, min, max, + meshdeform_ray_coords_func, meshdeform_ray_check_func); + + for(a=0; a<totface; a++, mface++) + RE_ray_tree_add_face(mdb->raytree, mface); + + RE_ray_tree_done(mdb->raytree); +} + +static void meshdeform_ray_tree_free(MeshDeformBind *mdb) +{ + MESHDEFORM_BIND= NULL; + RE_ray_tree_free(mdb->raytree); +} +#endif + +static int meshdeform_intersect(MeshDeformBind *mdb, Isect *isec) +{ + MFace *mface; + float face[4][3], co[3], uvw[3], len, nor[3]; + int f, hit, is= 0, totface; + + isec->labda= 1e10; + + mface= mdb->cagedm->getFaceArray(mdb->cagedm); + totface= mdb->cagedm->getNumFaces(mdb->cagedm); + + for(f=0; f<totface; f++, mface++) { + VECCOPY(face[0], mdb->cagecos[mface->v1]); + VECCOPY(face[1], mdb->cagecos[mface->v2]); + VECCOPY(face[2], mdb->cagecos[mface->v3]); + + if(mface->v4) { + VECCOPY(face[3], mdb->cagecos[mface->v4]); + hit= meshdeform_tri_intersect(isec->start, isec->end, face[0], face[1], face[2], co, uvw); + + if(hit) { + CalcNormFloat(face[0], face[1], face[2], nor); + } + else { + hit= meshdeform_tri_intersect(isec->start, isec->end, face[0], face[2], face[3], co, uvw); + CalcNormFloat(face[0], face[2], face[3], nor); + } + } + else { + hit= meshdeform_tri_intersect(isec->start, isec->end, face[0], face[1], face[2], co, uvw); + CalcNormFloat(face[0], face[1], face[2], nor); + } + + if(hit) { + len= VecLenf(isec->start, co)/VecLenf(isec->start, isec->end); + if(len < isec->labda) { + isec->labda= len; + isec->face= mface; + isec->isect= (INPR(isec->vec, nor) <= 0.0f); + is= 1; + } + } + } + + return is; +} + +static MDefBoundIsect *meshdeform_ray_tree_intersect(MeshDeformBind *mdb, float *co1, float *co2) +{ + MDefBoundIsect *isect; + Isect isec; + float (*cagecos)[3]; + MFace *mface; + float vert[4][3], len; + static float epsilon[3]= {0, 0, 0}; //1e-4, 1e-4, 1e-4}; + + /* setup isec */ + memset(&isec, 0, sizeof(isec)); + isec.mode= RE_RAY_MIRROR; /* we want the closest intersection */ + isec.lay= -1; + isec.face_last= NULL; + isec.faceorig= NULL; + isec.labda= 1e10f; + + VECADD(isec.start, co1, epsilon); + VECADD(isec.end, co2, epsilon); + VECSUB(isec.vec, isec.end, isec.start); + +#if 0 + /*if(RE_ray_tree_intersect(mdb->raytree, &isec)) {*/ +#endif + + if(meshdeform_intersect(mdb, &isec)) { + len= isec.labda; + mface= isec.face; + + /* create MDefBoundIsect */ + isect= BLI_memarena_alloc(mdb->memarena, sizeof(*isect)); + + /* compute intersection coordinate */ + isect->co[0]= co1[0] + isec.vec[0]*len; + isect->co[1]= co1[1] + isec.vec[1]*len; + isect->co[2]= co1[2] + isec.vec[2]*len; + + isect->len= VecLenf(co1, isect->co); + if(isect->len < MESHDEFORM_LEN_THRESHOLD) + isect->len= MESHDEFORM_LEN_THRESHOLD; + + isect->v[0]= mface->v1; + isect->v[1]= mface->v2; + isect->v[2]= mface->v3; + isect->v[3]= mface->v4; + isect->nvert= (mface->v4)? 4: 3; + + isect->facing= isec.isect; + + /* compute mean value coordinates for interpolation */ + cagecos= mdb->cagecos; + VECCOPY(vert[0], cagecos[mface->v1]); + VECCOPY(vert[1], cagecos[mface->v2]); + VECCOPY(vert[2], cagecos[mface->v3]); + if(mface->v4) VECCOPY(vert[3], cagecos[mface->v4]); + MeanValueWeights(vert, isect->nvert, isect->co, isect->uvw); + + return isect; + } + + return NULL; +} + +static int meshdeform_inside_cage(MeshDeformBind *mdb, float *co) +{ + MDefBoundIsect *isect; + float outside[3], start[3], dir[3]; + int i, counter; + + for(i=1; i<=6; i++) { + counter = 0; + + outside[0] = co[0] + (mdb->max[0] - mdb->min[0] + 1.0f)*MESHDEFORM_OFFSET[i][0]; + outside[1] = co[1] + (mdb->max[1] - mdb->min[1] + 1.0f)*MESHDEFORM_OFFSET[i][1]; + outside[2] = co[2] + (mdb->max[2] - mdb->min[2] + 1.0f)*MESHDEFORM_OFFSET[i][2]; + + VECCOPY(start, co); + VECSUB(dir, outside, start); + Normalize(dir); + + isect = meshdeform_ray_tree_intersect(mdb, start, outside); + if(isect && !isect->facing) + return 1; + } + + return 0; +} + +/* solving */ + +static int meshdeform_index(MeshDeformBind *mdb, int x, int y, int z, int n) +{ + int size= mdb->size; + + x += MESHDEFORM_OFFSET[n][0]; + y += MESHDEFORM_OFFSET[n][1]; + z += MESHDEFORM_OFFSET[n][2]; + + if(x < 0 || x >= mdb->size) + return -1; + if(y < 0 || y >= mdb->size) + return -1; + if(z < 0 || z >= mdb->size) + return -1; + + return x + y*size + z*size*size; +} + +static void meshdeform_cell_center(MeshDeformBind *mdb, int x, int y, int z, int n, float *center) +{ + x += MESHDEFORM_OFFSET[n][0]; + y += MESHDEFORM_OFFSET[n][1]; + z += MESHDEFORM_OFFSET[n][2]; + + center[0]= mdb->min[0] + x*mdb->width[0] + mdb->halfwidth[0]; + center[1]= mdb->min[1] + y*mdb->width[1] + mdb->halfwidth[1]; + center[2]= mdb->min[2] + z*mdb->width[2] + mdb->halfwidth[2]; +} + +static void meshdeform_add_intersections(MeshDeformBind *mdb, int x, int y, int z) +{ + MDefBoundIsect *isect; + float center[3], ncenter[3]; + int i, a; + + a= meshdeform_index(mdb, x, y, z, 0); + meshdeform_cell_center(mdb, x, y, z, 0, center); + + /* check each outgoing edge for intersection */ + for(i=1; i<=6; i++) { + if(meshdeform_index(mdb, x, y, z, i) == -1) + continue; + + meshdeform_cell_center(mdb, x, y, z, i, ncenter); + + isect= meshdeform_ray_tree_intersect(mdb, center, ncenter); + if(isect) { + mdb->boundisect[a][i-1]= isect; + mdb->tag[a]= MESHDEFORM_TAG_BOUNDARY; + } + } +} + +static void meshdeform_bind_floodfill(MeshDeformBind *mdb) +{ + int *stack, *tag= mdb->tag; + int a, b, i, xyz[3], stacksize, size= mdb->size; + + stack= MEM_callocN(sizeof(int)*mdb->size3, "MeshDeformBindStack"); + + /* we know lower left corner is EXTERIOR because of padding */ + tag[0]= MESHDEFORM_TAG_EXTERIOR; + stack[0]= 0; + stacksize= 1; + + /* floodfill exterior tag */ + while(stacksize > 0) { + a= stack[--stacksize]; + + xyz[2]= a/(size*size); + xyz[1]= (a - xyz[2]*size*size)/size; + xyz[0]= a - xyz[1]*size - xyz[2]*size*size; + + for(i=1; i<=6; i++) { + b= meshdeform_index(mdb, xyz[0], xyz[1], xyz[2], i); + + if(b != -1) { + if(tag[b] == MESHDEFORM_TAG_UNTYPED || + (tag[b] == MESHDEFORM_TAG_BOUNDARY && !mdb->boundisect[a][i-1])) { + tag[b]= MESHDEFORM_TAG_EXTERIOR; + stack[stacksize++]= b; + } + } + } + } + + /* other cells are interior */ + for(a=0; a<size*size*size; a++) + if(tag[a]==MESHDEFORM_TAG_UNTYPED) + tag[a]= MESHDEFORM_TAG_INTERIOR; + +#if 0 + { + int tb, ti, te, ts; + tb= ti= te= ts= 0; + for(a=0; a<size*size*size; a++) + if(tag[a]==MESHDEFORM_TAG_BOUNDARY) + tb++; + else if(tag[a]==MESHDEFORM_TAG_INTERIOR) + ti++; + else if(tag[a]==MESHDEFORM_TAG_EXTERIOR) { + te++; + + if(mdb->semibound[a]) + ts++; + } + + printf("interior %d exterior %d boundary %d semi-boundary %d\n", ti, te, tb, ts); + } +#endif + + MEM_freeN(stack); +} + +static float meshdeform_boundary_phi(MeshDeformBind *mdb, MDefBoundIsect *isect, int cagevert) +{ + int a; + + for(a=0; a<isect->nvert; a++) + if(isect->v[a] == cagevert) + return isect->uvw[a]; + + return 0.0f; +} + +static float meshdeform_interp_w(MeshDeformBind *mdb, float *gridvec, float *vec, int cagevert) +{ + float dvec[3], ivec[3], wx, wy, wz, result=0.0f; + float weight, totweight= 0.0f; + int i, a, x, y, z; + + for(i=0; i<3; i++) { + ivec[i]= (int)gridvec[i]; + dvec[i]= gridvec[i] - ivec[i]; + } + + for(i=0; i<8; i++) { + if(i & 1) { x= ivec[0]+1; wx= dvec[0]; } + else { x= ivec[0]; wx= 1.0f-dvec[0]; } + + if(i & 2) { y= ivec[1]+1; wy= dvec[1]; } + else { y= ivec[1]; wy= 1.0f-dvec[1]; } + + if(i & 4) { z= ivec[2]+1; wz= dvec[2]; } + else { z= ivec[2]; wz= 1.0f-dvec[2]; } + + CLAMP(x, 0, mdb->size-1); + CLAMP(y, 0, mdb->size-1); + CLAMP(z, 0, mdb->size-1); + + a= meshdeform_index(mdb, x, y, z, 0); + weight= wx*wy*wz; + result += weight*mdb->phi[a]; + totweight += weight; + } + + if(totweight > 0.0f) + result /= totweight; + + return result; +} + +static void meshdeform_check_semibound(MeshDeformBind *mdb, int x, int y, int z) +{ + int i, a; + + a= meshdeform_index(mdb, x, y, z, 0); + if(mdb->tag[a] != MESHDEFORM_TAG_EXTERIOR) + return; + + for(i=1; i<=6; i++) + if(mdb->boundisect[a][i-1]) + mdb->semibound[a]= 1; +} + +static float meshdeform_boundary_total_weight(MeshDeformBind *mdb, int x, int y, int z) +{ + float weight, totweight= 0.0f; + int i, a; + + a= meshdeform_index(mdb, x, y, z, 0); + + /* count weight for neighbour cells */ + for(i=1; i<=6; i++) { + if(meshdeform_index(mdb, x, y, z, i) == -1) + continue; + + if(mdb->boundisect[a][i-1]) + weight= 1.0f/mdb->boundisect[a][i-1]->len; + else if(!mdb->semibound[a]) + weight= 1.0f/mdb->width[0]; + else + weight= 0.0f; + + totweight += weight; + } + + return totweight; +} + +static void meshdeform_matrix_add_cell(MeshDeformBind *mdb, int x, int y, int z) +{ + MDefBoundIsect *isect; + float weight, totweight; + int i, a, acenter; + + acenter= meshdeform_index(mdb, x, y, z, 0); + if(mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR) + return; + + nlMatrixAdd(mdb->varidx[acenter], mdb->varidx[acenter], 1.0f); + + totweight= meshdeform_boundary_total_weight(mdb, x, y, z); + for(i=1; i<=6; i++) { + a= meshdeform_index(mdb, x, y, z, i); + if(a == -1 || mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR) + continue; + + isect= mdb->boundisect[acenter][i-1]; + if (!isect) { + weight= (1.0f/mdb->width[0])/totweight; + nlMatrixAdd(mdb->varidx[acenter], mdb->varidx[a], -weight); + } + } +} + +static void meshdeform_matrix_add_rhs(MeshDeformBind *mdb, int x, int y, int z, int cagevert) +{ + MDefBoundIsect *isect; + float rhs, weight, totweight; + int i, a, acenter; + + acenter= meshdeform_index(mdb, x, y, z, 0); + if(mdb->tag[acenter] == MESHDEFORM_TAG_EXTERIOR) + return; + + totweight= meshdeform_boundary_total_weight(mdb, x, y, z); + for(i=1; i<=6; i++) { + a= meshdeform_index(mdb, x, y, z, i); + if(a == -1) + continue; + + isect= mdb->boundisect[acenter][i-1]; + + if (isect) { + weight= (1.0f/isect->len)/totweight; + rhs= weight*meshdeform_boundary_phi(mdb, isect, cagevert); + nlRightHandSideAdd(0, mdb->varidx[acenter], rhs); + } + } +} + +static void meshdeform_matrix_add_semibound_phi(MeshDeformBind *mdb, int x, int y, int z, int cagevert) +{ + MDefBoundIsect *isect; + float rhs, weight, totweight; + int i, a; + + a= meshdeform_index(mdb, x, y, z, 0); + if(!mdb->semibound[a]) + return; + + mdb->phi[a]= 0.0f; + + totweight= meshdeform_boundary_total_weight(mdb, x, y, z); + for(i=1; i<=6; i++) { + isect= mdb->boundisect[a][i-1]; + + if (isect) { + weight= (1.0f/isect->len)/totweight; + rhs= weight*meshdeform_boundary_phi(mdb, isect, cagevert); + mdb->phi[a] += rhs; + } + } +} + +static void meshdeform_matrix_add_exterior_phi(MeshDeformBind *mdb, int x, int y, int z, int cagevert) +{ + float phi, totweight; + int i, a, acenter; + + acenter= meshdeform_index(mdb, x, y, z, 0); + if(mdb->tag[acenter] != MESHDEFORM_TAG_EXTERIOR || mdb->semibound[acenter]) + return; + + phi= 0.0f; + totweight= 0.0f; + for(i=1; i<=6; i++) { + a= meshdeform_index(mdb, x, y, z, i); + + if(a != -1 && mdb->semibound[a]) { + phi += mdb->phi[a]; + totweight += 1.0f; + } + } + + if(totweight != 0.0f) + mdb->phi[acenter]= phi/totweight; +} + +static void meshdeform_matrix_solve(MeshDeformBind *mdb) +{ + NLContext *context; + float vec[3], gridvec[3]; + int a, b, x, y, z, totvar; + char message[1024]; + + /* setup variable indices */ + mdb->varidx= MEM_callocN(sizeof(int)*mdb->size3, "MeshDeformDSvaridx"); + for(a=0, totvar=0; a<mdb->size3; a++) + mdb->varidx[a]= (mdb->tag[a] == MESHDEFORM_TAG_EXTERIOR)? -1: totvar++; + + if(totvar == 0) { + MEM_freeN(mdb->varidx); + return; + } + + progress_bar(0, "Starting mesh deform solve"); + + /* setup opennl solver */ + nlNewContext(); + context= nlGetCurrent(); + + nlSolverParameteri(NL_NB_VARIABLES, totvar); + nlSolverParameteri(NL_NB_ROWS, totvar); + nlSolverParameteri(NL_NB_RIGHT_HAND_SIDES, 1); + + nlBegin(NL_SYSTEM); + nlBegin(NL_MATRIX); + + /* build matrix */ + for(z=0; z<mdb->size; z++) + for(y=0; y<mdb->size; y++) + for(x=0; x<mdb->size; x++) + meshdeform_matrix_add_cell(mdb, x, y, z); + + /* solve for each cage vert */ + for(a=0; a<mdb->totcagevert; a++) { + if(a != 0) { + nlBegin(NL_SYSTEM); + nlBegin(NL_MATRIX); + } + + /* fill in right hand side and solve */ + for(z=0; z<mdb->size; z++) + for(y=0; y<mdb->size; y++) + for(x=0; x<mdb->size; x++) + meshdeform_matrix_add_rhs(mdb, x, y, z, a); + + nlEnd(NL_MATRIX); + nlEnd(NL_SYSTEM); + +#if 0 + nlPrintMatrix(); +#endif + + if(nlSolveAdvanced(NULL, NL_TRUE)) { + for(z=0; z<mdb->size; z++) + for(y=0; y<mdb->size; y++) + for(x=0; x<mdb->size; x++) + meshdeform_matrix_add_semibound_phi(mdb, x, y, z, a); + + for(z=0; z<mdb->size; z++) + for(y=0; y<mdb->size; y++) + for(x=0; x<mdb->size; x++) + meshdeform_matrix_add_exterior_phi(mdb, x, y, z, a); + + for(b=0; b<mdb->size3; b++) { + if(mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR) + mdb->phi[b]= nlGetVariable(0, mdb->varidx[b]); + mdb->totalphi[b] += mdb->phi[b]; + } + + if(mdb->weights) { + /* static bind : compute weights for each vertex */ + for(b=0; b<mdb->totvert; b++) { + if(mdb->inside[b]) { + VECCOPY(vec, mdb->vertexcos[b]); + Mat4MulVecfl(mdb->cagemat, vec); + gridvec[0]= (vec[0] - mdb->min[0] - mdb->halfwidth[0])/mdb->width[0]; + gridvec[1]= (vec[1] - mdb->min[1] - mdb->halfwidth[1])/mdb->width[1]; + gridvec[2]= (vec[2] - mdb->min[2] - mdb->halfwidth[2])/mdb->width[2]; + + mdb->weights[b*mdb->totcagevert + a]= meshdeform_interp_w(mdb, gridvec, vec, a); + } + } + } + else { + MDefBindInfluence *inf; + + /* dynamic bind */ + for(b=0; b<mdb->size3; b++) { + if(mdb->phi[b] >= MESHDEFORM_MIN_INFLUENCE) { + inf= BLI_memarena_alloc(mdb->memarena, sizeof(*inf)); + inf->vertex= a; + inf->weight= mdb->phi[b]; + inf->next= mdb->dyngrid[b]; + mdb->dyngrid[b]= inf; + } + } + } + } + else { + error("Mesh Deform: failed to find solution."); + break; + } + + sprintf(message, "Mesh deform solve %d / %d |||", a+1, mdb->totcagevert); + progress_bar((float)(a+1)/(float)(mdb->totcagevert), message); + } + +#if 0 + /* sanity check */ + for(b=0; b<mdb->size3; b++) + if(mdb->tag[b] != MESHDEFORM_TAG_EXTERIOR) + if(fabs(mdb->totalphi[b] - 1.0f) > 1e-4) + printf("totalphi deficiency [%s|%d] %d: %.10f\n", + (mdb->tag[b] == MESHDEFORM_TAG_INTERIOR)? "interior": "boundary", mdb->semibound[b], mdb->varidx[b], mdb->totalphi[b]); +#endif + + /* free */ + MEM_freeN(mdb->varidx); + + nlDeleteContext(context); +} + +void harmonic_coordinates_bind(Scene *scene, MeshDeformModifierData *mmd, float (*vertexcos)[3], int totvert, float cagemat[][4]) +{ + MeshDeformBind mdb; + MDefBindInfluence *inf; + MDefInfluence *mdinf; + MDefCell *cell; + MVert *mvert; + float center[3], vec[3], maxwidth, totweight; + int a, b, x, y, z, totinside, offset; + + waitcursor(1); + start_progress_bar(); + + memset(&mdb, 0, sizeof(MeshDeformBind)); + + /* get mesh and cage mesh */ + mdb.vertexcos= vertexcos; + mdb.totvert= totvert; + + mdb.cagedm= mesh_create_derived_no_deform(scene, mmd->object, NULL, CD_MASK_BAREMESH); + mdb.totcagevert= mdb.cagedm->getNumVerts(mdb.cagedm); + mdb.cagecos= MEM_callocN(sizeof(*mdb.cagecos)*mdb.totcagevert, "MeshDeformBindCos"); + Mat4CpyMat4(mdb.cagemat, cagemat); + + mvert= mdb.cagedm->getVertArray(mdb.cagedm); + for(a=0; a<mdb.totcagevert; a++) + VECCOPY(mdb.cagecos[a], mvert[a].co) + + /* compute bounding box of the cage mesh */ + INIT_MINMAX(mdb.min, mdb.max); + + for(a=0; a<mdb.totcagevert; a++) + DO_MINMAX(mdb.cagecos[a], mdb.min, mdb.max); + + /* allocate memory */ + mdb.size= (2<<(mmd->gridsize-1)) + 2; + mdb.size3= mdb.size*mdb.size*mdb.size; + mdb.tag= MEM_callocN(sizeof(int)*mdb.size3, "MeshDeformBindTag"); + mdb.phi= MEM_callocN(sizeof(float)*mdb.size3, "MeshDeformBindPhi"); + mdb.totalphi= MEM_callocN(sizeof(float)*mdb.size3, "MeshDeformBindTotalPhi"); + mdb.boundisect= MEM_callocN(sizeof(*mdb.boundisect)*mdb.size3, "MDefBoundIsect"); + mdb.semibound= MEM_callocN(sizeof(int)*mdb.size3, "MDefSemiBound"); + + mdb.inside= MEM_callocN(sizeof(int)*mdb.totvert, "MDefInside"); + + if(mmd->flag & MOD_MDEF_DYNAMIC_BIND) + mdb.dyngrid= MEM_callocN(sizeof(MDefBindInfluence*)*mdb.size3, "MDefDynGrid"); + else + mdb.weights= MEM_callocN(sizeof(float)*mdb.totvert*mdb.totcagevert, "MDefWeights"); + + mdb.memarena= BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE); + BLI_memarena_use_calloc(mdb.memarena); + + /* make bounding box equal size in all directions, add padding, and compute + * width of the cells */ + maxwidth = -1.0f; + for(a=0; a<3; a++) + if(mdb.max[a]-mdb.min[a] > maxwidth) + maxwidth= mdb.max[a]-mdb.min[a]; + + for(a=0; a<3; a++) { + center[a]= (mdb.min[a]+mdb.max[a])*0.5f; + mdb.min[a]= center[a] - maxwidth*0.5f; + mdb.max[a]= center[a] + maxwidth*0.5f; + + mdb.width[a]= (mdb.max[a]-mdb.min[a])/(mdb.size-4); + mdb.min[a] -= 2.1f*mdb.width[a]; + mdb.max[a] += 2.1f*mdb.width[a]; + + mdb.width[a]= (mdb.max[a]-mdb.min[a])/mdb.size; + mdb.halfwidth[a]= mdb.width[a]*0.5f; + } + + progress_bar(0, "Setting up mesh deform system"); + +#if 0 + /* create ray tree */ + meshdeform_ray_tree_create(&mdb); +#endif + + totinside= 0; + for(a=0; a<mdb.totvert; a++) { + VECCOPY(vec, mdb.vertexcos[a]); + Mat4MulVecfl(mdb.cagemat, vec); + mdb.inside[a]= meshdeform_inside_cage(&mdb, vec); + if(mdb.inside[a]) + totinside++; + } + + /* free temporary MDefBoundIsects */ + BLI_memarena_free(mdb.memarena); + mdb.memarena= BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE); + + /* start with all cells untyped */ + for(a=0; a<mdb.size3; a++) + mdb.tag[a]= MESHDEFORM_TAG_UNTYPED; + + /* detect intersections and tag boundary cells */ + for(z=0; z<mdb.size; z++) + for(y=0; y<mdb.size; y++) + for(x=0; x<mdb.size; x++) + meshdeform_add_intersections(&mdb, x, y, z); + +#if 0 + /* free ray tree */ + meshdeform_ray_tree_free(&mdb); +#endif + + /* compute exterior and interior tags */ + meshdeform_bind_floodfill(&mdb); + + for(z=0; z<mdb.size; z++) + for(y=0; y<mdb.size; y++) + for(x=0; x<mdb.size; x++) + meshdeform_check_semibound(&mdb, x, y, z); + + /* solve */ + meshdeform_matrix_solve(&mdb); + + /* assign results */ + mmd->bindcos= (float*)mdb.cagecos; + mmd->totvert= mdb.totvert; + mmd->totcagevert= mdb.totcagevert; + Mat4CpyMat4(mmd->bindmat, mmd->object->obmat); + + if(mmd->flag & MOD_MDEF_DYNAMIC_BIND) { + mmd->totinfluence= 0; + for(a=0; a<mdb.size3; a++) + for(inf=mdb.dyngrid[a]; inf; inf=inf->next) + mmd->totinfluence++; + + /* convert MDefBindInfluences to smaller MDefInfluences */ + mmd->dyngrid= MEM_callocN(sizeof(MDefCell)*mdb.size3, "MDefDynGrid"); + mmd->dyninfluences= MEM_callocN(sizeof(MDefInfluence)*mmd->totinfluence, "MDefInfluence"); + offset= 0; + for(a=0; a<mdb.size3; a++) { + cell= &mmd->dyngrid[a]; + cell->offset= offset; + + totweight= 0.0f; + mdinf= mmd->dyninfluences + cell->offset; + for(inf=mdb.dyngrid[a]; inf; inf=inf->next, mdinf++) { + mdinf->weight= inf->weight; + mdinf->vertex= inf->vertex; + totweight += mdinf->weight; + cell->totinfluence++; + } + + if(totweight > 0.0f) { + mdinf= mmd->dyninfluences + cell->offset; + for(b=0; b<cell->totinfluence; b++, mdinf++) + mdinf->weight /= totweight; + } + + offset += cell->totinfluence; + } + + mmd->dynverts= mdb.inside; + mmd->dyngridsize= mdb.size; + VECCOPY(mmd->dyncellmin, mdb.min); + mmd->dyncellwidth= mdb.width[0]; + MEM_freeN(mdb.dyngrid); + } + else { + mmd->bindweights= mdb.weights; + MEM_freeN(mdb.inside); + } + + /* transform bindcos to world space */ + for(a=0; a<mdb.totcagevert; a++) + Mat4MulVecfl(mmd->object->obmat, mmd->bindcos+a*3); + + /* free */ + mdb.cagedm->release(mdb.cagedm); + MEM_freeN(mdb.tag); + MEM_freeN(mdb.phi); + MEM_freeN(mdb.totalphi); + MEM_freeN(mdb.boundisect); + MEM_freeN(mdb.semibound); + BLI_memarena_free(mdb.memarena); + + end_progress_bar(); + waitcursor(0); +} + |