mirror of
https://github.com/lucaspalomodevelop/meshlab.git
synced 2026-03-18 02:24:38 +00:00
added for computing the curvature from sliding.
first version: works only with faces
This commit is contained in:
Fabio Ganovelli ganovelli
parent
319a412a6f
commit
07d2549f74
286
src/meshlabplugins/filter_meshing/curvature_from_sliding.h
Normal file
286
src/meshlabplugins/filter_meshing/curvature_from_sliding.h
Normal file
@ -0,0 +1,286 @@
|
||||
#include<vcg/complex/trimesh/update/position.h>
|
||||
#include<vcg/complex/trimesh/update/edges.h>
|
||||
#include <vcg/complex/trimesh/closest.h>
|
||||
#include <vcg/complex/trimesh/append.h>
|
||||
#include <vcg/complex/trimesh/clean.h>
|
||||
#include <vcg/complex/trimesh/base.h>
|
||||
#include <vcg/complex/trimesh/allocate.h>
|
||||
#include <vcg/space/index/grid_static_ptr.h>
|
||||
#include <lib/levmar/lm.h>
|
||||
|
||||
#include "frame.h"
|
||||
|
||||
template <class MeshType>
|
||||
class CurvatureFromSliding
|
||||
{
|
||||
|
||||
public:
|
||||
void Compute(MeshType & _mesh, typename MeshType::ScalarType rad ){
|
||||
|
||||
vcg::tri::UpdateFlags<MeshType>::FaceProjection(_mesh);
|
||||
|
||||
oPos = vcg::tri::Allocator<MeshType>::AddPerVertexAttribute<typename MeshType::CoordType>(patch,"oPpos");
|
||||
typename MeshType::VertexIterator vi;
|
||||
mesh = &_mesh;
|
||||
vcg::Box3<typename MeshType::ScalarType> 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<MeshType>::DeletePerVertexAttribute<typename MeshType::CoordType>(patch,oPos);
|
||||
}
|
||||
|
||||
|
||||
private:
|
||||
typedef vcg::GridStaticPtr <typename MeshType::FaceType, typename MeshType::ScalarType> MeshGridType;
|
||||
MeshGridType grid;
|
||||
MeshType patch,*mesh;
|
||||
typename MeshType::VertexPointer v;
|
||||
typename MeshType::ScalarType ballsize;
|
||||
typename MeshType::PerVertexAttributeHandle<typename MeshType::CoordType> 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<ScalarType> Point3x;
|
||||
typedef vcg::tri::FaceTmark<MeshType> MarkerFace;
|
||||
MarkerFace mf;
|
||||
mf.SetMesh(&mesh);
|
||||
vcg::face::PointDistanceBaseFunctor<ScalarType> 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<MeshType:: FaceType*> closests;
|
||||
std::vector<MeshType::ScalarType> distances;
|
||||
std::vector<MeshType::CoordType> points;
|
||||
MeshType::ScalarType distance;
|
||||
MeshType::CoordType point;
|
||||
|
||||
typedef vcg::tri::FaceTmark<MeshType> MarkerFace;
|
||||
MarkerFace mf;
|
||||
mf.SetMesh(mesh);
|
||||
typedef vcg::face::PointDistanceBaseFunctor<typename MeshType::ScalarType> FDistFunct;
|
||||
(this->grid.GetInSphere/*<FDistFunct,MarkerFace,OBJPTRCONTAINER,DISTCONTAINER,POINTCONTAINER>*/
|
||||
(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<vcg::Point3f > 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; ic<nc ;++ic){
|
||||
r =(rad/nc)*(ic+1);
|
||||
alpha = 0;
|
||||
dalpha = 2*M_PI / np[ic];
|
||||
for(int ip = 0; ip < np[ic];++ip){
|
||||
vcg::Point3f ps(cos(alpha)*r,sin(alpha)*r,0.0);
|
||||
vcg::Point3f pp = v->P() + uu*ps[0] + vv * ps[1];
|
||||
GetClosestFaceBase((*mesh),grid,pp,rad,distance,point);
|
||||
pts.push_back(point);
|
||||
alpha+=dalpha;
|
||||
}
|
||||
}
|
||||
|
||||
patch.Clear();
|
||||
|
||||
vcg::tri::Allocator<MeshType>::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<MeshType>::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<MeshType> * cfs = (CurvatureFromSliding<MeshType> *)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<typename MeshType::CoordType> &oPos=cfs->oPos;
|
||||
|
||||
typename MeshType::CoordType d1;
|
||||
vcg::Matrix44f rot,rt,rot90;
|
||||
rot.SetRotateRad(x[0],v->cN());
|
||||
d1 = rot * basedir;
|
||||
Frame<float>::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<MeshType> * cfs = (CurvatureFromSliding<MeshType> *)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<typename MeshType::CoordType> &oPos=cfs->oPos;
|
||||
|
||||
Frame<float>::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];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
170
src/meshlabplugins/filter_meshing/frame.h
Normal file
170
src/meshlabplugins/filter_meshing/frame.h
Normal file
@ -0,0 +1,170 @@
|
||||
#ifndef __FRAME__
|
||||
#define __FRAME__
|
||||
|
||||
#include <vcg/math/matrix44.h>
|
||||
#include <vcg/math/matrix33.h>
|
||||
|
||||
template <class S>
|
||||
struct Frame{
|
||||
S aX,aY,aZ; // euler angles
|
||||
vcg::Point3<S> orig; // origine del reference frame
|
||||
vcg::Matrix33<S> rot;
|
||||
S & Ang(const int & i){ return (i==0)?aX:( (i==1)? aY : aZ);}
|
||||
|
||||
void angles(vcg::Matrix33<S> mm) {
|
||||
rot = mm;
|
||||
vcg::Point3<S> 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<S> mp;
|
||||
|
||||
mp.RotateR(-aZ,vcg::Point3<S>(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<S>(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<S>(1,0,0));
|
||||
xp = mp * xp;
|
||||
yp = mp * yp;
|
||||
zp = mp * zp;
|
||||
}
|
||||
|
||||
Frame(){}
|
||||
Frame(vcg::Matrix33<S> m, vcg::Point3<S> O){
|
||||
angles(m);
|
||||
orig = O;
|
||||
}
|
||||
|
||||
void angles(vcg::Matrix44<S> mm) {
|
||||
vcg::Matrix33<S> 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<S> Apply(vcg::Point3<S> p){
|
||||
|
||||
vcg::Matrix33<S> mp;
|
||||
mp.RotateR(aX,vcg::Point3<S>(1,0,0));
|
||||
p = mp * p;
|
||||
mp.RotateR(aY,vcg::Point3<S>(0,1,0));
|
||||
p = mp * p;
|
||||
mp.RotateR(aZ,vcg::Point3<S>(0,0,1));
|
||||
p = mp * p;
|
||||
|
||||
p += orig;
|
||||
return p;
|
||||
}
|
||||
|
||||
void ToMatrix(vcg::Matrix44<S> &m){
|
||||
vcg::Matrix33<S> mp,tmp;
|
||||
mp.RotateR(aX,vcg::Point3<S>(1,0,0));
|
||||
tmp.RotateR(aY,vcg::Point3<S>(0,1,0));
|
||||
mp = mp * tmp;
|
||||
tmp.RotateR(aZ,vcg::Point3<S>(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<S> m4){
|
||||
vcg::Matrix33<S> 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<S> & point,
|
||||
const vcg::Point3<S> & norm,
|
||||
const vcg::Point3<S> & dir,
|
||||
const S & k,
|
||||
const S & delta,
|
||||
vcg::Matrix44<S> & rt){
|
||||
|
||||
vcg::Matrix44<S> 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<S> z = (dir ^ norm).Normalize();
|
||||
|
||||
// fw (from world) brings the canonical frame on the frame (dir,norm,z, point)
|
||||
vcg::Matrix44<S> 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
|
||||
Loading…
x
Reference in New Issue
Block a user