first release version
This commit is contained in:
parent
52648c58ad
commit
aa2d1751f7
|
|
@ -0,0 +1,351 @@
|
|||
#ifndef __MIQ__
|
||||
#define __MIQ__
|
||||
|
||||
|
||||
#define SIZEQUADS 512
|
||||
#define V_SIZE 1
|
||||
#define SIZEPARA 1024
|
||||
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
#include "mesh_type.h"
|
||||
#include "quadrangulator.h"
|
||||
#include "poisson_solver.h"
|
||||
#include "param_stats.h"
|
||||
#include "seams_initializer.h"
|
||||
#include "vertex_indexing.h"
|
||||
|
||||
#include <wrap/io_trimesh/import.h>
|
||||
#include <wrap/io_trimesh/export.h>
|
||||
|
||||
#define USECOMISO
|
||||
|
||||
template <class ScalarType>
|
||||
ScalarType Gauss(ScalarType &value)
|
||||
{
|
||||
const ScalarType E_NEPER=2.71828;
|
||||
const ScalarType den=sqrt(2.0*M_PI);
|
||||
ScalarType exponent=-0.5*pow(value,2);
|
||||
ScalarType res=((1.0/den)*pow(E_NEPER,exponent));
|
||||
return res;
|
||||
}
|
||||
|
||||
template <class MeshType,class PolyMesh>
|
||||
class MIQ{
|
||||
|
||||
public:
|
||||
typename MeshType::template PerFaceAttributeHandle<float> Handle_Stiffness;
|
||||
|
||||
// Different stiffening mode
|
||||
enum StiffMode{NO_STIFF = 0,GAUSSIAN = 1,ITERATIVE = 2};
|
||||
|
||||
// Init
|
||||
MIQ(MeshType &_mesh):mesh(_mesh),PSolver(mesh){};
|
||||
|
||||
// Load a mesh from file
|
||||
bool LoadMesh(const std::string PathMesh)
|
||||
{
|
||||
int position=PathMesh.find(".ply");
|
||||
int err;
|
||||
if (position==-1)
|
||||
{
|
||||
position=PathMesh.find(".obj");
|
||||
//vcg::tri::io::ImporterOBJ<CMesh>::Info objInf;
|
||||
int mask;
|
||||
vcg::tri::io::ImporterOBJ<CMesh>::LoadMask(PathMesh.c_str(),mask);
|
||||
err=vcg::tri::io::ImporterOBJ<CMesh>::Open(mesh,PathMesh.c_str(),mask);
|
||||
assert(position!=-1);
|
||||
if (err!=0)return false;
|
||||
}
|
||||
else
|
||||
{
|
||||
err=vcg::tri::io::ImporterPLY<CMesh>::Open(mesh,PathMesh.c_str());
|
||||
if (err!=ply::E_NOERROR)return false;
|
||||
}
|
||||
///UPDATE MESH STUFF
|
||||
vcg::tri::UpdateBounding<CMesh>::Box(mesh);
|
||||
vcg::tri::UpdateNormal<CMesh>::PerVertexNormalizedPerFace(mesh);
|
||||
vcg::tri::UpdateNormal<CMesh>::PerFaceNormalized(mesh);
|
||||
vcg::tri::UpdateTopology<CMesh>::FaceFace(mesh);
|
||||
vcg::tri::UpdateTopology<CMesh>::VertexFace(mesh);
|
||||
vcg::tri::UpdateFlags<CMesh>::FaceBorderFromFF(mesh);
|
||||
vcg::tri::UpdateFlags<CMesh>::VertexBorderFromFace(mesh);
|
||||
return true;
|
||||
}
|
||||
|
||||
// Load a field from file
|
||||
bool LoadField(const std::string PathFField,const std::string PathMesh=NULL)
|
||||
{
|
||||
SInit.Init(&mesh);
|
||||
|
||||
///FIELD LOADING
|
||||
int position=PathFField.find(".ffield");
|
||||
if (position==-1)
|
||||
{
|
||||
position=PathFField.find(".grad");
|
||||
assert(position!=-1);
|
||||
SInit.InitFromGradOBJ(PathFField,PathMesh);
|
||||
}
|
||||
else
|
||||
SInit.InitFromFField(PathFField);
|
||||
|
||||
VInd.Init(&mesh);
|
||||
//mesh.ScaleToUnitBox();
|
||||
VInd.InitMapping();
|
||||
VInd.InitFaceIntegerVal();
|
||||
VInd.InitSeamInfo();
|
||||
return true;
|
||||
}
|
||||
|
||||
// Load a mesh and a field from file
|
||||
bool LoadMeshField(const std::string PathMesh, const std::string PathFField)
|
||||
{
|
||||
bool loaded=LoadMesh(PathMesh);
|
||||
if (!loaded)return false;
|
||||
LoadField(PathFField,PathMesh);
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
// Parametrize the mesh
|
||||
void Solve(StiffMode stiffMode, double Stiffness = 5.0,
|
||||
double GradientSize = 30.0, bool DirectRound = false,
|
||||
int iter = 5, int localIter = 5, bool DoRound = true)
|
||||
{
|
||||
if (mesh.fn==0)return;
|
||||
InitDefaultStiffening();
|
||||
if (stiffMode==GAUSSIAN)
|
||||
{
|
||||
AddGaussStiffening(Stiffness);
|
||||
PSolver.SolvePoisson(GradientSize,1.f,DirectRound,localIter,DoRound);
|
||||
}
|
||||
else
|
||||
if (stiffMode==ITERATIVE)
|
||||
{
|
||||
for (int i=0;i<iter;i++)
|
||||
{
|
||||
PSolver.SolvePoisson(GradientSize,1.f,DirectRound,localIter,DoRound);
|
||||
int nflips=NumFlips(mesh);
|
||||
ScaleGLtoInt();
|
||||
bool folded=updateStiffening(GradientSize);
|
||||
ScaleInttoGL();
|
||||
printf("ITERATION %d FLIPS %d \n",i,nflips);
|
||||
if (!folded)break;
|
||||
}
|
||||
}
|
||||
else
|
||||
if (stiffMode==NO_STIFF)
|
||||
{
|
||||
PSolver.SolvePoisson(GradientSize,1.f,DirectRound,localIter,DoRound);
|
||||
}
|
||||
int nflips=NumFlips(mesh);
|
||||
printf("**** END OPTIMIZING #FLIPS %d ****\n",nflips);
|
||||
fflush(stdout);
|
||||
SelectFlippedFaces(mesh);
|
||||
UpdateUVBox();
|
||||
}
|
||||
|
||||
// Generate a quad mesh starting from the parametrization
|
||||
void Quadrangulate(PolyMesh &poly,double factor = 1)
|
||||
{
|
||||
Quad.Quadrangulate(mesh,poly,factor);
|
||||
}
|
||||
|
||||
// void removeDuplicateVertices(std::vector<std::vector<double> >& V, std::vector< std::vector<int > >& F);
|
||||
|
||||
// void exportQuadMesh(std::string out);
|
||||
|
||||
//bool LoadData(std::string PathMesh,
|
||||
// std::string PathFField);
|
||||
// Bounding box of the param domain
|
||||
vcg::Box2<double> UVBox;
|
||||
|
||||
void ScaleGLtoInt()
|
||||
{
|
||||
double factor=(double)SIZEPARA/(double)SIZEQUADS;
|
||||
for (unsigned int i=0;i<mesh.face.size();i++)
|
||||
{
|
||||
if (mesh.face[i].IsD())continue;
|
||||
mesh.face[i].WT(0).P()*=factor;
|
||||
mesh.face[i].WT(1).P()*=factor;
|
||||
mesh.face[i].WT(2).P()*=factor;
|
||||
}
|
||||
}
|
||||
|
||||
void ScaleInttoGL()
|
||||
{
|
||||
double factor=(double)SIZEQUADS/(double)SIZEPARA;
|
||||
for (unsigned int i=0;i<mesh.face.size();i++)
|
||||
{
|
||||
if (mesh.face[i].IsD())continue;
|
||||
mesh.face[i].WT(0).P()*=factor;
|
||||
mesh.face[i].WT(1).P()*=factor;
|
||||
mesh.face[i].WT(2).P()*=factor;
|
||||
}
|
||||
}
|
||||
|
||||
void UpdateUVBox()
|
||||
{
|
||||
UVBox.SetNull();
|
||||
for (unsigned int i=0;i<mesh.face.size();i++)
|
||||
{
|
||||
UVBox.Add((mesh.face[i].WT(0).P()));
|
||||
UVBox.Add((mesh.face[i].WT(1).P()));
|
||||
UVBox.Add((mesh.face[i].WT(2).P()));
|
||||
}
|
||||
}
|
||||
|
||||
// Quadrangulator
|
||||
Quadrangulator<CMesh,PolyMesh> Quad;
|
||||
|
||||
void colorByStiffening(typename MeshType::ScalarType MaxVal=16)
|
||||
{
|
||||
bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
|
||||
assert(hasStiffness);
|
||||
for (unsigned int i=0;i<mesh.face.size();i++)
|
||||
{
|
||||
//CMesh::ScalarType val=MaxVal-mesh.face[i].stiffening+1;
|
||||
CMesh::ScalarType val=MaxVal-Handle_Stiffness[i]+1;
|
||||
if (val<1)val=1;
|
||||
mesh.face[i].C()=vcg::Color4b::ColorRamp(1.0,MaxVal,val);
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
|
||||
void AddStiffening(typename MeshType::ScalarType C,int radius=4)
|
||||
{
|
||||
bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
|
||||
if(!hasStiffness)
|
||||
Handle_Stiffness=vcg::tri::Allocator<CMesh>::AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));
|
||||
|
||||
bool hasSingular = vcg::tri::HasPerVertexAttribute(mesh,std::string("Singular"));
|
||||
assert(hasSingular);
|
||||
CMesh::PerVertexAttributeHandle<bool> Handle_Singular;
|
||||
Handle_Singular=vcg::tri::Allocator<CMesh>::GetPerVertexAttribute<bool>(mesh,std::string("Singular"));
|
||||
|
||||
std::vector<CMesh::VertexType*> to_stiff;
|
||||
for(unsigned int i=0;i<mesh.vert.size();i++)
|
||||
{
|
||||
CMesh::VertexType *v=&mesh.vert[i];
|
||||
if (v->IsD())continue;
|
||||
//if (!v->IsSingular())continue;
|
||||
if (!Handle_Singular[v])continue;
|
||||
to_stiff.push_back(v);
|
||||
}
|
||||
for(unsigned int i=0;i<mesh.face.size();i++)
|
||||
{
|
||||
|
||||
CMesh::FaceType *f=&(mesh.face[i]);
|
||||
if (f->IsD())continue;
|
||||
if (!f->IsV())continue;
|
||||
to_stiff.push_back(f->V(0));
|
||||
to_stiff.push_back(f->V(1));
|
||||
to_stiff.push_back(f->V(2));
|
||||
}
|
||||
std::sort(to_stiff.begin(),to_stiff.end());
|
||||
std::vector<CMesh::VertexType*>::iterator new_end=std::unique(to_stiff.begin(),to_stiff.end());
|
||||
int dist=distance(to_stiff.begin(),new_end);
|
||||
to_stiff.resize(dist);
|
||||
for (unsigned int i=0;i<to_stiff.size();i++)
|
||||
{
|
||||
CMesh::VertexType *v=to_stiff[i];
|
||||
for (int r=0;r<radius;r++)
|
||||
{
|
||||
CMesh::ScalarType stiffVal=((CMesh::ScalarType)r)/(CMesh::ScalarType)radius;//((ScalarType)(radius-r))/(ScalarType)radius;
|
||||
stiffVal*=3.0;
|
||||
stiffVal=Gauss(stiffVal)/0.4;
|
||||
stiffVal=1+(stiffVal*C);
|
||||
std::vector<CMesh::FaceType*> ring;
|
||||
//mesh.GetConnectedFaces(v,r,ring);
|
||||
VFExtendedStarVF(v,r,ring);
|
||||
///then set stiffening
|
||||
for (unsigned int k=0;k<ring.size();k++)
|
||||
{
|
||||
CMesh::FaceType* f=ring[k];
|
||||
//if (f->stiffening<stiffVal)
|
||||
// f->stiffening=stiffVal;
|
||||
if (Handle_Stiffness[f]<stiffVal)
|
||||
Handle_Stiffness[f]=stiffVal;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
bool updateStiffening(typename MeshType::ScalarType grad_size)
|
||||
{
|
||||
bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
|
||||
if(!hasStiffness)
|
||||
Handle_Stiffness=vcg::tri::Allocator<CMesh>::AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));
|
||||
|
||||
bool flipped = NumFlips(mesh)>0;
|
||||
//if (h == 0.0)
|
||||
// return flipped;
|
||||
//
|
||||
//assert(h != 0.0);
|
||||
|
||||
if (!flipped)
|
||||
return false;
|
||||
CMesh::ScalarType maxL=0;
|
||||
CMesh::ScalarType maxD=0;
|
||||
if (flipped)
|
||||
{
|
||||
const double c = 1.0;
|
||||
const double d = 5.0;
|
||||
|
||||
for (unsigned int i = 0; i < mesh.face.size(); ++i)
|
||||
{
|
||||
CMesh::ScalarType dist=Distortion(mesh.face[i],grad_size);
|
||||
if (dist>maxD)maxD=dist;
|
||||
CMesh::ScalarType absLap=fabs(LaplaceDistortion(mesh.face[i], grad_size));
|
||||
if (absLap>maxL)maxL=absLap;
|
||||
CMesh::ScalarType stiffDelta = std::min(c * absLap, d);
|
||||
//mesh.face[i].stiffening+=stiffDelta;
|
||||
Handle_Stiffness[i]+=stiffDelta;
|
||||
}
|
||||
}
|
||||
printf("Maximum Distorsion %4.4f \n",maxD);
|
||||
printf("Maximum Laplacian %4.4f \n",maxL);
|
||||
return flipped;
|
||||
}
|
||||
|
||||
void InitDefaultStiffening()
|
||||
{
|
||||
bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
|
||||
if(!hasStiffness)
|
||||
Handle_Stiffness=vcg::tri::Allocator<CMesh>::AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));
|
||||
|
||||
for(unsigned int i=0;i<mesh.face.size();i++)
|
||||
{
|
||||
CMesh::FaceType *f=&(mesh.face[i]);
|
||||
//f->stiffening=1;
|
||||
Handle_Stiffness[f]=1;
|
||||
}
|
||||
}
|
||||
|
||||
void AddGaussStiffening(typename MeshType::ScalarType C)
|
||||
{
|
||||
int radius=floor(C);
|
||||
if (C<4)radius=4;
|
||||
AddStiffening(C,radius);
|
||||
}
|
||||
|
||||
// Mesh class
|
||||
MeshType &mesh;
|
||||
|
||||
// Quad mesh class
|
||||
//CMesh quadmesh;
|
||||
|
||||
///seams initializer
|
||||
SeamsInitializer<MeshType> SInit;
|
||||
|
||||
// Vertex indexing class used for the solver
|
||||
VertexIndexing<MeshType> VInd;
|
||||
|
||||
// Poisson solver
|
||||
PoissonSolver<MeshType> PSolver;
|
||||
|
||||
};
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,38 @@
|
|||
#ifndef AUXMATH_H
|
||||
#define AUXMATH_H
|
||||
#pragma once
|
||||
|
||||
#include <complex>
|
||||
#include <cmath>
|
||||
//#include "Claussen.h"
|
||||
|
||||
static const double pi = 3.1415926535897932384626433;
|
||||
|
||||
const double SQRT2_2 = 0.7071067811865; // sqrt(2)/2 = 1/sqrt(2)
|
||||
|
||||
const double EPSILON = 1.0e-8;
|
||||
|
||||
template <class T> T sqr( const T& a ) { return a*a; }
|
||||
|
||||
inline double l2s( const double angle ){
|
||||
const double s = 2*fabs(sin(angle));
|
||||
return s < 1e-40 ? -92.103403719761827360719658 : log(s);
|
||||
}
|
||||
//
|
||||
//inline double Lob( const double angle ){
|
||||
// return 0 <= angle && angle <= pi ? claussen(2*angle)/2 : -1e10;
|
||||
//}
|
||||
|
||||
inline double PropLength( const double l, const double an, const double ad ){
|
||||
return l*(sin(an)/sin(ad));
|
||||
}
|
||||
|
||||
inline double cot( const double x ){
|
||||
const double a = fabs(fmod(x+pi,2*pi)-pi)<1e-40 ? 1e-40 : x;
|
||||
return 1./tan(a);
|
||||
}
|
||||
|
||||
typedef std::complex<double> Cmplx;
|
||||
|
||||
#endif // AUXMATH_H
|
||||
|
||||
|
|
@ -0,0 +1,150 @@
|
|||
#ifndef MIQ_GL_UTILS
|
||||
#define MIQ_GL_UTILS
|
||||
//#include <wrap/gl/space.h>
|
||||
#include "vertex_indexing.h"
|
||||
#include "quadrangulator.h"
|
||||
|
||||
class Miq_Gl_Utils
|
||||
{
|
||||
public:
|
||||
|
||||
template <class VertexIndexingType>
|
||||
static void GLDrawVertexIndexing(VertexIndexingType &VI)
|
||||
{
|
||||
typedef typename VertexIndexingType::ScalarType ScalarType;
|
||||
|
||||
glPushAttrib(GL_ALL_ATTRIB_BITS);
|
||||
glEnable(GL_COLOR_MATERIAL);
|
||||
glDisable(GL_LIGHTING);
|
||||
glDepthRange(0,0.999);
|
||||
typename VertexIndexingType::ScalarType size=5;
|
||||
glPointSize(size);
|
||||
vcg::glColor(vcg::Color4b(0,255,0,255));
|
||||
glBegin(GL_POINTS);
|
||||
for (unsigned int i=0;i<VI.duplicated.size();i++)
|
||||
{
|
||||
assert(!VI.duplicated[i]->IsD());
|
||||
vcg::glVertex(VI.duplicated[i]->P());
|
||||
}
|
||||
glEnd();
|
||||
glPopAttrib();
|
||||
}
|
||||
|
||||
template <class MeshType>
|
||||
static void DrawFlippedFacesIfSelected(MeshType &Tmesh)
|
||||
{
|
||||
glPushAttrib(GL_ALL_ATTRIB_BITS);
|
||||
glDisable(GL_LIGHTING);
|
||||
glLineWidth(1.5);
|
||||
glDepthRange(0,0.998);
|
||||
vcg::glColor(vcg::Color4b(255,0,0,255));
|
||||
glBegin(GL_LINES);
|
||||
for (unsigned int i=0;i<Tmesh.face.size();i++)
|
||||
{
|
||||
if (!Tmesh.face[i].IsS())continue;
|
||||
for (int k=0;k<3;k++)
|
||||
{
|
||||
vcg::glVertex(Tmesh.face[i].P0(k));
|
||||
vcg::glVertex(Tmesh.face[i].P1(k));
|
||||
}
|
||||
}
|
||||
glEnd();
|
||||
glPopAttrib();
|
||||
}
|
||||
|
||||
template <class QuadrangulatorType>
|
||||
static void QuadGLDrawIntegerVertices(QuadrangulatorType &Quadr)
|
||||
{
|
||||
glPushAttrib(GL_ALL_ATTRIB_BITS);
|
||||
glDisable(GL_LIGHTING);
|
||||
glEnable(GL_COLOR_MATERIAL);
|
||||
glDisable(GL_TEXTURE_2D);
|
||||
glPointSize(8);
|
||||
|
||||
glDepthRange(0,0.997);
|
||||
/*glColor3d(1,0,0);*/
|
||||
glBegin(GL_POINTS);
|
||||
for (int i=0;i<Quadr.IntegerVertex.size();i++)
|
||||
{
|
||||
typename QuadrangulatorType::TriVertexType* v=Quadr.IntegerVertex[i];
|
||||
typename QuadrangulatorType::CoordType pos=v->P();
|
||||
if (v->IsV())
|
||||
glColor3d(1,0,0);
|
||||
else
|
||||
glColor3d(1,1,0);
|
||||
glVertex(pos);
|
||||
}
|
||||
glEnd();
|
||||
|
||||
glPopAttrib();
|
||||
}
|
||||
|
||||
template <class QuadrangulatorType>
|
||||
static void GLDrawIntegerLines(QuadrangulatorType &Quadr)
|
||||
{
|
||||
glPushAttrib(GL_ALL_ATTRIB_BITS);
|
||||
glDisable(GL_LIGHTING);
|
||||
glEnable(GL_COLOR_MATERIAL);
|
||||
glDisable(GL_TEXTURE_2D);
|
||||
glLineWidth(2);
|
||||
|
||||
glColor3d(0,1,0);
|
||||
glDepthRange(0,0.998);
|
||||
|
||||
for (int i=0;i<Quadr.IntegerLines.size();i++)
|
||||
{
|
||||
typename QuadrangulatorType::TriFaceType *f=Quadr.IntegerLines[i].first;
|
||||
int edge=Quadr.IntegerLines[i].second;
|
||||
typename QuadrangulatorType::TriVertexType* v0=f->V0(edge);
|
||||
typename QuadrangulatorType::TriVertexType* v1=f->V1(edge);
|
||||
glBegin(GL_LINES);
|
||||
glVertex(v0->P());
|
||||
glVertex(v1->P());
|
||||
glEnd();
|
||||
}
|
||||
|
||||
glPopAttrib();
|
||||
}
|
||||
|
||||
template <class QuadrangulatorType>
|
||||
static void GLDrawPolygons(QuadrangulatorType &Quadr)
|
||||
{
|
||||
glPushAttrib(GL_ALL_ATTRIB_BITS);
|
||||
glEnable(GL_LIGHTING);
|
||||
glEnable(GL_COLOR_MATERIAL);
|
||||
glDisable(GL_TEXTURE_2D);
|
||||
glColor3d(0.7,0.8,0.9);
|
||||
//glFrontFace(GL_CW);
|
||||
glDepthRange(0,0.998);
|
||||
for (unsigned int i=0;i<Quadr.polygons.size();i++)
|
||||
{
|
||||
glBegin(GL_POLYGON);
|
||||
for (unsigned int j=0;j<Quadr.polygons[i].size();j++)
|
||||
{
|
||||
typename QuadrangulatorType::TriVertexType* v=Quadr.polygons[i][j];
|
||||
glNormal(v->N());
|
||||
glVertex(v->P());
|
||||
}
|
||||
glEnd();
|
||||
}
|
||||
|
||||
glDepthRange(0,0.997);
|
||||
glDisable(GL_LIGHTING);
|
||||
glEnable(GL_COLOR_MATERIAL);
|
||||
glColor3d(0,0,0);
|
||||
for (unsigned int i=0;i<Quadr.polygons.size();i++)
|
||||
{
|
||||
glBegin(GL_LINE_LOOP);
|
||||
for (unsigned int j=0;j<Quadr.polygons[i].size();j++)
|
||||
{
|
||||
typename QuadrangulatorType::TriVertexType* v=Quadr.polygons[i][j];
|
||||
glVertex(v->P());
|
||||
}
|
||||
glEnd();
|
||||
}
|
||||
|
||||
glPopAttrib();
|
||||
}
|
||||
|
||||
};
|
||||
#endif
|
||||
|
|
@ -0,0 +1,199 @@
|
|||
#ifndef PARAM_STATS_H
|
||||
#define PARAM_STATS_H
|
||||
#ifdef MIQ_USE_ROBUST
|
||||
#include "predicates.h"
|
||||
#endif
|
||||
#include <vcg/space/triangle2.h>
|
||||
|
||||
template<typename CoordType2D>
|
||||
inline bool IsFlipped(const CoordType2D &uv0,
|
||||
const CoordType2D &uv1,
|
||||
const CoordType2D &uv2)
|
||||
{
|
||||
#ifdef USE_ROBUST
|
||||
double pa[2] = {uv0.X(), uv0.Y()};
|
||||
double pb[2] = {uv1.X(), uv1.Y()};
|
||||
double pc[2] = {uv2.X(), uv2.Y()};
|
||||
return (orient2d(pa, pb, pc) < 0);
|
||||
#else
|
||||
vcg::Triangle2<typename CoordType2D::ScalarType> t2(uv0,uv1,uv2);
|
||||
return (!t2.IsCCW());
|
||||
#endif
|
||||
}
|
||||
|
||||
template<typename FaceType>
|
||||
inline bool IsFlipped( FaceType &f)
|
||||
{
|
||||
|
||||
typedef typename FaceType::ScalarType ScalarType;
|
||||
typedef typename vcg::Point2<ScalarType> CoordType2D;
|
||||
CoordType2D uv0=f.WT(0).P();
|
||||
CoordType2D uv1=f.WT(1).P();
|
||||
CoordType2D uv2=f.WT(2).P();
|
||||
return (IsFlipped(uv0,uv1,uv2));
|
||||
|
||||
}
|
||||
|
||||
template< class MeshType>
|
||||
int NumFlips(MeshType &Tmesh)
|
||||
{
|
||||
int numFl=0;
|
||||
for (unsigned int i=0;i<Tmesh.face.size();i++)
|
||||
{
|
||||
if (Tmesh.face[i].IsD())continue;
|
||||
if (IsFlipped(Tmesh.face[i]))numFl++;
|
||||
}
|
||||
return numFl;
|
||||
}
|
||||
|
||||
template< class MeshType>
|
||||
void SelectFlippedFaces(MeshType &Tmesh)
|
||||
{
|
||||
typedef typename MeshType::ScalarType ScalarType;
|
||||
typedef typename vcg::Point2<ScalarType> CoordType2D;
|
||||
vcg::tri::UpdateFlags<MeshType>::FaceClearS(Tmesh);
|
||||
for (unsigned int i=0;i<Tmesh.face.size();i++)
|
||||
{
|
||||
CoordType2D uv0=Tmesh.face[i].WT(0).P();
|
||||
CoordType2D uv1=Tmesh.face[i].WT(1).P();
|
||||
CoordType2D uv2=Tmesh.face[i].WT(2).P();
|
||||
if (IsFlipped(uv0,uv1,uv2) ) Tmesh.face[i].SetS();
|
||||
}
|
||||
}
|
||||
|
||||
// Compute the lambda (distortion) of a triangle in the current
|
||||
// parameerization from the Mixed Integer paper.
|
||||
// f facet on which to compute distortion
|
||||
// h scaling factor applied to cross field
|
||||
// return the triangle's distortion
|
||||
template< class FaceType>
|
||||
typename FaceType::ScalarType Distortion(FaceType &f,typename FaceType::ScalarType h)
|
||||
{
|
||||
typedef typename FaceType::CoordType CoordType3D;
|
||||
typedef typename FaceType::ScalarType ScalarType;
|
||||
typedef typename vcg::Point2<ScalarType> CoordType2D;
|
||||
typedef typename FaceType::ScalarType ScalarType;
|
||||
|
||||
//Facet_const_handle self = f;
|
||||
assert(h > 0);
|
||||
|
||||
//TriangleIndex ftri = m_flattenedTriangles[f.index()];
|
||||
//TriangleIndex tri = m_triangles[f.index()];
|
||||
CoordType2D uv0 = f.WT(0).P();
|
||||
CoordType2D uv1 = f.WT(1).P();
|
||||
CoordType2D uv2 = f.WT(2).P();
|
||||
|
||||
CoordType3D p0 = f.P(0);
|
||||
CoordType3D p1 = f.P(1);
|
||||
CoordType3D p2 = f.P(2);
|
||||
|
||||
CoordType3D norm = (p1 - p0) ^ (p2 - p0);
|
||||
ScalarType area2 = (norm).Norm();
|
||||
ScalarType area2_inv = 1.0 / area2;
|
||||
norm *= area2_inv;
|
||||
|
||||
if (area2 > 0) {
|
||||
// Singular values of the Jacobian
|
||||
CoordType3D neg_t0 = norm ^ (p2 - p1);
|
||||
CoordType3D neg_t1 = norm ^ ( p0 - p2);
|
||||
CoordType3D neg_t2 = norm ^ ( p1 - p0);
|
||||
|
||||
/*CoordType3D diffu = area2_inv * (uv0.X() * neg_t0 +
|
||||
uv1.X() * neg_t1 + uv2.X() * neg_t2);
|
||||
CoordType3D diffv = area2_inv * (uv0.Y() * neg_t0 +
|
||||
uv1.Y() * neg_t1 + uv2.Y() * neg_t2);*/
|
||||
|
||||
CoordType3D diffu = (neg_t0 * uv0.X() +neg_t1 *uv1.X() + neg_t2 * uv2.X() )*area2_inv;
|
||||
CoordType3D diffv = (neg_t0 * uv0.Y() + neg_t1*uv1.Y() + neg_t2*uv2.Y() )*area2_inv;
|
||||
// first fundamental form
|
||||
ScalarType I00 = diffu*diffu; // guaranteed non-neg
|
||||
ScalarType I01 = diffu*diffv; // I01 = I10
|
||||
ScalarType I11 = diffv*diffv; // guaranteed non-neg
|
||||
// eigenvalues of a 2x2 matrix
|
||||
// [a00 a01]
|
||||
// [a10 a11]
|
||||
// 1/2 * [ (a00 + a11) +/- sqrt((a00 - a11)^2 + 4 a01 a10) ]
|
||||
ScalarType trI = I00 + I11; // guaranteed non-neg
|
||||
ScalarType diffDiag = I00 - I11; // guaranteed non-neg
|
||||
ScalarType sqrtDet = sqrt(std::max(0.0, diffDiag*diffDiag +
|
||||
4 * I01 * I01)); // guaranteed non-neg
|
||||
ScalarType sig1 = 0.5 * (trI + sqrtDet); // higher singular value
|
||||
ScalarType sig2 = 0.5 * (trI - sqrtDet); // lower singular value
|
||||
|
||||
// Avoid sig2 < 0 due to numerical error
|
||||
if (fabs(sig2) < 1.0e-8) sig2 = 0;
|
||||
assert(sig1 >= 0);
|
||||
assert(sig2 >= 0);
|
||||
|
||||
if (sig2 < 0) {
|
||||
printf("Distortion will be NaN! sig1^2 is negative (%lg)\n",
|
||||
sig2);
|
||||
}
|
||||
|
||||
// The singular values of the Jacobian are the sqrts of the
|
||||
// eigenvalues of the first fundamental form.
|
||||
sig1 = sqrt(sig1);
|
||||
sig2 = sqrt(sig2);
|
||||
|
||||
// distortion
|
||||
ScalarType tao = IsFlipped(f) ? -1 : 1;
|
||||
ScalarType factor = tao / h;
|
||||
ScalarType lam = fabs(factor * sig1 - 1) + fabs(factor * sig2 - 1);
|
||||
return lam;
|
||||
}
|
||||
else {
|
||||
return 10; // something "large"
|
||||
}
|
||||
}
|
||||
|
||||
////////////////////////////////////////////////////////////////////////////
|
||||
// Approximate the distortion laplacian using a uniform laplacian on
|
||||
// the dual mesh:
|
||||
// ___________
|
||||
// \-1 / \-1 /
|
||||
// \ / 3 \ /
|
||||
// \-----/
|
||||
// \-1 /
|
||||
// \ /
|
||||
//
|
||||
// @param[in] f facet on which to compute distortion laplacian
|
||||
// @param[in] h scaling factor applied to cross field
|
||||
// @return distortion laplacian for f
|
||||
///////////////////////////////////////////////////////////////////////////
|
||||
template< class FaceType>
|
||||
typename FaceType::ScalarType LaplaceDistortion(FaceType &f ,typename FaceType::ScalarType h)
|
||||
{
|
||||
typedef typename FaceType::ScalarType ScalarType;
|
||||
ScalarType mydist = Distortion(f, h);
|
||||
ScalarType lapl=0;
|
||||
for (int i=0;i<3;i++)
|
||||
lapl += (mydist- Distortion(*f.FFp(i), h));
|
||||
return lapl;
|
||||
}
|
||||
|
||||
template< class MeshType>
|
||||
void SetFaceQualityByDistortion(MeshType &Tmesh,
|
||||
typename MeshType::ScalarType h)
|
||||
{
|
||||
typedef typename MeshType::FaceType FaceType;
|
||||
typedef typename MeshType::ScalarType ScalarType;
|
||||
ScalarType minD=100;
|
||||
ScalarType maxD=0;
|
||||
|
||||
///evaluate min and max
|
||||
for (unsigned int i=0;i<Tmesh.face.size();i++)
|
||||
{
|
||||
ScalarType dist=Distortion<FaceType>(Tmesh.face[i],h);
|
||||
if (dist>maxD)maxD=dist;
|
||||
if (dist<minD)minD=dist;
|
||||
Tmesh.face[i].Q()= dist;
|
||||
}
|
||||
for (unsigned int i=0;i<Tmesh.face.size();i++)
|
||||
{
|
||||
Tmesh.face[i].Q()= maxD-Tmesh.face[i].Q();
|
||||
}
|
||||
printf("Distortion MIN %4.4f \n",(float)minD);
|
||||
printf("Distortion MAX %4.4f \n",(float)maxD);
|
||||
}
|
||||
|
||||
#endif // PARAM_STATS_H
|
||||
|
|
@ -0,0 +1,847 @@
|
|||
#ifndef MIQ_POISSON_SOLVER
|
||||
#define MIQ_POISSON_SOLVER
|
||||
|
||||
#include <gmm/gmm.h>
|
||||
#include <ConstrainedSolver.hh>
|
||||
#include <MISolver.hh>
|
||||
#include <GMM_Tools.hh>
|
||||
|
||||
#include "auxmath.h"
|
||||
#include "sparsesystemdata.h"
|
||||
#include "mesh_type.h"
|
||||
#include "vertex_indexing.h"
|
||||
|
||||
template <class MeshType>
|
||||
class PoissonSolver
|
||||
{
|
||||
typedef typename MeshType::ScalarType ScalarType;
|
||||
typedef typename MeshType::FaceType FaceType;
|
||||
typedef typename MeshType::VertexType VertexType;
|
||||
typedef typename MeshType::CoordType CoordType;
|
||||
//typedef VertexIndexing<MeshType> VertexIndexingType;
|
||||
//typedef typename VertexIndexingType::SeamInfo SeamInfo;
|
||||
|
||||
typename MeshType::template PerFaceAttributeHandle<vcg::Point3i> HandleS_Index;
|
||||
typename MeshType::template PerVertexAttributeHandle<bool> Handle_Singular;
|
||||
typename MeshType::template PerFaceAttributeHandle<float> Handle_Stiffness;
|
||||
///this handle for mesh
|
||||
typename MeshType::template PerMeshAttributeHandle<MeshSystemInfo> Handle_SystemInfo;
|
||||
|
||||
///range map data
|
||||
MeshType &mesh;
|
||||
|
||||
///vertex indexing structure
|
||||
//VertexIndexingType &VIndex;
|
||||
|
||||
///solver data
|
||||
SparseSystemData S;
|
||||
|
||||
///vector of unknowns
|
||||
std::vector< double > X;
|
||||
|
||||
////REAL PART
|
||||
///number of fixed vertex
|
||||
unsigned int n_fixed_vars;
|
||||
///the number of REAL variables for vertices
|
||||
unsigned int n_vert_vars;
|
||||
///total number of variables of the system,
|
||||
///do not consider constraints, but consider integer vars
|
||||
unsigned int num_total_vars;
|
||||
|
||||
///the number of scalar variables
|
||||
//unsigned int n_scalar_vars;
|
||||
|
||||
//////INTEGER PART
|
||||
///the total number of integer variables
|
||||
unsigned int n_integer_vars;
|
||||
|
||||
///CONSTRAINT PART
|
||||
///number of cuts constraints
|
||||
unsigned int num_cut_constraint;
|
||||
|
||||
///total number of constraints equations
|
||||
unsigned int num_constraint_equations;
|
||||
|
||||
///total size of the system including constraints
|
||||
unsigned int system_size;
|
||||
|
||||
///if you intend to make integer rotation
|
||||
///and translations
|
||||
bool integer_jumps_bary;
|
||||
|
||||
///vector of blocked vertices
|
||||
std::vector<VertexType*> Hard_constraints;
|
||||
|
||||
///vector of indexes to round
|
||||
std::vector<int> ids_to_round;
|
||||
|
||||
///boolean that is true if rounding to integer is needed
|
||||
bool integer_rounding;
|
||||
|
||||
///START SYSTEM ACCESS METHODS
|
||||
///add an entry to the LHS
|
||||
void AddValA(int Xindex,
|
||||
int Yindex,
|
||||
ScalarType val)
|
||||
{
|
||||
int size=(int)S.nrows();
|
||||
assert(0 <= Xindex && Xindex < size);
|
||||
assert(0 <= Yindex && Yindex < size);
|
||||
S.A().addEntryReal(Xindex,Yindex,val);
|
||||
}
|
||||
|
||||
///add a complex entry to the LHS
|
||||
void AddComplexA(int VarXindex,
|
||||
int VarYindex,
|
||||
Cmplx val)
|
||||
{
|
||||
int size=(int)S.nrows()/2;
|
||||
assert(0 <= VarXindex && VarXindex < size);
|
||||
assert(0 <= VarYindex && VarYindex < size);
|
||||
S.A().addEntryCmplx(VarXindex,VarYindex,val);
|
||||
}
|
||||
|
||||
///add a velue to the RHS
|
||||
void AddValB(int Xindex,
|
||||
ScalarType val)
|
||||
{
|
||||
int size=(int)S.nrows();
|
||||
assert(0 <= Xindex && Xindex < size);
|
||||
S.b()[Xindex] += val;
|
||||
}
|
||||
|
||||
///add the area term, scalefactor is used to sum up
|
||||
///and normalize on the overlap zones
|
||||
void AddAreaTerm(int index[3][3][2],ScalarType ScaleFactor)
|
||||
{
|
||||
const ScalarType entry=0.5*ScaleFactor;
|
||||
ScalarType val[3][3]= { {0, entry, -entry},
|
||||
{-entry, 0, entry},
|
||||
{entry, -entry, 0} };
|
||||
|
||||
for (int i=0;i<3;i++)
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
///add for both u and v
|
||||
int Xindex=index[i][j][0]*2;
|
||||
int Yindex=index[i][j][1]*2;
|
||||
|
||||
AddValA(Xindex+1,Yindex,-val[i][j]);
|
||||
AddValA(Xindex,Yindex+1,val[i][j]);
|
||||
}
|
||||
}
|
||||
|
||||
///set the diagonal of the matrix (which is zero at the beginning)
|
||||
///such that the sum of a row or a colums is zero
|
||||
void SetDiagonal(ScalarType val[3][3])
|
||||
{
|
||||
for (int i=0;i<3;i++)
|
||||
{
|
||||
ScalarType sum=0;
|
||||
for (int j=0;j<3;j++)
|
||||
sum+=val[i][j];
|
||||
val[i][i]=-sum;
|
||||
}
|
||||
}
|
||||
|
||||
///given a vector of scalar values and
|
||||
///a vector of indexes add such values
|
||||
///as specified by the indexes
|
||||
void AddRHS(ScalarType b[6],
|
||||
int index[3])
|
||||
{
|
||||
for (int i=0;i<3;i++)
|
||||
{
|
||||
ScalarType valU=b[i*2];
|
||||
ScalarType valV=b[(i*2)+1];
|
||||
AddValB((index[i]*2),valU);
|
||||
AddValB((index[i]*2)+1,valV);
|
||||
}
|
||||
}
|
||||
|
||||
///add a 3x3 block matrix to the system matrix...
|
||||
///indexes are specified in the 3x3 matrix of x,y pairs
|
||||
///indexes must be multiplied by 2 cause u and v
|
||||
void Add33Block(ScalarType val[3][3],int index[3][3][2])
|
||||
{
|
||||
for (int i=0;i<3;i++)
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
///add for both u and v
|
||||
int Xindex=index[i][j][0]*2;
|
||||
int Yindex=index[i][j][1]*2;
|
||||
assert((unsigned)Xindex<(n_vert_vars*2));
|
||||
assert((unsigned)Yindex<(n_vert_vars*2));
|
||||
AddValA(Xindex,Yindex,val[i][j]);
|
||||
AddValA(Xindex+1,Yindex+1,val[i][j]);
|
||||
}
|
||||
}
|
||||
|
||||
///add a 3x3 block matrix to the system matrix...
|
||||
///indexes are specified in the 3x3 matrix of x,y pairs
|
||||
///indexes must be multiplied by 2 cause u and v
|
||||
void Add44Block(ScalarType val[4][4],int index[4][4][2])
|
||||
{
|
||||
for (int i=0;i<4;i++)
|
||||
for (int j=0;j<4;j++)
|
||||
{
|
||||
///add for both u and v
|
||||
int Xindex=index[i][j][0]*2;
|
||||
int Yindex=index[i][j][1]*2;
|
||||
assert((unsigned)Xindex<(n_vert_vars*2));
|
||||
assert((unsigned)Yindex<(n_vert_vars*2));
|
||||
AddValA(Xindex,Yindex,val[i][j]);
|
||||
AddValA(Xindex+1,Yindex+1,val[i][j]);
|
||||
}
|
||||
}
|
||||
///END SYSTEM ACCESS METHODS
|
||||
|
||||
///START COMMON MATH FUNCTIONS
|
||||
///return the complex encoding the rotation
|
||||
///for a given missmatch interval
|
||||
Cmplx GetRotationComplex(int interval)
|
||||
{
|
||||
assert((interval>=0)&&(interval<4));
|
||||
|
||||
switch(interval)
|
||||
{
|
||||
case 0:return Cmplx(1,0);
|
||||
case 1:return Cmplx(0,1);
|
||||
case 2:return Cmplx(-1,0);
|
||||
default:return Cmplx(0,-1);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
vcg::Point2i Rotate(vcg::Point2i p,int interval)
|
||||
{
|
||||
assert(interval>=0);
|
||||
assert(interval<4);
|
||||
Cmplx rot=GetRotationComplex(interval);
|
||||
/*
|
||||
| real -imag| |p.x|
|
||||
| |*| |
|
||||
| imag real | |p.y|
|
||||
*/
|
||||
vcg::Point2i ret;
|
||||
ret.X()=rot.real()*p.X()-rot.imag()*p.Y();
|
||||
ret.Y()=rot.imag()*p.X()+rot.real()*p.Y();
|
||||
return (ret);
|
||||
}
|
||||
///END COMMON MATH FUNCTIONS
|
||||
|
||||
|
||||
///START ENERGY MINIMIZATION PART
|
||||
///initialize the LHS for a given face
|
||||
///for minimization of Dirichlet's energy
|
||||
void perElementLHS(FaceType *f,
|
||||
ScalarType val[3][3],
|
||||
int index[3][3][2])
|
||||
{
|
||||
///initialize to zero
|
||||
for (int x=0;x<3;x++)
|
||||
for (int y=0;y<3;y++)
|
||||
val[x][y]=0;
|
||||
|
||||
///get the vertices
|
||||
VertexType *v[3];
|
||||
v[0]=f->V(0);
|
||||
v[1]=f->V(1);
|
||||
v[2]=f->V(2);
|
||||
|
||||
///get the indexes of vertex instance (to consider cuts)
|
||||
///for the current face
|
||||
int Vindexes[3];
|
||||
//Vindexes[0]=f->syst_index[0];
|
||||
//Vindexes[1]=f->syst_index[1];
|
||||
//Vindexes[2]=f->syst_index[2];
|
||||
Vindexes[0]=HandleS_Index[f][0];
|
||||
Vindexes[1]=HandleS_Index[f][1];
|
||||
Vindexes[2]=HandleS_Index[f][2];
|
||||
//VIndex.getIndexInfo(f,Vindexes);
|
||||
|
||||
///initialize the indexes for the block
|
||||
for (int x=0;x<3;x++)
|
||||
for (int y=0;y<3;y++)
|
||||
{
|
||||
index[x][y][0]=Vindexes[x];
|
||||
index[x][y][1]=Vindexes[y];
|
||||
}
|
||||
|
||||
///initialize edges
|
||||
CoordType e[3];
|
||||
for (int k=0;k<3;k++)
|
||||
e[k]=v[(k+2)%3]->P()-v[(k+1)%3]->P();
|
||||
|
||||
///then consider area but also considering scale factor dur to overlaps
|
||||
ScalarType areaT=((f->P(1)-f->P(0))^(f->P(2)-f->P(0))).Norm()/2.0;
|
||||
//ScalarType ScaleFactor=f->area;
|
||||
for (int x=0;x<3;x++)
|
||||
for (int y=0;y<3;y++)
|
||||
if (x!=y)
|
||||
{
|
||||
ScalarType num=(e[x]*e[y]);
|
||||
val[x][y] =num/(4.0*areaT);//*ScaleFactor);//*(ScalarType)convex_size);
|
||||
val[x][y]*=Handle_Stiffness[f];//f->stiffening;
|
||||
//val[x][y]*=ScaleFactor;
|
||||
}
|
||||
|
||||
///set the matrix as diagonal
|
||||
SetDiagonal(val);
|
||||
}
|
||||
|
||||
///initialize the RHS for a given face
|
||||
///for minimization of Dirichlet's energy
|
||||
void perElementRHS(FaceType *f,
|
||||
ScalarType b[6],
|
||||
ScalarType vector_field_scale=1)
|
||||
{
|
||||
|
||||
/// then set the rhs
|
||||
CoordType scaled_Kreal;
|
||||
CoordType scaled_Kimag;
|
||||
CoordType fNorm=f->N();
|
||||
CoordType p[3];
|
||||
p[0]=f->P0(0);
|
||||
p[1]=f->P0(1);
|
||||
p[2]=f->P0(2);
|
||||
|
||||
CoordType neg_t[3];
|
||||
neg_t[0] = fNorm ^ (p[2] - p[1]);
|
||||
neg_t[1] = fNorm ^ (p[0] - p[2]);
|
||||
neg_t[2] = fNorm ^ (p[1] - p[0]);
|
||||
|
||||
CoordType K1,K2;
|
||||
//K1=CrossVector<FaceType>(*f,0);
|
||||
//K2=CrossVector<FaceType>(*f,1);
|
||||
K1=f->PD1();
|
||||
K2=f->PD2();
|
||||
|
||||
scaled_Kreal = K1*(vector_field_scale)/2;
|
||||
scaled_Kimag = K2*(vector_field_scale)/2;
|
||||
|
||||
ScalarType stiff_val=(ScalarType)Handle_Stiffness[f];
|
||||
b[0] = scaled_Kreal * neg_t[0]*stiff_val;
|
||||
b[1] = scaled_Kimag * neg_t[0]*stiff_val;
|
||||
b[2] = scaled_Kreal * neg_t[1]*stiff_val;
|
||||
b[3] = scaled_Kimag * neg_t[1]*stiff_val;
|
||||
b[4] = scaled_Kreal * neg_t[2]*stiff_val;
|
||||
b[5] = scaled_Kimag * neg_t[2]*stiff_val;
|
||||
}
|
||||
|
||||
///evaluate the LHS and RHS for a single face
|
||||
///for minimization of Dirichlet's energy
|
||||
void PerElementSystemReal(FaceType *f,
|
||||
ScalarType val[3][3],
|
||||
int index[3][3][2],
|
||||
ScalarType b[6],
|
||||
ScalarType vector_field_scale=1.0)
|
||||
{
|
||||
perElementLHS(f,val,index);
|
||||
perElementRHS(f,b,vector_field_scale);
|
||||
}
|
||||
///END ENERGY MINIMIZATION PART
|
||||
|
||||
///START FIXING VERTICES
|
||||
///set a given vertex as fixed
|
||||
void AddFixedVertex(VertexType *v)
|
||||
{
|
||||
n_fixed_vars++;
|
||||
Hard_constraints.push_back(v);
|
||||
//v->blocked=true;
|
||||
}
|
||||
|
||||
///find vertex to fix in case we're using
|
||||
///a vector field NB: multiple components not handled
|
||||
void FindFixedVertField()
|
||||
{
|
||||
Hard_constraints.clear();
|
||||
|
||||
n_fixed_vars=0;
|
||||
///fix the first singularity
|
||||
for (unsigned int j=0;j<mesh.vert.size();j++)
|
||||
{
|
||||
VertexType *v=&mesh.vert[j];
|
||||
if (v->IsD())continue;
|
||||
//if (v->IsSingular())
|
||||
if (Handle_Singular[v])
|
||||
{
|
||||
AddFixedVertex(v);
|
||||
v->T().P()=vcg::Point2<ScalarType>(0,0);
|
||||
return;
|
||||
}
|
||||
}
|
||||
///if anything fixed fix the first
|
||||
AddFixedVertex(&mesh.vert[0]);
|
||||
mesh.vert[0].T().P()=vcg::Point2<ScalarType>(0,0);
|
||||
}
|
||||
|
||||
///find hard constraint depending if using or not
|
||||
///a vector field
|
||||
void FindFixedVert()
|
||||
{
|
||||
Hard_constraints.clear();
|
||||
FindFixedVertField();
|
||||
}
|
||||
|
||||
int GetFirstVertexIndex(VertexType *v)
|
||||
{
|
||||
FaceType *f=v->VFp();
|
||||
int index=v->VFi();
|
||||
//return (f->syst_index[index]);
|
||||
return (HandleS_Index[f][index]);
|
||||
//int indexV[3];
|
||||
//VIndex.getIndexInfo(f,indexV);
|
||||
//return indexV[index];
|
||||
}
|
||||
|
||||
///fix the vertices which are flagged as fixed
|
||||
void FixBlockedVertex()
|
||||
{
|
||||
int offset_row=n_vert_vars+num_cut_constraint;
|
||||
|
||||
unsigned int constr_num = 0;
|
||||
for (unsigned int i=0;i<Hard_constraints.size();i++)
|
||||
{
|
||||
VertexType *v=Hard_constraints[i];
|
||||
assert(!v->IsD());
|
||||
///get first index of the vertex that must blocked
|
||||
//int index=v->vertex_index[0];
|
||||
int index=GetFirstVertexIndex(v);
|
||||
///multiply times 2 because of uv
|
||||
int indexvert=index*2;
|
||||
///find the first free row to add the constraint
|
||||
int indexRow=(offset_row+constr_num)*2;
|
||||
int indexCol=indexRow;
|
||||
|
||||
///add fixing constraint LHS
|
||||
AddValA(indexRow,indexvert,1);
|
||||
AddValA(indexRow+1,indexvert+1,1);
|
||||
|
||||
///add fixing constraint RHS
|
||||
//AddValB(indexCol,(ScalarType)v->range_index.X());
|
||||
//AddValB(indexCol+1,(ScalarType)v->range_index.Y());
|
||||
AddValB(indexCol,(ScalarType)v->T().P().X());
|
||||
AddValB(indexCol+1,(ScalarType)v->T().P().Y());
|
||||
//AddValB(indexCol,0);
|
||||
//AddValB(indexCol+1,0);
|
||||
constr_num++;
|
||||
}
|
||||
assert(constr_num==n_fixed_vars);
|
||||
}
|
||||
///END FIXING VERTICES
|
||||
|
||||
///HANDLING SINGULARITY
|
||||
//set the singularity round to integer location
|
||||
void AddSingularityRound()
|
||||
{
|
||||
for (unsigned int j=0;j<mesh.vert.size();j++)
|
||||
{
|
||||
VertexType *v=&mesh.vert[j];
|
||||
if (v->IsD())continue;
|
||||
//if (v->IsSingular())
|
||||
if (Handle_Singular[v])
|
||||
{
|
||||
//assert(v->vertex_index.size()==1);
|
||||
//ids_to_round.push_back(v->vertex_index[0]*2);
|
||||
//ids_to_round.push_back((v->vertex_index[0]*2)+1);
|
||||
int index0=GetFirstVertexIndex(v);
|
||||
ids_to_round.push_back(index0*2);
|
||||
ids_to_round.push_back((index0*2)+1);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
///START GENERIC SYSTEM FUNCTIONS
|
||||
//build the laplacian matrix cyclyng over all rangemaps
|
||||
//and over all faces
|
||||
void BuildLaplacianMatrix(double vfscale=1)
|
||||
{
|
||||
|
||||
///then for each face
|
||||
for (unsigned int j=0;j<mesh.face.size();j++)
|
||||
{
|
||||
|
||||
FaceType *f=&mesh.face[j];
|
||||
if (f->IsD())
|
||||
continue;
|
||||
|
||||
int var_idx[3];//vertex variable indices
|
||||
//for(int k = 0; k < 3; ++k)
|
||||
// var_idx[k] = f->syst_index[k];
|
||||
for(int k = 0; k < 3; ++k)
|
||||
var_idx[k] = HandleS_Index[f][k];
|
||||
|
||||
//VIndex.getIndexInfo(f,var_idx);
|
||||
|
||||
///block of variables
|
||||
ScalarType val[3][3];
|
||||
///block of vertex indexes
|
||||
int index[3][3][2];
|
||||
///righe hand side
|
||||
ScalarType b[6];
|
||||
///compute the system for the given face
|
||||
PerElementSystemReal(f, val,index, b, vfscale);
|
||||
|
||||
//Add the element to the matrix
|
||||
Add33Block(val,index);
|
||||
|
||||
///add area term.. to test if needed
|
||||
//if (!use_direction_field)
|
||||
// AddAreaTerm(index,1.0);//f->area);
|
||||
|
||||
///add right hand side
|
||||
//if (use_direction_field)
|
||||
AddRHS(b,var_idx);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
///find different sized of the system
|
||||
void FindSizes()
|
||||
{
|
||||
///find the vertex that need to be fixed
|
||||
FindFixedVert();
|
||||
|
||||
///REAL PART
|
||||
n_vert_vars=Handle_SystemInfo().num_vert_variables;//VIndex.NumVertexVariables();
|
||||
|
||||
///INTEGER PART
|
||||
///the total number of integer variables
|
||||
n_integer_vars=Handle_SystemInfo().num_integer_cuts;//VIndex.NumInteger();
|
||||
|
||||
///CONSTRAINT PART
|
||||
num_cut_constraint=Handle_SystemInfo().EdgeSeamInfo.size()*2;//VIndex.NumCutsConstraints();
|
||||
|
||||
|
||||
num_constraint_equations=num_cut_constraint+n_fixed_vars;
|
||||
|
||||
///total variable of the system
|
||||
num_total_vars=n_vert_vars+n_integer_vars;
|
||||
|
||||
///initialize matrix size
|
||||
//get max constraint-integer size
|
||||
int MaxAddSize=std::max(n_integer_vars,num_constraint_equations);
|
||||
|
||||
system_size = (n_vert_vars+MaxAddSize)*2;
|
||||
|
||||
printf("\n*** SYSTEM VARIABLES *** \n");
|
||||
printf("* NUM REAL VERTEX VARIABLES %d \n",n_vert_vars);
|
||||
|
||||
printf("\n*** SINGULARITY *** \n ");
|
||||
printf("* NUM SINGULARITY %d\n",(int)ids_to_round.size()/2);
|
||||
|
||||
printf("\n*** INTEGER VARIABLES *** \n");
|
||||
printf("* NUM INTEGER VARIABLES %d \n",(int)n_integer_vars);
|
||||
|
||||
printf("\n*** CONSTRAINTS *** \n ");
|
||||
printf("* NUM FIXED CONSTRAINTS %d\n",n_fixed_vars);
|
||||
printf("* NUM CUTS CONSTRAINTS %d\n",num_cut_constraint);
|
||||
|
||||
printf("\n*** TOTAL SIZE *** \n");
|
||||
printf("* TOTAL VARIABLE SIZE (WITH INTEGER TRASL) %d \n",num_total_vars);
|
||||
printf("* TOTAL CONSTRAINTS %d \n",num_constraint_equations);
|
||||
printf("* MATRIX SIZE %d \n",system_size);
|
||||
}
|
||||
|
||||
void AllocateSystem()
|
||||
{
|
||||
|
||||
S.initialize(system_size, system_size);
|
||||
printf("\n INITIALIZED SPARSE MATRIX OF %d x %d \n",system_size, system_size);
|
||||
// number of nonzero entries in the matrix
|
||||
int num_facet=mesh.fn;
|
||||
//unsigned int nentries_vert = 4*6*6*2*n_scalar_vars; // 6x6 for each facet
|
||||
unsigned int nentries_vert = 6*6*num_facet*2;
|
||||
unsigned int nentries_fixed = 2*n_fixed_vars; // 1 complex variable fixed in each constraint, i.e. 2 regular
|
||||
unsigned int nentries_transition = 2*6*(n_integer_vars); // 3 complex variables involved in each constraint = 8 regular
|
||||
unsigned int nentries_constraints = 4*4*6*(num_constraint_equations); //
|
||||
|
||||
|
||||
int total_reserve=nentries_vert
|
||||
+ nentries_fixed
|
||||
+ nentries_transition
|
||||
+nentries_constraints;
|
||||
|
||||
printf("\n*** SPACE ALLOCATION *** \n");
|
||||
printf("* number of reserved vertices variables %d \n",nentries_vert);
|
||||
printf("* number of reserved entries fixed %d \n",nentries_fixed);
|
||||
printf("* number of reserved entries integer %d \n",nentries_transition);
|
||||
printf("* number of reserved entries constraints %d \n",nentries_constraints);
|
||||
printf("* total number of reserved entries %d \n",total_reserve);
|
||||
S.A().reserve(2*total_reserve);
|
||||
|
||||
printf("\n*** ALLOCATED *** \n");
|
||||
|
||||
}
|
||||
|
||||
///intitialize the whole matrix
|
||||
void InitMatrix()
|
||||
{
|
||||
///find singularities that must be rounded
|
||||
//AddSingularityRound();
|
||||
FindSizes();
|
||||
AllocateSystem();
|
||||
}
|
||||
|
||||
///map back coordinates after that
|
||||
///the system has been solved
|
||||
void MapCoords()
|
||||
{
|
||||
///map coords to faces
|
||||
for (unsigned int j=0;j<mesh.face.size();j++)
|
||||
{
|
||||
FaceType *f=&mesh.face[j];
|
||||
if (f->IsD())continue;
|
||||
//int indexV[3];
|
||||
//VIndex.getIndexInfo(f,indexV);
|
||||
for (int k=0;k<3;k++)
|
||||
{
|
||||
//get the index of the variable in the system
|
||||
//int indexUV=indexV[k];//f->syst_index[k];
|
||||
int indexUV=HandleS_Index[f][k];
|
||||
///then get U and V coords
|
||||
double U=X[indexUV*2];
|
||||
double V=X[indexUV*2+1];
|
||||
vcg::Point2<ScalarType> uv=vcg::Point2<ScalarType>(U,V);
|
||||
///assing
|
||||
//f->realUV[k]=uv;
|
||||
ScalarType factor=(ScalarType)SIZEQUADS/(ScalarType)SIZEPARA;
|
||||
uv*=factor;
|
||||
///assing
|
||||
f->WT(k).P()=uv;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
///END GENERIC SYSTEM FUNCTIONS
|
||||
|
||||
///set the constraints for the inter-range cuts
|
||||
void BuildSeamConstraintsExplicitTranslation()
|
||||
{
|
||||
///add constraint(s) for every seam edge (not halfedge)
|
||||
int offset_row=n_vert_vars;
|
||||
///current constraint row
|
||||
int constr_row=offset_row;
|
||||
///current constraint
|
||||
unsigned int constr_num = 0;
|
||||
///get the info for seams
|
||||
//std::vector<SeamInfo> EdgeSeamInfo;
|
||||
//VIndex.GetSeamInfo(EdgeSeamInfo);
|
||||
for (unsigned int i=0;i<Handle_SystemInfo().EdgeSeamInfo.size();i++)
|
||||
{
|
||||
//if (i<25)continue;
|
||||
unsigned char interval=Handle_SystemInfo().EdgeSeamInfo[i].MMatch;
|
||||
if (interval==1)
|
||||
interval=3;
|
||||
else
|
||||
if(interval==3)
|
||||
interval=1;
|
||||
//printf("%d\n",interval);
|
||||
int p0 = Handle_SystemInfo().EdgeSeamInfo[i].v0;
|
||||
int p1 = Handle_SystemInfo().EdgeSeamInfo[i].v1;
|
||||
int p0p = Handle_SystemInfo().EdgeSeamInfo[i].v0p;
|
||||
int p1p = Handle_SystemInfo().EdgeSeamInfo[i].v1p;
|
||||
/*assert(p1!=p1p);
|
||||
assert(p0!=p0p);*/
|
||||
Cmplx rot=GetRotationComplex(interval);
|
||||
|
||||
///get the integer variable
|
||||
int integerVar=offset_row+Handle_SystemInfo().EdgeSeamInfo[i].integerVal;
|
||||
|
||||
//Cmplx rotInt=GetRotationComplex(interval);
|
||||
if (integer_rounding)
|
||||
{
|
||||
ids_to_round.push_back(integerVar*2);
|
||||
ids_to_round.push_back(integerVar*2+1);
|
||||
}
|
||||
|
||||
AddComplexA(constr_row, p0 , rot);
|
||||
AddComplexA(constr_row, p0p, -1);
|
||||
///then translation...considering the rotation
|
||||
///due to substitution
|
||||
AddComplexA(constr_row, integerVar, 1);
|
||||
|
||||
AddValB(2*constr_row,0);
|
||||
AddValB(2*constr_row+1,0);
|
||||
constr_row +=1;
|
||||
constr_num++;
|
||||
|
||||
AddComplexA(constr_row, p1, rot);
|
||||
AddComplexA(constr_row, p1p, -1);
|
||||
///other translation
|
||||
AddComplexA(constr_row, integerVar , 1);
|
||||
|
||||
AddValB(2*constr_row,0);
|
||||
AddValB(2*constr_row+1,0);
|
||||
|
||||
constr_row +=1;
|
||||
constr_num++;
|
||||
|
||||
//// r p1 - p1' + t = 0
|
||||
}
|
||||
//assert(constr_num==num_cut_constraint);
|
||||
}
|
||||
|
||||
|
||||
|
||||
///call of the mixed integer solver
|
||||
void MixedIntegerSolve(double cone_grid_res=1,
|
||||
bool direct_round=true,
|
||||
int localIter=0)
|
||||
{
|
||||
X=std::vector< double >((n_vert_vars+n_integer_vars)*2);
|
||||
|
||||
///variables part
|
||||
int ScalarSize=(n_vert_vars)*2;
|
||||
int SizeMatrix=(n_vert_vars+n_integer_vars)*2;
|
||||
printf("\n ALLOCATED X \n");
|
||||
|
||||
///matrix A
|
||||
gmm::col_matrix< gmm::wsvector< double > > A(SizeMatrix,SizeMatrix); // lhs matrix variables +
|
||||
///constraints part
|
||||
int CsizeX=num_constraint_equations*2;
|
||||
int CsizeY=SizeMatrix+1;
|
||||
gmm::row_matrix< gmm::wsvector< double > > C(CsizeX,CsizeY); // constraints
|
||||
printf("\n ALLOCATED QMM STRUCTURES \n");
|
||||
|
||||
std::vector< double > rhs(SizeMatrix,0); // rhs
|
||||
printf("\n ALLOCATED RHS STRUCTURES \n");
|
||||
//// copy LHS
|
||||
for(int i = 0; i < (int)S.A().nentries(); ++i)
|
||||
{
|
||||
int row = S.A().rowind()[i];
|
||||
int col = S.A().colind()[i];
|
||||
int size=(int)S.nrows();
|
||||
assert(0 <= row && row < size);
|
||||
assert(0 <= col && col < size);
|
||||
|
||||
// it's either part of the matrix
|
||||
if (row < ScalarSize) {
|
||||
A(row, col) += S.A().vals()[i];
|
||||
}
|
||||
// or it's a part of the constraint
|
||||
else
|
||||
{
|
||||
//if (row<(n_vert_vars+num_constraint_equations)*2)
|
||||
assert ((unsigned int)row<(n_vert_vars+num_constraint_equations)*2);
|
||||
//{
|
||||
int r = row - ScalarSize;
|
||||
assert(r<CsizeX);
|
||||
assert(col<CsizeY);
|
||||
C(r , col ) += S.A().vals()[i];
|
||||
//}
|
||||
}
|
||||
}
|
||||
printf("\n SET %d INTEGER VALUES \n",n_integer_vars);
|
||||
///add penalization term for integer variables
|
||||
double penalization=0.000001;
|
||||
int offline_index=ScalarSize;
|
||||
for(unsigned int i = 0; i < (n_integer_vars)*2; ++i)
|
||||
{
|
||||
int index=offline_index+i;
|
||||
A(index,index)=penalization;
|
||||
}
|
||||
printf("\n SET RHS \n");
|
||||
// copy RHS
|
||||
for(int i = 0; i < (int)ScalarSize; ++i)
|
||||
{
|
||||
rhs[i] = S.getRHSReal(i) * cone_grid_res;
|
||||
}
|
||||
// copy constraint RHS
|
||||
printf("\n SET %d CONSTRAINTS \n",num_constraint_equations);
|
||||
for(unsigned int i = 0; i < num_constraint_equations; ++i)
|
||||
{
|
||||
C(i, SizeMatrix) = -S.getRHSReal(ScalarSize + i) * cone_grid_res;
|
||||
}
|
||||
///copy values back into S
|
||||
COMISO::ConstrainedSolver solver;
|
||||
|
||||
solver.misolver().set_local_iters(localIter);
|
||||
solver.misolver().set_direct_rounding(direct_round);
|
||||
|
||||
std::sort(ids_to_round.begin(),ids_to_round.end());
|
||||
std::vector<int>::iterator new_end=std::unique(ids_to_round.begin(),ids_to_round.end());
|
||||
int dist=distance(ids_to_round.begin(),new_end);
|
||||
ids_to_round.resize(dist);
|
||||
solver.solve( C, A, X, rhs, ids_to_round, 0.0, true, true);
|
||||
}
|
||||
|
||||
void GetAttributes()
|
||||
{
|
||||
// you can query if an attribute is present or not
|
||||
bool hasSystIndex = vcg::tri::HasPerFaceAttribute(mesh,std::string("SystemIndex"));
|
||||
bool hasSingular = vcg::tri::HasPerVertexAttribute(mesh,std::string("Singular"));
|
||||
bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
|
||||
bool HasSystemInfo=vcg::tri::HasPerMeshAttribute(mesh,std::string("SystemInfo"));
|
||||
|
||||
assert(hasSystIndex);
|
||||
assert(hasSingular);
|
||||
assert(hasStiffness);
|
||||
assert(HasSystemInfo);
|
||||
|
||||
HandleS_Index = vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3i>(mesh,"SystemIndex");
|
||||
Handle_Singular=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<bool>(mesh,std::string("Singular"));
|
||||
Handle_Stiffness=vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<float>(mesh,std::string("Stiffness"));
|
||||
Handle_SystemInfo=vcg::tri::Allocator<MeshType>::template GetPerMeshAttribute<MeshSystemInfo>(mesh,"SystemInfo");
|
||||
}
|
||||
|
||||
public:
|
||||
|
||||
void SolvePoisson(ScalarType vector_field_scale=0.1f,
|
||||
ScalarType grid_res=1.f,
|
||||
bool direct_round=true,
|
||||
int localIter=0,
|
||||
bool _integer_rounding=true)
|
||||
{
|
||||
GetAttributes();
|
||||
//initialization of flags and data structures
|
||||
//ClearFlags();
|
||||
integer_rounding=_integer_rounding;
|
||||
|
||||
ids_to_round.clear();
|
||||
|
||||
///Initializing Matrix
|
||||
|
||||
int t0=clock();
|
||||
|
||||
///initialize the matrix ALLOCATING SPACE
|
||||
InitMatrix();
|
||||
printf("\n ALLOCATED THE MATRIX \n");
|
||||
|
||||
///build the laplacian system
|
||||
BuildLaplacianMatrix(vector_field_scale);
|
||||
|
||||
BuildSeamConstraintsExplicitTranslation();
|
||||
|
||||
////add the lagrange multiplier
|
||||
FixBlockedVertex();
|
||||
|
||||
printf("\n BUILT THE MATRIX \n");
|
||||
|
||||
if (integer_rounding)
|
||||
AddSingularityRound();
|
||||
|
||||
int t1=clock();
|
||||
printf("\n time:%d \n",t1-t0);
|
||||
printf("\n SOLVING \n");
|
||||
|
||||
MixedIntegerSolve(grid_res,direct_round,localIter);
|
||||
|
||||
int t2=clock();
|
||||
printf("\n time:%d \n",t2-t1);
|
||||
printf("\n ASSIGNING COORDS \n");
|
||||
MapCoords();
|
||||
int t3=clock();
|
||||
printf("\n time:%d \n",t3-t2);
|
||||
printf("\n FINISHED \n");
|
||||
//TagFoldedFaces();
|
||||
}
|
||||
|
||||
|
||||
PoissonSolver(MeshType &_mesh):mesh(_mesh)
|
||||
{}
|
||||
|
||||
};
|
||||
#endif
|
||||
|
|
@ -0,0 +1,517 @@
|
|||
#ifndef MIQ_QUADRANGULATOR_H
|
||||
#define MIQ_QUADRANGULATOR_H
|
||||
|
||||
|
||||
#include <vcg/simplex/face/pos.h>
|
||||
#include <vcg/simplex/face/jumping_pos.h>
|
||||
#include <vcg/complex/algorithms/attribute_seam.h>
|
||||
#include <vcg/complex/algorithms/refine.h>
|
||||
#include <vcg/complex/algorithms/smooth.h>
|
||||
#include <vcg/complex/algorithms/clean.h>
|
||||
|
||||
///return the list of faces that share a specified vertex
|
||||
///together with indexes those faces are sorted in ccw
|
||||
template <class FaceType>
|
||||
void SortedStar(const vcg::face::Pos<FaceType> &ep,
|
||||
std::vector< typename FaceType::VertexType*> &star)
|
||||
{
|
||||
typedef typename FaceType::VertexType VertexType;
|
||||
FaceType *f_init=ep.f;
|
||||
int edge_init=ep.z;
|
||||
vcg::face::JumpingPos<FaceType> VFI(f_init,edge_init);
|
||||
bool complete_turn=false;
|
||||
do
|
||||
{
|
||||
///take the current face
|
||||
FaceType *curr_f=VFI.F();
|
||||
int curr_edge=VFI.E();
|
||||
assert((curr_edge>=0)&&(curr_edge<=4));
|
||||
star.push_back(VFI.F()->V1(curr_edge));
|
||||
//go to next one
|
||||
VFI.NextFE();
|
||||
FaceType *next_f=VFI.F();
|
||||
///test the complete turn
|
||||
complete_turn=(next_f==f_init);
|
||||
}while (!complete_turn);
|
||||
}
|
||||
|
||||
template <class MeshType>
|
||||
inline void ExtractVertex(const MeshType & srcMesh,
|
||||
const typename MeshType::FaceType & f,
|
||||
int whichWedge,
|
||||
const MeshType &dstMesh,
|
||||
typename MeshType::VertexType & v)
|
||||
{
|
||||
(void)srcMesh;
|
||||
(void)dstMesh;
|
||||
|
||||
//v.P() = f.cP(whichWedge);
|
||||
v.ImportData(*f.cV(whichWedge));
|
||||
v.T() = f.cWT(whichWedge);
|
||||
}
|
||||
|
||||
template <class MeshType>
|
||||
inline bool CompareVertex(const MeshType & m,
|
||||
const typename MeshType::VertexType & vA,
|
||||
const typename MeshType::VertexType & vB)
|
||||
{
|
||||
(void)m;
|
||||
return ((vA.cT() == vB.cT())&&(vA.cP()==vB.cP()));
|
||||
}
|
||||
|
||||
|
||||
template <class TriMesh,class PolyMesh>
|
||||
class Quadrangulator
|
||||
{
|
||||
|
||||
public:
|
||||
typedef typename TriMesh::FaceType TriFaceType;
|
||||
typedef typename TriMesh::VertexType TriVertexType;
|
||||
typedef typename TriMesh::CoordType CoordType;
|
||||
typedef typename TriMesh::ScalarType ScalarType;
|
||||
|
||||
typedef typename PolyMesh::FaceType PolyFaceType;
|
||||
typedef typename PolyMesh::VertexType PolyVertexType;
|
||||
typedef typename PolyMesh::CoordType PolyCoordType;
|
||||
typedef typename PolyMesh::ScalarType PolyScalarType;
|
||||
|
||||
|
||||
///the set of all edges that belongs to integer lines
|
||||
std::set<std::pair<TriFaceType*,int> > IntegerEdges;
|
||||
|
||||
///the set of all integer vertices and the other vertices on integer lines which is connectes to
|
||||
std::map<TriVertexType*,std::vector<TriVertexType*> > IntegerLineAdj;
|
||||
///the set of integer vertices
|
||||
std::set<TriVertexType*> IntegerVertices;
|
||||
///temporary polygons
|
||||
std::vector<std::vector<TriVertexType *> > polygons;
|
||||
|
||||
///drawing debug structures
|
||||
std::vector<std::pair<TriFaceType*,int> > IntegerLines;
|
||||
std::vector<TriVertexType*> IntegerVertex;
|
||||
|
||||
static bool ToSplit(const vcg::Point2<ScalarType> &uv0,
|
||||
const vcg::Point2<ScalarType> &uv1,
|
||||
int Dir,
|
||||
int IntegerLine,
|
||||
ScalarType &alpha,
|
||||
ScalarType tolerance=0.0001)
|
||||
{
|
||||
ScalarType lineF=(ScalarType)IntegerLine;
|
||||
ScalarType val0=std::min(uv0.V(Dir),uv1.V(Dir));
|
||||
ScalarType val1=std::max(uv0.V(Dir),uv1.V(Dir));
|
||||
if (lineF<(val0+tolerance))return false;
|
||||
if (lineF>(val1-tolerance))return false;
|
||||
ScalarType dist=fabs(uv0.V(Dir)-uv1.V(Dir));
|
||||
if (dist<tolerance) return false;
|
||||
alpha=(1.0-fabs(uv0.V(Dir)-lineF)/dist);
|
||||
return true;
|
||||
}
|
||||
|
||||
///return true if the edge has to be splitted,
|
||||
///by considering the tolerance to the closest integer
|
||||
static bool ToSplit(const vcg::face::Pos<TriFaceType> &ep,
|
||||
ScalarType &alpha,
|
||||
ScalarType factor=1.0,
|
||||
ScalarType tolerance=0.0001)
|
||||
{
|
||||
//TriFaceType *f=ep.f;
|
||||
//int z=ep.z;
|
||||
TriVertexType* v0=ep.f->V(ep.z);
|
||||
TriVertexType* v1=ep.f->V1(ep.z);
|
||||
vcg::Point2<ScalarType> uv0=v0->T().P()*factor;
|
||||
vcg::Point2<ScalarType> uv1=v1->T().P()*factor;
|
||||
|
||||
///then test integer for each direction
|
||||
for (int dir=0;dir<2;dir++)
|
||||
{
|
||||
int Int0=std::min((int)uv0.V(dir),(int)uv1.V(dir));
|
||||
int Int1=std::max((int)uv0.V(dir),(int)uv1.V(dir));
|
||||
|
||||
for (int i=Int0;i<=Int1;i++)
|
||||
{
|
||||
bool to_split=ToSplit(uv0,uv1,dir,i,alpha,tolerance);
|
||||
if (to_split)return true;
|
||||
}
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
// Basic subdivision class
|
||||
// This class must provide methods for finding the position of the newly created vertices
|
||||
// In this implemenation we simply put the new vertex in the MidPoint position.
|
||||
// Color and TexCoords are interpolated accordingly.
|
||||
template<class MESH_TYPE>
|
||||
struct SplitMidPoint : public std::unary_function<vcg::face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType >
|
||||
{
|
||||
typedef typename MESH_TYPE::VertexType VertexType;
|
||||
typedef typename MESH_TYPE::FaceType FaceType;
|
||||
typedef typename MESH_TYPE::CoordType CoordType;
|
||||
|
||||
ScalarType factor;
|
||||
ScalarType tolerance;
|
||||
ScalarType alpha;
|
||||
|
||||
void operator()(typename MESH_TYPE::VertexType &nv,
|
||||
vcg::face::Pos<typename MESH_TYPE::FaceType> ep)
|
||||
{
|
||||
bool to_split=ToSplit(ep,alpha,factor,tolerance);
|
||||
assert(to_split);
|
||||
|
||||
///get the value on which the edge must be splitted
|
||||
VertexType* v0=ep.f->V(ep.z);
|
||||
VertexType* v1=ep.f->V1(ep.z);
|
||||
|
||||
nv.P()= v0->P()*alpha+v1->P()*(1.0-alpha);
|
||||
//nv.N()= v0->N()*alpha+v1->N()*(1.0-alpha);
|
||||
nv.T().P()=v0->T().P()*alpha+v1->T().P()*(1.0-alpha);
|
||||
}
|
||||
|
||||
vcg::TexCoord2<ScalarType> WedgeInterp(vcg::TexCoord2<ScalarType> &t0,
|
||||
vcg::TexCoord2<ScalarType> &t1)
|
||||
{
|
||||
vcg::TexCoord2<ScalarType> tmp;
|
||||
if (t0.n() != t1.n())
|
||||
cerr << "Failed assertion: Quadrangulator::WedgeInterp1" << endl;
|
||||
// assert(t0.n()== t1.n()); TODO put back
|
||||
tmp.n()=t0.n();
|
||||
// assert(alpha>=0); TODO put back
|
||||
if (alpha<0)
|
||||
cerr << "Failed assertion: Quadrangulator::WedgeInterp2" << endl;
|
||||
tmp.t()=(alpha*t0.t()+(1.0-alpha)*t1.t());
|
||||
return tmp;
|
||||
}
|
||||
|
||||
SplitMidPoint(){alpha=-1;}
|
||||
};
|
||||
|
||||
template <class MESH_TYPE>
|
||||
class EdgePredicate
|
||||
{
|
||||
typedef typename MESH_TYPE::VertexType VertexType;
|
||||
typedef typename MESH_TYPE::FaceType FaceType;
|
||||
typedef typename MESH_TYPE::ScalarType ScalarType;
|
||||
|
||||
public:
|
||||
ScalarType factor;
|
||||
ScalarType tolerance;
|
||||
|
||||
bool operator()(vcg::face::Pos<typename MESH_TYPE::FaceType> ep) const
|
||||
{
|
||||
ScalarType alpha;
|
||||
return(ToSplit(ep,alpha,factor,tolerance));
|
||||
}
|
||||
};
|
||||
|
||||
void SplitTris(TriMesh &to_split,
|
||||
ScalarType factor=1.0,
|
||||
ScalarType tolerance=0.0001)
|
||||
{
|
||||
bool done=true;
|
||||
SplitMidPoint<TriMesh> splMd;
|
||||
EdgePredicate<TriMesh> eP;
|
||||
|
||||
splMd.tolerance=tolerance;
|
||||
splMd.factor=factor;
|
||||
eP.tolerance=tolerance;
|
||||
eP.factor=factor;
|
||||
|
||||
while (done)
|
||||
done=vcg::tri::RefineE<TriMesh,SplitMidPoint<TriMesh>,EdgePredicate<TriMesh> >(to_split,splMd,eP);
|
||||
|
||||
for (unsigned int i=0;i<to_split.face.size();i++)
|
||||
for (int j=0;j<3;j++) to_split.face[i].WT(j).P()=to_split.face[i].V(j)->T().P();
|
||||
}
|
||||
|
||||
|
||||
bool IsOnIntegerLine(vcg::Point2<ScalarType> uv0,
|
||||
vcg::Point2<ScalarType> uv1,
|
||||
ScalarType tolerance=0.0001)
|
||||
{
|
||||
for (int dir=0;dir<2;dir++)
|
||||
{
|
||||
ScalarType val0=uv0.V(dir);
|
||||
ScalarType val1=uv1.V(dir);
|
||||
int integer0=floor(uv0.V(dir)+0.5);
|
||||
int integer1=floor(uv1.V(dir)+0.5);
|
||||
if (integer0!=integer1)continue;
|
||||
if ((fabs(val0-(ScalarType)integer0))>tolerance)continue;
|
||||
if ((fabs(val1-(ScalarType)integer1))>tolerance)continue;
|
||||
return true;
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
bool IsOnIntegerVertex(vcg::Point2<ScalarType> uv,
|
||||
ScalarType tolerance=0.0001)
|
||||
{
|
||||
for (int dir=0;dir<2;dir++)
|
||||
{
|
||||
ScalarType val0=uv.V(dir);
|
||||
int integer0=floor(val0+0.5);
|
||||
if ((fabs(val0-(ScalarType)integer0))>tolerance)return false;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
void InitIntegerVectors()
|
||||
{
|
||||
IntegerLines=std::vector<std::pair<TriFaceType*,int> >(IntegerEdges.begin(),IntegerEdges.end());
|
||||
IntegerVertex=std::vector<TriVertexType* > (IntegerVertices.begin(),IntegerVertices.end());
|
||||
}
|
||||
|
||||
void EraseIntegerEdge(const vcg::face::Pos<TriFaceType> &ep)
|
||||
{
|
||||
std::pair<TriFaceType*,int> edge(ep.F(),ep.E());
|
||||
assert(IntegerEdges.count(edge)!=0);
|
||||
IntegerEdges.erase(edge);
|
||||
}
|
||||
|
||||
void EraseIntegerEdge(const std::vector<std::pair<TriFaceType*,int> > &to_erase)
|
||||
{
|
||||
for (unsigned int i=0;i<to_erase.size();i++)
|
||||
IntegerEdges.erase(to_erase[i]);
|
||||
}
|
||||
|
||||
void TestIntegerEdges()
|
||||
{
|
||||
typedef typename std::pair<TriFaceType*,int> pair_type;
|
||||
typedef typename std::vector< pair_type > vect_type;
|
||||
typename vect_type::iterator IteIntl;
|
||||
for (IteIntl=IntegerLines.begin();
|
||||
IteIntl!=IntegerLines.end();
|
||||
IteIntl++)
|
||||
{
|
||||
int E=(*IteIntl).second;
|
||||
TriFaceType *F=(*IteIntl).first;
|
||||
TriFaceType *F1=F->FFp(E);
|
||||
if (F==F1) continue;
|
||||
int E1=F->FFi(E);
|
||||
std::pair<TriFaceType*,int> curr_edge(F1,E1);
|
||||
assert(IntegerEdges.count(curr_edge)!=0);
|
||||
}
|
||||
}
|
||||
|
||||
void InitIntegerEdgesVert(TriMesh &Tmesh,
|
||||
ScalarType factor=1.0,
|
||||
ScalarType tolerance=0.0001)
|
||||
{
|
||||
IntegerEdges.clear();
|
||||
for (unsigned int i=0;i<Tmesh.face.size();i++)
|
||||
{
|
||||
TriFaceType *f=&Tmesh.face[i];
|
||||
if (f->IsD())continue;
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
TriFaceType *f1=f->FFp(j);
|
||||
int e1=f->FFi(j);
|
||||
bool IsBorder=f->IsB(j);
|
||||
TriVertexType *v0=f->V0(j);
|
||||
TriVertexType *v1=f->V1(j);
|
||||
vcg::Point2<ScalarType> uv0=f->WT(j).P()*factor;
|
||||
vcg::Point2<ScalarType> uv1=f->WT((j+1)%3).P()*factor;
|
||||
if (IsOnIntegerLine(uv0,uv1,tolerance)||IsBorder)
|
||||
{
|
||||
//IntegerEdges.insert(std::pair<TriVertexType*,TriVertexType*>(v0,v1));
|
||||
IntegerEdges.insert(std::pair<TriFaceType*,int>(f,j));
|
||||
if (!IsBorder)
|
||||
IntegerEdges.insert(std::pair<TriFaceType*,int>(f1,e1));
|
||||
else
|
||||
{
|
||||
IntegerVertices.insert(v0);
|
||||
IntegerVertices.insert(v1);
|
||||
}
|
||||
}
|
||||
if (IsOnIntegerVertex(uv0))
|
||||
IntegerVertices.insert(v0);
|
||||
|
||||
if (IsOnIntegerVertex(uv1))
|
||||
IntegerVertices.insert(v1);
|
||||
|
||||
}
|
||||
}
|
||||
//InitIntegerNeigh(Tmesh);
|
||||
InitIntegerVectors();
|
||||
TestIntegerEdges();
|
||||
}
|
||||
|
||||
|
||||
///return the first and the last edge
|
||||
///following an integer line
|
||||
///until if reach anothe integer edge
|
||||
bool OneIntegerStep(vcg::face::Pos<TriFaceType> &ep)
|
||||
{
|
||||
TriFaceType *f_init=ep.f;
|
||||
TriFaceType *currF=f_init;
|
||||
//int edge_init=ep.z;
|
||||
//ep.V()=f_init->V(edge_init);
|
||||
TriVertexType* v_init=ep.V();
|
||||
bool complete_turn=false;
|
||||
do
|
||||
{
|
||||
ep.FlipE();
|
||||
///see if found an integer vert
|
||||
currF=ep.F();
|
||||
int currE=ep.E();
|
||||
assert((currE>=0)&&(currE<=4));
|
||||
|
||||
std::pair<TriFaceType*,int> curr_edge(currF,currE);
|
||||
|
||||
if (IntegerEdges.count(curr_edge)!=0)
|
||||
{
|
||||
///go to the other side
|
||||
ep.FlipV();
|
||||
assert(ep.V()!=v_init);
|
||||
return true;
|
||||
}
|
||||
ep.FlipF();
|
||||
///see if there's a border
|
||||
bool jumped=(currF==ep.F());
|
||||
if (jumped)
|
||||
return false;
|
||||
///test the complete turn
|
||||
complete_turn=(ep.F()==f_init);
|
||||
}while (!complete_turn);
|
||||
return false;
|
||||
}
|
||||
|
||||
///find a polygon starting from half edge ep, return true if found
|
||||
bool FindPolygon(vcg::face::Pos<TriFaceType> &ep,
|
||||
std::vector<TriVertexType*> &poly)
|
||||
{
|
||||
|
||||
poly.clear();
|
||||
TriVertexType* v_init=ep.V();
|
||||
///it must start from an integer vert
|
||||
assert(IntegerVertices.count(v_init)!=0);
|
||||
poly.push_back(v_init);
|
||||
std::vector<std::pair<TriFaceType*,int> > to_erase;
|
||||
to_erase.push_back(std::pair<TriFaceType*,int>(ep.F(),ep.E()));
|
||||
do
|
||||
{
|
||||
bool done=OneIntegerStep(ep);
|
||||
if (!done)
|
||||
{
|
||||
EraseIntegerEdge(to_erase);
|
||||
return false;
|
||||
}
|
||||
to_erase.push_back(std::pair<TriFaceType*,int>(ep.F(),ep.E()));
|
||||
TriVertexType* v_curr=ep.V();
|
||||
if ((IntegerVertices.count(v_curr)!=0)&&
|
||||
(v_curr!=v_init))
|
||||
poly.push_back(v_curr);
|
||||
}while(ep.V()!=v_init);
|
||||
EraseIntegerEdge(to_erase);
|
||||
return true;
|
||||
}
|
||||
|
||||
void FindPolygons(TriMesh &Tmesh,
|
||||
std::vector<std::vector<TriVertexType *> > &polygons)
|
||||
{
|
||||
//int limit=2;
|
||||
for (unsigned int i=0;i<Tmesh.face.size();i++)
|
||||
{
|
||||
TriFaceType * f=&Tmesh.face[i];
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
TriVertexType* v0=f->V0(j);
|
||||
//TriVertexType* v1=f->V1(j);
|
||||
|
||||
//std::pair<TriVertexType*,TriVertexType*> edge(v0,v1);*/
|
||||
std::pair<TriFaceType*,int> edge(f,j);
|
||||
if (IntegerEdges.count(edge)==0)continue;///edge already used or not integer
|
||||
if (IntegerVertices.count(v0)==0)continue; ///must start from integer vert
|
||||
///create the pos
|
||||
vcg::face::Pos<TriFaceType> ep(f,j);
|
||||
std::vector<TriVertexType *> poly;
|
||||
|
||||
bool found=FindPolygon(ep,poly);
|
||||
if (found)
|
||||
{
|
||||
std::reverse(poly.begin(),poly.end());///REVERSE ORDER
|
||||
polygons.push_back(poly);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void InitVertexQuadMesh(TriMesh &Tmesh)
|
||||
{
|
||||
FindPolygons(Tmesh,polygons);
|
||||
}
|
||||
|
||||
public:
|
||||
|
||||
|
||||
void TestIsProper(TriMesh &Tmesh)
|
||||
{
|
||||
///test manifoldness
|
||||
int test=vcg::tri::Clean<TriMesh>::CountNonManifoldVertexFF(Tmesh);
|
||||
//assert(test==0);
|
||||
if (test != 0)
|
||||
cerr << "Assertion failed: TestIsProper NonManifoldVertices!" << endl;
|
||||
|
||||
test=vcg::tri::Clean<TriMesh>::CountNonManifoldEdgeFF(Tmesh);
|
||||
//assert(test==0);
|
||||
if (test != 0)
|
||||
cerr << "Assertion failed: TestIsProper NonManifoldEdges" << endl;
|
||||
|
||||
for (unsigned int i=0;i<Tmesh.face.size();i++)
|
||||
{
|
||||
TriFaceType *f=&Tmesh.face[i];
|
||||
assert (!f->IsD());
|
||||
for (int z=0;z<3;z++)
|
||||
{
|
||||
//int indexOpp=f->FFi(z);
|
||||
TriFaceType *Fopp=f->FFp(z);
|
||||
if (Fopp==f) continue;
|
||||
//assert( f->FFp(z)->FFp(f->FFi(z))==f );
|
||||
if (f->FFp(z)->FFp(f->FFi(z))!=f)
|
||||
cerr << "Assertion failed: TestIsProper f->FFp(z)->FFp(f->FFi(z))!=f " << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void Quadrangulate(TriMesh &Tmesh,
|
||||
PolyMesh &Pmesh,
|
||||
ScalarType factor=1.0,
|
||||
ScalarType tolerance=0.000001)
|
||||
{
|
||||
TestIsProper(Tmesh);
|
||||
vcg::tri::AttributeSeam::SplitVertex(Tmesh, ExtractVertex<TriMesh>, CompareVertex<TriMesh>);
|
||||
vcg::tri::Allocator<TriMesh>::CompactVertexVector(Tmesh);
|
||||
vcg::tri::Allocator<TriMesh>::CompactFaceVector(Tmesh);
|
||||
vcg::tri::UpdateTopology<TriMesh>::FaceFace(Tmesh);
|
||||
(void)Pmesh;
|
||||
TestIsProper(Tmesh);
|
||||
|
||||
///then split the tris
|
||||
SplitTris(Tmesh,factor,tolerance);
|
||||
///join the vertices back!
|
||||
ScalarType EPS=(ScalarType)0.00000001;
|
||||
vcg::tri::Clean<TriMesh>::MergeCloseVertex(Tmesh,EPS);
|
||||
|
||||
vcg::tri::UpdateNormal<TriMesh>::PerFaceNormalized(Tmesh); // update Normals
|
||||
vcg::tri::UpdateNormal<TriMesh>::PerVertexNormalized(Tmesh);// update Normals
|
||||
///compact the mesh
|
||||
vcg::tri::Allocator<TriMesh>::CompactVertexVector(Tmesh);
|
||||
vcg::tri::Allocator<TriMesh>::CompactFaceVector(Tmesh);
|
||||
vcg::tri::UpdateTopology<TriMesh>::FaceFace(Tmesh); // update Topology
|
||||
vcg::tri::UpdateTopology<TriMesh>::TestFaceFace(Tmesh); //and test it
|
||||
///set flags
|
||||
vcg::tri::UpdateFlags<TriMesh>::VertexClearV(Tmesh);
|
||||
vcg::tri::UpdateFlags<TriMesh>::FaceBorderFromFF(Tmesh);
|
||||
vcg::tri::UpdateFlags<TriMesh>::VertexBorderFromFace(Tmesh);
|
||||
///test manifoldness
|
||||
TestIsProper(Tmesh);
|
||||
|
||||
vcg::tri::UpdateFlags<TriMesh>::VertexClearV(Tmesh);
|
||||
|
||||
InitIntegerEdgesVert(Tmesh,factor,tolerance);
|
||||
InitVertexQuadMesh(Tmesh);
|
||||
}
|
||||
};
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,325 @@
|
|||
#ifndef MIQ_SEAMS_INTIALIZER
|
||||
#define MIQ_SEAMS_INTIALIZER
|
||||
#include <vcg/complex/complex.h>
|
||||
#include <vcg/simplex/face/pos.h>
|
||||
#include <vcg/simplex/face/jumping_pos.h>
|
||||
#include <vcg/simplex/face/topology.h>
|
||||
|
||||
|
||||
template <class MeshType>
|
||||
class SeamsInitializer
|
||||
{
|
||||
|
||||
private:
|
||||
typedef typename MeshType::ScalarType ScalarType;
|
||||
typedef typename MeshType::FaceType FaceType;
|
||||
typedef typename MeshType::VertexType VertexType;
|
||||
typedef typename MeshType::CoordType CoordType;
|
||||
typedef typename MeshType::FaceIterator FaceIterator;
|
||||
|
||||
MeshType *mesh;
|
||||
|
||||
///per face per edge of mmatch in the solver
|
||||
typename MeshType::template PerFaceAttributeHandle<vcg::Point3i> Handle_MMatch;
|
||||
///per vertex singular or not
|
||||
typename MeshType::template PerVertexAttributeHandle<bool> Handle_Singular;
|
||||
///per vertex degree of a singularity
|
||||
typename MeshType::template PerVertexAttributeHandle<int> Handle_SingularDegree;
|
||||
///seam per face
|
||||
typename MeshType::template PerFaceAttributeHandle<vcg::Point3<bool> > Handle_Seams;
|
||||
|
||||
bool IsRotSeam(const FaceType *f0,const int edge)
|
||||
{
|
||||
unsigned char MM=Handle_MMatch[f0][edge];//MissMatch(f0,edge);
|
||||
return (MM!=0);
|
||||
}
|
||||
|
||||
///return true if a vertex is singluar by looking at initialized missmatches
|
||||
bool IsSingularByMMatch(const CVertex &v,int &missmatch)
|
||||
{
|
||||
///check that is on border..
|
||||
if (v.IsB())return false;
|
||||
|
||||
std::vector<CFace*> faces;
|
||||
std::vector<int> edges;
|
||||
//SortedFaces(v,faces);
|
||||
vcg::face::VFOrderedStarVF_FF<CFace>(v,faces,edges);
|
||||
|
||||
missmatch=0;
|
||||
for (unsigned int i=0;i<faces.size();i++)
|
||||
{
|
||||
CFace *curr_f=faces[i];
|
||||
int currMM=Handle_MMatch[curr_f][edges[i]];
|
||||
missmatch+=currMM;
|
||||
}
|
||||
missmatch=missmatch%4;
|
||||
return(missmatch!=0);
|
||||
}
|
||||
|
||||
///initialized mapping structures if are not already initialized
|
||||
void AddAttributesIfNeeded()
|
||||
{
|
||||
|
||||
bool HasHandleMMatch=vcg::tri::HasPerFaceAttribute(*mesh,std::string("MissMatch"));
|
||||
if (!HasHandleMMatch)
|
||||
Handle_MMatch = vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<vcg::Point3i>(*mesh,std::string("MissMatch"));
|
||||
else
|
||||
Handle_MMatch = vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3i>(*mesh,std::string("MissMatch"));
|
||||
|
||||
bool HasHandleSingularity=vcg::tri::HasPerVertexAttribute(*mesh,std::string("Singular"));
|
||||
if (!HasHandleSingularity)
|
||||
Handle_Singular=vcg::tri::Allocator<MeshType>::template AddPerVertexAttribute<bool>(*mesh,std::string("Singular"));
|
||||
else
|
||||
Handle_Singular=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<bool>(*mesh,std::string("Singular"));
|
||||
|
||||
bool HasHandleSingularityDegree=vcg::tri::HasPerVertexAttribute(*mesh,std::string("SingularityDegree"));
|
||||
if (!HasHandleSingularityDegree)
|
||||
Handle_SingularDegree=vcg::tri::Allocator<MeshType>::template AddPerVertexAttribute<int>(*mesh,std::string("SingularityDegree"));
|
||||
else
|
||||
Handle_SingularDegree=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<int>(*mesh,std::string("SingularityDegree"));
|
||||
|
||||
bool HasHandleSeams=vcg::tri::HasPerFaceAttribute(*mesh,std::string("Seams"));
|
||||
if (!HasHandleSeams)
|
||||
Handle_Seams=vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<vcg::Point3<bool> >(*mesh,std::string("Seams"));
|
||||
else
|
||||
Handle_Seams=vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3<bool> >(*mesh,std::string("Seams"));
|
||||
}
|
||||
|
||||
|
||||
void FloodFill(FaceType* start)
|
||||
{
|
||||
std::deque<FaceType*> d;
|
||||
///clean the visited flag
|
||||
start->SetV();
|
||||
d.push_back(start);
|
||||
|
||||
while (!d.empty()){
|
||||
FaceType *f = d.at(0); d.pop_front();
|
||||
for (int s = 0; s<3; s++)
|
||||
{
|
||||
FaceType *g = f->FFp(s);
|
||||
int j = f->FFi(s);
|
||||
if ((!(IsRotSeam(f,s))) && (!(IsRotSeam(g,j))) && (!g->IsV()) )
|
||||
{
|
||||
//f->seam[s] = false;
|
||||
Handle_Seams[f][s]=false;
|
||||
//g->seam[j] = false; // dissolve seam
|
||||
Handle_Seams[g][j]=false;
|
||||
g->SetV();
|
||||
d.push_back(g);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void Retract(){
|
||||
std::vector<int> e(mesh->vert.size(),0); // number of edges per vert
|
||||
VertexType *vb = &(mesh->vert[0]);
|
||||
for (FaceIterator f = mesh->face.begin(); f!=mesh->face.end(); f++) if (!f->IsD()){
|
||||
for (int s = 0; s<3; s++){
|
||||
//if (f->seam[s])
|
||||
if (Handle_Seams[f][s])
|
||||
if (f->FFp(s)<=&*f) {
|
||||
e[ f->V(s) - vb ] ++;
|
||||
e[ f->V1(s) - vb ] ++;
|
||||
}
|
||||
}
|
||||
}
|
||||
bool over=true;
|
||||
int guard = 0;
|
||||
do {
|
||||
over = true;
|
||||
for (FaceIterator f = mesh->face.begin(); f!=mesh->face.end(); f++) if (!f->IsD()){
|
||||
for (int s = 0; s<3; s++){
|
||||
//if (f->seam[s])
|
||||
if (Handle_Seams[f][s])
|
||||
if (!(IsRotSeam(&(*f),s))) // never retract rot seams
|
||||
//if (f->FFp(s)<=&*f)
|
||||
{
|
||||
if (e[ f->V(s) - vb ] == 1) {
|
||||
// dissolve seam
|
||||
//f->seam[s] = false;
|
||||
Handle_Seams[f][s]=false;
|
||||
//f->FFp(s)->seam[(int)f->FFi(s)] = false;
|
||||
Handle_Seams[f->FFp(s)][(int)f->FFi(s)]=false;
|
||||
e[ f->V(s) - vb ] --;
|
||||
e[ f->V1(s) - vb ] --;
|
||||
over = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (guard++>10000) over = true;
|
||||
|
||||
} while (!over);
|
||||
}
|
||||
|
||||
bool LoadSeamsMMFromOBJ(std::string PathOBJ)
|
||||
{
|
||||
FILE *f = fopen(PathOBJ.c_str(),"rt");
|
||||
if (!f)
|
||||
return false;
|
||||
|
||||
for (unsigned int i=0;i<mesh->face.size();i++)
|
||||
{
|
||||
for (int j=0;j<3;j++)
|
||||
Handle_Seams[i][j]=false;
|
||||
}
|
||||
|
||||
while (!feof(f))
|
||||
{
|
||||
|
||||
int f_int,v_int,rot;
|
||||
int readed=fscanf(f,"sm %d %d %d\n",&f_int,&v_int,&rot);
|
||||
///skip lines
|
||||
if (readed==0)
|
||||
{
|
||||
char buff[200];
|
||||
fscanf(f,"%s\n",&buff[0]);
|
||||
}
|
||||
else ///add the actual seams
|
||||
{
|
||||
VertexType *v=&mesh->vert[v_int-1];
|
||||
FaceType *f0=&mesh->face[f_int-1];
|
||||
int e0=-1;
|
||||
if (f0->V(0)==v)e0=0;
|
||||
if (f0->V(1)==v)e0=1;
|
||||
if (f0->V(2)==v)e0=2;
|
||||
e0=(e0+2)%3;
|
||||
assert(e0!=-1);
|
||||
FaceType *f1;
|
||||
int e1;
|
||||
f1=f0->FFp(e0);
|
||||
e1=f0->FFi(e0);
|
||||
Handle_Seams[f0][e0]=true;
|
||||
Handle_Seams[f1][e1]=true;
|
||||
|
||||
Handle_MMatch[f0][e0]=rot;
|
||||
int rot1;
|
||||
if (rot==0)rot1=0;
|
||||
if (rot==1)rot1=3;
|
||||
if (rot==2)rot1=2;
|
||||
if (rot==3)rot1=1;
|
||||
Handle_MMatch[f1][e1]=rot1;
|
||||
}
|
||||
}
|
||||
//printf("NEED %d LINES\n",i);
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
|
||||
void AddSeamsByMM()
|
||||
{
|
||||
for (unsigned int i=0;i<mesh->face.size();i++)
|
||||
{
|
||||
FaceType *f=&mesh->face[i];
|
||||
if (f->IsD())continue;
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
if (IsRotSeam(f,j))
|
||||
Handle_Seams[f][j]=true;
|
||||
//f->SetSeam(j);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void SelectSingularityByMM()
|
||||
{
|
||||
for (unsigned int i=0;i<mesh->vert.size();i++)
|
||||
{
|
||||
if (mesh->vert[i].IsD())continue;
|
||||
int missmatch;
|
||||
bool isSing=IsSingularByMMatch(mesh->vert[i],missmatch);
|
||||
if (isSing)
|
||||
{
|
||||
//mesh->vert[i].SetS();
|
||||
Handle_Singular[i]=true;
|
||||
Handle_SingularDegree[i]=missmatch;
|
||||
}
|
||||
else
|
||||
{
|
||||
//mesh->vert[i].ClearS();
|
||||
Handle_Singular[i]=false;
|
||||
Handle_SingularDegree[i]=0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
int InitTopologycalCuts(){
|
||||
vcg::tri::UpdateFlags<CMesh>::FaceClearV(*mesh);
|
||||
|
||||
for (FaceIterator f = mesh->face.begin(); f!=mesh->face.end(); f++)
|
||||
if (!f->IsD())
|
||||
{
|
||||
Handle_Seams[f][0]=true;
|
||||
Handle_Seams[f][1]=true;
|
||||
Handle_Seams[f][2]=true;
|
||||
}
|
||||
|
||||
int index=0;
|
||||
for (FaceIterator f = mesh->face.begin(); f!=mesh->face.end(); f++)
|
||||
if (!f->IsD())
|
||||
{
|
||||
if (!f->IsV())
|
||||
{
|
||||
index++;
|
||||
FloodFill(&*f);
|
||||
}
|
||||
}
|
||||
|
||||
Retract();
|
||||
return index;
|
||||
}
|
||||
|
||||
void InitMMatch()
|
||||
{
|
||||
for (unsigned int i=0;i<mesh->face.size();i++)
|
||||
{
|
||||
CFace *curr_f=&mesh->face[i];
|
||||
for (int j=0;j<3;j++)
|
||||
{
|
||||
CFace *opp_f=curr_f->FFp(j);
|
||||
if (curr_f==opp_f)
|
||||
Handle_MMatch[curr_f][j]=0;
|
||||
else
|
||||
Handle_MMatch[curr_f][j]=vcg::tri::CrossField<CMesh>::MissMatchByCross(*curr_f,*opp_f);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
public:
|
||||
|
||||
bool InitFromGradOBJ(const std::string &PathGrad,
|
||||
const std::string &PathObj)
|
||||
{
|
||||
AddAttributesIfNeeded();
|
||||
///OPEN THE GRAD FILE
|
||||
bool field_loaded=vcg::tri::CrossField<MeshType>::LoadGrad(mesh,PathGrad.c_str());
|
||||
if (!field_loaded)return false;
|
||||
LoadSeamsMMFromOBJ(PathObj);
|
||||
SelectSingularityByMM();
|
||||
return true;
|
||||
}
|
||||
|
||||
bool InitFromFField(const std::string &PathFField)
|
||||
{
|
||||
AddAttributesIfNeeded();
|
||||
bool field_loaded=vcg::tri::CrossField<MeshType>::LoadFFIELD(mesh,PathFField.c_str());
|
||||
if (!field_loaded)return false;
|
||||
vcg::tri::CrossField<MeshType>::MakeDirectionFaceCoherent(*mesh);
|
||||
InitMMatch();
|
||||
SelectSingularityByMM();
|
||||
InitTopologycalCuts();
|
||||
AddSeamsByMM();
|
||||
return true;
|
||||
}
|
||||
|
||||
void Init(MeshType *_mesh){mesh=_mesh;}//AllocateMappingStructures();}
|
||||
|
||||
SeamsInitializer(){mesh=NULL;}
|
||||
};
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,256 @@
|
|||
#ifndef __SPARSE_SYSTEM_DATA_H__
|
||||
#define __SPARSE_SYSTEM_DATA_H__
|
||||
|
||||
#include "auxmath.h"
|
||||
//#include "MatlabInterface.h"
|
||||
|
||||
class SparseMatrixData{
|
||||
protected:
|
||||
unsigned int m_nrows;
|
||||
unsigned int m_ncols;
|
||||
unsigned int m_nentries; // total number of (including repeated) (i,j) entries
|
||||
unsigned int m_reserved_entries;
|
||||
unsigned int *m_rowind;
|
||||
unsigned int *m_colind;
|
||||
double *m_vals;
|
||||
unsigned int *m_nz; // nonzeros per row
|
||||
|
||||
|
||||
public:
|
||||
unsigned int nrows() { return m_nrows ; }
|
||||
unsigned int ncols() { return m_ncols ; }
|
||||
unsigned int nentries() { return m_nentries; }
|
||||
unsigned int* nz() { return m_nz ; }
|
||||
unsigned int& nz(unsigned int i) { assert(i < m_nrows); return m_nz[i]; }
|
||||
unsigned int* rowind() { return m_rowind ; }
|
||||
unsigned int* colind() { return m_colind ; }
|
||||
double* vals() { return m_vals ; }
|
||||
|
||||
// create an empty matrix with a fixed number of rows
|
||||
SparseMatrixData(): m_nrows(0), m_ncols(0), m_nentries(0), m_reserved_entries(0), m_rowind(0), m_colind(0), m_vals(0),
|
||||
m_nz(NULL){ }
|
||||
|
||||
// create an empty matrix with a fixed number of rows
|
||||
void initialize(int nr, int nc) {
|
||||
assert(nr >= 0 && nc >=0);
|
||||
m_nrows = nr;
|
||||
m_ncols = nc;
|
||||
m_nentries = 0;
|
||||
m_reserved_entries = 0;
|
||||
m_rowind = NULL;
|
||||
m_colind = NULL;
|
||||
m_vals = NULL;
|
||||
if(nr > 0) {
|
||||
m_nz = new unsigned int[m_nrows];
|
||||
assert(m_nz);
|
||||
} else m_nz = 0;
|
||||
std::fill( m_nz, m_nz+m_nrows, 0 );
|
||||
}
|
||||
// allocate space for nnz nonzero entries
|
||||
void reserve(int nnz){
|
||||
assert(nnz >= 0);
|
||||
cleanup();
|
||||
if(nnz > 0) {
|
||||
m_rowind = new unsigned int[nnz];
|
||||
m_colind = new unsigned int[nnz];
|
||||
m_vals = new double[nnz];
|
||||
assert( m_rowind && m_colind && m_vals);
|
||||
} else { m_rowind = 0; m_colind = 0; m_vals = 0; }
|
||||
m_reserved_entries = nnz;
|
||||
}
|
||||
// extend space for entries
|
||||
void extend_entries(int extra_nnz){
|
||||
assert(m_nrows > 0);
|
||||
if( extra_nnz <=0) return;
|
||||
unsigned int* old_rowind = m_rowind;
|
||||
unsigned int* old_colind = m_colind;
|
||||
double* old_vals = m_vals;
|
||||
m_rowind = new unsigned int[m_reserved_entries+extra_nnz];
|
||||
m_colind = new unsigned int[m_reserved_entries+extra_nnz];
|
||||
m_vals = new double[m_reserved_entries+extra_nnz];
|
||||
assert( m_rowind && m_colind && m_vals);
|
||||
if(old_rowind) { std::copy(old_rowind, old_rowind+m_reserved_entries, m_rowind); delete[] old_rowind;}
|
||||
if(old_colind) { std::copy(old_colind, old_colind+m_reserved_entries, m_colind); delete[] old_colind;}
|
||||
if(old_vals) { std::copy(old_vals, old_vals+m_reserved_entries, m_vals); delete[] old_vals; }
|
||||
m_reserved_entries += extra_nnz;
|
||||
}
|
||||
|
||||
virtual void extend_rows(int extra_rows){
|
||||
if(extra_rows <= 0) return;
|
||||
unsigned int* old_nz = m_nz;
|
||||
m_nz = new unsigned int[m_nrows+extra_rows];
|
||||
if(old_nz) { std::copy(old_nz, old_nz+m_nrows, m_nz); delete[] old_nz;}
|
||||
std::fill( m_nz+m_nrows, m_nz+m_nrows+extra_rows, 0);
|
||||
m_nrows += extra_rows;
|
||||
}
|
||||
|
||||
// reset the matrix to empty, deallocate space
|
||||
void cleanup(){
|
||||
if(m_rowind) delete[] m_rowind; m_rowind = 0;
|
||||
if(m_colind) delete[] m_colind; m_colind = 0;
|
||||
if(m_vals) delete[] m_vals; m_vals = 0;
|
||||
m_nentries = 0;
|
||||
m_reserved_entries = 0;
|
||||
}
|
||||
|
||||
// add a nonzero entry to the matrix
|
||||
// no checks are done for coinciding entries
|
||||
// the interpretation of the repeated entries (replace or add)
|
||||
// depends on how the actual sparse matrix datastructure is constructed
|
||||
|
||||
void addEntryCmplx(unsigned int i, unsigned int j, Cmplx val) {
|
||||
assert(m_nentries < m_reserved_entries-3);
|
||||
m_rowind[m_nentries ] = 2*i; m_colind[m_nentries ] = 2*j; m_vals[m_nentries] = val.real();
|
||||
m_rowind[m_nentries+1] = 2*i; m_colind[m_nentries+1] = 2*j+1; m_vals[m_nentries+1] = -val.imag();
|
||||
m_rowind[m_nentries+2] = 2*i+1; m_colind[m_nentries+2] = 2*j; m_vals[m_nentries+2] = val.imag();
|
||||
m_rowind[m_nentries+3] = 2*i+1; m_colind[m_nentries+3] = 2*j+1; m_vals[m_nentries+3] = val.real();
|
||||
m_nentries += 4;
|
||||
}
|
||||
|
||||
void addEntryReal(unsigned int i, unsigned int j, double val) {
|
||||
assert(m_nentries < m_reserved_entries);
|
||||
m_rowind[m_nentries] = i; m_colind[m_nentries] = j; m_vals[m_nentries] = val;
|
||||
m_nentries++;
|
||||
}
|
||||
|
||||
|
||||
void symmetrize(unsigned int startind, unsigned int endind) {
|
||||
assert( startind <= m_nentries && endind <= m_nentries);
|
||||
for( unsigned int i = startind; i < endind; i++) {
|
||||
addEntryReal( m_colind[i], m_rowind[i], m_vals[i]);
|
||||
}
|
||||
}
|
||||
|
||||
/*void sendToMatlab(const char* name) {
|
||||
MatlabInterface matlab;
|
||||
if( m_rowind && m_colind && m_vals)
|
||||
matlab.SetEngineSparseRealMatrix(name, m_nentries, &m_rowind[0], &m_colind[0], &m_vals[0], m_nrows, m_ncols);
|
||||
else
|
||||
matlab.SetEngineSparseRealMatrix(name, 0, (const int*)0, (const int*)0, (const double*)0, m_nrows, m_ncols);
|
||||
}*/
|
||||
|
||||
|
||||
/* void getFromMatlab(const char* name) {
|
||||
MatlabInterface matlab;
|
||||
cleanup();
|
||||
if (m_nz) delete[] m_nz;
|
||||
matlab.GetEncodedSparseRealMatrix(name,m_rowind,m_colind,m_vals, m_nentries);
|
||||
m_reserved_entries = m_nentries;
|
||||
m_nrows = 0;
|
||||
m_ncols = 0;
|
||||
for(int i = 0; i < (int)m_nentries; i++) {
|
||||
m_nrows = std::max(m_rowind[i]+1,m_nrows);
|
||||
m_ncols = std::max(m_colind[i]+1,m_ncols);
|
||||
}
|
||||
m_nz = new unsigned int[m_nrows];
|
||||
std::fill(m_nz, m_nz+m_nrows,0);
|
||||
for(int i = 0; i < (int)m_nentries; i++) {
|
||||
m_nz[m_rowind[i]]++;
|
||||
}
|
||||
}*/
|
||||
|
||||
virtual ~SparseMatrixData() {
|
||||
cleanup();
|
||||
delete [] m_nz;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
// a small class to manage storage for matrix data
|
||||
// not using stl vectors: want to make all memory management
|
||||
// explicit to avoid hidden automatic reallocation
|
||||
// TODO: redo with STL vectors but with explicit mem. management
|
||||
|
||||
class SparseSystemData {
|
||||
private:
|
||||
// matrix representation, A[rowind[i],colind[i]] = vals[i]
|
||||
// right-hand side
|
||||
SparseMatrixData m_A;
|
||||
double *m_b;
|
||||
double *m_x;
|
||||
|
||||
public:
|
||||
SparseMatrixData& A() { return m_A; }
|
||||
double* b() { return m_b ; }
|
||||
double* x() { return m_x ; }
|
||||
unsigned int nrows() { return m_A.nrows(); }
|
||||
|
||||
public:
|
||||
|
||||
SparseSystemData(): m_A(), m_b(NULL), m_x(NULL){ }
|
||||
|
||||
void initialize(unsigned int nr, unsigned int nc) {
|
||||
m_A.initialize(nr,nc);
|
||||
m_b = new double[nr];
|
||||
m_x = new double[nr];
|
||||
assert(m_b);
|
||||
std::fill( m_b, m_b+nr, 0.);
|
||||
}
|
||||
|
||||
void addRHSCmplx(unsigned int i, Cmplx val) {
|
||||
assert( 2*i+1 < m_A.nrows());
|
||||
m_b[2*i] += val.real(); m_b[2*i+1] += val.imag();
|
||||
}
|
||||
|
||||
void setRHSCmplx(unsigned int i, Cmplx val) {
|
||||
assert( 2*i+1 < m_A.nrows());
|
||||
m_b[2*i] = val.real(); m_b[2*i+1] = val.imag();
|
||||
}
|
||||
|
||||
Cmplx getRHSCmplx(unsigned int i) {
|
||||
assert( 2*i+1 < m_A.nrows());
|
||||
return Cmplx( m_b[2*i], m_b[2*i+1]);
|
||||
}
|
||||
|
||||
double getRHSReal(unsigned int i) {
|
||||
assert( i < m_A.nrows());
|
||||
return m_b[i];
|
||||
}
|
||||
|
||||
Cmplx getXCmplx(unsigned int i) {
|
||||
assert( 2*i+1 < m_A.nrows());
|
||||
return Cmplx( m_x[2*i], m_x[2*i+1]);
|
||||
}
|
||||
|
||||
virtual void extend_rows(int extra_rows){
|
||||
if(extra_rows <= 0) return;
|
||||
|
||||
double* old_b = m_b;
|
||||
m_b = new double[m_A.nrows()+extra_rows];
|
||||
if(old_b) { std::copy(old_b, old_b+m_A.nrows(), m_b); delete[] old_b;}
|
||||
std::fill( m_b+m_A.nrows(), m_b+m_A.nrows()+extra_rows, 0.);
|
||||
double* old_x = m_x;
|
||||
m_x = new double[m_A.nrows()+extra_rows];
|
||||
if(old_x) { std::copy(old_x, old_x+m_A.nrows(), m_x); delete[] old_x;}
|
||||
m_A.extend_rows(extra_rows);
|
||||
}
|
||||
|
||||
|
||||
/*void getFromMatlab(const char* nameA, const char* nameb, const char* namex, unsigned int n_vars, bool getsol=false) {
|
||||
MatlabInterface matlab;
|
||||
m_A.getFromMatlab(nameA);
|
||||
assert(n_vars <= m_A.nrows());
|
||||
if (m_b) delete [] m_b;
|
||||
if (m_x) delete [] m_x;
|
||||
m_b = new double[m_A.nrows()];
|
||||
m_x = new double[n_vars];
|
||||
matlab.GetEngineRealMatrix(nameb, m_A.nrows(),1, m_b);
|
||||
if(getsol) {
|
||||
if (m_x) delete [] m_x;
|
||||
matlab.GetEngineRealMatrix(nameb, n_vars,1, m_x);
|
||||
}
|
||||
}*/
|
||||
|
||||
void cleanMem() {
|
||||
m_A.cleanup();
|
||||
delete [] m_b;
|
||||
delete [] m_x;
|
||||
}
|
||||
|
||||
virtual ~SparseSystemData() {
|
||||
delete [] m_b;
|
||||
delete [] m_x;
|
||||
}
|
||||
};
|
||||
|
||||
#endif
|
||||
|
|
@ -0,0 +1,462 @@
|
|||
#ifndef MIQ_VERTEX_INDEXING
|
||||
#define MIQ_VERTEX_INDEXING
|
||||
#include <vcg/complex/complex.h>
|
||||
#include <vcg/simplex/face/pos.h>
|
||||
#include <vcg/simplex/face/jumping_pos.h>
|
||||
#include <vcg/simplex/face/topology.h>
|
||||
|
||||
struct SeamInfo
|
||||
{
|
||||
int v0,v0p,v1,v1p;
|
||||
int integerVal;
|
||||
//unsigned char RotInt;
|
||||
unsigned char MMatch;
|
||||
|
||||
SeamInfo(int _v0,
|
||||
int _v1,
|
||||
int _v0p,
|
||||
int _v1p,
|
||||
int _MMatch,
|
||||
int _integerVal)
|
||||
{
|
||||
v0=_v0;
|
||||
v1=_v1;
|
||||
v0p=_v0p;
|
||||
v1p=_v1p;
|
||||
integerVal=_integerVal;
|
||||
MMatch=_MMatch;
|
||||
}
|
||||
|
||||
SeamInfo(const SeamInfo &S1)
|
||||
{
|
||||
v0=S1.v0;
|
||||
v1=S1.v1;
|
||||
v0p=S1.v0p;
|
||||
v1p=S1.v1p;
|
||||
integerVal=S1.integerVal;
|
||||
MMatch=S1.MMatch;
|
||||
}
|
||||
};
|
||||
|
||||
struct MeshSystemInfo
|
||||
{
|
||||
///total number of scalar variables
|
||||
int num_scalar_variables;
|
||||
////number of vertices variables
|
||||
int num_vert_variables;
|
||||
///num of integer for cuts
|
||||
int num_integer_cuts;
|
||||
///this are used for drawing purposes
|
||||
std::vector<SeamInfo> EdgeSeamInfo;
|
||||
};
|
||||
|
||||
template <class MeshType>
|
||||
class VertexIndexing
|
||||
{
|
||||
|
||||
private:
|
||||
typedef typename MeshType::ScalarType ScalarType;
|
||||
typedef typename MeshType::FaceType FaceType;
|
||||
typedef typename MeshType::VertexType VertexType;
|
||||
typedef typename MeshType::CoordType CoordType;
|
||||
typedef typename MeshType::FaceIterator FaceIterator;
|
||||
|
||||
MeshType *mesh;
|
||||
|
||||
///this maps back index to vertices
|
||||
std::vector<VertexType*> IndexToVert;
|
||||
///this maps the integer for edge
|
||||
typename MeshType::template PerFaceAttributeHandle<vcg::Point3i> Handle_Integer;
|
||||
///per face indexes of vertex in the solver
|
||||
typename MeshType::template PerFaceAttributeHandle<vcg::Point3i> HandleS_Index;
|
||||
///per face per edge of mmatch in the solver
|
||||
typename MeshType::template PerFaceAttributeHandle<vcg::Point3i> Handle_MMatch;
|
||||
///per vertex variable indexes
|
||||
typename MeshType::template PerVertexAttributeHandle<std::vector<int> > HandleV_Integer;
|
||||
///per vertex singular or not
|
||||
typename MeshType::template PerVertexAttributeHandle<bool> Handle_Singular;
|
||||
///per vertex degree of a singularity
|
||||
typename MeshType::template PerVertexAttributeHandle<int> Handle_SingularDegree;
|
||||
///seam per face
|
||||
typename MeshType::template PerFaceAttributeHandle<vcg::Point3<bool> > Handle_Seams;
|
||||
///this handle for mesh
|
||||
typename MeshType::template PerMeshAttributeHandle<MeshSystemInfo> Handle_SystemInfo;
|
||||
|
||||
///this are used for drawing purposes
|
||||
std::vector<VertexType*> duplicated;
|
||||
|
||||
void FirstPos(const VertexType* v,FaceType *&f,int &edge)
|
||||
{
|
||||
f=v->cVFp();
|
||||
edge=v->cVFi();
|
||||
}
|
||||
|
||||
|
||||
int AddNewIndex(VertexType* v0)
|
||||
{
|
||||
Handle_SystemInfo().num_scalar_variables++;
|
||||
HandleV_Integer[v0].push_back(Handle_SystemInfo().num_scalar_variables);
|
||||
IndexToVert.push_back(v0);
|
||||
return Handle_SystemInfo().num_scalar_variables;
|
||||
}
|
||||
|
||||
bool HasIndex(int indexVert,int indexVar)
|
||||
{
|
||||
for (unsigned int i=0;i<HandleV_Integer[indexVert].size();i++)
|
||||
if (HandleV_Integer[indexVert][i]==indexVar)return true;
|
||||
return false;
|
||||
}
|
||||
|
||||
bool HasIndex(VertexType* v0,int indexVar)
|
||||
{
|
||||
for (unsigned int i=0;i<HandleV_Integer[v0].size();i++)
|
||||
if (HandleV_Integer[v0][i]==indexVar)return true;
|
||||
return false;
|
||||
}
|
||||
|
||||
void GetSeamInfo(const FaceType *f0,
|
||||
const FaceType *f1,
|
||||
const int indexE,
|
||||
int &v0,int &v1,
|
||||
int &v0p,int &v1p,
|
||||
unsigned char &_MMatch,
|
||||
int &integerVar)
|
||||
{
|
||||
int edgef0=indexE;
|
||||
v0=HandleS_Index[f0][edgef0];
|
||||
v1=HandleS_Index[f0][(edgef0+1)%3];
|
||||
////get the index on opposite side
|
||||
assert(f0->FFp(edgef0)==f1);
|
||||
int edgef1=f0->cFFi(edgef0);
|
||||
v1p=HandleS_Index[f1][edgef1];
|
||||
v0p=HandleS_Index[f1][(edgef1+1)%3];
|
||||
|
||||
integerVar=Handle_Integer[f0][edgef0];
|
||||
_MMatch=Handle_MMatch[f0][edgef0];
|
||||
assert(f0->cV0(edgef0)==f1->cV1(edgef1));
|
||||
assert(f0->cV1(edgef0)==f1->cV0(edgef1));
|
||||
}
|
||||
|
||||
bool IsSeam(CFace *f0,CFace *f1)
|
||||
{
|
||||
for (int i=0;i<3;i++)
|
||||
{
|
||||
FaceType *f_clos=f0->FFp(i);
|
||||
assert(!f_clos->IsD()); ///check not deleted
|
||||
if (f_clos==f0)continue;///border
|
||||
if(f_clos==f1)
|
||||
return(Handle_Seams[f0][i]);
|
||||
}
|
||||
assert(0);
|
||||
return false;
|
||||
}
|
||||
|
||||
///find initial position of the pos to
|
||||
// assing face to vert inxex correctly
|
||||
void FindInitialPos(const VertexType * vert,
|
||||
int &edge,
|
||||
FaceType *&face)
|
||||
{
|
||||
FaceType *f_init;
|
||||
int edge_init;
|
||||
FirstPos(vert,f_init,edge_init);
|
||||
vcg::face::JumpingPos<FaceType> VFI(f_init,edge_init);
|
||||
bool vertexB=vert->IsB();
|
||||
bool possible_split=false;
|
||||
bool complete_turn=false;
|
||||
do
|
||||
{
|
||||
FaceType *curr_f=VFI.F();
|
||||
int curr_edge=VFI.E();
|
||||
VFI.NextFE();
|
||||
FaceType *next_f=VFI.F();
|
||||
///test if I've just crossed a border
|
||||
bool on_border=(curr_f->FFp(curr_edge)==curr_f);
|
||||
//bool mismatch=false;
|
||||
bool seam=false;
|
||||
|
||||
///or if I've just crossed a seam
|
||||
///if I'm on a border I MUST start from the one next t othe border
|
||||
if (!vertexB)
|
||||
//seam=curr_f->IsSeam(next_f);
|
||||
seam=IsSeam(curr_f,next_f);
|
||||
if (vertexB)
|
||||
assert(!Handle_Singular[vert]);
|
||||
;
|
||||
//assert(!vert->IsSingular());
|
||||
possible_split=((on_border)||(seam));
|
||||
complete_turn=(next_f==f_init);
|
||||
}while ((!possible_split)&&(!complete_turn));
|
||||
face=VFI.F();
|
||||
edge=VFI.E();
|
||||
///test that is not on a border
|
||||
//assert(face->FFp(edge)!=face);
|
||||
}
|
||||
|
||||
///intialize the mapping given an initial pos
|
||||
///whih must be initialized with FindInitialPos
|
||||
void MapIndexes(VertexType * vert,
|
||||
int &edge_init,
|
||||
FaceType *&f_init)
|
||||
{
|
||||
///check that is not on border..
|
||||
///in such case maybe it's non manyfold
|
||||
///insert an initial index
|
||||
int curr_index=AddNewIndex(vert);
|
||||
///and initialize the jumping pos
|
||||
vcg::face::JumpingPos<FaceType> VFI(f_init,edge_init);
|
||||
bool complete_turn=false;
|
||||
do
|
||||
{
|
||||
FaceType *curr_f=VFI.F();
|
||||
int curr_edge=VFI.E();
|
||||
///assing the current index
|
||||
HandleS_Index[curr_f][curr_edge]=curr_index;
|
||||
VFI.NextFE();
|
||||
FaceType *next_f=VFI.F();
|
||||
///test if I've finiseh with the face exploration
|
||||
complete_turn=(next_f==f_init);
|
||||
///or if I've just crossed a mismatch
|
||||
if (!complete_turn)
|
||||
{
|
||||
bool seam=false;
|
||||
//seam=curr_f->IsSeam(next_f);
|
||||
seam=IsSeam(curr_f,next_f);
|
||||
if (seam)
|
||||
{
|
||||
///then add a new index
|
||||
curr_index=AddNewIndex(vert);
|
||||
}
|
||||
}
|
||||
}while (!complete_turn);
|
||||
}
|
||||
|
||||
///intialize the mapping for a given vertex
|
||||
void InitMappingSeam(VertexType *vert)
|
||||
{
|
||||
///first rotate until find the first pos after a mismatch
|
||||
///or a border or return to the first position...
|
||||
FaceType *f_init=vert->VFp();
|
||||
int indexE=vert->VFi();
|
||||
vcg::face::JumpingPos<FaceType> VFI(f_init,indexE);
|
||||
|
||||
int edge_init;
|
||||
FaceType *face_init;
|
||||
FindInitialPos(vert,edge_init,face_init);
|
||||
MapIndexes(vert,edge_init,face_init);
|
||||
}
|
||||
|
||||
///intialize the mapping for a given sampled mesh
|
||||
void InitMappingSeam()
|
||||
{
|
||||
//num_scalar_variables=-1;
|
||||
Handle_SystemInfo().num_scalar_variables=-1;
|
||||
for (unsigned int i=0;i<mesh->vert.size();i++)
|
||||
if (!mesh->vert[i].IsD())
|
||||
InitMappingSeam(&mesh->vert[i]);
|
||||
|
||||
for (unsigned int j=0;j<mesh->vert.size();j++)
|
||||
{
|
||||
VertexType *v= &(mesh->vert[j]);
|
||||
if (!v->IsD()){
|
||||
assert(HandleV_Integer[j].size()>0);
|
||||
if (HandleV_Integer[j].size()>1)
|
||||
duplicated.push_back(v);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
///initialized mapping structures if are not already initialized
|
||||
void AddAttributesIfNeeded()
|
||||
{
|
||||
IndexToVert.clear();
|
||||
|
||||
///other per mesh attributes
|
||||
///see if already exists otherwise it creates it
|
||||
bool HasSystemInfo=vcg::tri::HasPerMeshAttribute(*mesh,std::string("SystemInfo"));
|
||||
if (!HasSystemInfo)
|
||||
Handle_SystemInfo=vcg::tri::Allocator<MeshType>::template AddPerMeshAttribute<MeshSystemInfo>(*mesh,std::string("SystemInfo"));
|
||||
else
|
||||
Handle_SystemInfo=vcg::tri::Allocator<MeshType>::template GetPerMeshAttribute<MeshSystemInfo>(*mesh,"SystemInfo");
|
||||
|
||||
Handle_SystemInfo().num_scalar_variables=0;
|
||||
Handle_SystemInfo().num_vert_variables=0;
|
||||
Handle_SystemInfo().num_integer_cuts=0;
|
||||
|
||||
duplicated.clear();
|
||||
|
||||
bool HasSystemIndex=vcg::tri::HasPerFaceAttribute(*mesh,std::string("SystemIndex"));
|
||||
if (!HasSystemIndex)
|
||||
HandleS_Index = vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<vcg::Point3i>(*mesh,std::string("SystemIndex"));
|
||||
else
|
||||
HandleS_Index = vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3i>(*mesh,std::string("SystemIndex"));
|
||||
|
||||
bool HasHandleInteger=vcg::tri::HasPerFaceAttribute(*mesh,std::string("Integer"));
|
||||
if (!HasHandleInteger)
|
||||
Handle_Integer= vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<vcg::Point3i>(*mesh,std::string("Integer"));
|
||||
else
|
||||
Handle_Integer= vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3i>(*mesh,std::string("Integer"));
|
||||
|
||||
bool HasHandleVertexInteger=vcg::tri::HasPerVertexAttribute(*mesh,std::string("VertInteger"));
|
||||
if (!HasHandleVertexInteger)
|
||||
HandleV_Integer=vcg::tri::Allocator<MeshType>::template AddPerVertexAttribute<std::vector<int> >(*mesh,std::string("VertInteger"));
|
||||
else
|
||||
HandleV_Integer=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<std::vector<int> >(*mesh,std::string("VertInteger"));
|
||||
|
||||
|
||||
bool HasHandleMMatch=vcg::tri::HasPerFaceAttribute(*mesh,std::string("MissMatch"));
|
||||
assert(HasHandleMMatch);
|
||||
Handle_MMatch = vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3i>(*mesh,std::string("MissMatch"));
|
||||
|
||||
bool HasHandleSingularity=vcg::tri::HasPerVertexAttribute(*mesh,std::string("Singular"));
|
||||
assert(HasHandleSingularity);
|
||||
Handle_Singular=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<bool>(*mesh,std::string("Singular"));
|
||||
|
||||
bool HasHandleSingularityDegree=vcg::tri::HasPerVertexAttribute(*mesh,std::string("SingularityDegree"));
|
||||
assert(HasHandleSingularityDegree);
|
||||
Handle_SingularDegree=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<int>(*mesh,std::string("SingularityDegree"));
|
||||
|
||||
bool HasHandleSeams=vcg::tri::HasPerFaceAttribute(*mesh,std::string("Seams"));
|
||||
assert(HasHandleSeams);
|
||||
Handle_Seams=vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3<bool> >(*mesh,std::string("Seams"));
|
||||
|
||||
}
|
||||
|
||||
///test consistency of face variables per vert mapping
|
||||
void TestSeamMapping(FaceType *f)
|
||||
{
|
||||
for (int k=0;k<3;k++)
|
||||
{
|
||||
int indexV=HandleS_Index[f][k];
|
||||
VertexType *v=f->V(k);
|
||||
bool has_index=HasIndex(v,indexV);
|
||||
assert(has_index);
|
||||
}
|
||||
}
|
||||
|
||||
///test consistency of face variables per vert mapping
|
||||
void TestSeamMapping(int indexVert)
|
||||
{
|
||||
VertexType *v=&mesh->vert[indexVert];
|
||||
for (unsigned int k=0;k<HandleV_Integer[indexVert].size();k++)
|
||||
{
|
||||
int indexV=HandleV_Integer[indexVert][k];
|
||||
|
||||
///get faces sharing vertex
|
||||
std::vector<FaceType*> faces;
|
||||
std::vector<int> indexes;
|
||||
|
||||
vcg::face::VFStarVF(v,faces,indexes);
|
||||
for (unsigned int j=0;j<faces.size();j++)
|
||||
{
|
||||
FaceType *f=faces[j];
|
||||
int index=indexes[j];
|
||||
assert(!f->IsD());
|
||||
assert(f->V(index)==v);
|
||||
assert((index>=0)&&(index<3));
|
||||
|
||||
if (HandleS_Index[f][index]==indexV)
|
||||
return;
|
||||
}
|
||||
}
|
||||
assert(0);
|
||||
}
|
||||
|
||||
|
||||
///check consistency of variable mapping across seams
|
||||
void TestSeamMapping()
|
||||
{
|
||||
printf("\n TESTING SEAM INDEXES \n");
|
||||
///test F-V mapping
|
||||
for (unsigned int j=0;j<mesh->face.size();j++)
|
||||
{
|
||||
FaceType *f=&mesh->face[j];
|
||||
if (f->IsD()) continue;
|
||||
TestSeamMapping(f);
|
||||
}
|
||||
|
||||
///TEST V-F MAPPING
|
||||
for (unsigned int j=0;j<mesh->vert.size();j++)
|
||||
{
|
||||
VertexType *v=&mesh->vert[j];
|
||||
if (v->IsD()) continue;
|
||||
TestSeamMapping(j);
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
public:
|
||||
|
||||
///vertex to variable mapping
|
||||
void InitMapping()
|
||||
{
|
||||
//use_direction_field=_use_direction_field;
|
||||
|
||||
IndexToVert.clear();
|
||||
duplicated.clear();
|
||||
|
||||
//ClearMapping();
|
||||
|
||||
InitMappingSeam();
|
||||
|
||||
Handle_SystemInfo().num_vert_variables=Handle_SystemInfo().num_scalar_variables+1;
|
||||
|
||||
///end testing...
|
||||
TestSeamMapping();
|
||||
}
|
||||
|
||||
void InitFaceIntegerVal()
|
||||
{
|
||||
Handle_SystemInfo().num_integer_cuts=0;
|
||||
for (unsigned int j=0;j<mesh->face.size();j++)
|
||||
{
|
||||
FaceType *f=&mesh->face[j];
|
||||
if (f->IsD())continue;
|
||||
//f->IntegerVar.clear();
|
||||
for (int k=0;k<3;k++)
|
||||
{
|
||||
//if (f->IsSeam(k))
|
||||
if (Handle_Seams[f][k])
|
||||
{
|
||||
//f->IntegerVar.push_back(num_integer_cuts);
|
||||
Handle_Integer[j][k]=Handle_SystemInfo().num_integer_cuts;
|
||||
Handle_SystemInfo().num_integer_cuts++;
|
||||
}
|
||||
else
|
||||
Handle_Integer[j][k]=-1;
|
||||
//f->IntegerVar.push_back(-1);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void InitSeamInfo()
|
||||
{
|
||||
Handle_SystemInfo().EdgeSeamInfo.clear();
|
||||
for (unsigned int j=0;j<mesh->face.size();j++)
|
||||
{
|
||||
FaceType *f0=&mesh->face[j];
|
||||
if (f0->IsD())continue;
|
||||
for (int k=0;k<3;k++)
|
||||
{
|
||||
FaceType *f1=f0->FFp(k);
|
||||
if (f1==f0)continue;
|
||||
bool seam=Handle_Seams[f0][k];//f0->IsSeam(k);
|
||||
if ((seam) &&(f0<f1))
|
||||
{
|
||||
int v0,v0p,v1,v1p;
|
||||
unsigned char MM;
|
||||
int integerVar;
|
||||
GetSeamInfo(f0,f1,k,v0,v1,v0p,v1p,MM,integerVar);
|
||||
Handle_SystemInfo().EdgeSeamInfo.push_back(SeamInfo(v0,v1,v0p,v1p,MM,integerVar));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void Init(MeshType *_mesh){mesh=_mesh;AddAttributesIfNeeded();}
|
||||
|
||||
VertexIndexing(){mesh=NULL;}
|
||||
};
|
||||
|
||||
#endif
|
||||
Loading…
Reference in New Issue