lucaspalomodevelop/meshlab
1
0
Fork 0
mirror of https://github.com/lucaspalomodevelop/meshlab.git synced 2026-03-20 03:16:10 +00:00
Code Issues Packages Projects Releases Wiki Activity

Corrections of include paths to comply the new folder arrangement of the VCG library (a second set of the many filters; removing useless includes)

Browse Source
This commit is contained in:
Paolo Cignoni cignoni 2011-04-03 23:51:26 +00:00
parent 3a2303e66d
commit 1128b129ea
25 changed files with 9204 additions and 9216 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1124
src/meshlabplugins/filter_isoparametrization/diam_parametrization.h
View File

File diff suppressed because it is too large Load Diff

348
src/meshlabplugins/filter_isoparametrization/diam_topology.h
View File

@ -1,175 +1,175 @@
#include <algorithm>
#include <time.h>
#include <vcg/complex/trimesh/refine.h>
#include <vcg/simplex/face/pos.h>
#include <vcg/space/color4.h>
#include <vcg/space/intersection2.h>
#ifndef _DIAMONDTOPOLOGY
#define _DIAMONDTOPOLOGY
class DiamTopology
{
typedef struct DiamondPatch
{
int FFp[4];
int FFi[4];
AbstractFace *f;
int edgeInd;
DiamondPatch(){
FFp[0]=-1;
FFp[1]=-1;
FFp[2]=-1;
FFp[3]=-1;
FFi[0]=-1;
FFi[1]=-1;
FFi[2]=-1;
FFi[3]=-1;
}
};
std::vector<DiamondPatch> patches;
IsoParametrization *isoParam;
typedef std::pair<AbstractFace*,int> TriEdge;
std::vector<DiamondPatch> DDAdiacency;
std::map<TriEdge,int> edgeMap;
void InsertEdges()
{
//insert all the edges
AbstractMesh *domain=isoParam->AbsMesh();
int index=0;
for (int i=0;i<domain->fn;i++)
{
AbstractFace *f0=&(domain->face[i]);
if (!f0->IsD())
{
for(int j=0;j<3;j++)
{
AbstractFace * f1=f0->FFp(j);
if (f1>f0)
{
DDAdiacency.push_back(DiamondPatch());
DDAdiacency[index].f=f0;
DDAdiacency[index].edgeInd=j;
TriEdge key=TriEdge(f0,j);
edgeMap.insert(std::pair<TriEdge,int>(key,index));
index++;
}
}
}
}
}
void SetTopology()
{
///data used to calculate topology
typedef std::pair<AbstractVertex*,AbstractFace*> DiamEdge;
struct InfoAdiacency
{
int index_global;
int index_local;
};
std::map<DiamEdge,InfoAdiacency> EdgeDiamMap;
for (unsigned int i=0;i<DDAdiacency.size();i++)
{
int indexEdge=DDAdiacency[i].edgeInd;
AbstractFace* face=DDAdiacency[i].f;
AbstractVertex *v0=face->V0(indexEdge);
AbstractVertex *v1=face->V1(indexEdge);
AbstractFace *f0=face;
AbstractFace *f1=face->FFp(indexEdge);
///set right order
if (v1<v0)
std::swap(v0,v1);
if (f1<f0)
std::swap(f0,f1);
///then put in the map the variuos edge
DiamEdge de[4];
de[0]=DiamEdge(v0,f0);
de[1]=DiamEdge(v1,f0);
de[2]=DiamEdge(v1,f1);
de[3]=DiamEdge(v0,f1);
for (int j=0;j<4;j++)
{
///find if already exist
std::map<DiamEdge,InfoAdiacency>::iterator ItMap=EdgeDiamMap.find(de[j]);
///first time insert in the table
if (ItMap==EdgeDiamMap.end())
{
InfoAdiacency InfAd;
InfAd.index_global=i;
InfAd.index_local=j;
EdgeDiamMap.insert(std::pair<DiamEdge,InfoAdiacency>(de[j],InfAd));
}
else
{
int index_global=(*ItMap).second.index_global;
int index_local =(*ItMap).second.index_local;
DDAdiacency[index_global].FFp[index_local]=i;
DDAdiacency[index_global].FFi[index_local]=j;
DDAdiacency[i].FFp[j]=index_global;
DDAdiacency[i].FFi[j]=index_local;
}
}
}
}
bool TestTopology()
{
for (unsigned int i=0;i<DDAdiacency.size();i++)
{
for (int j=0;j<4;j++)
{
int index_global=DDAdiacency[i].FFp[j];
int index_local=DDAdiacency[i].FFi[j];
if ((index_global==-1)||(index_local==-1))
return false;
int opposite0=DDAdiacency[index_global].FFp[index_local];
int index_local1=DDAdiacency[index_global].FFi[index_local];
///int opposite1=DDAdiacency[i].FFp[index_local];
if (opposite0!=i)
return false;
if (index_local1!=j)
return false;
}
}
return true;
}
public:
///initialize the parameterization
void Init(IsoParametrization *_isoParam)
{
isoParam=_isoParam;
///initialize edges
InsertEdges();
///then set topology
SetTopology();
assert(TestTopology());
}
};
#include <algorithm>
#include <time.h>
#include <vcg/complex/algorithms/refine.h>
#include <vcg/simplex/face/pos.h>
#include <vcg/space/color4.h>
#include <vcg/space/intersection2.h>
#ifndef _DIAMONDTOPOLOGY
#define _DIAMONDTOPOLOGY
class DiamTopology
{
typedef struct DiamondPatch
{
int FFp[4];
int FFi[4];
AbstractFace *f;
int edgeInd;
DiamondPatch(){
FFp[0]=-1;
FFp[1]=-1;
FFp[2]=-1;
FFp[3]=-1;
FFi[0]=-1;
FFi[1]=-1;
FFi[2]=-1;
FFi[3]=-1;
}
};
std::vector<DiamondPatch> patches;
IsoParametrization *isoParam;
typedef std::pair<AbstractFace*,int> TriEdge;
std::vector<DiamondPatch> DDAdiacency;
std::map<TriEdge,int> edgeMap;
void InsertEdges()
{
//insert all the edges
AbstractMesh *domain=isoParam->AbsMesh();
int index=0;
for (int i=0;i<domain->fn;i++)
{
AbstractFace *f0=&(domain->face[i]);
if (!f0->IsD())
{
for(int j=0;j<3;j++)
{
AbstractFace * f1=f0->FFp(j);
if (f1>f0)
{
DDAdiacency.push_back(DiamondPatch());
DDAdiacency[index].f=f0;
DDAdiacency[index].edgeInd=j;
TriEdge key=TriEdge(f0,j);
edgeMap.insert(std::pair<TriEdge,int>(key,index));
index++;
}
}
}
}
}
void SetTopology()
{
///data used to calculate topology
typedef std::pair<AbstractVertex*,AbstractFace*> DiamEdge;
struct InfoAdiacency
{
int index_global;
int index_local;
};
std::map<DiamEdge,InfoAdiacency> EdgeDiamMap;
for (unsigned int i=0;i<DDAdiacency.size();i++)
{
int indexEdge=DDAdiacency[i].edgeInd;
AbstractFace* face=DDAdiacency[i].f;
AbstractVertex *v0=face->V0(indexEdge);
AbstractVertex *v1=face->V1(indexEdge);
AbstractFace *f0=face;
AbstractFace *f1=face->FFp(indexEdge);
///set right order
if (v1<v0)
std::swap(v0,v1);
if (f1<f0)
std::swap(f0,f1);
///then put in the map the variuos edge
DiamEdge de[4];
de[0]=DiamEdge(v0,f0);
de[1]=DiamEdge(v1,f0);
de[2]=DiamEdge(v1,f1);
de[3]=DiamEdge(v0,f1);
for (int j=0;j<4;j++)
{
///find if already exist
std::map<DiamEdge,InfoAdiacency>::iterator ItMap=EdgeDiamMap.find(de[j]);
///first time insert in the table
if (ItMap==EdgeDiamMap.end())
{
InfoAdiacency InfAd;
InfAd.index_global=i;
InfAd.index_local=j;
EdgeDiamMap.insert(std::pair<DiamEdge,InfoAdiacency>(de[j],InfAd));
}
else
{
int index_global=(*ItMap).second.index_global;
int index_local =(*ItMap).second.index_local;
DDAdiacency[index_global].FFp[index_local]=i;
DDAdiacency[index_global].FFi[index_local]=j;
DDAdiacency[i].FFp[j]=index_global;
DDAdiacency[i].FFi[j]=index_local;
}
}
}
}
bool TestTopology()
{
for (unsigned int i=0;i<DDAdiacency.size();i++)
{
for (int j=0;j<4;j++)
{
int index_global=DDAdiacency[i].FFp[j];
int index_local=DDAdiacency[i].FFi[j];
if ((index_global==-1)||(index_local==-1))
return false;
int opposite0=DDAdiacency[index_global].FFp[index_local];
int index_local1=DDAdiacency[index_global].FFi[index_local];
///int opposite1=DDAdiacency[i].FFp[index_local];
if (opposite0!=i)
return false;
if (index_local1!=j)
return false;
}
}
return true;
}
public:
///initialize the parameterization
void Init(IsoParametrization *_isoParam)
{
isoParam=_isoParam;
///initialize edges
InsertEdges();
///then set topology
SetTopology();
assert(TestTopology());
}
};
#endif

472
src/meshlabplugins/filter_isoparametrization/diamond_sampler.h
View File

@ -1,237 +1,237 @@
#ifndef _DIAMSAMPLER
#define _DIAMSAMPLER
#include <stat_remeshing.h>
#include <vcg/complex/trimesh/clean.h>
class DiamSampler{
typedef IsoParametrization::CoordType CoordType;
typedef IsoParametrization::ScalarType ScalarType;
std::vector<std::vector<std::vector<CoordType> > > SampledPos;
IsoParametrization *isoParam;
unsigned int sampleSize;
void DeAllocatePos()
{
///positions
for (unsigned int i=0;i<SampledPos.size();i++)
{
for (unsigned int j=0;j<SampledPos[i].size();j++)
SampledPos[i][j].clear();
SampledPos[i].clear();
}
SampledPos.clear();
}
void AllocatePos(const int & sizeSampl)
{
///positions
AbstractMesh *domain=isoParam->AbsMesh();
int num_edges=0;
for (unsigned int i=0;i<domain->face.size();i++)
{
AbstractFace *f=&domain->face[i];
for (int j=0;j<3;j++)
if (f->FFp(j)>f)
num_edges++;
}
SampledPos.resize(num_edges);
for (unsigned int i=0;i<SampledPos.size();i++)
{
SampledPos[i].resize(sizeSampl);
for (unsigned int j=0;j<SampledPos[i].size();j++)
SampledPos[i][j].resize(sizeSampl);
}
}
CoordType AveragePos(const std::vector<ParamFace*> &faces,
std::vector<CoordType> &barys)
{
CoordType pos=CoordType(0,0,0);
for (unsigned int i=0;i<faces.size();i++)
{
pos+=(faces[i]->V(0)->P()*barys[i].X()+
faces[i]->V(1)->P()*barys[i].Y()+
faces[i]->V(2)->P()*barys[i].Z());
}
pos/=(ScalarType)faces.size();
return pos;
}
int n_diamonds;
int inFace;
int inEdge;
int inStar;
int n_merged;
public:
///initialize the sampler
void Init(IsoParametrization *_isoParam)
{
isoParam=_isoParam;
}
template <class OutputMesh>
void GetMesh(OutputMesh &SaveMesh)
{
typedef typename OutputMesh::FaceType MyFace;
typedef typename OutputMesh::VertexType MyVertex;
SaveMesh.Clear();
SaveMesh.vert.reserve(SampledPos.size()*
sampleSize*
sampleSize);
SaveMesh.face.reserve(SampledPos.size()*
(sampleSize-1)*
(sampleSize-1)*2);
SaveMesh.vn=0;
SaveMesh.fn=0;
///suposed to be the same everywhere
///allocate space
std::vector<std::vector<MyVertex*> > vertMatrix;
vertMatrix.resize(sampleSize);
for (unsigned int i=0;i<sampleSize;i++)
vertMatrix[i].resize(sampleSize);
///add vertices
for (unsigned int i=0;i<SampledPos.size();i++)
//for (unsigned int i=11;i<12;i++)
{
///add vertices & positions
for (unsigned int j=0;j<sampleSize;j++)
for (unsigned int k=0;k<sampleSize;k++)
{
SaveMesh.vert.push_back(MyVertex());
SaveMesh.vert.back().P()=SampledPos[i][j][k];
vertMatrix[j][k]=&SaveMesh.vert.back();
SaveMesh.vn++;
}
/*printf("samppling %d\n",i);*/
///add faces
for (unsigned int j=0;j<sampleSize-1;j++)
for (unsigned int k=0;k<sampleSize-1;k++)
{
MyVertex *v0=vertMatrix[j][k];
MyVertex *v1=vertMatrix[j+1][k];
MyVertex *v2=vertMatrix[j+1][k+1];
MyVertex *v3=vertMatrix[j][k+1];
SaveMesh.face.push_back(MyFace());
SaveMesh.face.back().V(0)=v0;
SaveMesh.face.back().V(1)=v1;
SaveMesh.face.back().V(2)=v3;
SaveMesh.face.push_back(MyFace());
SaveMesh.face.back().V(0)=v1;
SaveMesh.face.back().V(1)=v2;
SaveMesh.face.back().V(2)=v3;
SaveMesh.fn+=2;
}
}
///get minimum edge size
ScalarType minE,maxE;
MaxMinEdge<OutputMesh>(SaveMesh,minE,maxE);
/*int num_tri=SampledPos.size()*sampleSize*sampleSize*2;
ScalarType Area_mesh=Area<OutputMesh>(SaveMesh);
ScalarType Edge=sqrt((((Area_mesh/(ScalarType)num_tri)*4.0)/(ScalarType)sqrt(3.0)));*/
n_merged=vcg::tri::Clean<OutputMesh>::MergeCloseVertex(SaveMesh,(ScalarType)minE*(ScalarType)0.9);
vcg::tri::UpdateNormals<OutputMesh>::PerVertexNormalized(SaveMesh);
/*Log("Merged %d vertices\n",merged);*/
}
//typedef enum SampleAttr{SMNormal, SMColor, SMPosition};
///sample the parametrization
bool SamplePos(const int &size)
{
if (size<2)
return false;
sampleSize=size;
DeAllocatePos();//clear old data
AllocatePos(size); ///allocate for new one
inFace=0;
inEdge=0;
inStar=0;
int global=0;
///start sampling values
/*Log("Num Diamonds: %d \n",SampledPos.size());*/
for (unsigned int diam=0;diam<SampledPos.size();diam++)
for (unsigned int j=0;j<sampleSize;j++)
for (unsigned int k=0;k<sampleSize;k++)
{
vcg::Point2<ScalarType> UVQuad,UV;
UVQuad.X()=(ScalarType)j/(ScalarType)(sampleSize-1);
UVQuad.Y()=(ScalarType)k/(ScalarType)(sampleSize-1);
int I;
//printf("Quad: %d,%f,%f \n",diam,UV.X(),UV.Y());
///get coordinate in parametric space
isoParam->inv_GE1Quad(diam,UVQuad,I,UV);
//printf("Alfabeta: %d,%f,%f \n",I,UV.X(),UV.Y());
///and sample
std::vector<ParamFace*> faces;
std::vector<CoordType> barys;
int domain=isoParam->Theta(I,UV,faces,barys);
if (domain==0)
inFace++;
else
if (domain==1)
inEdge++;
else
if (domain==2)
inStar++;
global++;
//printf("Find in domain: %d \n",domain);
///store value
CoordType val=AveragePos(faces,barys);
/*if (domain==2)
val=CoordType(0,0,0);*/
SampledPos[diam][j][k]=val;
}
return true;
/*#ifndef _MESHLAB
printf("In Face: %f \n",(ScalarType)inFace/(ScalarType)global);
printf("In Diamond: %f \n",(ScalarType)inEdge/(ScalarType)global);
printf("In Star: %f \n",(ScalarType)inStar/(ScalarType)global);
#endif*/
/*Log("In Face: %f \n",(ScalarType)inFace/(ScalarType)global);
Log("In Diamond: %f \n",(ScalarType)inEdge/(ScalarType)global);
Log("In Star: %f \n",(ScalarType)inStar/(ScalarType)global);*/
}
void getResData(int &_n_diamonds,int &_inFace,
int &_inEdge,int &_inStar,int &_n_merged)
{
_n_diamonds=n_diamonds;
_inFace=inFace;
_inEdge=inEdge;
_inStar=inStar;
_n_merged=n_merged;
}
void GetCoords(std::vector<CoordType> &positions)
{
for (unsigned int diam=0;diam<SampledPos.size();diam++)
for (unsigned int j=0;j<sampleSize;j++)
for (unsigned int k=0;k<sampleSize;k++)
positions.push_back(SampledPos[diam][j][k]);
}
};
#ifndef _DIAMSAMPLER
#define _DIAMSAMPLER
#include <stat_remeshing.h>
#include <vcg/complex/algorithms/clean.h>
class DiamSampler{
typedef IsoParametrization::CoordType CoordType;
typedef IsoParametrization::ScalarType ScalarType;
std::vector<std::vector<std::vector<CoordType> > > SampledPos;
IsoParametrization *isoParam;
unsigned int sampleSize;
void DeAllocatePos()
{
///positions
for (unsigned int i=0;i<SampledPos.size();i++)
{
for (unsigned int j=0;j<SampledPos[i].size();j++)
SampledPos[i][j].clear();
SampledPos[i].clear();
}
SampledPos.clear();
}
void AllocatePos(const int & sizeSampl)
{
///positions
AbstractMesh *domain=isoParam->AbsMesh();
int num_edges=0;
for (unsigned int i=0;i<domain->face.size();i++)
{
AbstractFace *f=&domain->face[i];
for (int j=0;j<3;j++)
if (f->FFp(j)>f)
num_edges++;
}
SampledPos.resize(num_edges);
for (unsigned int i=0;i<SampledPos.size();i++)
{
SampledPos[i].resize(sizeSampl);
for (unsigned int j=0;j<SampledPos[i].size();j++)
SampledPos[i][j].resize(sizeSampl);
}
}
CoordType AveragePos(const std::vector<ParamFace*> &faces,
std::vector<CoordType> &barys)
{
CoordType pos=CoordType(0,0,0);
for (unsigned int i=0;i<faces.size();i++)
{
pos+=(faces[i]->V(0)->P()*barys[i].X()+
faces[i]->V(1)->P()*barys[i].Y()+
faces[i]->V(2)->P()*barys[i].Z());
}
pos/=(ScalarType)faces.size();
return pos;
}
int n_diamonds;
int inFace;
int inEdge;
int inStar;
int n_merged;
public:
///initialize the sampler
void Init(IsoParametrization *_isoParam)
{
isoParam=_isoParam;
}
template <class OutputMesh>
void GetMesh(OutputMesh &SaveMesh)
{
typedef typename OutputMesh::FaceType MyFace;
typedef typename OutputMesh::VertexType MyVertex;
SaveMesh.Clear();
SaveMesh.vert.reserve(SampledPos.size()*
sampleSize*
sampleSize);
SaveMesh.face.reserve(SampledPos.size()*
(sampleSize-1)*
(sampleSize-1)*2);
SaveMesh.vn=0;
SaveMesh.fn=0;
///suposed to be the same everywhere
///allocate space
std::vector<std::vector<MyVertex*> > vertMatrix;
vertMatrix.resize(sampleSize);
for (unsigned int i=0;i<sampleSize;i++)
vertMatrix[i].resize(sampleSize);
///add vertices
for (unsigned int i=0;i<SampledPos.size();i++)
//for (unsigned int i=11;i<12;i++)
{
///add vertices & positions
for (unsigned int j=0;j<sampleSize;j++)
for (unsigned int k=0;k<sampleSize;k++)
{
SaveMesh.vert.push_back(MyVertex());
SaveMesh.vert.back().P()=SampledPos[i][j][k];
vertMatrix[j][k]=&SaveMesh.vert.back();
SaveMesh.vn++;
}
/*printf("samppling %d\n",i);*/
///add faces
for (unsigned int j=0;j<sampleSize-1;j++)
for (unsigned int k=0;k<sampleSize-1;k++)
{
MyVertex *v0=vertMatrix[j][k];
MyVertex *v1=vertMatrix[j+1][k];
MyVertex *v2=vertMatrix[j+1][k+1];
MyVertex *v3=vertMatrix[j][k+1];
SaveMesh.face.push_back(MyFace());
SaveMesh.face.back().V(0)=v0;
SaveMesh.face.back().V(1)=v1;
SaveMesh.face.back().V(2)=v3;
SaveMesh.face.push_back(MyFace());
SaveMesh.face.back().V(0)=v1;
SaveMesh.face.back().V(1)=v2;
SaveMesh.face.back().V(2)=v3;
SaveMesh.fn+=2;
}
}
///get minimum edge size
ScalarType minE,maxE;
MaxMinEdge<OutputMesh>(SaveMesh,minE,maxE);
/*int num_tri=SampledPos.size()*sampleSize*sampleSize*2;
ScalarType Area_mesh=Area<OutputMesh>(SaveMesh);
ScalarType Edge=sqrt((((Area_mesh/(ScalarType)num_tri)*4.0)/(ScalarType)sqrt(3.0)));*/
n_merged=vcg::tri::Clean<OutputMesh>::MergeCloseVertex(SaveMesh,(ScalarType)minE*(ScalarType)0.9);
vcg::tri::UpdateNormals<OutputMesh>::PerVertexNormalized(SaveMesh);
/*Log("Merged %d vertices\n",merged);*/
}
//typedef enum SampleAttr{SMNormal, SMColor, SMPosition};
///sample the parametrization
bool SamplePos(const int &size)
{
if (size<2)
return false;
sampleSize=size;
DeAllocatePos();//clear old data
AllocatePos(size); ///allocate for new one
inFace=0;
inEdge=0;
inStar=0;
int global=0;
///start sampling values
/*Log("Num Diamonds: %d \n",SampledPos.size());*/
for (unsigned int diam=0;diam<SampledPos.size();diam++)
for (unsigned int j=0;j<sampleSize;j++)
for (unsigned int k=0;k<sampleSize;k++)
{
vcg::Point2<ScalarType> UVQuad,UV;
UVQuad.X()=(ScalarType)j/(ScalarType)(sampleSize-1);
UVQuad.Y()=(ScalarType)k/(ScalarType)(sampleSize-1);
int I;
//printf("Quad: %d,%f,%f \n",diam,UV.X(),UV.Y());
///get coordinate in parametric space
isoParam->inv_GE1Quad(diam,UVQuad,I,UV);
//printf("Alfabeta: %d,%f,%f \n",I,UV.X(),UV.Y());
///and sample
std::vector<ParamFace*> faces;
std::vector<CoordType> barys;
int domain=isoParam->Theta(I,UV,faces,barys);
if (domain==0)
inFace++;
else
if (domain==1)
inEdge++;
else
if (domain==2)
inStar++;
global++;
//printf("Find in domain: %d \n",domain);
///store value
CoordType val=AveragePos(faces,barys);
/*if (domain==2)
val=CoordType(0,0,0);*/
SampledPos[diam][j][k]=val;
}
return true;
/*#ifndef _MESHLAB
printf("In Face: %f \n",(ScalarType)inFace/(ScalarType)global);
printf("In Diamond: %f \n",(ScalarType)inEdge/(ScalarType)global);
printf("In Star: %f \n",(ScalarType)inStar/(ScalarType)global);
#endif*/
/*Log("In Face: %f \n",(ScalarType)inFace/(ScalarType)global);
Log("In Diamond: %f \n",(ScalarType)inEdge/(ScalarType)global);
Log("In Star: %f \n",(ScalarType)inStar/(ScalarType)global);*/
}
void getResData(int &_n_diamonds,int &_inFace,
int &_inEdge,int &_inStar,int &_n_merged)
{
_n_diamonds=n_diamonds;
_inFace=inFace;
_inEdge=inEdge;
_inStar=inStar;
_n_merged=n_merged;
}
void GetCoords(std::vector<CoordType> &positions)
{
for (unsigned int diam=0;diam<SampledPos.size();diam++)
for (unsigned int j=0;j<sampleSize;j++)
for (unsigned int k=0;k<sampleSize;k++)
positions.push_back(SampledPos[diam][j][k]);
}
};
#endif

1063
src/meshlabplugins/filter_isoparametrization/filter_isoparametrization.cpp
View File

File diff suppressed because it is too large Load Diff

4030
src/meshlabplugins/filter_isoparametrization/iso_parametrization.h
View File

File diff suppressed because it is too large Load Diff

270
src/meshlabplugins/filter_isoparametrization/iso_transfer.h
View File

@ -1,135 +1,135 @@
#include <iso_parametrization.h>
#include<vcg/simplex/edge/base.h>
#include<vcg/simplex/vertex/base.h>
#include<vcg/simplex/face/base.h>
#include <vcg/complex/trimesh/base.h>
#include <vcg/complex/trimesh/update/topology.h>
#include <vcg/complex/trimesh/update/edges.h>
#include <vcg/complex/trimesh/update/bounding.h>
#include <vcg/complex/trimesh/update/flag.h>
#include <vcg/space/index/grid_static_ptr.h>
#include <vcg/complex/trimesh/closest.h>
#ifndef _ISO_TRANSFER
#define _ISO_TRANSFER
class IsoTransfer
{
typedef vcg::GridStaticPtr<ParamMesh::FaceType, ParamMesh::ScalarType> TriMeshGrid;
typedef ParamMesh::CoordType CoordType;
typedef ParamMesh::ScalarType ScalarType;
TriMeshGrid TRGrid;
void Clamp(CoordType &bary)
{
/* float eps=0.01;*/
float sum=0;
int bigger=0;
int lower=0;
/* for (int i=0;i<3;i++)
{
if ((bary.V(i)<eps)&&(bary.V(i)>-eps))
bary.V(i)=0;
if ((bary.V(i)<1+eps)&&(bary.V(i)>1-eps))
bary.V(i)=1;
sum+=bary.V(i);
if (bary.V(i)>bary.V(bigger))
bigger=i;
if (bary.V(i)<bary.V(lower))
lower=i;
}
assert(bigger!=lower);
if (sum>(1.0+eps))
{
float diff=sum-1.0;
bary.V(bigger)-=diff;
}
else
if (sum<(1.0-eps))
{
float diff=1.0-sum;
bary.V(lower)+=diff;
}*/
for (int i=0;i<3;i++)
{
if (bary.V(i)<0)
bary.V(i)=0;
if (bary.V(i)>1)
bary.V(i)=1;
sum+=bary.V(i);
if (bary.V(i)>bary.V(bigger))
bigger=i;
if (bary.V(i)<bary.V(lower))
lower=i;
}
//assert(bigger!=lower);
if (sum>(1.0))
{
float diff=sum-1.0;
bary.V(bigger)-=diff;
}
else
if (sum<(1.0))
{
float diff=1.0-sum;
bary.V(lower)+=diff;
}
}
public:
template <class MeshType>
void Transfer(IsoParametrization &IsoParam,
MeshType &to_assing)
{
///put the mesh in the grid
typedef typename MeshType::ScalarType ScalarType;
vcg::tri::UpdateBounding<ParamMesh>::Box(*IsoParam.ParaMesh());
vcg::tri::UpdateNormals<ParamMesh>::PerFaceNormalized(*IsoParam.ParaMesh());
vcg::tri::UpdateNormals<ParamMesh>::PerVertexAngleWeighted(*IsoParam.ParaMesh());
vcg::tri::UpdateNormals<ParamMesh>::NormalizeVertex(*IsoParam.ParaMesh());
vcg::tri::UpdateEdges<ParamMesh>::Set(*IsoParam.ParaMesh());
vcg::tri::UpdateFlags<ParamMesh>::FaceProjection(*IsoParam.ParaMesh());
TRGrid.Set(IsoParam.ParaMesh()->face.begin(),IsoParam.ParaMesh()->face.end());
ScalarType maxDist=IsoParam.ParaMesh()->bbox.Diag();
///then for each vertex find the closest
for (size_t i=0;i<to_assing.vert.size();i++)
{
typename MeshType::VertexType *vert=&to_assing.vert[i];
if (!vert->IsD())
{
ScalarType dist;
CoordType closest,norm,bary;
ParamMesh::FaceType * f=NULL;
f=GetClosestFace(*IsoParam.ParaMesh(),TRGrid,vert->P(), maxDist,dist,closest,norm,bary);
assert(f!=NULL);
///then find back the coordinates
if (!((bary.X()>=0)&&(bary.X()<=1)&&
(bary.Y()>=0)&&(bary.Y()<=1)&&
(bary.Z()>=0)&&(bary.Z()<=1)))
{
printf("%i,%3.3f,%3.3f,%3.3f",int(i),bary.X(),bary.Y(),bary.Z());
system("pause");
}
Clamp(bary);
int I;
vcg::Point2<ScalarType> UV;
IsoParam.Phi(f,bary,I,UV);
///and finally set to the vertex
assert(I>=0);
vert->T().P()=UV;
vert->T().N()=I;
vert->Q()=(typename MeshType::ScalarType)I;
}
}
}
};
#endif
#include <iso_parametrization.h>
#include<vcg/simplex/edge/base.h>
#include<vcg/simplex/vertex/base.h>
#include<vcg/simplex/face/base.h>
#include <vcg/complex/complex.h>
#include <vcg/complex/algorithms/update/topology.h>
#include <vcg/complex/algorithms/update/edges.h>
#include <vcg/complex/algorithms/update/bounding.h>
#include <vcg/complex/algorithms/update/flag.h>
#include <vcg/space/index/grid_static_ptr.h>
#include <vcg/complex/algorithms/closest.h>
#ifndef _ISO_TRANSFER
#define _ISO_TRANSFER
class IsoTransfer
{
typedef vcg::GridStaticPtr<ParamMesh::FaceType, ParamMesh::ScalarType> TriMeshGrid;
typedef ParamMesh::CoordType CoordType;
typedef ParamMesh::ScalarType ScalarType;
TriMeshGrid TRGrid;
void Clamp(CoordType &bary)
{
/* float eps=0.01;*/
float sum=0;
int bigger=0;
int lower=0;
/* for (int i=0;i<3;i++)
{
if ((bary.V(i)<eps)&&(bary.V(i)>-eps))
bary.V(i)=0;
if ((bary.V(i)<1+eps)&&(bary.V(i)>1-eps))
bary.V(i)=1;
sum+=bary.V(i);
if (bary.V(i)>bary.V(bigger))
bigger=i;
if (bary.V(i)<bary.V(lower))
lower=i;
}
assert(bigger!=lower);
if (sum>(1.0+eps))
{
float diff=sum-1.0;
bary.V(bigger)-=diff;
}
else
if (sum<(1.0-eps))
{
float diff=1.0-sum;
bary.V(lower)+=diff;
}*/
for (int i=0;i<3;i++)
{
if (bary.V(i)<0)
bary.V(i)=0;
if (bary.V(i)>1)
bary.V(i)=1;
sum+=bary.V(i);
if (bary.V(i)>bary.V(bigger))
bigger=i;
if (bary.V(i)<bary.V(lower))
lower=i;
}
//assert(bigger!=lower);
if (sum>(1.0))
{
float diff=sum-1.0;
bary.V(bigger)-=diff;
}
else
if (sum<(1.0))
{
float diff=1.0-sum;
bary.V(lower)+=diff;
}
}
public:
template <class MeshType>
void Transfer(IsoParametrization &IsoParam,
MeshType &to_assing)
{
///put the mesh in the grid
typedef typename MeshType::ScalarType ScalarType;
vcg::tri::UpdateBounding<ParamMesh>::Box(*IsoParam.ParaMesh());
vcg::tri::UpdateNormals<ParamMesh>::PerFaceNormalized(*IsoParam.ParaMesh());
vcg::tri::UpdateNormals<ParamMesh>::PerVertexAngleWeighted(*IsoParam.ParaMesh());
vcg::tri::UpdateNormals<ParamMesh>::NormalizeVertex(*IsoParam.ParaMesh());
vcg::tri::UpdateEdges<ParamMesh>::Set(*IsoParam.ParaMesh());
vcg::tri::UpdateFlags<ParamMesh>::FaceProjection(*IsoParam.ParaMesh());
TRGrid.Set(IsoParam.ParaMesh()->face.begin(),IsoParam.ParaMesh()->face.end());
ScalarType maxDist=IsoParam.ParaMesh()->bbox.Diag();
///then for each vertex find the closest
for (size_t i=0;i<to_assing.vert.size();i++)
{
typename MeshType::VertexType *vert=&to_assing.vert[i];
if (!vert->IsD())
{
ScalarType dist;
CoordType closest,norm,bary;
ParamMesh::FaceType * f=NULL;
f=GetClosestFace(*IsoParam.ParaMesh(),TRGrid,vert->P(), maxDist,dist,closest,norm,bary);
assert(f!=NULL);
///then find back the coordinates
if (!((bary.X()>=0)&&(bary.X()<=1)&&
(bary.Y()>=0)&&(bary.Y()<=1)&&
(bary.Z()>=0)&&(bary.Z()<=1)))
{
printf("%i,%3.3f,%3.3f,%3.3f",int(i),bary.X(),bary.Y(),bary.Z());
system("pause");
}
Clamp(bary);
int I;
vcg::Point2<ScalarType> UV;
IsoParam.Phi(f,bary,I,UV);
///and finally set to the vertex
assert(I>=0);
vert->T().P()=UV;
vert->T().N()=I;
vert->Q()=(typename MeshType::ScalarType)I;
}
}
}
};
#endif

2756
src/meshlabplugins/filter_isoparametrization/local_parametrization.h
View File

File diff suppressed because it is too large Load Diff

10
src/meshlabplugins/filter_isoparametrization/mesh_operators.h
View File

@ -1,11 +1,11 @@
#ifndef MESH_OPERATORS
#define MESH_OPERATORS
#include <vcg/simplex/face/pos.h>
#include <vcg/complex/trimesh/update/normal.h>
#include <vcg/complex/trimesh/update/topology.h>
#include <vcg/complex/trimesh/update/flag.h>
#include <vcg/complex/trimesh/update/bounding.h>
#include <vcg/complex/trimesh/update/edges.h>
#include <vcg/complex/algorithms/update/normal.h>
#include <vcg/complex/algorithms/update/topology.h>
#include <vcg/complex/algorithms/update/flag.h>
#include <vcg/complex/algorithms/update/bounding.h>
#include <vcg/complex/algorithms/update/edges.h>
#include <vector>
#include <map>

1404
src/meshlabplugins/filter_isoparametrization/param_collapse.h
View File

File diff suppressed because it is too large Load Diff

630
src/meshlabplugins/filter_isoparametrization/param_flip.h
View File

@ -1,315 +1,315 @@
#ifndef PARAM_FLIP
#define PARAM_FLIP
#include <vcg/complex/local_optimization/tri_edge_flip.h>
///Flip function
template <class BaseMesh>
class ParamEdgeFlip : public vcg::tri::PlanarEdgeFlip<BaseMesh, ParamEdgeFlip<BaseMesh> >
{
typedef typename BaseMesh::VertexType::EdgeType EdgeType;
typedef typename BaseMesh::VertexType BaseVertex;
typedef typename BaseMesh::VertexType VertexType;
typedef typename BaseMesh::FaceType BaseFace;
typedef typename BaseMesh::FaceType FaceType;
typedef typename BaseMesh::CoordType CoordType;
typedef typename BaseMesh::ScalarType ScalarType;
typedef vcg::tri::PlanarEdgeFlip<BaseMesh, ParamEdgeFlip<BaseMesh> > Super;
ScalarType diff;
public:
static EnergyType &EType(){static EnergyType E;return E;};
bool savedomain;
bool IsFeasible()
{
if(!vcg::face::CheckFlipEdge(*this->_pos.F(), this->_pos.E()))
return false;
return (this->_priority>0);
}
inline ParamEdgeFlip() {}
/*!
* Constructor with <I>pos</I> type
*/
inline ParamEdgeFlip(const typename Super::PosType pos, int mark)
{
this->_pos = pos;
this->_localMark = mark;
this->_priority = this->ComputePriority();
savedomain=false;
}
///do the effective flip
void ExecuteFlip(FaceType &f, const int &edge, BaseMesh *base_domain=NULL)
{
std::vector<FaceType*> faces;
faces.push_back(&f);
faces.push_back(f.FFp(edge));
std::vector<VertexType*> HresVert;
getHresVertex<FaceType>(faces,HresVert);
///parametrize H_res mesh respect to diamond
for (unsigned int i=0;i<HresVert.size();i++)
{
VertexType* v=HresVert[i];
///get father & bary
FaceType* father=v->father;
CoordType bary=v->Bary;
assert((father==faces[0])||(father==faces[1]));
vcg::Point2<ScalarType> t0=father->V(0)->T().P();
vcg::Point2<ScalarType> t1=father->V(1)->T().P();
vcg::Point2<ScalarType> t2=father->V(2)->T().P();
//assert(testBaryCoords(bary));
if(!testBaryCoords(bary))
{
printf("BAry0 :%lf,%lf,%lf",bary.X(),bary.Y(),bary.Z());
//system("pause");
}
GetUV<BaseMesh>(father,bary,v->T().U(),v->T().V());
}
///update VF topology
FaceType *f1=f.FFp(edge);
FaceType *f0=&f;
vcg::face::VFDetach(*f1,0);
vcg::face::VFDetach(*f1,1);
vcg::face::VFDetach(*f1,2);
vcg::face::VFDetach(*f0,0);
vcg::face::VFDetach(*f0,1);
vcg::face::VFDetach(*f0,2);
///then do the effective flip
vcg::face::FlipEdge(f,edge);
vcg::face::VFAppend(f1,0);
vcg::face::VFAppend(f1,1);
vcg::face::VFAppend(f1,2);
vcg::face::VFAppend(f0,0);
vcg::face::VFAppend(f0,1);
vcg::face::VFAppend(f0,2);
///edh updating topology
///set son->father new link
for (unsigned int i=0;i<HresVert.size();i++)
{
VertexType* v=HresVert[i];
ScalarType U=v->T().U();
ScalarType V=v->T().V();
CoordType bary;
int index;
bool found=GetBaryFaceFromUV(faces,U,V,bary,index);
if (!found)
{
printf("\n U : %lf; V : %lf \n",U,V);
//system("pause");
}
//assert(found);
assert(testBaryCoords(bary));
if (base_domain!=NULL)
AssingFather(*v,faces[index],bary,*base_domain);
else
{
v->father=faces[index];
assert(!faces[index]->IsD());
v->Bary=bary;
}
}
///set father->son new link
for (unsigned int i=0;i<faces.size();i++)
faces[i]->vertices_bary.clear();
for (unsigned int i=0;i<HresVert.size();i++)
{
VertexType *son=HresVert[i];
FaceType *father=son->father;
CoordType bary=son->Bary;
father->vertices_bary.push_back(std::pair<BaseVertex*,vcg::Point3f>(son,bary));
}
}
ScalarType EdgeDiff()
{
/*
1
/|\
/ | \
2 f0|f1 3
\ | /
\|/
0
*/
VertexType *v0, *v1, *v2, *v3;
int edge0 = this->_pos.E();
v0 = this->_pos.F()->V0(edge0);
v1 = this->_pos.F()->V1(edge0);
v2 = this->_pos.F()->V2(edge0);
v3 = this->_pos.F()->FFp(edge0)->V2(this->_pos.F()->FFi(edge0));
int edge1=this->_pos.F()->FFi(edge0);
FaceType* f0=this->_pos.F();
FaceType* f1=this->_pos.F()->FFp(edge0);
///parametrize all possible diamonds
///diam0 & diam1
///make a copy of the mesh
std::vector<FaceType*> OrdFace;
OrdFace.push_back(f0);
OrdFace.push_back(f1);
BaseMesh Diam;
BaseMesh DiamHres;
///create a copy of the domain and of the H resolution
CopySubMeshLevels(OrdFace,Diam,DiamHres);
///parametrize domains
ParametrizeDiamondEquilateral(Diam,edge0,edge1);
///copy parametrization on original mesh
FaceType* on_edge[2];
on_edge[0]=&Diam.face[0];
on_edge[1]=&Diam.face[1];
assert(Diam.face[0].FFp(edge0)==&Diam.face[1]);///test
assert(Diam.face[1].FFp(edge1)==&Diam.face[0]);///test
///Evaluate lenght of shared edge
ScalarType L0=EstimateLenghtByParam<BaseMesh>(Diam.face[0].V(edge0),Diam.face[0].V((edge0+1)%3),on_edge);
///do the flip on the copied mesh do not affect the original mesh
ExecuteFlip(Diam.face[0],edge0);
UpdateTopologies(&Diam);
///get the non border edge of face0
int NB_edge=-1;
if (!Diam.face[0].IsB(0))
NB_edge=0;
else
if (!Diam.face[0].IsB(1))
NB_edge=1;
else
if (!Diam.face[0].IsB(2))
NB_edge=2;
assert(NB_edge!=-1);
ScalarType L1=EstimateLenghtByParam<BaseMesh>(Diam.face[0].V(NB_edge),Diam.face[0].V((NB_edge+1)%3),on_edge);
ScalarType value=L0-L1;
diff=value;
this->_priority = 1.0/value;
return (this->_priority);
}
void Execute(BaseMesh &m)
{
assert(this->_priority>0);
/*
1
/|\
/ | \
2 f0|f1 3
\ | /
\|/
0
*/
VertexType *v0, *v1, *v2, *v3;
int edge0 = this->_pos.E();
v0 = this->_pos.F()->V0(edge0);
v1 = this->_pos.F()->V1(edge0);
v2 = this->_pos.F()->V2(edge0);
v3 = this->_pos.F()->FFp(edge0)->V2(this->_pos.F()->FFi(edge0));
///assing texcoords
ScalarType h=(sqrt((ScalarType)3.0)/(ScalarType)2.0);
v0->T().P()=vcg::Point2<ScalarType>(0,(ScalarType)-0.5);
v1->T().P()=vcg::Point2<ScalarType>(0,(ScalarType)0.5);
v2->T().P()=vcg::Point2<ScalarType>(-h,0);
v3->T().P()=vcg::Point2<ScalarType>(h,0);
#ifndef _MESHLAB
///save domain if need for demos
if (savedomain)
{
BaseMesh hlev_mesh;
std::vector<FaceType*> faces;
FaceType* f=this->_pos.F();
int edge=this->_pos.E();
faces.push_back(f);
faces.push_back(f->FFp(edge));
std::vector<VertexType*> HresVert;
getHresVertex<FaceType>(faces,HresVert);
///parametrize H_res mesh respect to diamond
std::vector<VertexType*> OrderedVertices;
std::vector<FaceType*> OrderedFaces;
CopyMeshFromVertices(HresVert,OrderedVertices,OrderedFaces,hlev_mesh);
//for (int i=0;i<hlev_mesh.vert.size();i++)
//hlev_mesh.vert[i].C()=hlev_mesh.vert[i].OriginalCol;
vcg::tri::io::ExporterPLY<BaseMesh>::Save(hlev_mesh,"c:/export_submeshes/FLIPHlev3D.ply",vcg::tri::io::Mask::IOM_VERTCOLOR);
for (unsigned int i=0;i<hlev_mesh.vert.size();i++)
{
hlev_mesh.vert[i].P().X()=hlev_mesh.vert[i].T().U();
hlev_mesh.vert[i].P().Y()=hlev_mesh.vert[i].T().V();
hlev_mesh.vert[i].P().Z()=0;
}
vcg::tri::io::ExporterPLY<BaseMesh>::Save(hlev_mesh,"c:/export_submeshes/FLIPHlevUV.ply",vcg::tri::io::Mask::IOM_VERTCOLOR);
}
#endif
ExecuteFlip(*this->_pos.F(),this->_pos.E(),&m);
UpdateTopologies(&m);
///stars optimization
/*int t0=clock();*/
/*OptimizeStar<BaseMesh>(v0);
OptimizeStar<BaseMesh>(v1);
OptimizeStar<BaseMesh>(v2);
OptimizeStar<BaseMesh>(v3);*/
SmartOptimizeStar<BaseMesh>(v0,m,Accuracy(),EType());
SmartOptimizeStar<BaseMesh>(v1,m,Accuracy(),EType());
SmartOptimizeStar<BaseMesh>(v2,m,Accuracy(),EType());
SmartOptimizeStar<BaseMesh>(v3,m,Accuracy(),EType());
/*int t1=clock();
time_opt+=(t1-t0);*/
}
ScalarType ComputePriority()
{
this->_priority=EdgeDiff();
return this->_priority;
}
BaseFace *getF()
{return this->_pos.F();}
int getE()
{return this->_pos.E();}
public:
static int &Accuracy()
{
static int _acc;
return _acc;
}
};
#endif
#ifndef PARAM_FLIP
#define PARAM_FLIP
#include <vcg/complex/algorithms/local_optimization/tri_edge_flip.h>
///Flip function
template <class BaseMesh>
class ParamEdgeFlip : public vcg::tri::PlanarEdgeFlip<BaseMesh, ParamEdgeFlip<BaseMesh> >
{
typedef typename BaseMesh::VertexType::EdgeType EdgeType;
typedef typename BaseMesh::VertexType BaseVertex;
typedef typename BaseMesh::VertexType VertexType;
typedef typename BaseMesh::FaceType BaseFace;
typedef typename BaseMesh::FaceType FaceType;
typedef typename BaseMesh::CoordType CoordType;
typedef typename BaseMesh::ScalarType ScalarType;
typedef vcg::tri::PlanarEdgeFlip<BaseMesh, ParamEdgeFlip<BaseMesh> > Super;
ScalarType diff;
public:
static EnergyType &EType(){static EnergyType E;return E;};
bool savedomain;
bool IsFeasible()
{
if(!vcg::face::CheckFlipEdge(*this->_pos.F(), this->_pos.E()))
return false;
return (this->_priority>0);
}
inline ParamEdgeFlip() {}
/*!
* Constructor with <I>pos</I> type
*/
inline ParamEdgeFlip(const typename Super::PosType pos, int mark)
{
this->_pos = pos;
this->_localMark = mark;
this->_priority = this->ComputePriority();
savedomain=false;
}
///do the effective flip
void ExecuteFlip(FaceType &f, const int &edge, BaseMesh *base_domain=NULL)
{
std::vector<FaceType*> faces;
faces.push_back(&f);
faces.push_back(f.FFp(edge));
std::vector<VertexType*> HresVert;
getHresVertex<FaceType>(faces,HresVert);
///parametrize H_res mesh respect to diamond
for (unsigned int i=0;i<HresVert.size();i++)
{
VertexType* v=HresVert[i];
///get father & bary
FaceType* father=v->father;
CoordType bary=v->Bary;
assert((father==faces[0])||(father==faces[1]));
vcg::Point2<ScalarType> t0=father->V(0)->T().P();
vcg::Point2<ScalarType> t1=father->V(1)->T().P();
vcg::Point2<ScalarType> t2=father->V(2)->T().P();
//assert(testBaryCoords(bary));
if(!testBaryCoords(bary))
{
printf("BAry0 :%lf,%lf,%lf",bary.X(),bary.Y(),bary.Z());
//system("pause");
}
GetUV<BaseMesh>(father,bary,v->T().U(),v->T().V());
}
///update VF topology
FaceType *f1=f.FFp(edge);
FaceType *f0=&f;
vcg::face::VFDetach(*f1,0);
vcg::face::VFDetach(*f1,1);
vcg::face::VFDetach(*f1,2);
vcg::face::VFDetach(*f0,0);
vcg::face::VFDetach(*f0,1);
vcg::face::VFDetach(*f0,2);
///then do the effective flip
vcg::face::FlipEdge(f,edge);
vcg::face::VFAppend(f1,0);
vcg::face::VFAppend(f1,1);
vcg::face::VFAppend(f1,2);
vcg::face::VFAppend(f0,0);
vcg::face::VFAppend(f0,1);
vcg::face::VFAppend(f0,2);
///edh updating topology
///set son->father new link
for (unsigned int i=0;i<HresVert.size();i++)
{
VertexType* v=HresVert[i];
ScalarType U=v->T().U();
ScalarType V=v->T().V();
CoordType bary;
int index;
bool found=GetBaryFaceFromUV(faces,U,V,bary,index);
if (!found)
{
printf("\n U : %lf; V : %lf \n",U,V);
//system("pause");
}
//assert(found);
assert(testBaryCoords(bary));
if (base_domain!=NULL)
AssingFather(*v,faces[index],bary,*base_domain);
else
{
v->father=faces[index];
assert(!faces[index]->IsD());
v->Bary=bary;
}
}
///set father->son new link
for (unsigned int i=0;i<faces.size();i++)
faces[i]->vertices_bary.clear();
for (unsigned int i=0;i<HresVert.size();i++)
{
VertexType *son=HresVert[i];
FaceType *father=son->father;
CoordType bary=son->Bary;
father->vertices_bary.push_back(std::pair<BaseVertex*,vcg::Point3f>(son,bary));
}
}
ScalarType EdgeDiff()
{
/*
1
/|\
/ | \
2 f0|f1 3
\ | /
\|/
0
*/
VertexType *v0, *v1, *v2, *v3;
int edge0 = this->_pos.E();
v0 = this->_pos.F()->V0(edge0);
v1 = this->_pos.F()->V1(edge0);
v2 = this->_pos.F()->V2(edge0);
v3 = this->_pos.F()->FFp(edge0)->V2(this->_pos.F()->FFi(edge0));
int edge1=this->_pos.F()->FFi(edge0);
FaceType* f0=this->_pos.F();
FaceType* f1=this->_pos.F()->FFp(edge0);
///parametrize all possible diamonds
///diam0 & diam1
///make a copy of the mesh
std::vector<FaceType*> OrdFace;
OrdFace.push_back(f0);
OrdFace.push_back(f1);
BaseMesh Diam;
BaseMesh DiamHres;
///create a copy of the domain and of the H resolution
CopySubMeshLevels(OrdFace,Diam,DiamHres);
///parametrize domains
ParametrizeDiamondEquilateral(Diam,edge0,edge1);
///copy parametrization on original mesh
FaceType* on_edge[2];
on_edge[0]=&Diam.face[0];
on_edge[1]=&Diam.face[1];
assert(Diam.face[0].FFp(edge0)==&Diam.face[1]);///test
assert(Diam.face[1].FFp(edge1)==&Diam.face[0]);///test
///Evaluate lenght of shared edge
ScalarType L0=EstimateLenghtByParam<BaseMesh>(Diam.face[0].V(edge0),Diam.face[0].V((edge0+1)%3),on_edge);
///do the flip on the copied mesh do not affect the original mesh
ExecuteFlip(Diam.face[0],edge0);
UpdateTopologies(&Diam);
///get the non border edge of face0
int NB_edge=-1;
if (!Diam.face[0].IsB(0))
NB_edge=0;
else
if (!Diam.face[0].IsB(1))
NB_edge=1;
else
if (!Diam.face[0].IsB(2))
NB_edge=2;
assert(NB_edge!=-1);
ScalarType L1=EstimateLenghtByParam<BaseMesh>(Diam.face[0].V(NB_edge),Diam.face[0].V((NB_edge+1)%3),on_edge);
ScalarType value=L0-L1;
diff=value;
this->_priority = 1.0/value;
return (this->_priority);
}
void Execute(BaseMesh &m)
{
assert(this->_priority>0);
/*
1
/|\
/ | \
2 f0|f1 3
\ | /
\|/
0
*/
VertexType *v0, *v1, *v2, *v3;
int edge0 = this->_pos.E();
v0 = this->_pos.F()->V0(edge0);
v1 = this->_pos.F()->V1(edge0);
v2 = this->_pos.F()->V2(edge0);
v3 = this->_pos.F()->FFp(edge0)->V2(this->_pos.F()->FFi(edge0));
///assing texcoords
ScalarType h=(sqrt((ScalarType)3.0)/(ScalarType)2.0);
v0->T().P()=vcg::Point2<ScalarType>(0,(ScalarType)-0.5);
v1->T().P()=vcg::Point2<ScalarType>(0,(ScalarType)0.5);
v2->T().P()=vcg::Point2<ScalarType>(-h,0);
v3->T().P()=vcg::Point2<ScalarType>(h,0);
#ifndef _MESHLAB
///save domain if need for demos
if (savedomain)
{
BaseMesh hlev_mesh;
std::vector<FaceType*> faces;
FaceType* f=this->_pos.F();
int edge=this->_pos.E();
faces.push_back(f);
faces.push_back(f->FFp(edge));
std::vector<VertexType*> HresVert;
getHresVertex<FaceType>(faces,HresVert);
///parametrize H_res mesh respect to diamond
std::vector<VertexType*> OrderedVertices;
std::vector<FaceType*> OrderedFaces;
CopyMeshFromVertices(HresVert,OrderedVertices,OrderedFaces,hlev_mesh);
//for (int i=0;i<hlev_mesh.vert.size();i++)
//hlev_mesh.vert[i].C()=hlev_mesh.vert[i].OriginalCol;
vcg::tri::io::ExporterPLY<BaseMesh>::Save(hlev_mesh,"c:/export_submeshes/FLIPHlev3D.ply",vcg::tri::io::Mask::IOM_VERTCOLOR);
for (unsigned int i=0;i<hlev_mesh.vert.size();i++)
{
hlev_mesh.vert[i].P().X()=hlev_mesh.vert[i].T().U();
hlev_mesh.vert[i].P().Y()=hlev_mesh.vert[i].T().V();
hlev_mesh.vert[i].P().Z()=0;
}
vcg::tri::io::ExporterPLY<BaseMesh>::Save(hlev_mesh,"c:/export_submeshes/FLIPHlevUV.ply",vcg::tri::io::Mask::IOM_VERTCOLOR);
}
#endif
ExecuteFlip(*this->_pos.F(),this->_pos.E(),&m);
UpdateTopologies(&m);
///stars optimization
/*int t0=clock();*/
/*OptimizeStar<BaseMesh>(v0);
OptimizeStar<BaseMesh>(v1);
OptimizeStar<BaseMesh>(v2);
OptimizeStar<BaseMesh>(v3);*/
SmartOptimizeStar<BaseMesh>(v0,m,Accuracy(),EType());
SmartOptimizeStar<BaseMesh>(v1,m,Accuracy(),EType());
SmartOptimizeStar<BaseMesh>(v2,m,Accuracy(),EType());
SmartOptimizeStar<BaseMesh>(v3,m,Accuracy(),EType());
/*int t1=clock();
time_opt+=(t1-t0);*/
}
ScalarType ComputePriority()
{
this->_priority=EdgeDiff();
return this->_priority;
}
BaseFace *getF()
{return this->_pos.F();}
int getE()
{return this->_pos.E();}
public:
static int &Accuracy()
{
static int _acc;
return _acc;
}
};
#endif

368
src/meshlabplugins/filter_isoparametrization/param_mesh.h
View File

@ -1,184 +1,184 @@
#ifndef PARAM_MESH
#define PARAM_MESH
// stuff to define the mesh
#include <vcg/simplex/face/component_rt.h>
#include <vcg/complex/trimesh/base.h>
#include "iso_parametrization.h"
class BaseVertex;
class BaseEdge;
class BaseFace;
class BaseUsedTypes: public vcg::UsedTypes < vcg::Use<BaseVertex>::AsVertexType,
vcg::Use<BaseEdge >::AsEdgeType,
vcg::Use<BaseFace >::AsFaceType >{};
///AUXILIARY STRUCTURES USED FOR PARAMETRIZATION COMPUTATION
class BaseVertex : public vcg::Vertex< BaseUsedTypes,
vcg::vertex::VFAdj,
vcg::vertex::Coord3f,
vcg::vertex::Normal3f,
vcg::vertex::Mark,
vcg::vertex::BitFlags,
vcg::vertex::Color4b,
vcg::vertex::TexCoord2f>
{
public:
ScalarType area;
CoordType RPos;
BaseVertex *brother;
BaseFace *father;
CoordType Bary;
//ScalarType Damp;
BaseVertex()
{
brother=NULL;
//Damp=1;
//num_collapse=1;
}
static const bool Has_Auxiliary(){return true;}
vcg::Point2<ScalarType> RestUV;
vcg::Color4b OriginalCol;
//int num_collapse;
void ImportData(const BaseVertex & left )
{
vcg::Vertex< BaseUsedTypes,
vcg::vertex::VFAdj,
vcg::vertex::Coord3f,
vcg::vertex::Normal3f,
vcg::vertex::Mark,
vcg::vertex::BitFlags,
vcg::vertex::Color4b,
vcg::vertex::TexCoord2f>::ImportData(left);
this->area=left.area;
this->RPos=left.RPos;
this->brother=left.brother;
this->father=left.father;
this->Bary=left.Bary;
}
template < class LeftV>
void ImportData(const LeftV & left )
{
vcg::Vertex< BaseUsedTypes,
vcg::vertex::VFAdj,
vcg::vertex::Coord3f,
vcg::vertex::Normal3f,
vcg::vertex::Mark,
vcg::vertex::BitFlags,
vcg::vertex::Color4b,
vcg::vertex::TexCoord2f>::ImportData(left);
}
};
///class maintaing additional auxiliary data used during the parameterization
class AuxiliaryVertData{
public:
BaseVertex::ScalarType area;
BaseVertex::CoordType RPos;
BaseVertex *brother;
BaseFace *father;
BaseVertex::CoordType Bary;
vcg::Point2<BaseVertex::ScalarType> RestUV;
vcg::Color4b OriginalCol;
AuxiliaryVertData()
{brother=NULL;}
};
class AuxiliaryFaceData
{
public:
std::vector<std::pair<BaseVertex*,vcg::Point3f> > vertices_bary;
typedef std::vector<std::pair<BaseVertex*,vcg::Point3f> >::iterator IteVBary;
vcg::Color4b group;
BaseVertex::ScalarType areadelta;
vcg::Color4b colorDivision;
};
class BaseFace;
class BaseEdge : public vcg::Edge<BaseUsedTypes,vcg::edge::EVAdj> {
public:
inline BaseEdge() {};
inline BaseEdge( BaseVertex * v0, BaseVertex * v1)
{
V(0)=v0;
V(1)=v1;
}
//EdgeBase<BaseVertex,BaseEdge,BaseFace>(v0,v1){};
static inline BaseEdge OrderedEdge(BaseVertex* v0,BaseVertex* v1){
if(v0<v1) return BaseEdge(v0,v1);
else return BaseEdge(v1,v0);
}
};
class BaseFace : public vcg::Face < BaseUsedTypes,
vcg::face::VFAdj,
vcg::face::FFAdj,
vcg::face::VertexRef,
vcg::face::BitFlags,
vcg::face::EdgePlane,
vcg::face::Mark,
vcg::face::Normal3f,
vcg::face::Color4b>
{
public:
std::vector<std::pair<BaseVertex*,vcg::Point3f> > vertices_bary;
typedef std::vector<std::pair<BaseVertex*,vcg::Point3f> >::iterator IteVBary;
vcg::Color4b group;
ScalarType areadelta;
vcg::Color4b colorDivision;
template < class LeftV>
void ImportData(const LeftV & left )
{
vcg::Face < BaseUsedTypes,
vcg::face::VFAdj,
vcg::face::FFAdj,
vcg::face::VertexRef,
vcg::face::BitFlags,
vcg::face::EdgePlane,
vcg::face::Mark,
vcg::face::Normal3f,
vcg::face::Color4b>::ImportData(left);
}
void ImportData(const BaseFace & left )
{
vcg::Face < BaseUsedTypes,
vcg::face::VFAdj,
vcg::face::FFAdj,
vcg::face::VertexRef,
vcg::face::BitFlags,
vcg::face::EdgePlane,
vcg::face::Mark,
vcg::face::Normal3f,
vcg::face::Color4b>::ImportData(left);
this->vertices_bary = std::vector<std::pair<BaseVertex*,vcg::Point3f> > (left.vertices_bary);
this->group=left.group;
this->areadelta=left.areadelta;
this->colorDivision=left.colorDivision;
}
};
/// the main mesh class
class BaseMesh: public vcg::tri::TriMesh<std::vector<BaseVertex>, std::vector<BaseFace> > {};
#endif
#ifndef PARAM_MESH
#define PARAM_MESH
// stuff to define the mesh
#include <vcg/simplex/face/component_rt.h>
#include <vcg/complex/complex.h>
#include "iso_parametrization.h"
class BaseVertex;
class BaseEdge;
class BaseFace;
class BaseUsedTypes: public vcg::UsedTypes < vcg::Use<BaseVertex>::AsVertexType,
vcg::Use<BaseEdge >::AsEdgeType,
vcg::Use<BaseFace >::AsFaceType >{};
///AUXILIARY STRUCTURES USED FOR PARAMETRIZATION COMPUTATION
class BaseVertex : public vcg::Vertex< BaseUsedTypes,
vcg::vertex::VFAdj,
vcg::vertex::Coord3f,
vcg::vertex::Normal3f,
vcg::vertex::Mark,
vcg::vertex::BitFlags,
vcg::vertex::Color4b,
vcg::vertex::TexCoord2f>
{
public:
ScalarType area;
CoordType RPos;
BaseVertex *brother;
BaseFace *father;
CoordType Bary;
//ScalarType Damp;
BaseVertex()
{
brother=NULL;
//Damp=1;
//num_collapse=1;
}
static const bool Has_Auxiliary(){return true;}
vcg::Point2<ScalarType> RestUV;
vcg::Color4b OriginalCol;
//int num_collapse;
void ImportData(const BaseVertex & left )
{
vcg::Vertex< BaseUsedTypes,
vcg::vertex::VFAdj,
vcg::vertex::Coord3f,
vcg::vertex::Normal3f,
vcg::vertex::Mark,
vcg::vertex::BitFlags,
vcg::vertex::Color4b,
vcg::vertex::TexCoord2f>::ImportData(left);
this->area=left.area;
this->RPos=left.RPos;
this->brother=left.brother;
this->father=left.father;
this->Bary=left.Bary;
}
template < class LeftV>
void ImportData(const LeftV & left )
{
vcg::Vertex< BaseUsedTypes,
vcg::vertex::VFAdj,
vcg::vertex::Coord3f,
vcg::vertex::Normal3f,
vcg::vertex::Mark,
vcg::vertex::BitFlags,
vcg::vertex::Color4b,
vcg::vertex::TexCoord2f>::ImportData(left);
}
};
///class maintaing additional auxiliary data used during the parameterization
class AuxiliaryVertData{
public:
BaseVertex::ScalarType area;
BaseVertex::CoordType RPos;
BaseVertex *brother;
BaseFace *father;
BaseVertex::CoordType Bary;
vcg::Point2<BaseVertex::ScalarType> RestUV;
vcg::Color4b OriginalCol;
AuxiliaryVertData()
{brother=NULL;}
};
class AuxiliaryFaceData
{
public:
std::vector<std::pair<BaseVertex*,vcg::Point3f> > vertices_bary;
typedef std::vector<std::pair<BaseVertex*,vcg::Point3f> >::iterator IteVBary;
vcg::Color4b group;
BaseVertex::ScalarType areadelta;
vcg::Color4b colorDivision;
};
class BaseFace;
class BaseEdge : public vcg::Edge<BaseUsedTypes,vcg::edge::EVAdj> {
public:
inline BaseEdge() {};
inline BaseEdge( BaseVertex * v0, BaseVertex * v1)
{
V(0)=v0;
V(1)=v1;
}
//EdgeBase<BaseVertex,BaseEdge,BaseFace>(v0,v1){};
static inline BaseEdge OrderedEdge(BaseVertex* v0,BaseVertex* v1){
if(v0<v1) return BaseEdge(v0,v1);
else return BaseEdge(v1,v0);
}
};
class BaseFace : public vcg::Face < BaseUsedTypes,
vcg::face::VFAdj,
vcg::face::FFAdj,
vcg::face::VertexRef,
vcg::face::BitFlags,
vcg::face::EdgePlane,
vcg::face::Mark,
vcg::face::Normal3f,
vcg::face::Color4b>
{
public:
std::vector<std::pair<BaseVertex*,vcg::Point3f> > vertices_bary;
typedef std::vector<std::pair<BaseVertex*,vcg::Point3f> >::iterator IteVBary;
vcg::Color4b group;
ScalarType areadelta;
vcg::Color4b colorDivision;
template < class LeftV>
void ImportData(const LeftV & left )
{
vcg::Face < BaseUsedTypes,
vcg::face::VFAdj,
vcg::face::FFAdj,
vcg::face::VertexRef,
vcg::face::BitFlags,
vcg::face::EdgePlane,
vcg::face::Mark,
vcg::face::Normal3f,
vcg::face::Color4b>::ImportData(left);
}
void ImportData(const BaseFace & left )
{
vcg::Face < BaseUsedTypes,
vcg::face::VFAdj,
vcg::face::FFAdj,
vcg::face::VertexRef,
vcg::face::BitFlags,
vcg::face::EdgePlane,
vcg::face::Mark,
vcg::face::Normal3f,
vcg::face::Color4b>::ImportData(left);
this->vertices_bary = std::vector<std::pair<BaseVertex*,vcg::Point3f> > (left.vertices_bary);
this->group=left.group;
this->areadelta=left.areadelta;
this->colorDivision=left.colorDivision;
}
};
/// the main mesh class
class BaseMesh: public vcg::tri::TriMesh<std::vector<BaseVertex>, std::vector<BaseFace> > {};
#endif

2554
src/meshlabplugins/filter_isoparametrization/parametrizator.h
View File

File diff suppressed because it is too large Load Diff

686
src/meshlabplugins/filter_isoparametrization/stat_remeshing.h
View File

@ -1,343 +1,343 @@
#include <vcg/complex/trimesh/update/topology.h>
#include <vcg/complex/trimesh/update/flag.h>
#include <vcg/space/triangle3.h>
#include <wrap/io_trimesh/import_ply.h>
#include <wrap/io_trimesh/export_ply.h>
#include <vcg/simplex/face/pos.h>
#include <vcg/math/histogram.h>
#ifndef STAT_REMESHING
#define STAT_REMESHING
#define PI 3.14159265
//return moltiplication of aspect ratio of faces of the mesh
///to see if add + normalize is better
template <class MeshType>
typename MeshType::ScalarType MinimumAspectRatio(const MeshType &mesh)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
ScalarType res=1.f;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
ScalarType test=vcg::QualityRadii((*Fi).P(0),(*Fi).P(1),(*Fi).P(2));
if (test<res)
res=test;
}
return (res);
}
//return moltiplication of aspect ratio of faces of the mesh
///to see if add + normalize is better
template <class MeshType>
typename MeshType::ScalarType MinimumArea(const MeshType &mesh)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
ScalarType res=10000.f;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
ScalarType test=vcg::DoubleArea(*Fi)/2.0;
if (test<res)
res=test;
}
return (res);
}
//return moltiplication of aspect ratio of faces of the mesh
///to see if add + normalize is better
template <class MeshType>
typename MeshType::ScalarType MaximumArea(const MeshType &mesh)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
ScalarType res=0.f;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
ScalarType test=vcg::DoubleArea(*Fi)/2.0;
if (test>res)
res=test;
}
return (res);
}
//return moltiplication of aspect ratio of faces of the mesh
///to see if add + normalize is better
template <class FaceType>
typename FaceType::ScalarType MinAngleFace(const FaceType &f)
{
typedef typename FaceType::ScalarType ScalarType;
typedef typename FaceType::CoordType CoordType;
ScalarType res=360.0;
for (int i=0;i<3;i++)
{
CoordType v0=f.P((i+1)%3)-f.P(i);
CoordType v1=f.P((i+2)%3)-f.P(i);
v0.Normalize();
v1.Normalize();
ScalarType angle=acos(v0*v1)* 180.0 / PI;
if (angle<res)
res=angle;
}
assert(res<=60.0);
return (res);
}
//return moltiplication of aspect ratio of faces of the mesh
///to see if add + normalize is better
template <class FaceType>
typename FaceType::ScalarType MaxAngleFace(const FaceType &f)
{
typedef typename FaceType::ScalarType ScalarType;
typedef typename FaceType::CoordType CoordType;
ScalarType res=0;
for (int i=0;i<3;i++)
{
CoordType v0=f.P((i+1)%3)-f.P(i);
CoordType v1=f.P((i+2)%3)-f.P(i);
v0.Normalize();
v1.Normalize();
ScalarType angle=acos(v0*v1)* 180.0 / PI;
if (angle>res)
res=angle;
}
return (res);
}
template <class MeshType>
typename MeshType::ScalarType MinAngle(const MeshType &mesh)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
ScalarType res=360.0;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
ScalarType testAngle=MinAngleFace(*Fi);
if (testAngle<res)
res=testAngle;
}
return (res);
}
template <class MeshType>
typename MeshType::ScalarType MaxAngle(const MeshType &mesh)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
ScalarType res=0;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
ScalarType testAngle=MaxAngleFace(*Fi);
if (testAngle>res)
res=testAngle;
}
return (res);
}
template <class MeshType>
void MaxMinEdge(const MeshType &mesh,typename MeshType::ScalarType &min,
typename MeshType::ScalarType &max)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
min=10000.0;
max=0.0;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
for (int i=0;i<3;i++)
{
typename MeshType::VertexType* v0=(*Fi).V(i);
typename MeshType::VertexType* v1=(*Fi).V((i+1)%3);
if (v0>v1)
{
ScalarType dist=(v0->P()-v1->P()).Norm();
if (dist<min)
min=dist;
if (dist>max)
max=dist;
}
}
}
}
////return moltiplication of aspect ratio of faces of the mesh
/////to see if add + normalize is better
//template <class MeshType>
//ScalarType AverageMinAngle(const MeshType &mesh)
//{
// ScalarType res=0;
// MeshType::ConstFaceIterator Fi;
// for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
// if ((!(*Fi).IsD()))
// {
// ScalarType testAngle=MinAngleFace(*Fi);
// res+=testAngle;
// }
// return (res/((ScalarType)mesh.fn));
//}
template <class MeshType>
void StatAspectRatio(MeshType &mesh,
typename MeshType::ScalarType &min,
typename MeshType::ScalarType &average,
typename MeshType::ScalarType &stand_dev)
{
typedef typename MeshType::ScalarType ScalarType;
vcg::Histogram<ScalarType> HAspectRatio;
ScalarType minRatio=MinimumAspectRatio(mesh);
ScalarType maxRatio=1.f;
HAspectRatio.SetRange(minRatio,maxRatio,500);
typename MeshType::FaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
HAspectRatio.Add(vcg::QualityRadii((*Fi).P(0),(*Fi).P(1),(*Fi).P(2)));
average=HAspectRatio.Avg();
stand_dev=HAspectRatio.StandardDeviation();
min=minRatio;
}
template <class MeshType>
void StatArea(MeshType &mesh,
typename MeshType::ScalarType &min,
typename MeshType::ScalarType &max,
typename MeshType::ScalarType &average,
typename MeshType::ScalarType &stand_dev)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
vcg::Histogram<ScalarType> HArea;
ScalarType minArea=MinimumArea(mesh);
ScalarType maxArea=MaximumArea(mesh);
HArea.SetRange(minArea,maxArea,500);
typename MeshType::FaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
{
CoordType p0=(*Fi).P0(0);
CoordType p1=(*Fi).P1(0);
CoordType p2=(*Fi).P2(0);
ScalarType area=((p1-p0)^(p2-p0)).Norm()/2.0;
HArea.Add(area);
}
average=HArea.Avg();
stand_dev=HArea.StandardDeviation();
min=minArea;
max=maxArea;
}
template <class MeshType>
void StatAngle(MeshType &mesh,
typename MeshType::ScalarType &min,
typename MeshType::ScalarType &max,
typename MeshType::ScalarType &average,
typename MeshType::ScalarType &stand_dev)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexType VertexType;
vcg::Histogram<ScalarType> HAngle;
ScalarType minAngle=MinAngle(mesh);
ScalarType maxAngle=MaxAngle(mesh);
HAngle.SetRange(minAngle,maxAngle,500);
typename MeshType::FaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
HAngle.Add(MinAngleFace((*Fi)));
average=HAngle.Avg();
stand_dev=HAngle.StandardDeviation();
min=minAngle;
max=maxAngle;
}
template <class MeshType>
void StatEdge(MeshType &mesh,
typename MeshType::ScalarType &min,
typename MeshType::ScalarType &max,
typename MeshType::ScalarType &average,
typename MeshType::ScalarType &stand_dev)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexType VertexType;
vcg::Histogram<ScalarType> HEdge;
ScalarType minEdge;
ScalarType maxEdge;
MaxMinEdge(mesh,minEdge,maxEdge);
HEdge.SetRange(minEdge,maxEdge,500);
typedef typename MeshType::ScalarType ScalarType;
typename MeshType::FaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
{
for (int i=0;i<3;i++)
{
VertexType* v0=(*Fi).V(i);
VertexType* v1=(*Fi).V((i+1)%3);
if ((v0>v1)||((*Fi).FFp(i)==&(*Fi)))
{
ScalarType dist=(v0->P()-v1->P()).Norm();
HEdge.Add(dist);
}
}
}
average=HEdge.Avg();
stand_dev=HEdge.StandardDeviation();
min=minEdge;
max=maxEdge;
}
template <class MeshType>
int NumRegular(MeshType &mesh)
{
typedef typename MeshType::FaceType FaceType;
///update topology
vcg::tri::UpdateTopology<MeshType>::VertexFace(mesh);
typename MeshType::VertexIterator Vi;
int irregular=0;
for (Vi=mesh.vert.begin();Vi!=mesh.vert.end();Vi++)
{
if ((!(*Vi).IsD())&&(!(*Vi).IsB()))
{
vcg::face::VFIterator<FaceType> VFI(&(*Vi));
int num=0;
while (!(VFI.End()))
{
num++;
++VFI;
}
if (num!=6)
irregular++;
}
}
return irregular;
}
#endif
#include <vcg/complex/algorithms/update/topology.h>
#include <vcg/complex/algorithms/update/flag.h>
#include <vcg/space/triangle3.h>
#include <wrap/io_trimesh/import_ply.h>
#include <wrap/io_trimesh/export_ply.h>
#include <vcg/simplex/face/pos.h>
#include <vcg/math/histogram.h>
#ifndef STAT_REMESHING
#define STAT_REMESHING
#define PI 3.14159265
//return moltiplication of aspect ratio of faces of the mesh
///to see if add + normalize is better
template <class MeshType>
typename MeshType::ScalarType MinimumAspectRatio(const MeshType &mesh)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
ScalarType res=1.f;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
ScalarType test=vcg::QualityRadii((*Fi).P(0),(*Fi).P(1),(*Fi).P(2));
if (test<res)
res=test;
}
return (res);
}
//return moltiplication of aspect ratio of faces of the mesh
///to see if add + normalize is better
template <class MeshType>
typename MeshType::ScalarType MinimumArea(const MeshType &mesh)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
ScalarType res=10000.f;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
ScalarType test=vcg::DoubleArea(*Fi)/2.0;
if (test<res)
res=test;
}
return (res);
}
//return moltiplication of aspect ratio of faces of the mesh
///to see if add + normalize is better
template <class MeshType>
typename MeshType::ScalarType MaximumArea(const MeshType &mesh)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
ScalarType res=0.f;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
ScalarType test=vcg::DoubleArea(*Fi)/2.0;
if (test>res)
res=test;
}
return (res);
}
//return moltiplication of aspect ratio of faces of the mesh
///to see if add + normalize is better
template <class FaceType>
typename FaceType::ScalarType MinAngleFace(const FaceType &f)
{
typedef typename FaceType::ScalarType ScalarType;
typedef typename FaceType::CoordType CoordType;
ScalarType res=360.0;
for (int i=0;i<3;i++)
{
CoordType v0=f.P((i+1)%3)-f.P(i);
CoordType v1=f.P((i+2)%3)-f.P(i);
v0.Normalize();
v1.Normalize();
ScalarType angle=acos(v0*v1)* 180.0 / PI;
if (angle<res)
res=angle;
}
assert(res<=60.0);
return (res);
}
//return moltiplication of aspect ratio of faces of the mesh
///to see if add + normalize is better
template <class FaceType>
typename FaceType::ScalarType MaxAngleFace(const FaceType &f)
{
typedef typename FaceType::ScalarType ScalarType;
typedef typename FaceType::CoordType CoordType;
ScalarType res=0;
for (int i=0;i<3;i++)
{
CoordType v0=f.P((i+1)%3)-f.P(i);
CoordType v1=f.P((i+2)%3)-f.P(i);
v0.Normalize();
v1.Normalize();
ScalarType angle=acos(v0*v1)* 180.0 / PI;
if (angle>res)
res=angle;
}
return (res);
}
template <class MeshType>
typename MeshType::ScalarType MinAngle(const MeshType &mesh)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
ScalarType res=360.0;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
ScalarType testAngle=MinAngleFace(*Fi);
if (testAngle<res)
res=testAngle;
}
return (res);
}
template <class MeshType>
typename MeshType::ScalarType MaxAngle(const MeshType &mesh)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
ScalarType res=0;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
ScalarType testAngle=MaxAngleFace(*Fi);
if (testAngle>res)
res=testAngle;
}
return (res);
}
template <class MeshType>
void MaxMinEdge(const MeshType &mesh,typename MeshType::ScalarType &min,
typename MeshType::ScalarType &max)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
min=10000.0;
max=0.0;
typename MeshType::ConstFaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
if ((!(*Fi).IsD()))
{
for (int i=0;i<3;i++)
{
typename MeshType::VertexType* v0=(*Fi).V(i);
typename MeshType::VertexType* v1=(*Fi).V((i+1)%3);
if (v0>v1)
{
ScalarType dist=(v0->P()-v1->P()).Norm();
if (dist<min)
min=dist;
if (dist>max)
max=dist;
}
}
}
}
////return moltiplication of aspect ratio of faces of the mesh
/////to see if add + normalize is better
//template <class MeshType>
//ScalarType AverageMinAngle(const MeshType &mesh)
//{
// ScalarType res=0;
// MeshType::ConstFaceIterator Fi;
// for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
// if ((!(*Fi).IsD()))
// {
// ScalarType testAngle=MinAngleFace(*Fi);
// res+=testAngle;
// }
// return (res/((ScalarType)mesh.fn));
//}
template <class MeshType>
void StatAspectRatio(MeshType &mesh,
typename MeshType::ScalarType &min,
typename MeshType::ScalarType &average,
typename MeshType::ScalarType &stand_dev)
{
typedef typename MeshType::ScalarType ScalarType;
vcg::Histogram<ScalarType> HAspectRatio;
ScalarType minRatio=MinimumAspectRatio(mesh);
ScalarType maxRatio=1.f;
HAspectRatio.SetRange(minRatio,maxRatio,500);
typename MeshType::FaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
HAspectRatio.Add(vcg::QualityRadii((*Fi).P(0),(*Fi).P(1),(*Fi).P(2)));
average=HAspectRatio.Avg();
stand_dev=HAspectRatio.StandardDeviation();
min=minRatio;
}
template <class MeshType>
void StatArea(MeshType &mesh,
typename MeshType::ScalarType &min,
typename MeshType::ScalarType &max,
typename MeshType::ScalarType &average,
typename MeshType::ScalarType &stand_dev)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::CoordType CoordType;
vcg::Histogram<ScalarType> HArea;
ScalarType minArea=MinimumArea(mesh);
ScalarType maxArea=MaximumArea(mesh);
HArea.SetRange(minArea,maxArea,500);
typename MeshType::FaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
{
CoordType p0=(*Fi).P0(0);
CoordType p1=(*Fi).P1(0);
CoordType p2=(*Fi).P2(0);
ScalarType area=((p1-p0)^(p2-p0)).Norm()/2.0;
HArea.Add(area);
}
average=HArea.Avg();
stand_dev=HArea.StandardDeviation();
min=minArea;
max=maxArea;
}
template <class MeshType>
void StatAngle(MeshType &mesh,
typename MeshType::ScalarType &min,
typename MeshType::ScalarType &max,
typename MeshType::ScalarType &average,
typename MeshType::ScalarType &stand_dev)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexType VertexType;
vcg::Histogram<ScalarType> HAngle;
ScalarType minAngle=MinAngle(mesh);
ScalarType maxAngle=MaxAngle(mesh);
HAngle.SetRange(minAngle,maxAngle,500);
typename MeshType::FaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
HAngle.Add(MinAngleFace((*Fi)));
average=HAngle.Avg();
stand_dev=HAngle.StandardDeviation();
min=minAngle;
max=maxAngle;
}
template <class MeshType>
void StatEdge(MeshType &mesh,
typename MeshType::ScalarType &min,
typename MeshType::ScalarType &max,
typename MeshType::ScalarType &average,
typename MeshType::ScalarType &stand_dev)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexType VertexType;
vcg::Histogram<ScalarType> HEdge;
ScalarType minEdge;
ScalarType maxEdge;
MaxMinEdge(mesh,minEdge,maxEdge);
HEdge.SetRange(minEdge,maxEdge,500);
typedef typename MeshType::ScalarType ScalarType;
typename MeshType::FaceIterator Fi;
for (Fi=mesh.face.begin();Fi!=mesh.face.end();Fi++)
{
for (int i=0;i<3;i++)
{
VertexType* v0=(*Fi).V(i);
VertexType* v1=(*Fi).V((i+1)%3);
if ((v0>v1)||((*Fi).FFp(i)==&(*Fi)))
{
ScalarType dist=(v0->P()-v1->P()).Norm();
HEdge.Add(dist);
}
}
}
average=HEdge.Avg();
stand_dev=HEdge.StandardDeviation();
min=minEdge;
max=maxEdge;
}
template <class MeshType>
int NumRegular(MeshType &mesh)
{
typedef typename MeshType::FaceType FaceType;
///update topology
vcg::tri::UpdateTopology<MeshType>::VertexFace(mesh);
typename MeshType::VertexIterator Vi;
int irregular=0;
for (Vi=mesh.vert.begin();Vi!=mesh.vert.end();Vi++)
{
if ((!(*Vi).IsD())&&(!(*Vi).IsB()))
{
vcg::face::VFIterator<FaceType> VFI(&(*Vi));
int num=0;
while (!(VFI.End()))
{
num++;
++VFI;
}
if (num!=6)
irregular++;
}
}
return irregular;
}
#endif

996
src/meshlabplugins/filter_isoparametrization/tangent_space.h
View File

@ -1,498 +1,498 @@
#ifndef _ISO_TANGENTSPACE
#define _ISO_TANGENTSPACE
#include "iso_parametrization.h"
#include <vcg/complex/trimesh/update/curvature.h>
#include <vcg/complex/trimesh/update/normal.h>
class TangentSpace{
typedef IsoParametrization::CoordType CoordType;
typedef IsoParametrization::ScalarType ScalarType;
vcg::SimpleTempData<typename ParamMesh::VertContainer,vcg::Matrix33<ScalarType> > *ProjMatrix;
public:
IsoParametrization *isoParam;
bool isoParamTheta(int i, vcg::Point2<ScalarType> p, vcg::Point3<ScalarType> &res) const{
return isoParam->Theta(i,p,res);
}
//
void Theta(int i,
const vcg::Point2<ScalarType> &UV,
CoordType &pos3D){
isoParamTheta(i,UV,pos3D);
}
bool Theta(const int &I,
const vcg::Point2<ScalarType> &alpha_beta, // alphaBeta
std::vector<ParamFace*> &face,
std::vector<CoordType> &baryVal)
{
int ret=isoParam->Theta(I,alpha_beta,face,baryVal);
return (ret!=-1);
}
///initialize the sampler
void Init(IsoParametrization *_isoParam,ScalarType radius=(ScalarType)0.1)
{
isoParam=_isoParam;
ProjMatrix = new vcg::SimpleTempData<typename ParamMesh::VertContainer,vcg::Matrix33<ScalarType> > (isoParam->ParaMesh()->vert);
InitProjectionMatrix(radius);
vcg::tri::UpdateNormals<ParamMesh>::PerFaceNormalized(*isoParam->ParaMesh());
vcg::tri::UpdateCurvature<ParamMesh>::PrincipalDirectionsNormalCycles(*isoParam->ParaMesh());
}
///given an initial position in parametric space (I0,bary0)
///and a 2D vector (vect) expressed in parametric space modify the final
///position (I1,bary1) abd return true if everithing was ok, false otherwise
bool Sum(const int &I0,const vcg::Point2<ScalarType> &bary0,
const vcg::Point2<ScalarType> &vect,
int &I1,vcg::Point2<ScalarType> &bary1,int &domain) const
{
vcg::Point2<ScalarType> dest=bary0+vect;
//vect[0]*Xaxis + vect[1]*Yaxis;
ScalarType alpha=dest.X();
ScalarType beta=dest.Y();
///point inside the face
if ((alpha>=0)&&(alpha<=1)&&(beta>=0)&&(beta<=1)&&((alpha+beta)<=1))
{
bary1=dest;
I1=I0;
domain=0;
return true;
}
///control edges
int edge=-1;
if ((alpha<=1)&&(beta<=1)&&((alpha+beta)>=1))
edge=0;
else
if ((alpha<=0)&&(beta<=1)&&((alpha+beta)>=0))
edge=1;
else
if ((alpha<=1)&&(beta<=0)&&((alpha+beta)>=0))
edge=2;
if (edge!=-1)
{
int DiamIndex=isoParam->GetDiamond(I0,edge);
vcg::Point2<ScalarType> UVDiam;
///transform to diamond coordinates
isoParam->GE1(I0,dest,DiamIndex,UVDiam);
///trasform back to I,alpha,beta
isoParam->inv_GE1(DiamIndex,UVDiam,I1,bary1);
domain=1;
return true;
}
int star=-1;
ScalarType gamma=(1-alpha-beta);
if ((alpha>beta)&&(alpha>gamma))
star=0;
else
if ((beta>alpha)&&(beta>gamma))
star=1;
else
star=2;
///get the index of star
int StarIndex=isoParam->GetStarIndex(I0,star);
vcg::Point2<ScalarType> UVHstar;
///transform to UV
bool found=isoParam->GE0(I0,dest,StarIndex,UVHstar);
///trasform back to I,alpha,beta
if (!found)
return false;
found=isoParam->inv_GE0(StarIndex,UVHstar,I1,bary1);
/*AbstractFace* f0=&isoParam->AbsMesh()->face[I0];
AbstractFace* f1=&isoParam->AbsMesh()->face[I1];
AbstractVertex *v0=f0->V(0);
AbstractVertex *v1=f0->V(1);
AbstractVertex *v2=f0->V(2);
AbstractVertex *v3=f1->V(0);
AbstractVertex *v4=f1->V(1);
AbstractVertex *v5=f1->V(2);*/
domain=2;
return found;
}
///given two positions in parametric space (I0,bary0) and (I1,bary1)
///modify the 2D vector (vect) and return true if everithing was ok, false otherwise
bool Sub(const int &I0,const vcg::Point2<ScalarType> &bary0,
const int &I1,const vcg::Point2<ScalarType> &bary1,
vcg::Point2<ScalarType> &vect,int &num) const
{
int IndexDomain;
num=isoParam->InterpolationSpace(I0,I1,IndexDomain);
///is a face
if (num==0)
{
//printf("F");
assert(I0==I1);
vect=bary0-bary1;
return true;
}
else
///a diamond
if (num==1)
{
//printf("D");
///tranform in UV space
vcg::Point2<ScalarType> UVDiam;
isoParam->GE1(I1,bary1,IndexDomain,UVDiam);
///then find bary coords wich respect to the first face
vcg::Point2<ScalarType> bary2;
isoParam->inv_GE1_fixedI(IndexDomain,UVDiam,I0,bary2);
vect=bary0-bary2;
return true;
}
else
///a star
if (num==2)
{
//printf("S");
///tranform in UV space
vcg::Point2<ScalarType> UVStar;
isoParam->GE0(I1,bary1,IndexDomain,UVStar);
///then find bary coords wich respect to the first face
vcg::Point2<ScalarType> bary2;
isoParam->inv_GE0_fixedI(IndexDomain,UVStar,I0,bary2);
vect=bary0-bary2;
return true;
}
else
return false;
}
bool TangentDir(const int &I,const vcg::Point2<ScalarType> &bary,
CoordType &XAxis,CoordType &YAxis,ScalarType radius=(ScalarType)0.5) const
{
///3d position of origin
//CoordType origin;
//isoParam->Theta(I,bary,origin);
// two axis in alpha-beta space that will be orthogonal in UV-Space
const ScalarType h=(ScalarType)2.0/sqrt((ScalarType)3.0);
vcg::Point2<ScalarType> XP=vcg::Point2<ScalarType>(1,0);
vcg::Point2<ScalarType> YP=vcg::Point2<ScalarType>(-0.5,h);
XP*=radius;
YP*=radius;
vcg::Point2<ScalarType> XM=-XP;
vcg::Point2<ScalarType> YM=-YP;
//Yaxis.Normalize();
vcg::Point2<ScalarType> bary0,bary1,bary2,bary3;
int I0,I1,I2,I3;
CoordType X0,X1,Y0,Y1;
///find paraCoords of four neighbors
int domain;
bool done=true;
done&=Sum(I,bary,XP,I0,bary0,domain);
done&=Sum(I,bary,XM,I1,bary1,domain);
done&=Sum(I,bary,YP,I2,bary2,domain);
done&=Sum(I,bary,YM,I3,bary3,domain);
if (!done)
return false;
//if (I0!=I1 || I1!=I2 || I2!=I3) return false;
///get 3d position
isoParamTheta(I0,bary0,X0);
isoParamTheta(I1,bary1,X1);
isoParamTheta(I2,bary2,Y0);
isoParamTheta(I3,bary3,Y1);
///average .. considering one is opposite respect to the other
XAxis=(X0-X1)/(ScalarType)2.0;
YAxis=(Y0-Y1)/(ScalarType)2.0;
///final scaling
XAxis/=radius;
YAxis/=radius;
return true;
}
void Test(int Sample=100,int Ite=100)
{
int max=isoParam->AbsMesh()->face.size();
for (int I=0;I<max;I++)
{
printf("\n TESTING %d \n",I);
for (int i=0;i<Sample;i++)
{
for (int j=0;j<Sample;j++)
{
if ((i+j)<Sample)
{
ScalarType alpha=(ScalarType)i/(ScalarType)Sample;
ScalarType beta=(ScalarType)j/(ScalarType)Sample;
assert((alpha+beta)<=1.0);
vcg::Point2<ScalarType> bary=vcg::Point2<ScalarType>(alpha,beta);
assert((bary.X()>=0)&&(bary.X()<=1));
assert((bary.Y()>=0)&&(bary.Y()<=1));
assert((bary.X()+bary.Y()<=1));
for (int k=0;k<Ite;k++)
{
ScalarType d0=((ScalarType)(rand()%1000));
ScalarType d1=((ScalarType)(rand()%1000));
vcg::Point2<ScalarType> vect=vcg::Point2<ScalarType>(d0,d1);
vect.Normalize();
ScalarType norm=(ScalarType)0.05;//((ScalarType)(rand()%1000))/(ScalarType)2000;
assert(norm<1);
vect*=norm;
int I1;
vcg::Point2<ScalarType> bary1;
vcg::Point2<ScalarType> vect1;
int domain0;
bool b1=Sum(I,bary,vect,I1,bary1,domain0);
assert(b1);
assert((bary1.X()>=0)&&(bary1.X()<=1));
assert((bary1.Y()>=0)&&(bary1.Y()<=1));
if(!((bary1.X()+bary1.Y())<=1.00001))
{
printf("\n SUM %.4f \n",bary1.X()+bary1.Y());
assert(0);
}
assert(I1<max);
int domain;
bool b2=Sub(I,bary,I1,bary1,vect1,domain);
assert(b2);
ScalarType diff=(vect1+vect).Norm();
if (domain0==0)
assert(domain==0);
if ((domain0==1)&&(domain==2))
assert(0);
/*if (domain0!=domain)
assert(0);*/
if (diff>0.001)
{
printf("\n DIFF %.4f domain SUM %d domain SUB %d \n",diff,domain0,domain);
//assert(0);
}
//assert(fabs(vect1.X()-vect.X())<0.0001);
//assert(fabs(vect1.Y()-vect.Y())<0.0001);
}
}
}
}
}
}
void InitProjectionMatrix(ScalarType radius=(ScalarType)0.1)
{
for (int i=0;i<isoParam->ParaMesh()->vert.size();i++)
{
int I=isoParam->ParaMesh()->vert[i].T().N();
vcg::Point2<ScalarType> bary=isoParam->ParaMesh()->vert[i].T().P();
CoordType origin;
isoParamTheta(I,bary,origin);
CoordType XAxis,YAxis; // tangent axis in 3D Object space
///get tangent directions
bool done=TangentDir(I,bary,XAxis,YAxis,(ScalarType)0.1);
if (!done)
{
(*ProjMatrix)[i].SetIdentity();
continue;
}
// must compute res2d such that: res2d.X() * XAxis + res2d.Y() * YAxis + dontCare * ZAxis = vect3d
CoordType ZAxis = -(XAxis^YAxis).Normalize()*XAxis.Norm();
(*ProjMatrix)[i].SetColumn(0,XAxis);
(*ProjMatrix)[i].SetColumn(1,YAxis);
(*ProjMatrix)[i].SetColumn(2,ZAxis);
vcg::Invert((*ProjMatrix)[i]);
}
}
void GetCurvature(const int &I,const vcg::Point2<ScalarType> &alpha_beta,
CoordType &d1,CoordType &d2,ScalarType &k1,ScalarType &k2)
{
ParamFace* face=NULL;
CoordType baryVal;
isoParam->Theta(I,alpha_beta,face,baryVal);
ParamVertex *v0=face->V(0);
ParamVertex *v1=face->V(1);
ParamVertex *v2=face->V(2);
d1=v0->PD1()*baryVal.X()+v1->PD1()*baryVal.Y()+v2->PD1()*baryVal.Z();
d2=v0->PD2()*baryVal.X()+v1->PD2()*baryVal.Y()+v2->PD2()*baryVal.Z();
k1=v0->K1()*baryVal.X()+v1->K1()*baryVal.Y()+v2->K1()*baryVal.Z();
k2=v0->K2()*baryVal.X()+v1->K2()*baryVal.Y()+v2->K2()*baryVal.Z();
}
void GetProjectionMatrix(const int &I,const vcg::Point2<ScalarType> &alpha_beta,
vcg::Matrix33<ScalarType> &projMatr) const
{
ParamFace* face=NULL;
CoordType baryVal;
int dom=isoParam->Theta(I,alpha_beta,face,baryVal);
if (dom==-1)
{
projMatr.SetIdentity();
return;
}
ParamVertex *v0=face->V(0);
ParamVertex *v1=face->V(1);
ParamVertex *v2=face->V(2);
projMatr=(*ProjMatrix)[v0]*baryVal.X();
projMatr+=(*ProjMatrix)[v1]*baryVal.Y();
projMatr+=(*ProjMatrix)[v2]*baryVal.Z();
}
bool Project(const int &I,const vcg::Point2<ScalarType> &bary,
const CoordType &vect3d, // in object space
vcg::Point2<ScalarType> &res2d) const
{
vcg::Matrix33f m;
GetProjectionMatrix(I,bary,m);
ScalarType deltaX=vect3d*m.GetRow(0);
ScalarType deltaY=vect3d*m.GetRow(1);
// two axis in alpha-beta space that will be orthogonal in UV-Space
const ScalarType h=(ScalarType)2.0/sqrt((ScalarType)3.0);
vcg::Point2<ScalarType> XP=vcg::Point2<ScalarType>(1,0);
vcg::Point2<ScalarType> YP=vcg::Point2<ScalarType>(-0.5,h);
res2d = XP*deltaX + YP*deltaY;
//res2d=vcg::Point2<ScalarType>(deltaX,deltaY);
return true;
}
// WEIGHTED INTERPOLATION OF POINTS IN TANGENT SPACE
/// weights MUST be normalized
bool Interpolate(const int &I0,const vcg::Point2<ScalarType> &alpha_beta0,
const int &I1,const vcg::Point2<ScalarType> &alpha_beta1,
ScalarType weight,
int &I_res,
vcg::Point2<ScalarType> &alpha_beta_res)
{
int IndexDomain;
int kind=isoParam->InterpolationSpace(I0,I1,IndexDomain);
if (kind==-1)
return false;
vcg::Point2<ScalarType> transformed0,transformed1;
///interpolate in a face
if (kind==0)
{
isoParam->GE2(IndexDomain,alpha_beta0,transformed0);
isoParam->GE2(IndexDomain,alpha_beta1,transformed1);
}
else
///interpolate in a diamond
if (kind==1)
{
isoParam->GE1(I0,alpha_beta0,IndexDomain,transformed0);
isoParam->GE1(I1,alpha_beta1,IndexDomain,transformed1);
}
else
{
isoParam->GE0(I0,alpha_beta0,IndexDomain,transformed0);
isoParam->GE0(I1,alpha_beta1,IndexDomain,transformed1);
}
vcg::Point2<ScalarType> UV_interp=transformed0*weight+transformed1*(1.0-weight);
///FINALLY......
///transform back to alpha,beta,I
if (kind==0)
{
isoParam->inv_GE2(IndexDomain,UV_interp,alpha_beta_res);
I_res=IndexDomain;
}
else
if (kind==1)
{
isoParam->inv_GE1(IndexDomain,UV_interp,I_res,alpha_beta_res);
}
else
isoParam->inv_GE0(IndexDomain,UV_interp,I_res,alpha_beta_res);
return true;
}
ScalarType AbstractArea()
{
ScalarType Cnum=sqrt(3.0)/4.0;
return(isoParam->AbsMesh()->fn*Cnum);
}
// WEIGHTED INTERPOLATION OF POINTS IN TANGENT SPACE
/// weights MUST be normalized
bool Interpolate(const std::vector<int> &I,
const std::vector<vcg::Point2<ScalarType> > &alpha_beta,
const std::vector<ScalarType> &weights,
int &I_res,
vcg::Point2<ScalarType> &alpha_beta_res,
int *num=NULL)
{
int size;
if (num==NULL)
size=alpha_beta.size();
else
size=*num;
int IndexDomain;
int kind=isoParam->InterpolationSpace(I,IndexDomain,num);
if (kind==-1)
return false;
std::vector<vcg::Point2<ScalarType> > transformed;
transformed.resize(size);
///interpolate in a face
if (kind==0)
{
for (int i=0;i<size;i++)
isoParam->GE2(IndexDomain,alpha_beta[i],transformed[i]);
}
else
///interpolate in a diamond
if (kind==1)
{
for (int i=0;i<size;i++)
isoParam->GE1(I[i],alpha_beta[i],IndexDomain,transformed[i]);
}
else
{
for (int i=0;i<size;i++)
bool b0=isoParam->GE0(I[i],alpha_beta[i],IndexDomain,transformed[i]);
}
/// do the interpolation
vcg::Point2<ScalarType> UV_interp=vcg::Point2<ScalarType>(0,0);
for (int i=0;i<size;i++)
UV_interp+=(transformed[i]*weights[i]);
///FINALLY......
///transform back to alpha,beta,I
if (kind==0)
{
isoParam->inv_GE2(IndexDomain,UV_interp,alpha_beta_res);
I_res=IndexDomain;
}
else
if (kind==1)
{
isoParam->inv_GE1(IndexDomain,UV_interp,I_res,alpha_beta_res);
}
else
isoParam->inv_GE0(IndexDomain,UV_interp,I_res,alpha_beta_res);
return true;
}
};
#endif
#ifndef _ISO_TANGENTSPACE
#define _ISO_TANGENTSPACE
#include "iso_parametrization.h"
#include <vcg/complex/algorithms/update/curvature.h>
#include <vcg/complex/algorithms/update/normal.h>
class TangentSpace{
typedef IsoParametrization::CoordType CoordType;
typedef IsoParametrization::ScalarType ScalarType;
vcg::SimpleTempData<typename ParamMesh::VertContainer,vcg::Matrix33<ScalarType> > *ProjMatrix;
public:
IsoParametrization *isoParam;
bool isoParamTheta(int i, vcg::Point2<ScalarType> p, vcg::Point3<ScalarType> &res) const{
return isoParam->Theta(i,p,res);
}
//
void Theta(int i,
const vcg::Point2<ScalarType> &UV,
CoordType &pos3D){
isoParamTheta(i,UV,pos3D);
}
bool Theta(const int &I,
const vcg::Point2<ScalarType> &alpha_beta, // alphaBeta
std::vector<ParamFace*> &face,
std::vector<CoordType> &baryVal)
{
int ret=isoParam->Theta(I,alpha_beta,face,baryVal);
return (ret!=-1);
}
///initialize the sampler
void Init(IsoParametrization *_isoParam,ScalarType radius=(ScalarType)0.1)
{
isoParam=_isoParam;
ProjMatrix = new vcg::SimpleTempData<typename ParamMesh::VertContainer,vcg::Matrix33<ScalarType> > (isoParam->ParaMesh()->vert);
InitProjectionMatrix(radius);
vcg::tri::UpdateNormals<ParamMesh>::PerFaceNormalized(*isoParam->ParaMesh());
vcg::tri::UpdateCurvature<ParamMesh>::PrincipalDirectionsNormalCycles(*isoParam->ParaMesh());
}
///given an initial position in parametric space (I0,bary0)
///and a 2D vector (vect) expressed in parametric space modify the final
///position (I1,bary1) abd return true if everithing was ok, false otherwise
bool Sum(const int &I0,const vcg::Point2<ScalarType> &bary0,
const vcg::Point2<ScalarType> &vect,
int &I1,vcg::Point2<ScalarType> &bary1,int &domain) const
{
vcg::Point2<ScalarType> dest=bary0+vect;
//vect[0]*Xaxis + vect[1]*Yaxis;
ScalarType alpha=dest.X();
ScalarType beta=dest.Y();
///point inside the face
if ((alpha>=0)&&(alpha<=1)&&(beta>=0)&&(beta<=1)&&((alpha+beta)<=1))
{
bary1=dest;
I1=I0;
domain=0;
return true;
}
///control edges
int edge=-1;
if ((alpha<=1)&&(beta<=1)&&((alpha+beta)>=1))
edge=0;
else
if ((alpha<=0)&&(beta<=1)&&((alpha+beta)>=0))
edge=1;
else
if ((alpha<=1)&&(beta<=0)&&((alpha+beta)>=0))
edge=2;
if (edge!=-1)
{
int DiamIndex=isoParam->GetDiamond(I0,edge);
vcg::Point2<ScalarType> UVDiam;
///transform to diamond coordinates
isoParam->GE1(I0,dest,DiamIndex,UVDiam);
///trasform back to I,alpha,beta
isoParam->inv_GE1(DiamIndex,UVDiam,I1,bary1);
domain=1;
return true;
}
int star=-1;
ScalarType gamma=(1-alpha-beta);
if ((alpha>beta)&&(alpha>gamma))
star=0;
else
if ((beta>alpha)&&(beta>gamma))
star=1;
else
star=2;
///get the index of star
int StarIndex=isoParam->GetStarIndex(I0,star);
vcg::Point2<ScalarType> UVHstar;
///transform to UV
bool found=isoParam->GE0(I0,dest,StarIndex,UVHstar);
///trasform back to I,alpha,beta
if (!found)
return false;
found=isoParam->inv_GE0(StarIndex,UVHstar,I1,bary1);
/*AbstractFace* f0=&isoParam->AbsMesh()->face[I0];
AbstractFace* f1=&isoParam->AbsMesh()->face[I1];
AbstractVertex *v0=f0->V(0);
AbstractVertex *v1=f0->V(1);
AbstractVertex *v2=f0->V(2);
AbstractVertex *v3=f1->V(0);
AbstractVertex *v4=f1->V(1);
AbstractVertex *v5=f1->V(2);*/
domain=2;
return found;
}
///given two positions in parametric space (I0,bary0) and (I1,bary1)
///modify the 2D vector (vect) and return true if everithing was ok, false otherwise
bool Sub(const int &I0,const vcg::Point2<ScalarType> &bary0,
const int &I1,const vcg::Point2<ScalarType> &bary1,
vcg::Point2<ScalarType> &vect,int &num) const
{
int IndexDomain;
num=isoParam->InterpolationSpace(I0,I1,IndexDomain);
///is a face
if (num==0)
{
//printf("F");
assert(I0==I1);
vect=bary0-bary1;
return true;
}
else
///a diamond
if (num==1)
{
//printf("D");
///tranform in UV space
vcg::Point2<ScalarType> UVDiam;
isoParam->GE1(I1,bary1,IndexDomain,UVDiam);
///then find bary coords wich respect to the first face
vcg::Point2<ScalarType> bary2;
isoParam->inv_GE1_fixedI(IndexDomain,UVDiam,I0,bary2);
vect=bary0-bary2;
return true;
}
else
///a star
if (num==2)
{
//printf("S");
///tranform in UV space
vcg::Point2<ScalarType> UVStar;
isoParam->GE0(I1,bary1,IndexDomain,UVStar);
///then find bary coords wich respect to the first face
vcg::Point2<ScalarType> bary2;
isoParam->inv_GE0_fixedI(IndexDomain,UVStar,I0,bary2);
vect=bary0-bary2;
return true;
}
else
return false;
}
bool TangentDir(const int &I,const vcg::Point2<ScalarType> &bary,
CoordType &XAxis,CoordType &YAxis,ScalarType radius=(ScalarType)0.5) const
{
///3d position of origin
//CoordType origin;
//isoParam->Theta(I,bary,origin);
// two axis in alpha-beta space that will be orthogonal in UV-Space
const ScalarType h=(ScalarType)2.0/sqrt((ScalarType)3.0);
vcg::Point2<ScalarType> XP=vcg::Point2<ScalarType>(1,0);
vcg::Point2<ScalarType> YP=vcg::Point2<ScalarType>(-0.5,h);
XP*=radius;
YP*=radius;
vcg::Point2<ScalarType> XM=-XP;
vcg::Point2<ScalarType> YM=-YP;
//Yaxis.Normalize();
vcg::Point2<ScalarType> bary0,bary1,bary2,bary3;
int I0,I1,I2,I3;
CoordType X0,X1,Y0,Y1;
///find paraCoords of four neighbors
int domain;
bool done=true;
done&=Sum(I,bary,XP,I0,bary0,domain);
done&=Sum(I,bary,XM,I1,bary1,domain);
done&=Sum(I,bary,YP,I2,bary2,domain);
done&=Sum(I,bary,YM,I3,bary3,domain);
if (!done)
return false;
//if (I0!=I1 || I1!=I2 || I2!=I3) return false;
///get 3d position
isoParamTheta(I0,bary0,X0);
isoParamTheta(I1,bary1,X1);
isoParamTheta(I2,bary2,Y0);
isoParamTheta(I3,bary3,Y1);
///average .. considering one is opposite respect to the other
XAxis=(X0-X1)/(ScalarType)2.0;
YAxis=(Y0-Y1)/(ScalarType)2.0;
///final scaling
XAxis/=radius;
YAxis/=radius;
return true;
}
void Test(int Sample=100,int Ite=100)
{
int max=isoParam->AbsMesh()->face.size();
for (int I=0;I<max;I++)
{
printf("\n TESTING %d \n",I);
for (int i=0;i<Sample;i++)
{
for (int j=0;j<Sample;j++)
{
if ((i+j)<Sample)
{
ScalarType alpha=(ScalarType)i/(ScalarType)Sample;
ScalarType beta=(ScalarType)j/(ScalarType)Sample;
assert((alpha+beta)<=1.0);
vcg::Point2<ScalarType> bary=vcg::Point2<ScalarType>(alpha,beta);
assert((bary.X()>=0)&&(bary.X()<=1));
assert((bary.Y()>=0)&&(bary.Y()<=1));
assert((bary.X()+bary.Y()<=1));
for (int k=0;k<Ite;k++)
{
ScalarType d0=((ScalarType)(rand()%1000));
ScalarType d1=((ScalarType)(rand()%1000));
vcg::Point2<ScalarType> vect=vcg::Point2<ScalarType>(d0,d1);
vect.Normalize();
ScalarType norm=(ScalarType)0.05;//((ScalarType)(rand()%1000))/(ScalarType)2000;
assert(norm<1);
vect*=norm;
int I1;
vcg::Point2<ScalarType> bary1;
vcg::Point2<ScalarType> vect1;
int domain0;
bool b1=Sum(I,bary,vect,I1,bary1,domain0);
assert(b1);
assert((bary1.X()>=0)&&(bary1.X()<=1));
assert((bary1.Y()>=0)&&(bary1.Y()<=1));
if(!((bary1.X()+bary1.Y())<=1.00001))
{
printf("\n SUM %.4f \n",bary1.X()+bary1.Y());
assert(0);
}
assert(I1<max);
int domain;
bool b2=Sub(I,bary,I1,bary1,vect1,domain);
assert(b2);
ScalarType diff=(vect1+vect).Norm();
if (domain0==0)
assert(domain==0);
if ((domain0==1)&&(domain==2))
assert(0);
/*if (domain0!=domain)
assert(0);*/
if (diff>0.001)
{
printf("\n DIFF %.4f domain SUM %d domain SUB %d \n",diff,domain0,domain);
//assert(0);
}
//assert(fabs(vect1.X()-vect.X())<0.0001);
//assert(fabs(vect1.Y()-vect.Y())<0.0001);
}
}
}
}
}
}
void InitProjectionMatrix(ScalarType radius=(ScalarType)0.1)
{
for (int i=0;i<isoParam->ParaMesh()->vert.size();i++)
{
int I=isoParam->ParaMesh()->vert[i].T().N();
vcg::Point2<ScalarType> bary=isoParam->ParaMesh()->vert[i].T().P();
CoordType origin;
isoParamTheta(I,bary,origin);
CoordType XAxis,YAxis; // tangent axis in 3D Object space
///get tangent directions
bool done=TangentDir(I,bary,XAxis,YAxis,(ScalarType)0.1);
if (!done)
{
(*ProjMatrix)[i].SetIdentity();
continue;
}
// must compute res2d such that: res2d.X() * XAxis + res2d.Y() * YAxis + dontCare * ZAxis = vect3d
CoordType ZAxis = -(XAxis^YAxis).Normalize()*XAxis.Norm();
(*ProjMatrix)[i].SetColumn(0,XAxis);
(*ProjMatrix)[i].SetColumn(1,YAxis);
(*ProjMatrix)[i].SetColumn(2,ZAxis);
vcg::Invert((*ProjMatrix)[i]);
}
}
void GetCurvature(const int &I,const vcg::Point2<ScalarType> &alpha_beta,
CoordType &d1,CoordType &d2,ScalarType &k1,ScalarType &k2)
{
ParamFace* face=NULL;
CoordType baryVal;
isoParam->Theta(I,alpha_beta,face,baryVal);
ParamVertex *v0=face->V(0);
ParamVertex *v1=face->V(1);
ParamVertex *v2=face->V(2);
d1=v0->PD1()*baryVal.X()+v1->PD1()*baryVal.Y()+v2->PD1()*baryVal.Z();
d2=v0->PD2()*baryVal.X()+v1->PD2()*baryVal.Y()+v2->PD2()*baryVal.Z();
k1=v0->K1()*baryVal.X()+v1->K1()*baryVal.Y()+v2->K1()*baryVal.Z();
k2=v0->K2()*baryVal.X()+v1->K2()*baryVal.Y()+v2->K2()*baryVal.Z();
}
void GetProjectionMatrix(const int &I,const vcg::Point2<ScalarType> &alpha_beta,
vcg::Matrix33<ScalarType> &projMatr) const
{
ParamFace* face=NULL;
CoordType baryVal;
int dom=isoParam->Theta(I,alpha_beta,face,baryVal);
if (dom==-1)
{
projMatr.SetIdentity();
return;
}
ParamVertex *v0=face->V(0);
ParamVertex *v1=face->V(1);
ParamVertex *v2=face->V(2);
projMatr=(*ProjMatrix)[v0]*baryVal.X();
projMatr+=(*ProjMatrix)[v1]*baryVal.Y();
projMatr+=(*ProjMatrix)[v2]*baryVal.Z();
}
bool Project(const int &I,const vcg::Point2<ScalarType> &bary,
const CoordType &vect3d, // in object space
vcg::Point2<ScalarType> &res2d) const
{
vcg::Matrix33f m;
GetProjectionMatrix(I,bary,m);
ScalarType deltaX=vect3d*m.GetRow(0);
ScalarType deltaY=vect3d*m.GetRow(1);
// two axis in alpha-beta space that will be orthogonal in UV-Space
const ScalarType h=(ScalarType)2.0/sqrt((ScalarType)3.0);
vcg::Point2<ScalarType> XP=vcg::Point2<ScalarType>(1,0);
vcg::Point2<ScalarType> YP=vcg::Point2<ScalarType>(-0.5,h);
res2d = XP*deltaX + YP*deltaY;
//res2d=vcg::Point2<ScalarType>(deltaX,deltaY);
return true;
}
// WEIGHTED INTERPOLATION OF POINTS IN TANGENT SPACE
/// weights MUST be normalized
bool Interpolate(const int &I0,const vcg::Point2<ScalarType> &alpha_beta0,
const int &I1,const vcg::Point2<ScalarType> &alpha_beta1,
ScalarType weight,
int &I_res,
vcg::Point2<ScalarType> &alpha_beta_res)
{
int IndexDomain;
int kind=isoParam->InterpolationSpace(I0,I1,IndexDomain);
if (kind==-1)
return false;
vcg::Point2<ScalarType> transformed0,transformed1;
///interpolate in a face
if (kind==0)
{
isoParam->GE2(IndexDomain,alpha_beta0,transformed0);
isoParam->GE2(IndexDomain,alpha_beta1,transformed1);
}
else
///interpolate in a diamond
if (kind==1)
{
isoParam->GE1(I0,alpha_beta0,IndexDomain,transformed0);
isoParam->GE1(I1,alpha_beta1,IndexDomain,transformed1);
}
else
{
isoParam->GE0(I0,alpha_beta0,IndexDomain,transformed0);
isoParam->GE0(I1,alpha_beta1,IndexDomain,transformed1);
}
vcg::Point2<ScalarType> UV_interp=transformed0*weight+transformed1*(1.0-weight);
///FINALLY......
///transform back to alpha,beta,I
if (kind==0)
{
isoParam->inv_GE2(IndexDomain,UV_interp,alpha_beta_res);
I_res=IndexDomain;
}
else
if (kind==1)
{
isoParam->inv_GE1(IndexDomain,UV_interp,I_res,alpha_beta_res);
}
else
isoParam->inv_GE0(IndexDomain,UV_interp,I_res,alpha_beta_res);
return true;
}
ScalarType AbstractArea()
{
ScalarType Cnum=sqrt(3.0)/4.0;
return(isoParam->AbsMesh()->fn*Cnum);
}
// WEIGHTED INTERPOLATION OF POINTS IN TANGENT SPACE
/// weights MUST be normalized
bool Interpolate(const std::vector<int> &I,
const std::vector<vcg::Point2<ScalarType> > &alpha_beta,
const std::vector<ScalarType> &weights,
int &I_res,
vcg::Point2<ScalarType> &alpha_beta_res,
int *num=NULL)
{
int size;
if (num==NULL)
size=alpha_beta.size();
else
size=*num;
int IndexDomain;
int kind=isoParam->InterpolationSpace(I,IndexDomain,num);
if (kind==-1)
return false;
std::vector<vcg::Point2<ScalarType> > transformed;
transformed.resize(size);
///interpolate in a face
if (kind==0)
{
for (int i=0;i<size;i++)
isoParam->GE2(IndexDomain,alpha_beta[i],transformed[i]);
}
else
///interpolate in a diamond
if (kind==1)
{
for (int i=0;i<size;i++)
isoParam->GE1(I[i],alpha_beta[i],IndexDomain,transformed[i]);
}
else
{
for (int i=0;i<size;i++)
bool b0=isoParam->GE0(I[i],alpha_beta[i],IndexDomain,transformed[i]);
}
/// do the interpolation
vcg::Point2<ScalarType> UV_interp=vcg::Point2<ScalarType>(0,0);
for (int i=0;i<size;i++)
UV_interp+=(transformed[i]*weights[i]);
///FINALLY......
///transform back to alpha,beta,I
if (kind==0)
{
isoParam->inv_GE2(IndexDomain,UV_interp,alpha_beta_res);
I_res=IndexDomain;
}
else
if (kind==1)
{
isoParam->inv_GE1(IndexDomain,UV_interp,I_res,alpha_beta_res);
}
else
isoParam->inv_GE0(IndexDomain,UV_interp,I_res,alpha_beta_res);
return true;
}
};
#endif

3
src/meshlabplugins/filter_layer/filter_layer.cpp
View File

@ -29,8 +29,7 @@
#include "filter_layer.h"
#include <vcg/complex/trimesh/clean.h>
#include<vcg/complex/trimesh/append.h>
#include<vcg/complex/append.h>

14
src/meshlabplugins/filter_measure/filter_measure.cpp
View File

@ -26,15 +26,15 @@
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <vcg/complex/trimesh/clean.h>
#include <vcg/complex/trimesh/inertia.h>
#include <vcg/complex/trimesh/stat.h>
#include <vcg/complex/algorithms/clean.h>
#include <vcg/complex/algorithms/inertia.h>
#include <vcg/complex/algorithms/stat.h>
#include <vcg/complex/trimesh/update/selection.h>
#include<vcg/complex/trimesh/append.h>
#include <vcg/complex/algorithms/update/selection.h>
#include<vcg/complex/append.h>
#include<vcg/simplex/face/pos.h>
#include<vcg/complex/trimesh/bitquad_support.h>
#include<vcg/complex/trimesh/bitquad_optimization.h>
#include<vcg/complex/algorithms/bitquad_support.h>
#include<vcg/complex/algorithms/bitquad_optimization.h>
#include "filter_measure.h"
using namespace std;

2
src/meshlabplugins/filter_meshing/filter_meshing.pro
View File

@ -1,6 +1,6 @@
include (../../shared.pri)
HEADERS += $$VCGDIR/vcg/complex/trimesh/clean.h\
HEADERS += $$VCGDIR/vcg/complex/algorithms/clean.h\
quadric_simp.h \
quadric_tex_simp.h \
meshfilter.h

22
src/meshlabplugins/filter_meshing/meshfilter.cpp
View File

@ -21,17 +21,17 @@
****************************************************************************/
#include "meshfilter.h"
#include <vcg/complex/trimesh/clean.h>
#include <vcg/complex/trimesh/stat.h>
#include <vcg/complex/trimesh/smooth.h>
#include <vcg/complex/trimesh/hole.h>
#include <vcg/complex/trimesh/refine_loop.h>
#include <vcg/complex/trimesh/bitquad_support.h>
#include <vcg/complex/trimesh/bitquad_creation.h>
#include <vcg/complex/trimesh/clustering.h>
#include <vcg/complex/trimesh/attribute_seam.h>
#include <vcg/complex/trimesh/update/curvature.h>
#include <vcg/complex/trimesh/update/curvature_fitting.h>
#include <vcg/complex/algorithms/clean.h>
#include <vcg/complex/algorithms/stat.h>
#include <vcg/complex/algorithms/smooth.h>
#include <vcg/complex/algorithms/hole.h>
#include <vcg/complex/algorithms/refine_loop.h>
#include <vcg/complex/algorithms/bitquad_support.h>
#include <vcg/complex/algorithms/bitquad_creation.h>
#include <vcg/complex/algorithms/clustering.h>
#include <vcg/complex/algorithms/attribute_seam.h>
#include <vcg/complex/algorithms/update/curvature.h>
#include <vcg/complex/algorithms/update/curvature_fitting.h>
#include <vcg/space/normal_extrapolation.h>
#include <vcg/space/fitting3.h>
#include "quadric_tex_simp.h"

3
src/meshlabplugins/filter_meshing/quadric_simp.cpp
View File

@ -20,9 +20,6 @@
* *
****************************************************************************/
#include "meshfilter.h"
#include <vcg/complex/local_optimization.h>
#include <vcg/complex/local_optimization/tri_edge_collapse_quadric.h>
#include <vcg/container/simple_temporary_data.h>
#include "quadric_simp.h"
using namespace vcg;

10
src/meshlabplugins/filter_meshing/quadric_simp.h
View File

@ -55,11 +55,11 @@ Added remove non manifold and quadric simplification filter.
#include <limits>
#include "meshfilter.h"
#include <vcg/complex/trimesh/update/position.h>
#include <vcg/complex/trimesh/update/bounding.h>
#include <vcg/complex/trimesh/update/selection.h>
#include <vcg/complex/local_optimization.h>
#include <vcg/complex/local_optimization/tri_edge_collapse_quadric.h>
#include <vcg/complex/algorithms/update/position.h>
#include <vcg/complex/algorithms/update/bounding.h>
#include <vcg/complex/algorithms/update/selection.h>
#include <vcg/complex/algorithms/local_optimization.h>
#include <vcg/complex/algorithms/local_optimization/tri_edge_collapse_quadric.h>
#include <vcg/container/simple_temporary_data.h>
namespace vcg {
namespace tri {

4
src/meshlabplugins/filter_meshing/quadric_tex_simp.h
View File

@ -23,8 +23,8 @@
#ifndef __QUADRICTEXSIMP_H
#define __QUADRICTEXSIMP_H
#include <vcg/complex/local_optimization.h>
#include <vcg/complex/local_optimization/tri_edge_collapse_quadric.h>
#include <vcg/complex/algorithms/local_optimization.h>
#include <vcg/complex/algorithms/local_optimization/tri_edge_collapse_quadric.h>
#include <vcg/container/simple_temporary_data.h>
#include "algebra5.h"
#include "quadric5.h"

12
src/meshlabplugins/filter_mls/mlsplugin.cpp
View File

@ -30,12 +30,12 @@
#include <common/interfaces.h>
#include <vcg/complex/trimesh/clean.h>
#include <vcg/complex/trimesh/refine.h>
#include <vcg/complex/trimesh/refine_loop.h>
#include <vcg/complex/trimesh/append.h>
#include <vcg/complex/trimesh/create/advancing_front.h>
#include <vcg/complex/trimesh/create/marching_cubes.h>
#include <vcg/complex/algorithms/clean.h>
#include <vcg/complex/algorithms/refine.h>
#include <vcg/complex/algorithms/refine_loop.h>
#include <vcg/complex/append.h>
#include <vcg/complex/algorithms/create/advancing_front.h>
#include <vcg/complex/algorithms/create/marching_cubes.h>
#include "mlsmarchingcube.h"
#include "mlsplugin.h"

714
src/meshlabplugins/filter_mutualinfo/alignset.cpp
View File

@ -1,357 +1,357 @@
#include <iostream>
#include <GL/glew.h>
#include <QGLContext>
#include <QDomNode>
#include <QDir>
#include <QFileInfo>
#include <QFile>
#include <QTextStream>
#include <QGLFramebufferObject>
#include "alignset.h"
#include <wrap/gl/shot.h>
#include <wrap/gl/camera.h>
//#include <vcg/math/shot.h>
#include <vcg/complex/trimesh/update/normal.h>
#include <vcg/complex/trimesh/update/bounding.h>
#include <wrap/io_trimesh/import_ply.h>
#include "shutils.h"
using namespace std;
AlignSet::AlignSet(): mode(COMBINE),
target(NULL), render(NULL),
vbo(0), nbo(0), cbo(0), ibo(0), error(0){
box.SetNull();
correspList = new QList<PointCorrespondence*>();
imageRatio=1;
}
AlignSet::~AlignSet() {
if(target) delete []target;
if(render) delete []render;
delete correspList;
}
void AlignSet::initializeGL() {
programs[COLOR] = createShaders("varying vec4 color; void main() { gl_Position = ftransform(); color = gl_Color; }",
"varying vec4 color; void main() { gl_FragColor = color; }");
programs[NORMALMAP] = createShaders("varying vec3 normal; void main() { normal = gl_NormalMatrix * gl_Normal; gl_Position = ftransform(); }",
"varying vec3 normal; void main() { "
"vec3 color = normalize(normal); color = color * 0.5 + 0.5; gl_FragColor = vec4(color, 1.0); }");
programs[COMBINE] = createShaders("varying vec3 normal; varying vec4 color; void main() { "
"normal = gl_NormalMatrix * gl_Normal; gl_Position = ftransform(); color = gl_Color; }",
"varying vec3 normal; varying vec4 color; void main() { "
"vec3 ncolor = normalize(normal); ncolor = ncolor * 0.5 + 0.5; "
"float t = color.x*color.x; gl_FragColor = (1-t)*color + t*(vec4(ncolor, 1.0)); }");
programs[SPECULAR] = createShaders("varying vec3 reflection; void main() { "
"vec3 normal = normalize(gl_NormalMatrix * gl_Normal); vec4 position = gl_ModelViewMatrix * gl_Vertex; "
"reflection = reflect(position.xyz, normal); gl_Position = ftransform(); }",
"varying vec3 reflection; varying vec4 color; void main() { "
"vec4 ncolor; ncolor.xyz = normalize(reflection); ncolor.w = 1.0; gl_FragColor = ncolor * 0.5 + 0.5; }");
programs[SILHOUETTE] = createShaders("varying vec4 color; void main() { gl_Position = ftransform(); color = gl_Color; }",
"varying vec4 color; void main() { gl_FragColor = color; }");
programs[SPECAMB] = createShaders("varying vec3 reflection; varying vec4 color; void main() { "
"vec3 normal = normalize(gl_NormalMatrix * gl_Normal); vec4 position = gl_ModelViewMatrix * gl_Vertex; "
"reflection = reflect(position.xyz, normal); gl_Position = ftransform(); color = gl_Color; }",
"varying vec3 reflection; varying vec4 color; void main() { "
"vec3 ncolor = normalize(reflection); ncolor = ncolor * 0.5 + 0.5; "
"float t = color.x*color.x; gl_FragColor = (1-t)*color + t*(vec4(ncolor, 1.0)); }");
// generate a new VBO and get the associated ID
glGenBuffersARB(1, &vbo);
glGenBuffersARB(1, &nbo);
glGenBuffersARB(1, &cbo);
glGenBuffersARB(1, &ibo);
}
//resample image IF too big.
void AlignSet::resize(int max_side) {
int w = image->width();
int h = image->height();
if(image->isNull()) {
w = 1024;
h = 768;
}
if(w > max_side) {
h = h*max_side/w;
w = max_side;
}
if(h > max_side) {
w = w*max_side/h;
h = max_side;
}
wt=w;
ht=h;
if(target) delete []target;
if(render) delete []render;
target = new unsigned char[w*h];
render = new unsigned char[w*h];
if(image->isNull()) return;
//resize image and store values into render
QImage im;
if(w != image->width() || h != image->height())
im = image->scaled(w, h, Qt::IgnoreAspectRatio); //Qt::KeepAspectRatio);
else im = *image;
im.save("image.jpg");
assert(w == im.width());
assert(h == im.height());
QColor color;
int offset = 0;
//equalize image
int histo[256];
memset(histo, 0, 256*sizeof(int));
for (int y = h-1; y >= 0; y--) {
for (int x = 0; x < w; x++) {
color.setRgb(im.pixel(x, y));
unsigned char c = (unsigned char)(color.red() * 0.3f + color.green() * 0.59f + color.blue() * 0.11f);
target[offset] = c;
histo[c]++;
offset++;
}
}
#ifdef RESCALE_HISTO
int cumulative[256];
cumulative[0] = histo[0];
for(int i = 1; i < 256; i++)
cumulative[i] = cumulative[i-1] + histo[i];
int min = 0;
int max = 255;
for(int i = 0; i < 256; i++) {
if(cumulative[i] > 20) break;
min = i;
}
//invert cumulative..
cumulative[255] = histo[255];
for(int i = 254; i >= 0; i--)
cumulative[i] = cumulative[i+1] + histo[i];
for(int i = 255; i >= 0; i--) {
if(cumulative[i] > 20) break;
max = i;
}
assert(max > min);
//rescale between min and max (should use bresenham but I am lazy
unsigned char equa[256];
for(int i = 0; i < 256; i++) {
if(i < min) equa[i] = 0;
if(i > max) equa[i] = 255;
equa[i] = (255*(i - min))/(max - min);
}
for(int i = 0; i < w*h; i++)
target[i] = equa[target[i]];
#endif
}
void AlignSet::renderScene(vcg::Shot<float> &view, int component) {
QSize fbosize(wt,ht);
QGLFramebufferObjectFormat frmt;
frmt.setInternalTextureFormat(GL_RGBA);
frmt.setAttachment(QGLFramebufferObject::Depth);
QGLFramebufferObject fbo(fbosize,frmt);
float _near, _far;
_near=0.1;
_far=10000;
GlShot< vcg::Shot<float> >::GetNearFarPlanes(view, mesh->bbox, _near, _far);
//assert(_near <= _far);
if(_near <= 0) _near = 0.1;
if(_far < _near) _far = 1000;
//GLenum err = glGetError();
//render to FBO
fbo.bind();
glViewport(0, 0, wt, ht);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
GlShot< vcg::Shot<float> >::SetView(shot, 0.5*_near, 2*_far);
// err = glGetError();
bool use_colors = false;
bool use_normals = false;
int program = programs[mode]; //standard pipeline
switch(mode){
case COLOR:
use_colors = true;
break;
case NORMALMAP:
case SPECULAR:
use_normals = true;
break;
case COMBINE:
case SPECAMB:
use_colors = true;
use_normals = true;
break;
case SILHOUETTE:
break;
default: assert(0);
}
glDisable(GL_LIGHTING);
//bind indices
glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, ibo);
//bind vertices
glEnable(GL_COLOR_MATERIAL);
glBindBufferARB(GL_ARRAY_BUFFER_ARB, vbo);
glEnableClientState(GL_VERTEX_ARRAY); // activate vertex coords array
glVertexPointer(3, GL_FLOAT, 0, 0); // last param is offset, not ptr
// err = glGetError();
glUseProgram(program);
// err = glGetError();
if(use_colors) {
glBindBufferARB(GL_ARRAY_BUFFER_ARB, cbo);
glEnableClientState(GL_COLOR_ARRAY);
glColorPointer(4, GL_UNSIGNED_BYTE, 0, 0);
}
if(use_normals) {
glBindBufferARB(GL_ARRAY_BUFFER_ARB, nbo);
glEnableClientState(GL_NORMAL_ARRAY);
glNormalPointer(GL_FLOAT, 0, 0);
}
// err = glGetError();
int start = 0;
int tot = 30000;
if (mesh->fn>0)
{
while(start < mesh->fn) {
glDrawElements(GL_TRIANGLES, tot*3, GL_UNSIGNED_INT, (void *)(start*3*sizeof(int)));
start += tot;
if(start + tot > mesh->fn)
tot = mesh->fn - start;
}
}
else glDrawArrays(GL_POINTS, 0, mesh->vn);
render = new unsigned char[wt*ht];
glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
glPixelStorei(GL_PACK_ALIGNMENT, 1);
switch(component) {
case 0: glReadPixels( 0, 0, wt, ht, GL_RED, GL_UNSIGNED_BYTE, render); break;
case 1: glReadPixels( 0, 0, wt, ht, GL_GREEN, GL_UNSIGNED_BYTE, render); break;
case 2: glReadPixels( 0, 0, wt, ht, GL_BLUE, GL_UNSIGNED_BYTE, render); break;
case 3: glReadPixels( 0, 0, wt, ht, GL_ALPHA, GL_UNSIGNED_BYTE, render); break;
case 4: break;
}
//err = glGetError();
glDisableClientState(GL_VERTEX_ARRAY); // deactivate vertex array
if(use_colors) glDisableClientState(GL_COLOR_ARRAY);
if(use_normals) glDisableClientState(GL_NORMAL_ARRAY);
//err = glGetError();
// bind with 0, so, switch back to normal pointer operation
glBindBufferARB(GL_ARRAY_BUFFER_ARB, 0);
glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
switch(mode) {
case SILHOUETTE:
case COLOR:
case COMBINE:
case NORMALMAP: glEnable(GL_LIGHTING); break;
default: break;
}
// standard opengl pipeline is re-activated
glUseProgram(0);
GlShot< vcg::Shot<float> >::UnsetView();
glFinish();
/*QImage l=fbo.toImage();
l.save("rendering.jpg");*/
fbo.release();
}
void AlignSet::readRender(int component) {
QSize fbosize(wt,ht);
QGLFramebufferObjectFormat frmt;
frmt.setInternalTextureFormat(GL_RGBA);
frmt.setAttachment(QGLFramebufferObject::Depth);
QGLFramebufferObject fbo(fbosize,frmt);
fbo.bind();
glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
glPixelStorei(GL_PACK_ALIGNMENT, 1);
switch(component) {
case 0: glReadPixels( 0, 0, width(), height(), GL_RED, GL_UNSIGNED_BYTE, render); break;
case 1: glReadPixels( 0, 0, width(), height(), GL_GREEN, GL_UNSIGNED_BYTE, render); break;
case 2: glReadPixels( 0, 0, width(), height(), GL_BLUE, GL_UNSIGNED_BYTE, render); break;
case 3: glReadPixels( 0, 0, width(), height(), GL_ALPHA, GL_UNSIGNED_BYTE, render); break;
}
QImage l=fbo.toImage();
l.save("puppo.jpg");
fbo.release();
}
GLuint AlignSet::createShaderFromFiles(QString name) {
QString vert = "shaders/" + name + ".vert";
QString frag = "shaders/" + name + ".frag";
const char *vs_src = ShaderUtils::importShaders(vert.toAscii().data());
if(!vs_src) {
cerr << "Could not load shader: " << qPrintable(vert) << endl;
return 0;
}
const char *fs_src = ShaderUtils::importShaders(frag.toAscii().data());
if(!fs_src) {
cerr << "Could not load shader: " << qPrintable(frag) << endl;
return 0;
}
return createShaders(vs_src, fs_src);
}
GLuint AlignSet::createShaders(const char *vert, const char *frag) {
GLuint vs = glCreateShader(GL_VERTEX_SHADER);
glShaderSource(vs, 1, &vert, NULL);
ShaderUtils::compileShader(vs);
GLuint fs = glCreateShader(GL_FRAGMENT_SHADER);
glShaderSource(fs, 1, &frag, NULL);
ShaderUtils::compileShader(fs);
GLuint prog = glCreateProgram();
glAttachShader(prog, vs);
glAttachShader(prog, fs);
ShaderUtils::linkShaderProgram(prog);
return prog;
}
#include <iostream>
#include <GL/glew.h>
#include <QGLContext>
#include <QDomNode>
#include <QDir>
#include <QFileInfo>
#include <QFile>
#include <QTextStream>
#include <QGLFramebufferObject>
#include "alignset.h"
#include <wrap/gl/shot.h>
#include <wrap/gl/camera.h>
//#include <vcg/math/shot.h>
#include <vcg/complex/algorithms/update/normal.h>
#include <vcg/complex/algorithms/update/bounding.h>
#include <wrap/io_trimesh/import_ply.h>
#include "shutils.h"
using namespace std;
AlignSet::AlignSet(): mode(COMBINE),
target(NULL), render(NULL),
vbo(0), nbo(0), cbo(0), ibo(0), error(0){
box.SetNull();
correspList = new QList<PointCorrespondence*>();
imageRatio=1;
}
AlignSet::~AlignSet() {
if(target) delete []target;
if(render) delete []render;
delete correspList;
}
void AlignSet::initializeGL() {
programs[COLOR] = createShaders("varying vec4 color; void main() { gl_Position = ftransform(); color = gl_Color; }",
"varying vec4 color; void main() { gl_FragColor = color; }");
programs[NORMALMAP] = createShaders("varying vec3 normal; void main() { normal = gl_NormalMatrix * gl_Normal; gl_Position = ftransform(); }",
"varying vec3 normal; void main() { "
"vec3 color = normalize(normal); color = color * 0.5 + 0.5; gl_FragColor = vec4(color, 1.0); }");
programs[COMBINE] = createShaders("varying vec3 normal; varying vec4 color; void main() { "
"normal = gl_NormalMatrix * gl_Normal; gl_Position = ftransform(); color = gl_Color; }",
"varying vec3 normal; varying vec4 color; void main() { "
"vec3 ncolor = normalize(normal); ncolor = ncolor * 0.5 + 0.5; "
"float t = color.x*color.x; gl_FragColor = (1-t)*color + t*(vec4(ncolor, 1.0)); }");
programs[SPECULAR] = createShaders("varying vec3 reflection; void main() { "
"vec3 normal = normalize(gl_NormalMatrix * gl_Normal); vec4 position = gl_ModelViewMatrix * gl_Vertex; "
"reflection = reflect(position.xyz, normal); gl_Position = ftransform(); }",
"varying vec3 reflection; varying vec4 color; void main() { "
"vec4 ncolor; ncolor.xyz = normalize(reflection); ncolor.w = 1.0; gl_FragColor = ncolor * 0.5 + 0.5; }");
programs[SILHOUETTE] = createShaders("varying vec4 color; void main() { gl_Position = ftransform(); color = gl_Color; }",
"varying vec4 color; void main() { gl_FragColor = color; }");
programs[SPECAMB] = createShaders("varying vec3 reflection; varying vec4 color; void main() { "
"vec3 normal = normalize(gl_NormalMatrix * gl_Normal); vec4 position = gl_ModelViewMatrix * gl_Vertex; "
"reflection = reflect(position.xyz, normal); gl_Position = ftransform(); color = gl_Color; }",
"varying vec3 reflection; varying vec4 color; void main() { "
"vec3 ncolor = normalize(reflection); ncolor = ncolor * 0.5 + 0.5; "
"float t = color.x*color.x; gl_FragColor = (1-t)*color + t*(vec4(ncolor, 1.0)); }");
// generate a new VBO and get the associated ID
glGenBuffersARB(1, &vbo);
glGenBuffersARB(1, &nbo);
glGenBuffersARB(1, &cbo);
glGenBuffersARB(1, &ibo);
}
//resample image IF too big.
void AlignSet::resize(int max_side) {
int w = image->width();
int h = image->height();
if(image->isNull()) {
w = 1024;
h = 768;
}
if(w > max_side) {
h = h*max_side/w;
w = max_side;
}
if(h > max_side) {
w = w*max_side/h;
h = max_side;
}
wt=w;
ht=h;
if(target) delete []target;
if(render) delete []render;
target = new unsigned char[w*h];
render = new unsigned char[w*h];
if(image->isNull()) return;
//resize image and store values into render
QImage im;
if(w != image->width() || h != image->height())
im = image->scaled(w, h, Qt::IgnoreAspectRatio); //Qt::KeepAspectRatio);
else im = *image;
im.save("image.jpg");
assert(w == im.width());
assert(h == im.height());
QColor color;
int offset = 0;
//equalize image
int histo[256];
memset(histo, 0, 256*sizeof(int));
for (int y = h-1; y >= 0; y--) {
for (int x = 0; x < w; x++) {
color.setRgb(im.pixel(x, y));
unsigned char c = (unsigned char)(color.red() * 0.3f + color.green() * 0.59f + color.blue() * 0.11f);
target[offset] = c;
histo[c]++;
offset++;
}
}
#ifdef RESCALE_HISTO
int cumulative[256];
cumulative[0] = histo[0];
for(int i = 1; i < 256; i++)
cumulative[i] = cumulative[i-1] + histo[i];
int min = 0;
int max = 255;
for(int i = 0; i < 256; i++) {
if(cumulative[i] > 20) break;
min = i;
}
//invert cumulative..
cumulative[255] = histo[255];
for(int i = 254; i >= 0; i--)
cumulative[i] = cumulative[i+1] + histo[i];
for(int i = 255; i >= 0; i--) {
if(cumulative[i] > 20) break;
max = i;
}
assert(max > min);
//rescale between min and max (should use bresenham but I am lazy
unsigned char equa[256];
for(int i = 0; i < 256; i++) {
if(i < min) equa[i] = 0;
if(i > max) equa[i] = 255;
equa[i] = (255*(i - min))/(max - min);
}
for(int i = 0; i < w*h; i++)
target[i] = equa[target[i]];
#endif
}
void AlignSet::renderScene(vcg::Shot<float> &view, int component) {
QSize fbosize(wt,ht);
QGLFramebufferObjectFormat frmt;
frmt.setInternalTextureFormat(GL_RGBA);
frmt.setAttachment(QGLFramebufferObject::Depth);
QGLFramebufferObject fbo(fbosize,frmt);
float _near, _far;
_near=0.1;
_far=10000;
GlShot< vcg::Shot<float> >::GetNearFarPlanes(view, mesh->bbox, _near, _far);
//assert(_near <= _far);
if(_near <= 0) _near = 0.1;
if(_far < _near) _far = 1000;
//GLenum err = glGetError();
//render to FBO
fbo.bind();
glViewport(0, 0, wt, ht);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
GlShot< vcg::Shot<float> >::SetView(shot, 0.5*_near, 2*_far);
// err = glGetError();
bool use_colors = false;
bool use_normals = false;
int program = programs[mode]; //standard pipeline
switch(mode){
case COLOR:
use_colors = true;
break;
case NORMALMAP:
case SPECULAR:
use_normals = true;
break;
case COMBINE:
case SPECAMB:
use_colors = true;
use_normals = true;
break;
case SILHOUETTE:
break;
default: assert(0);
}
glDisable(GL_LIGHTING);
//bind indices
glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, ibo);
//bind vertices
glEnable(GL_COLOR_MATERIAL);
glBindBufferARB(GL_ARRAY_BUFFER_ARB, vbo);
glEnableClientState(GL_VERTEX_ARRAY); // activate vertex coords array
glVertexPointer(3, GL_FLOAT, 0, 0); // last param is offset, not ptr
// err = glGetError();
glUseProgram(program);
// err = glGetError();
if(use_colors) {
glBindBufferARB(GL_ARRAY_BUFFER_ARB, cbo);
glEnableClientState(GL_COLOR_ARRAY);
glColorPointer(4, GL_UNSIGNED_BYTE, 0, 0);
}
if(use_normals) {
glBindBufferARB(GL_ARRAY_BUFFER_ARB, nbo);
glEnableClientState(GL_NORMAL_ARRAY);
glNormalPointer(GL_FLOAT, 0, 0);
}
// err = glGetError();
int start = 0;
int tot = 30000;
if (mesh->fn>0)
{
while(start < mesh->fn) {
glDrawElements(GL_TRIANGLES, tot*3, GL_UNSIGNED_INT, (void *)(start*3*sizeof(int)));
start += tot;
if(start + tot > mesh->fn)
tot = mesh->fn - start;
}
}
else glDrawArrays(GL_POINTS, 0, mesh->vn);
render = new unsigned char[wt*ht];
glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
glPixelStorei(GL_PACK_ALIGNMENT, 1);
switch(component) {
case 0: glReadPixels( 0, 0, wt, ht, GL_RED, GL_UNSIGNED_BYTE, render); break;
case 1: glReadPixels( 0, 0, wt, ht, GL_GREEN, GL_UNSIGNED_BYTE, render); break;
case 2: glReadPixels( 0, 0, wt, ht, GL_BLUE, GL_UNSIGNED_BYTE, render); break;
case 3: glReadPixels( 0, 0, wt, ht, GL_ALPHA, GL_UNSIGNED_BYTE, render); break;
case 4: break;
}
//err = glGetError();
glDisableClientState(GL_VERTEX_ARRAY); // deactivate vertex array
if(use_colors) glDisableClientState(GL_COLOR_ARRAY);
if(use_normals) glDisableClientState(GL_NORMAL_ARRAY);
//err = glGetError();
// bind with 0, so, switch back to normal pointer operation
glBindBufferARB(GL_ARRAY_BUFFER_ARB, 0);
glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
switch(mode) {
case SILHOUETTE:
case COLOR:
case COMBINE:
case NORMALMAP: glEnable(GL_LIGHTING); break;
default: break;
}
// standard opengl pipeline is re-activated
glUseProgram(0);
GlShot< vcg::Shot<float> >::UnsetView();
glFinish();
/*QImage l=fbo.toImage();
l.save("rendering.jpg");*/
fbo.release();
}
void AlignSet::readRender(int component) {
QSize fbosize(wt,ht);
QGLFramebufferObjectFormat frmt;
frmt.setInternalTextureFormat(GL_RGBA);
frmt.setAttachment(QGLFramebufferObject::Depth);
QGLFramebufferObject fbo(fbosize,frmt);
fbo.bind();
glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
glPixelStorei(GL_PACK_ALIGNMENT, 1);
switch(component) {
case 0: glReadPixels( 0, 0, width(), height(), GL_RED, GL_UNSIGNED_BYTE, render); break;
case 1: glReadPixels( 0, 0, width(), height(), GL_GREEN, GL_UNSIGNED_BYTE, render); break;
case 2: glReadPixels( 0, 0, width(), height(), GL_BLUE, GL_UNSIGNED_BYTE, render); break;
case 3: glReadPixels( 0, 0, width(), height(), GL_ALPHA, GL_UNSIGNED_BYTE, render); break;
}
QImage l=fbo.toImage();
l.save("puppo.jpg");
fbo.release();
}
GLuint AlignSet::createShaderFromFiles(QString name) {
QString vert = "shaders/" + name + ".vert";
QString frag = "shaders/" + name + ".frag";
const char *vs_src = ShaderUtils::importShaders(vert.toAscii().data());
if(!vs_src) {
cerr << "Could not load shader: " << qPrintable(vert) << endl;
return 0;
}
const char *fs_src = ShaderUtils::importShaders(frag.toAscii().data());
if(!fs_src) {
cerr << "Could not load shader: " << qPrintable(frag) << endl;
return 0;
}
return createShaders(vs_src, fs_src);
}
GLuint AlignSet::createShaders(const char *vert, const char *frag) {
GLuint vs = glCreateShader(GL_VERTEX_SHADER);
glShaderSource(vs, 1, &vert, NULL);
ShaderUtils::compileShader(vs);
GLuint fs = glCreateShader(GL_FRAGMENT_SHADER);
glShaderSource(fs, 1, &frag, NULL);
ShaderUtils::compileShader(fs);
GLuint prog = glCreateProgram();
glAttachShader(prog, vs);
glAttachShader(prog, fs);
ShaderUtils::linkShaderProgram(prog);
return prog;
}

441
src/meshlabplugins/filter_perceptualmetric/filter_perceptualmetric.cpp
View File

@ -1,224 +1,217 @@
/****************************************************************************
* MeshLab o o *
* A versatile mesh processing toolbox o o *
* _ O _ *
* Copyright(C) 2005 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* 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 (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
/****************************************************************************
History
$Log: filter_perceptualmetric.cpp,v $
****************************************************************************/
#include <QtGui>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <meshlab/meshmodel.h>
#include <meshlab/interfaces.h>
#include <vcg/complex/trimesh/clean.h>
#include <vcg/complex/trimesh/update/position.h>
#include <vcg/complex/trimesh/update/normal.h>
#include <vcg/complex/trimesh/update/flag.h>
#include <vcg/complex/trimesh/update/bounding.h>
#include <vcg/complex/trimesh/update/quality.h>
#include <vcg/complex/trimesh/point_sampling.h>
#include <vcg/complex/trimesh/create/resampler.h>
#include <vcg/simplex/face/distance.h>
#include <vcg/complex/trimesh/update/color.h>
#include <vcg/complex/trimesh/geodesic.h>
#include <vcg/space/index/grid_static_ptr.h>
// filter interface
#include "filter_perceptualmetric.h"
// perceptual metrics core
#include "perceptualmetrics.h"
using namespace vcg;
using namespace std;
// Constructor usually performs only two simple tasks of filling the two lists
// - typeList: with all the possible id of the filtering actions
// - actionList with the corresponding actions. If you want to add icons to your filtering actions you can do here by construction the QActions accordingly
FilterPerceptualMetric::FilterPerceptualMetric()
{
typeList
<< FP_ROUGHNESS_MULTISCALE
<< FP_ROUGHNESS_SMOOTHING
<< FP_STRAIN_ENERGY
;
foreach(FilterIDType tt , types())
actionList << new QAction(filterName(tt), this);
}
// ST() must return the very short string describing each filtering action
// (this string is used also to define the menu entry)
QString FilterPerceptualMetric::filterName(FilterIDType filterId)
{
switch(filterId)
{
case FP_ROUGHNESS_MULTISCALE : return QString("Roughness-multiscale perceptual metric");
case FP_ROUGHNESS_SMOOTHING : return QString("Roughness-smoothing perceptual metric");
case FP_STRAIN_ENERGY : return QString("Energy Strain-based perceptual metric");
default : assert(0); return QString("unknown filter!!!!");
}
}
// Info() must return the longer string describing each filtering action
// (this string is used in the About plugin dialog)
QString FilterPerceptualMetric::filterInfo(FilterIDType filterId)
{
switch(filterId)
{
case FP_ROUGHNESS_MULTISCALE : return QString("Evaluate the perceptual difference between two given meshes; this perceptual metric is based on roughness measure, the roughness is evaluated at multiple-scale using dihedral angles.");
case FP_ROUGHNESS_SMOOTHING : return QString("Evaluate the perceptual difference between two given meshes; this perceptual metric is based on roughness measure, the roughness is evaluated using a smoothed version of the 3D models.");
case FP_STRAIN_ENERGY : return QString("Evaluate the perceptual difference between two given meshes; this perceptual metric is based on strain energy.");
default : assert(0); return QString("unknown filter!!!!");
}
}
int FilterPerceptualMetric::getRequirements(QAction *action)
{
switch(ID(action))
{
case FP_ROUGHNESS_MULTISCALE : return MeshModel::MM_VERTFACETOPO | MeshModel::MM_FACENORMAL;
case FP_ROUGHNESS_SMOOTHING : return MeshModel::MM_VERTFACETOPO;
case FP_STRAIN_ENERGY : return MeshModel::MM_VERTFACETOPO;
default: assert(0);
}
return 0;
}
// This function define the needed parameters for each filter. Return true if the filter has some parameters
// it is called every time, so you can set the default value of parameters according to the mesh
// For each parameter you need to define,
// - the name of the parameter,
// - the string shown in the dialog
// - the default value
// - a possibly long string describing the meaning of that parameter (shown as a popup help in the dialog)
void FilterPerceptualMetric::initParameterSet(QAction *action, MeshDocument & md, FilterParameterSet & parlst)
{
switch(ID(action))
{
case FP_ROUGHNESS_MULTISCALE :
{
MeshModel *refmesh = md.mm();
foreach (refmesh, md.meshList)
if (refmesh != md.mm()) break;
parlst.addMesh("ReferenceMesh", refmesh, "Reference Mesh",
"The original mesh.");
parlst.addMesh("InputMesh", md.mm(), "Mesh",
"The mesh where the perceptual impairment of the processing is evaluated.");
} break;
case FP_ROUGHNESS_SMOOTHING :
{
MeshModel *refmesh = md.mm();
foreach (refmesh, md.meshList)
if (refmesh != md.mm()) break;
parlst.addMesh("ReferenceMesh", refmesh, "Reference Mesh",
"The original mesh.");
parlst.addMesh("InputMesh", md.mm(), "Mesh",
"The mesh where the perceptual impairment of the processing is evaluated.");
} break;
case FP_STRAIN_ENERGY :
{
MeshModel *refmesh = md.mm();
foreach (refmesh, md.meshList)
if (refmesh != md.mm()) break;
parlst.addMesh("ReferenceMesh", refmesh, "Reference Mesh",
"The original mesh.");
parlst.addMesh("InputMesh", md.mm(), "Mesh",
"The mesh where the perceptual impairment of the processing is evaluated.");
} break;
default : assert(0);
}
}
bool FilterPerceptualMetric::applyFilter(QAction *action, MeshDocument &md, FilterParameterSet & par, vcg::CallBackPos *cb)
{
switch(ID(action))
{
case FP_ROUGHNESS_MULTISCALE :
{
MeshModel* mesh0 = par.getMesh("ReferenceMesh"); // reference mesh
MeshModel* mesh1 = par.getMesh("InputMesh"); // the processed mesh
double globalimpact = PerceptualMetrics<CMeshO>::roughnessMultiscale(mesh0->cm, mesh1->cm);
Log(0,"This metric is not implemented yet!!");
}
break;
case FP_ROUGHNESS_SMOOTHING :
{
MeshModel* mesh0 = par.getMesh("ReferenceMesh"); // reference mesh
MeshModel* mesh1 = par.getMesh("InputMesh"); // the processed mesh
double globalimpact = PerceptualMetrics<CMeshO>::roughnessSmoothing(mesh0->cm, mesh1->cm);
Log(0,"This metric is not implemented yet!!");
}
break;
case FP_STRAIN_ENERGY :
{
MeshModel* mesh0 = par.getMesh("ReferenceMesh"); // reference mesh
MeshModel* mesh1 = par.getMesh("InputMesh"); // the processed mesh
double globalimpact = PerceptualMetrics<CMeshO>::strainEnergy(mesh0->cm, mesh1->cm);
Log(0,"Perceptual Distance: %f",globalimpact);
}
break;
default : assert(0);
}
return true;
}
MeshFilterInterface::FilterClass FilterPerceptualMetric::getClass(QAction *action)
{
switch(ID(action))
{
case FP_ROUGHNESS_MULTISCALE : return Quality;
case FP_ROUGHNESS_SMOOTHING : return Quality;
case FP_STRAIN_ENERGY : return Quality;
default: assert(0);
}
return FilterClass(0);
}
Q_EXPORT_PLUGIN(FilterPerceptualMetric)
/****************************************************************************
* MeshLab o o *
* A versatile mesh processing toolbox o o *
* _ O _ *
* Copyright(C) 2005 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* 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 (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
/****************************************************************************
History
$Log: filter_perceptualmetric.cpp,v $
****************************************************************************/
#include <QtGui>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <meshlab/meshmodel.h>
#include <meshlab/interfaces.h>
#include <vcg/complex/algorithms/point_sampling.h>
#include <vcg/complex/algorithms/create/resampler.h>
#include <vcg/simplex/face/distance.h>
#include <vcg/complex/algorithms/geodesic.h>
#include <vcg/space/index/grid_static_ptr.h>
// filter interface
#include "filter_perceptualmetric.h"
// perceptual metrics core
#include "perceptualmetrics.h"
using namespace vcg;
using namespace std;
// Constructor usually performs only two simple tasks of filling the two lists
// - typeList: with all the possible id of the filtering actions
// - actionList with the corresponding actions. If you want to add icons to your filtering actions you can do here by construction the QActions accordingly
FilterPerceptualMetric::FilterPerceptualMetric()
{
typeList
<< FP_ROUGHNESS_MULTISCALE
<< FP_ROUGHNESS_SMOOTHING
<< FP_STRAIN_ENERGY
;
foreach(FilterIDType tt , types())
actionList << new QAction(filterName(tt), this);
}
// ST() must return the very short string describing each filtering action
// (this string is used also to define the menu entry)
QString FilterPerceptualMetric::filterName(FilterIDType filterId)
{
switch(filterId)
{
case FP_ROUGHNESS_MULTISCALE : return QString("Roughness-multiscale perceptual metric");
case FP_ROUGHNESS_SMOOTHING : return QString("Roughness-smoothing perceptual metric");
case FP_STRAIN_ENERGY : return QString("Energy Strain-based perceptual metric");
default : assert(0); return QString("unknown filter!!!!");
}
}
// Info() must return the longer string describing each filtering action
// (this string is used in the About plugin dialog)
QString FilterPerceptualMetric::filterInfo(FilterIDType filterId)
{
switch(filterId)
{
case FP_ROUGHNESS_MULTISCALE : return QString("Evaluate the perceptual difference between two given meshes; this perceptual metric is based on roughness measure, the roughness is evaluated at multiple-scale using dihedral angles.");
case FP_ROUGHNESS_SMOOTHING : return QString("Evaluate the perceptual difference between two given meshes; this perceptual metric is based on roughness measure, the roughness is evaluated using a smoothed version of the 3D models.");
case FP_STRAIN_ENERGY : return QString("Evaluate the perceptual difference between two given meshes; this perceptual metric is based on strain energy.");
default : assert(0); return QString("unknown filter!!!!");
}
}
int FilterPerceptualMetric::getRequirements(QAction *action)
{
switch(ID(action))
{
case FP_ROUGHNESS_MULTISCALE : return MeshModel::MM_VERTFACETOPO | MeshModel::MM_FACENORMAL;
case FP_ROUGHNESS_SMOOTHING : return MeshModel::MM_VERTFACETOPO;
case FP_STRAIN_ENERGY : return MeshModel::MM_VERTFACETOPO;
default: assert(0);
}
return 0;
}
// This function define the needed parameters for each filter. Return true if the filter has some parameters
// it is called every time, so you can set the default value of parameters according to the mesh
// For each parameter you need to define,
// - the name of the parameter,
// - the string shown in the dialog
// - the default value
// - a possibly long string describing the meaning of that parameter (shown as a popup help in the dialog)
void FilterPerceptualMetric::initParameterSet(QAction *action, MeshDocument & md, FilterParameterSet & parlst)
{
switch(ID(action))
{
case FP_ROUGHNESS_MULTISCALE :
{
MeshModel *refmesh = md.mm();
foreach (refmesh, md.meshList)
if (refmesh != md.mm()) break;
parlst.addMesh("ReferenceMesh", refmesh, "Reference Mesh",
"The original mesh.");
parlst.addMesh("InputMesh", md.mm(), "Mesh",
"The mesh where the perceptual impairment of the processing is evaluated.");
} break;
case FP_ROUGHNESS_SMOOTHING :
{
MeshModel *refmesh = md.mm();
foreach (refmesh, md.meshList)
if (refmesh != md.mm()) break;
parlst.addMesh("ReferenceMesh", refmesh, "Reference Mesh",
"The original mesh.");
parlst.addMesh("InputMesh", md.mm(), "Mesh",
"The mesh where the perceptual impairment of the processing is evaluated.");
} break;
case FP_STRAIN_ENERGY :
{
MeshModel *refmesh = md.mm();
foreach (refmesh, md.meshList)
if (refmesh != md.mm()) break;
parlst.addMesh("ReferenceMesh", refmesh, "Reference Mesh",
"The original mesh.");
parlst.addMesh("InputMesh", md.mm(), "Mesh",
"The mesh where the perceptual impairment of the processing is evaluated.");
} break;
default : assert(0);
}
}
bool FilterPerceptualMetric::applyFilter(QAction *action, MeshDocument &md, FilterParameterSet & par, vcg::CallBackPos *cb)
{
switch(ID(action))
{
case FP_ROUGHNESS_MULTISCALE :
{
MeshModel* mesh0 = par.getMesh("ReferenceMesh"); // reference mesh
MeshModel* mesh1 = par.getMesh("InputMesh"); // the processed mesh
double globalimpact = PerceptualMetrics<CMeshO>::roughnessMultiscale(mesh0->cm, mesh1->cm);
Log(0,"This metric is not implemented yet!!");
}
break;
case FP_ROUGHNESS_SMOOTHING :
{
MeshModel* mesh0 = par.getMesh("ReferenceMesh"); // reference mesh
MeshModel* mesh1 = par.getMesh("InputMesh"); // the processed mesh
double globalimpact = PerceptualMetrics<CMeshO>::roughnessSmoothing(mesh0->cm, mesh1->cm);
Log(0,"This metric is not implemented yet!!");
}
break;
case FP_STRAIN_ENERGY :
{
MeshModel* mesh0 = par.getMesh("ReferenceMesh"); // reference mesh
MeshModel* mesh1 = par.getMesh("InputMesh"); // the processed mesh
double globalimpact = PerceptualMetrics<CMeshO>::strainEnergy(mesh0->cm, mesh1->cm);
Log(0,"Perceptual Distance: %f",globalimpact);
}
break;
default : assert(0);
}
return true;
}
MeshFilterInterface::FilterClass FilterPerceptualMetric::getClass(QAction *action)
{
switch(ID(action))
{
case FP_ROUGHNESS_MULTISCALE : return Quality;
case FP_ROUGHNESS_SMOOTHING : return Quality;
case FP_STRAIN_ENERGY : return Quality;
default: assert(0);
}
return FilterClass(0);
}
Q_EXPORT_PLUGIN(FilterPerceptualMetric)

484
src/meshlabplugins/filter_perceptualmetric/perceptualmetrics.h
View File

@ -1,242 +1,242 @@
/****************************************************************************
* MeshLab o o *
* A versatile mesh processing toolbox o o *
* _ O _ *
* Copyright(C) 2005 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* 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 (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
/****************************************************************************
History
$Log: perceptualmetrics.h,v $
****************************************************************************/
#ifndef PERCEPTUALMETRICS_H
#define PERCEPTUALMETRICS_H
#include <QObject>
#include <common/interfaces.h>
#include <vcg/complex/trimesh/stat.h>
template <class MeshType>
class PerceptualMetrics
{
// definitions
public:
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FacePointer FacePointer;
typedef typename MeshType::FaceIterator FaceIterator;
typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::FaceContainer FaceContainer;
// private methods
private:
static void planeProjection(double a, double b, double c, double d,
const CoordType &axis1, const CoordType &axis2, const CoordType & O,
const CoordType &p, double &xcoor, double &ycoor)
{
double u = p[0];
double v = p[1];
double w = p[2];
CoordType Pproj;
double num = a*u + b*v + c*w + d;
// den = a*a + b*b + c*c that is assumed one in this case (!)
Pproj[0] = u - a * num; // num/den
Pproj[1] = v - b * num; // num/den
Pproj[2] = w - c * num; // num/den
Pproj -= O;
xcoor = axis1 * Pproj;
ycoor = axis2 * Pproj;
}
static void computingDisplacementFunctionCoefficients(FacePointer f0, FacePointer f1,
double & a1u, double & a2u, double & a3u, double & a1v, double & a2v, double & a3v,
double & area)
{
CoordType pD = f0->V(0)->P();
CoordType pE = f0->V(1)->P();
CoordType pF = f0->V(2)->P();
CoordType pDp = f1->V(0)->P();
CoordType pEp = f1->V(1)->P();
CoordType pFp = f1->V(2)->P();
CoordType axis1 = (pE-pD);
CoordType axis2 = (pF-pD);
CoordType N = axis1 ^ axis2;
// axis adjustment
axis2 = N ^ axis1;
axis1.Normalize();
axis2.Normalize();
N.Normalize();
double a = N[0];
double b = N[1];
double c = N[2];
double d = -(a * pD[0] + b * pD[1] + c * pD[2]);
// triangle in the local reference system
double xD,yD,xE,yE,xF,yF;
planeProjection(a,b,c,d, axis1, axis2, pD, pD, xD, yD);
planeProjection(a,b,c,d, axis1, axis2, pD, pE, xE, yE);
planeProjection(a,b,c,d, axis1, axis2, pD, pF, xF, yF);
// triangle in the local reference system after deformations
double xDp,yDp,xEp,yEp,xFp,yFp;
planeProjection(a,b,c,d, axis1, axis2, pD, pDp, xDp, yDp);
planeProjection(a,b,c,d, axis1, axis2, pD, pEp, xEp, yEp);
planeProjection(a,b,c,d, axis1, axis2, pD, pFp, xFp, yFp);
// deformation
double uD = xDp - xD;
double vD = yDp - yD;
double uE = xEp - xE;
double vE = yEp - yE;
double uF = xFp - xF;
double vF = yFp - yF;
// parameterization
double AA = (yE-yF) * xD + (yF-yD) * xE + (yD-yE) * xF; // AA = 2A
double factor = 1.0 / (AA);
a1u = factor * ((xE*yF - xF*yE)*uD + (xF*yD - xD*yF)*uE + (xD*yE - xE*yD)*uF);
a2u = factor * ((yE-yF)*uD + (yF-yD)*uE + (yD-yE)*uF);
a3u = factor * ((xF-xE)*uD + (xD-xF)*uE + (xE-xD)*uF);
a1v = factor * ((xE*yF - xF*yE)*vD + (xF*yD - xD*yF)*vE + (xD*yE - xE*yD)*vF);
a2v = factor * ((yE-yF)*vD + (yF-yD)*vE + (yD-yE)*vF);
a3v = factor * ((xF-xE)*vD + (xD-xF)*vE + (xE-xD)*vF);
area = 0.5 * AA;
}
// public methods
public:
static double roughnessMultiscale(MeshType & refmesh, MeshType & mesh)
{
//...TODO...
return 0.0;
}
static double roughnessSmoothing(MeshType & refmesh, MeshType & mesh)
{
//...TODO...
return 0.0;
}
static double strainEnergy(MeshType & refmesh, MeshType & mesh)
{
double epsilonx, epsilony, epsilonz;
double gammaxy, gammayz, gammazx;
double epsilonx_prime, epsilony_prime, epsilonz_prime;
double epsilon_ii_squared;
double epsilon_ij_prime_epsilon_ij_prime;
double Sdelta; // Area of each triangle
double ni = 0.0; // Poisson's ratio
double E = 1.0; // Young's modulus
double lambda; // Lame's first parameter
double G;
double Warea, Wdistortion, W;
double u,v;
double a1u,a2u,a3u;
double a1v,a2v,a3v;
double x,y;
double D;
FacePointer f0;
FacePointer f1;
W = 0.0;
for (int i = 0; i < refmesh.fn; i++)
{
f0 = &refmesh.face[i];
f1 = &mesh.face[i];
// displacement function are:
// u = a1u + x a2u + y a3u
// v = a1v + x a2v + y a3v
computingDisplacementFunctionCoefficients(f0, f1, a1u, a2u, a3u, a1v, a2v, a3v, Sdelta);
// epsilonx = du / dx
epsilonx = a2u;
// epsilony = dv / dy
epsilony = a3v;
// epsilonz
epsilonz = (ni / (ni - 1.0)) * (epsilonx + epsilony);
// gammaxy = dv / dx + du / dy
gammaxy = a2v + a3u;
gammayz = 0.0;
gammazx = 0.0;
// strain energy for a single triangle
///////////////////////////////////////////////////////////////////////////////////////
epsilonx_prime = epsilonx - (epsilonx + epsilony + epsilonz) / 3.0;
epsilony_prime = epsilony - (epsilonx + epsilony + epsilonz) / 3.0;
epsilonz_prime = epsilonz - (epsilonx + epsilony + epsilonz) / 3.0;
epsilon_ii_squared = (epsilonx + epsilony + epsilonz) * (epsilonx + epsilony + epsilonz);
epsilon_ij_prime_epsilon_ij_prime = epsilonx_prime * epsilonx_prime + epsilony_prime * epsilony_prime +
epsilonz_prime * epsilonz_prime + 0.5 * (gammaxy*gammaxy + gammayz*gammayz + gammazx*gammazx);
lambda = (E * ni) / ((1.0 + ni) * (1 - 2.0 * ni));
G = E / (2.0 * (1.0 + ni));
Warea = 0.5 * (lambda + 0.666666666666 * G) * epsilon_ii_squared * Sdelta;
Wdistortion = G * epsilon_ij_prime_epsilon_ij_prime * Sdelta;
W += Warea + Wdistortion;
}
// Average Strain Energy (ASE)
//////////////////////////////////////////////////////////////////////////
double area1 = vcg::tri::Stat<MeshType>::ComputeMeshArea(refmesh);
//double area2 = vcg::tri::Stat<MeshType>::ComputeMeshArea(mesh);
// the area if the deformed mesh is considered the same of the original one,
// since the deformation is assumed to be small
return W / (area1);
}
};
#endif /* PERCEPTUALMETRICS_H */
/****************************************************************************
* MeshLab o o *
* A versatile mesh processing toolbox o o *
* _ O _ *
* Copyright(C) 2005 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* 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 (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
/****************************************************************************
History
$Log: perceptualmetrics.h,v $
****************************************************************************/
#ifndef PERCEPTUALMETRICS_H
#define PERCEPTUALMETRICS_H
#include <QObject>
#include <common/interfaces.h>
#include <vcg/complex/algorithms/stat.h>
template <class MeshType>
class PerceptualMetrics
{
// definitions
public:
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::FacePointer FacePointer;
typedef typename MeshType::FaceIterator FaceIterator;
typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::FaceContainer FaceContainer;
// private methods
private:
static void planeProjection(double a, double b, double c, double d,
const CoordType &axis1, const CoordType &axis2, const CoordType & O,
const CoordType &p, double &xcoor, double &ycoor)
{
double u = p[0];
double v = p[1];
double w = p[2];
CoordType Pproj;
double num = a*u + b*v + c*w + d;
// den = a*a + b*b + c*c that is assumed one in this case (!)
Pproj[0] = u - a * num; // num/den
Pproj[1] = v - b * num; // num/den
Pproj[2] = w - c * num; // num/den
Pproj -= O;
xcoor = axis1 * Pproj;
ycoor = axis2 * Pproj;
}
static void computingDisplacementFunctionCoefficients(FacePointer f0, FacePointer f1,
double & a1u, double & a2u, double & a3u, double & a1v, double & a2v, double & a3v,
double & area)
{
CoordType pD = f0->V(0)->P();
CoordType pE = f0->V(1)->P();
CoordType pF = f0->V(2)->P();
CoordType pDp = f1->V(0)->P();
CoordType pEp = f1->V(1)->P();
CoordType pFp = f1->V(2)->P();
CoordType axis1 = (pE-pD);
CoordType axis2 = (pF-pD);
CoordType N = axis1 ^ axis2;
// axis adjustment
axis2 = N ^ axis1;
axis1.Normalize();
axis2.Normalize();
N.Normalize();
double a = N[0];
double b = N[1];
double c = N[2];
double d = -(a * pD[0] + b * pD[1] + c * pD[2]);
// triangle in the local reference system
double xD,yD,xE,yE,xF,yF;
planeProjection(a,b,c,d, axis1, axis2, pD, pD, xD, yD);
planeProjection(a,b,c,d, axis1, axis2, pD, pE, xE, yE);
planeProjection(a,b,c,d, axis1, axis2, pD, pF, xF, yF);
// triangle in the local reference system after deformations
double xDp,yDp,xEp,yEp,xFp,yFp;
planeProjection(a,b,c,d, axis1, axis2, pD, pDp, xDp, yDp);
planeProjection(a,b,c,d, axis1, axis2, pD, pEp, xEp, yEp);
planeProjection(a,b,c,d, axis1, axis2, pD, pFp, xFp, yFp);
// deformation
double uD = xDp - xD;
double vD = yDp - yD;
double uE = xEp - xE;
double vE = yEp - yE;
double uF = xFp - xF;
double vF = yFp - yF;
// parameterization
double AA = (yE-yF) * xD + (yF-yD) * xE + (yD-yE) * xF; // AA = 2A
double factor = 1.0 / (AA);
a1u = factor * ((xE*yF - xF*yE)*uD + (xF*yD - xD*yF)*uE + (xD*yE - xE*yD)*uF);
a2u = factor * ((yE-yF)*uD + (yF-yD)*uE + (yD-yE)*uF);
a3u = factor * ((xF-xE)*uD + (xD-xF)*uE + (xE-xD)*uF);
a1v = factor * ((xE*yF - xF*yE)*vD + (xF*yD - xD*yF)*vE + (xD*yE - xE*yD)*vF);
a2v = factor * ((yE-yF)*vD + (yF-yD)*vE + (yD-yE)*vF);
a3v = factor * ((xF-xE)*vD + (xD-xF)*vE + (xE-xD)*vF);
area = 0.5 * AA;
}
// public methods
public:
static double roughnessMultiscale(MeshType & refmesh, MeshType & mesh)
{
//...TODO...
return 0.0;
}
static double roughnessSmoothing(MeshType & refmesh, MeshType & mesh)
{
//...TODO...
return 0.0;
}
static double strainEnergy(MeshType & refmesh, MeshType & mesh)
{
double epsilonx, epsilony, epsilonz;
double gammaxy, gammayz, gammazx;
double epsilonx_prime, epsilony_prime, epsilonz_prime;
double epsilon_ii_squared;
double epsilon_ij_prime_epsilon_ij_prime;
double Sdelta; // Area of each triangle
double ni = 0.0; // Poisson's ratio
double E = 1.0; // Young's modulus
double lambda; // Lame's first parameter
double G;
double Warea, Wdistortion, W;
double u,v;
double a1u,a2u,a3u;
double a1v,a2v,a3v;
double x,y;
double D;
FacePointer f0;
FacePointer f1;
W = 0.0;
for (int i = 0; i < refmesh.fn; i++)
{
f0 = &refmesh.face[i];
f1 = &mesh.face[i];
// displacement function are:
// u = a1u + x a2u + y a3u
// v = a1v + x a2v + y a3v
computingDisplacementFunctionCoefficients(f0, f1, a1u, a2u, a3u, a1v, a2v, a3v, Sdelta);
// epsilonx = du / dx
epsilonx = a2u;
// epsilony = dv / dy
epsilony = a3v;
// epsilonz
epsilonz = (ni / (ni - 1.0)) * (epsilonx + epsilony);
// gammaxy = dv / dx + du / dy
gammaxy = a2v + a3u;
gammayz = 0.0;
gammazx = 0.0;
// strain energy for a single triangle
///////////////////////////////////////////////////////////////////////////////////////
epsilonx_prime = epsilonx - (epsilonx + epsilony + epsilonz) / 3.0;
epsilony_prime = epsilony - (epsilonx + epsilony + epsilonz) / 3.0;
epsilonz_prime = epsilonz - (epsilonx + epsilony + epsilonz) / 3.0;
epsilon_ii_squared = (epsilonx + epsilony + epsilonz) * (epsilonx + epsilony + epsilonz);
epsilon_ij_prime_epsilon_ij_prime = epsilonx_prime * epsilonx_prime + epsilony_prime * epsilony_prime +
epsilonz_prime * epsilonz_prime + 0.5 * (gammaxy*gammaxy + gammayz*gammayz + gammazx*gammazx);
lambda = (E * ni) / ((1.0 + ni) * (1 - 2.0 * ni));
G = E / (2.0 * (1.0 + ni));
Warea = 0.5 * (lambda + 0.666666666666 * G) * epsilon_ii_squared * Sdelta;
Wdistortion = G * epsilon_ij_prime_epsilon_ij_prime * Sdelta;
W += Warea + Wdistortion;
}
// Average Strain Energy (ASE)
//////////////////////////////////////////////////////////////////////////
double area1 = vcg::tri::Stat<MeshType>::ComputeMeshArea(refmesh);
//double area2 = vcg::tri::Stat<MeshType>::ComputeMeshArea(mesh);
// the area if the deformed mesh is considered the same of the original one,
// since the deformation is assumed to be small
return W / (area1);
}
};
#endif /* PERCEPTUALMETRICS_H */