From 07d2549f74ba8a74f0be07b2f8da098d94b8b3ff Mon Sep 17 00:00:00 2001 From: Fabio Ganovelli ganovelli Date: Wed, 14 Jan 2009 17:36:23 +0000 Subject: [PATCH] added for computing the curvature from sliding. first version: works only with faces --- .../filter_meshing/curvature_from_sliding.h | 286 ++++++++++++++++++ src/meshlabplugins/filter_meshing/frame.h | 170 +++++++++++ 2 files changed, 456 insertions(+) create mode 100644 src/meshlabplugins/filter_meshing/curvature_from_sliding.h create mode 100644 src/meshlabplugins/filter_meshing/frame.h diff --git a/src/meshlabplugins/filter_meshing/curvature_from_sliding.h b/src/meshlabplugins/filter_meshing/curvature_from_sliding.h new file mode 100644 index 000000000..ab115e328 --- /dev/null +++ b/src/meshlabplugins/filter_meshing/curvature_from_sliding.h @@ -0,0 +1,286 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "frame.h" + +template +class CurvatureFromSliding +{ + +public: + void Compute(MeshType & _mesh, typename MeshType::ScalarType rad ){ + + vcg::tri::UpdateFlags::FaceProjection(_mesh); + + oPos = vcg::tri::Allocator::AddPerVertexAttribute(patch,"oPpos"); + typename MeshType::VertexIterator vi; + mesh = &_mesh; + vcg::Box3 bbox = mesh->bbox; + typename MeshType::ScalarType infl = bbox.Diag()*0.1; + typename MeshType::CoordType bb(infl,infl,infl); + bbox.min-=bb; + bbox.max+=bb; + + grid.Set((*mesh).face.begin(),(*mesh).face.end(),bbox); + for(vi = mesh->vert.begin(); vi!= mesh->vert.end(); ++vi){ + v = &(*vi); + Sample( rad); + OnVertex(); + } + vcg::tri::Allocator::DeletePerVertexAttribute(patch,oPos); + } + + +private: + typedef vcg::GridStaticPtr MeshGridType; + MeshGridType grid; + MeshType patch,*mesh; + typename MeshType::VertexPointer v; + typename MeshType::ScalarType ballsize; + typename MeshType::PerVertexAttributeHandle oPos; + typename MeshType::CoordType basedir ; + + static typename MeshType::FaceType * GetClosestFaceBase( MeshType & mesh,MeshGridType & gr,const typename MeshGridType::CoordType & _p, + const typename MeshGridType::ScalarType & _maxDist,typename MeshGridType::ScalarType & _minDist, + typename MeshGridType::CoordType &_closestPt) + { + typedef typename MeshGridType::ScalarType ScalarType; + typedef Point3 Point3x; + typedef vcg::tri::FaceTmark MarkerFace; + MarkerFace mf; + mf.SetMesh(&mesh); + vcg::face::PointDistanceBaseFunctor PDistFunct; + _minDist=_maxDist; + return (gr.GetClosest(PDistFunct,mf,_p,_maxDist,_minDist,_closestPt)); + } + + + void Sample( typename MeshType::ScalarType rad){ + + MeshType::FaceType * f =0; + MeshType::ScalarType dist; + + std::vector closests; + std::vector distances; + std::vector points; + MeshType::ScalarType distance; + MeshType::CoordType point; + + typedef vcg::tri::FaceTmark MarkerFace; + MarkerFace mf; + mf.SetMesh(mesh); + typedef vcg::face::PointDistanceBaseFunctor FDistFunct; + (this->grid.GetInSphere/**/ + (FDistFunct(),mf, this->v->P(),rad,closests,distances,points)); + + patch.Clear(); + vcg::tri::SubSet(patch,closests); + + // sampling points + vcg::Point3f uu,vv; + vcg::GetUV( v->N(),uu,vv); + + std::vector pts; // collection of sampling points + + int ip = 0; + const int nc = 2; + int np[3]={ 4,8 }; + float alpha,dalpha,r; + alpha = 0;r = rad/nc; + pts.push_back( v->P()); + for(int ic = 0; icP() + uu*ps[0] + vv * ps[1]; + GetClosestFaceBase((*mesh),grid,pp,rad,distance,point); + pts.push_back(point); + alpha+=dalpha; + } + } + + patch.Clear(); + + vcg::tri::Allocator::AddVertices(patch,pts.size()); + for(int i = 0; i < patch.vert.size(); ++i) { + patch.vert[i].P() = oPos[patch.vert[i]] = pts[i]; + } + + vcg::tri::UpdateBounding::Box(patch); + this->ballsize = patch.bbox.Diag(); + +} + + + void OnVertex(){ + + static int m = 0; + static double *fvec = NULL; + + const int n = 2; + + if(patch.vert.size()>m){ + if (fvec!=NULL) {free(fvec); fvec = NULL;} + fvec = (double*)malloc(patch.vert.size()*sizeof(double)); + } + m = patch.vert.size(); + + if(m==0) return; + + double x[2]; + memset(&fvec[0],0,m*sizeof(double)); + memset(x,0,n*sizeof(double)); + + x[0] = x[1] = 0.0; + + MeshType::CoordType d1,d2; + vcg::GetUV(v->N(),basedir,d2); + vcg::Matrix44f tmp; + + double info[9]; + /* I: opts[0-4] = minim. options [\mu, \epsilon1, \epsilon2, \epsilon3, \delta]. Respectively the + * scale factor for initial \mu, stopping thresholds for ||J^T e||_inf, ||Dp||_2 and ||e||_2 and + * the step used in difference approximation to the jacobian. Set to NULL for defaults to be used. + * If \delta<0, the jacobian is approximated with central differences which are more accurate + * (but slower!) compared to the forward differences employed by default. + */ + double opts[5] = { + LM_INIT_MU + ,LM_STOP_THRESH + ,LM_STOP_THRESH + ,LM_STOP_THRESH + ,-0.01//LM_DIFF_DELTA// + }; + + dlevmar_dif(&CurvatureFromSliding::eval,x,&fvec[0],n,m,50,opts,info,NULL,NULL,this); + + tmp.SetRotateRad(x[0],v->N()); + v->PD1() = tmp * basedir; + v->K1() = x[1]; + + //printf("las err %f\n",maxerr); + + + float c = x[1]; + + x[0] = x[0] + M_PI * 0.5; + tmp.SetRotateRad(x[0],v->N()); + v->PD2() = tmp * basedir; + memset(&fvec[0],0,m*sizeof(double)); + + dlevmar_dif(&CurvatureFromSliding::eval_only_k,&x[1],&fvec[0],1,m,50,NULL,info,NULL,NULL,this); + + tmp.SetRotateRad(x[0],v->N()); + + v->K2() = x[1]; + + printf("%f %f \n",v->K1(),v->K2()); + if(fabs((float)v->K1())> fabs((float)v->K2())){ + std::swap(v->PD2(),v->PD1()); + std::swap(v->K2(),v->K1()); + } + +} +static void eval( double *x,double *p, int n,int m, void *data){ +// x[0] --> angle +// x[1] --> curvature + CurvatureFromSliding * cfs = (CurvatureFromSliding *)data; + typename MeshType::VertexPointer v = cfs->v; + MeshType & patch = cfs->patch; + MeshType * mesh = cfs->mesh; + MeshGridType & grid = cfs->grid; + typename MeshType::CoordType basedir =cfs->basedir; + typename MeshType::ScalarType ballsize = cfs->ballsize; + typename MeshType::PerVertexAttributeHandle &oPos=cfs->oPos; + + typename MeshType::CoordType d1; + vcg::Matrix44f rot,rt,rot90; + rot.SetRotateRad(x[0],v->cN()); + d1 = rot * basedir; + Frame::BuildRotoTranslation(v->cP(),v->cN(),d1,x[1], ballsize*0.1,rt); + + float tot_dist = 0.0; + // trasform the mesh + for(int i = 0; i < patch.vert.size(); ++i) + patch.vert[i].P() = rt * oPos[ patch.vert[i]]; + + typename MeshType::FaceType * f = 0; + float dist_upper_bound = patch.bbox.Diag(); + + float dist; + vcg::Point3f normf,bestq,ip; + // sampling + for(int i = 0 ; i < patch.vert.size(); ++i){ + f=NULL; + if(grid.bbox.IsIn(patch.vert[i].P())){ + dist = dist_upper_bound; + f = GetClosestFaceBase(*mesh, grid, patch.vert[i].P(), dist_upper_bound, dist, bestq); + } + else f=NULL; + if(f==NULL) + p[i] = 0; + else + p[i] = dist; + tot_dist+=p[i]; + } + } + + +static void eval_only_k( double *x,double *p, int n,int m, void *data){ +// *x --> curvature + vcg::Matrix44f rt; + typename MeshType::CoordType d1; + CurvatureFromSliding * cfs = (CurvatureFromSliding *)data; + typename MeshType::VertexPointer v = cfs->v; + MeshType & patch = cfs->patch; + MeshType * mesh = cfs->mesh; + MeshGridType & grid = cfs->grid; + typename MeshType::ScalarType ballsize = cfs->ballsize; + typename MeshType::CoordType basedir =cfs->basedir; + typename MeshType::PerVertexAttributeHandle &oPos=cfs->oPos; + + Frame::BuildRotoTranslation( v->cP(), + v->cN(), + v->PD2(), + *x, + ballsize*0.1, + rt); + + float tot_dist = 0.0; + // trasforma la mesh + for(int i = 0; i < patch.vert.size(); ++i) + patch.vert[i].P() = rt * oPos[patch.vert[i]]; + + typename MeshType::FaceType * f = 0; + float dist_upper_bound = patch.bbox.Diag(); + + float dist; + vcg::Point3f normf,bestq,ip; + // sampling + for(int i = 0 ; i < patch.vert.size(); ++i){ + f=NULL; + if(grid.bbox.IsIn(patch.vert[i].P())){ + dist = dist_upper_bound; + f = GetClosestFaceBase(* mesh , grid, patch.vert[i].P(), dist_upper_bound, dist, bestq); + } + else f=NULL; + if(f==NULL) + p[i] = dist_upper_bound; + else + p[i] = dist; + tot_dist+=p[i]; + } + + } + + +}; \ No newline at end of file diff --git a/src/meshlabplugins/filter_meshing/frame.h b/src/meshlabplugins/filter_meshing/frame.h new file mode 100644 index 000000000..6d25a89b6 --- /dev/null +++ b/src/meshlabplugins/filter_meshing/frame.h @@ -0,0 +1,170 @@ +#ifndef __FRAME__ +#define __FRAME__ + +#include +#include + +template +struct Frame{ + S aX,aY,aZ; // euler angles + vcg::Point3 orig; // origine del reference frame + vcg::Matrix33 rot; + S & Ang(const int & i){ return (i==0)?aX:( (i==1)? aY : aZ);} + + void angles(vcg::Matrix33 mm) { + rot = mm; + vcg::Point3 x = mm.GetRow(0),y = mm.GetRow(1),z = mm.GetRow(2),xp,yp,zp; + + // rotation around z that brings x on XZ + float norm = vcg::math::Sqrt(x[0]*x[0]+x[1]*x[1]); + float cosa = (norm < 0.0001)?1.0: x[0]/norm; + aZ = vcg::math::Acos(cosa); + if(x[1] < 0) aZ= -aZ; + + vcg::Matrix33 mp; + + mp.RotateR(-aZ,vcg::Point3(0,0,1)); + xp = mp * x; + yp = mp * y; + zp = mp * z; + + // rotation around y that brings (x su XY) su X + norm = vcg::math::Sqrt(xp[0]*xp[0]+xp[2]*xp[2]); + cosa = (norm < 0.0001)?1.0:xp[0]/norm; + aY = vcg::math::Acos(cosa); + if(xp[2] > 0 ) aY= -aY; + + mp.RotateR(-aY,vcg::Point3(0,1,0)); + xp = mp * xp; + yp = mp * yp; + zp = mp * zp; + + norm = vcg::math::Sqrt(yp[1]*yp[1]+yp[2]*yp[2]); + cosa = (norm < 0.0001)?1.0:yp[1]/norm; + aX = vcg::math::Acos(cosa); + if(yp[2] < 0) aX= -aX; + + mp.RotateR(-aX,vcg::Point3(1,0,0)); + xp = mp * xp; + yp = mp * yp; + zp = mp * zp; + } + + Frame(){} + Frame(vcg::Matrix33 m, vcg::Point3 O){ + angles(m); + orig = O; + } + + void angles(vcg::Matrix44 mm) { + vcg::Matrix33 m; + memcpy(&m[0][0],&mm[0][0],3*sizeof(S)); + memcpy(&m[1][0],&mm[1][0],3*sizeof(S)); + memcpy(&m[2][0],&mm[2][0],3*sizeof(S)); + angles(m); + } + + vcg::Point3 Apply(vcg::Point3 p){ + + vcg::Matrix33 mp; + mp.RotateR(aX,vcg::Point3(1,0,0)); + p = mp * p; + mp.RotateR(aY,vcg::Point3(0,1,0)); + p = mp * p; + mp.RotateR(aZ,vcg::Point3(0,0,1)); + p = mp * p; + + p += orig; + return p; + } + + void ToMatrix(vcg::Matrix44 &m){ + vcg::Matrix33 mp,tmp; + mp.RotateR(aX,vcg::Point3(1,0,0)); + tmp.RotateR(aY,vcg::Point3(0,1,0)); + mp = mp * tmp; + tmp.RotateR(aZ,vcg::Point3(0,0,1)); + mp = mp * tmp; + memcpy(&m[0][0],&mp[0][0],3*sizeof(S)); + memcpy(&m[1][0],&mp[1][0],3*sizeof(S)); + memcpy(&m[2][0],&mp[2][0],3*sizeof(S)); + m[0][3] = orig[0]; + m[1][3] = orig[1]; + m[2][3] = orig[2]; + } + + Frame(vcg::Matrix44 m4){ + vcg::Matrix33 m; + memcpy(&m[0][0],&m4[0][0],3*sizeof(S)); + memcpy(&m[1][0],&m4[1][0],3*sizeof(S)); + memcpy(&m[2][0],&m4[2][0],3*sizeof(S)); + + m.Transpose(); + angles(m); + orig[0] = m4[0][3]; + orig[1] = m4[1][3]; + orig[2] = m4[2][3]; + } + + static void BuildRotoTranslation(const vcg::Point3 & point, + const vcg::Point3 & norm, + const vcg::Point3 & dir, + const S & k, + const S & delta, + vcg::Matrix44 & rt){ + + vcg::Matrix44 rt2d; + rt2d.SetIdentity(); + + // rt2d is a rototranslation around the axis parallel to z + // and passing though (0,-1/k,0) + rt2d[0][0] = rt2d[1][1] = cos(k*delta); + rt2d[0][1] = -sin(k*delta); + rt2d[1][0] = -rt2d[0][1]; + + + // taylor expansion for sinx/x and (cosx-1)/x + S delta2 = delta*delta; + S delta3 = delta2*delta; + S delta4 = delta3*delta; + S delta5 = delta4*delta; + S delta6 = delta5*delta; + S delta7 = delta6*delta; + S k2 = k*k; + S k3 = k2*k; + S k4 = k3*k; + S k5 = k4*k; + S k6 = k5*k; + + float sinkd = sin ( k* delta) / (-k); + float coskd = (1.0 - cos ( k* delta)) / (-k); -1.0 + cos (k*delta)/k; + rt2d[0][3] = -(delta - k2 * delta3/6.0 + k4*delta5/120.0 -k6*delta7/5040); + rt2d[1][3] = - k * delta2/2.0 + k3*delta4/24 -k5*delta6/720; + + // build a frame with y = norm e x = dir + vcg::Point3 z = (dir ^ norm).Normalize(); + + // fw (from world) brings the canonical frame on the frame (dir,norm,z, point) + vcg::Matrix44 tw,fw; + fw.SetIdentity(); + memcpy(&fw[0][0],&dir[0] ,3*sizeof(S)); + memcpy(&fw[1][0],&norm[0] ,3*sizeof(S)); + memcpy(&fw[2][0],&z[0] ,3*sizeof(S)); + fw[0][3] = - dir * point; + fw[1][3] = - norm * point; + fw[2][3] = - z * point; + + + tw = fw; + vcg::Transpose(tw); + memset(&tw[3][0],0 ,3*sizeof(S)); + tw[0][3] = point[0]; + tw[1][3] = point[1]; + tw[2][3] = point[2]; + + rt = tw * rt2d * fw;// write the result +} + +}; + +#endif \ No newline at end of file