From c38743aaffc09f893a461d09dd871f7300e38155 Mon Sep 17 00:00:00 2001 From: nicopietroni Date: Wed, 12 Nov 2014 15:42:43 +0000 Subject: [PATCH] changed to version on igl linked on wrap --- wrap/miq/MIQ.h | 99 ---- wrap/miq/core/auxmath.h | 38 -- wrap/miq/core/glUtils.h | 264 --------- wrap/miq/core/param_stats.h | 218 -------- wrap/miq/core/poisson_solver.h | 863 ------------------------------ wrap/miq/core/seams_initializer.h | 367 ------------- wrap/miq/core/sparsesystemdata.h | 256 --------- wrap/miq/core/stiffening.h | 170 ------ wrap/miq/core/vertex_indexing.h | 443 --------------- wrap/miq/quadrangulator.h | 740 ------------------------- 10 files changed, 3458 deletions(-) delete mode 100644 wrap/miq/MIQ.h delete mode 100644 wrap/miq/core/auxmath.h delete mode 100644 wrap/miq/core/glUtils.h delete mode 100644 wrap/miq/core/param_stats.h delete mode 100644 wrap/miq/core/poisson_solver.h delete mode 100644 wrap/miq/core/seams_initializer.h delete mode 100644 wrap/miq/core/sparsesystemdata.h delete mode 100644 wrap/miq/core/stiffening.h delete mode 100644 wrap/miq/core/vertex_indexing.h delete mode 100644 wrap/miq/quadrangulator.h diff --git a/wrap/miq/MIQ.h b/wrap/miq/MIQ.h deleted file mode 100644 index 79a3bff4..00000000 --- a/wrap/miq/MIQ.h +++ /dev/null @@ -1,99 +0,0 @@ -#ifndef __MIQ__ -#define __MIQ__ - -#include -#include "quadrangulator.h" -#include "core/poisson_solver.h" -#include "core/param_stats.h" -#include "core/seams_initializer.h" -#include "core/vertex_indexing.h" -#include "core/stiffening.h" -#include - -#define USECOMISO - -template -class MIQ_parametrization{ - -public: - typename MeshType::template PerFaceAttributeHandle Handle_Stiffness; - - // Different stiffening mode - enum StiffMode{NO_STIFF = 0,GAUSSIAN = 1,ITERATIVE = 2}; - - // Parametrize the mesh - static void DoParameterize(MeshType &mesh,StiffMode stiffMode, - double Stiffness = 5.0,double GradientSize = 30.0, - bool DirectRound = false,int iter = 5, - int localIter = 5, bool DoRound = true) - { - vcg::PoissonSolver PSolver(mesh); - if (mesh.fn==0)return; - StiffeningInitializer::InitDefaultStiffening(mesh); - if (stiffMode==GAUSSIAN) - { - StiffeningInitializer::AddGaussStiffening(mesh,Stiffness); - PSolver.SolvePoisson(GradientSize,1.f,DirectRound,localIter,DoRound); - } - else - if (stiffMode==ITERATIVE) - { - for (int i=0;i::updateStiffeningJacobianDistorsion(mesh,GradientSize); - 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); - - } - -public: - - static bool IsValid(MeshType &mesh) - { - int n_comp=vcg::tri::Clean::CountConnectedComponents(mesh); - int non_manifE=vcg::tri::Clean::CountNonManifoldEdgeFF(mesh); - int non_manifV=vcg::tri::Clean::CountNonManifoldVertexFF(mesh); - return ((n_comp==1)&&(non_manifE==0)&&(non_manifV==0)); - } - static void InitSeamsSing(MeshType &mesh, - bool orient_globally, - bool initMM, - bool initCuts) - { - SeamsInitializer SInit; - SInit.Init(&mesh,orient_globally,initMM,initCuts); - } - - static void Parametrize(MeshType &mesh,StiffMode stiffMode, - double Stiffness = 5.0, - double GradientSize = 30.0, - bool DirectRound = false,int iter = 5, - int localIter = 5, bool DoRound = true) - { - VertexIndexing VInd; - - VInd.Init(&mesh); - VInd.InitMapping(); - VInd.InitFaceIntegerVal(); - VInd.InitSeamInfo(); - - DoParameterize(mesh,stiffMode,Stiffness,GradientSize,DirectRound,iter,localIter , DoRound); - } - - -}; - -#endif diff --git a/wrap/miq/core/auxmath.h b/wrap/miq/core/auxmath.h deleted file mode 100644 index cbbc032a..00000000 --- a/wrap/miq/core/auxmath.h +++ /dev/null @@ -1,38 +0,0 @@ -#ifndef AUXMATH_H -#define AUXMATH_H -#pragma once - -#include -#include -//#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 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 Cmplx; - -#endif // AUXMATH_H - diff --git a/wrap/miq/core/glUtils.h b/wrap/miq/core/glUtils.h deleted file mode 100644 index 3eb05ea8..00000000 --- a/wrap/miq/core/glUtils.h +++ /dev/null @@ -1,264 +0,0 @@ -#ifndef MIQ_GL_UTILS -#define MIQ_GL_UTILS -//#include -#include "vertex_indexing.h" - -class Miq_Gl_Utils -{ -public: - - ///singular vertices should be selected - template - static void GLDrawSingularities(MeshType &mesh, - float size=8, - bool DrawUV=false) - { - bool hasSingular = vcg::tri::HasPerVertexAttribute(mesh,std::string("Singular")); - bool hasSingularDegree = vcg::tri::HasPerVertexAttribute(mesh,std::string("SingularityDegree")); - if (!hasSingular)return; - - typename MeshType::template PerVertexAttributeHandle Handle_Singular; - typename MeshType::template PerVertexAttributeHandle Handle_SingularDegree; - - Handle_Singular=vcg::tri::Allocator::template GetPerVertexAttribute(mesh,std::string("Singular")); - - Handle_SingularDegree=vcg::tri::Allocator::template GetPerVertexAttribute(mesh,std::string("SingularityDegree")); - - glPushAttrib(GL_ALL_ATTRIB_BITS); - glEnable(GL_COLOR_MATERIAL); - glDisable(GL_LIGHTING); - glDepthRange(0,0.999); - glPointSize(size); - glBegin(GL_POINTS); - glColor4d(0,1,1,0.7); - for (unsigned int j=0;jIsD())continue; - for (int i=0;i<3;i++) - { - CVertex *v=f->V(i); - if (!Handle_Singular[v])continue; - //int mmatch=3; - if (hasSingularDegree) - { - int mmatch=Handle_SingularDegree[v]; - if (mmatch==1)vcg::glColor(vcg::Color4b(255,0,0,255));///valence 5 - if (mmatch==2)vcg::glColor(vcg::Color4b(0,255,0,255));///valence 6 - if (mmatch==3)vcg::glColor(vcg::Color4b(0,0,255,255));///valence 3 - } - if (!DrawUV) - vcg::glVertex(v->P()); - else - glVertex3d(f->WT(i).P().X(),f->WT(i).P().Y(),0); - } - } - - glEnd(); - glPopAttrib(); - } - - template - static void GLDrawFaceSeams(const FaceType &f, - vcg::Point3 seams, - vcg::Color4b seamCol[3], - float size=3, - bool UV=false) - { - typedef typename FaceType::CoordType CoordType; - typedef typename FaceType::ScalarType ScalarType; - glLineWidth(size); - glBegin(GL_LINES); - for (int i=0;i<3;i++) - { - if (!seams[i])continue; - vcg::glColor(seamCol[i]); - if (!UV) - { - glVertex(f.cV0(i)->P()); - glVertex(f.cV1(i)->P()); - } - else - { - int i0=i; - int i1=(i0+1)%3; - CoordType p0=CoordType(f.cWT(i0).P().X(),f.cWT(i0).P().Y(),0); - CoordType p1=CoordType(f.cWT(i1).P().X(),f.cWT(i1).P().Y(),0); - glVertex(p0); - glVertex(p1); - } - } - glEnd(); - } - - template - static void GLDrawSeams(MeshType &mesh, - float size=3, - bool UV=false, - int numCuts=400) - { - bool hasSeam = vcg::tri::HasPerFaceAttribute(mesh,std::string("Seams")); - if(!hasSeam)return; - bool HasSeamIndex=vcg::tri::HasPerFaceAttribute(mesh,std::string("SeamsIndex")); - - typedef typename MeshType::template PerFaceAttributeHandle > SeamsHandleType; - typedef typename MeshType::template PerFaceAttributeHandle SeamsIndexHandleType; - - typedef typename vcg::tri::Allocator SeamsAllocator; - - SeamsHandleType Handle_Seam; - Handle_Seam=SeamsAllocator::template GetPerFaceAttribute >(mesh,std::string("Seams")); - - SeamsIndexHandleType Handle_SeamIndex; - if (HasSeamIndex) - Handle_SeamIndex=SeamsAllocator::template GetPerFaceAttribute(mesh,std::string("SeamsIndex")); - - glPushAttrib(GL_ALL_ATTRIB_BITS); - glEnable(GL_COLOR_MATERIAL); - glDisable(GL_LIGHTING); - - glDepthRange(0,0.999); - for (unsigned int i=0;i seams=Handle_Seam[i]; - vcg::Color4b seamCol[3]; - for (int j=0;j<3;j++) - { - seamCol[j]=vcg::Color4b(0,255,0,255); - if (HasSeamIndex) - { - int index=Handle_SeamIndex[i][j]; - //assert(index>0); - if (index>=0) - seamCol[j]=vcg::Color4b(255,0,0,255); - //seamCol[j]=vcg::Color4b::Scatter(numCuts,index); - } - } - - GLDrawFaceSeams(mesh.face[i],seams,seamCol,size,UV); - - } - glPopAttrib(); - } - - ///this is useful to debug the output - ///of the vertex indexing class - template - 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;iIsD()); - vcg::glVertex(VI.duplicated[i]->P()); - } - glEnd(); - glPopAttrib(); - } - - - template - 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;iP(); - if (v->IsV()) - glColor3d(1,0,0); - else - glColor3d(1,1,0); - glVertex(pos); - } - glEnd(); - - glPopAttrib(); - } - - template - 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;iV0(edge); - typename QuadrangulatorType::TriVertexType* v1=f->V1(edge); - glBegin(GL_LINES); - glVertex(v0->P()); - glVertex(v1->P()); - glEnd(); - } - - glPopAttrib(); - } - - template - 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;iN()); - glVertex(v->P()); - } - glEnd(); - } - - glDepthRange(0,0.997); - glDisable(GL_LIGHTING); - glEnable(GL_COLOR_MATERIAL); - glColor3d(0,0,0); - for (unsigned int i=0;iP()); - } - glEnd(); - } - - glPopAttrib(); - } -}; -#endif diff --git a/wrap/miq/core/param_stats.h b/wrap/miq/core/param_stats.h deleted file mode 100644 index 45162dc7..00000000 --- a/wrap/miq/core/param_stats.h +++ /dev/null @@ -1,218 +0,0 @@ -#ifndef PARAM_STATS_H -#define PARAM_STATS_H -#ifdef MIQ_USE_ROBUST -#include -#endif -#include - -template -inline bool IsFlipped(const CoordType2D &uv0, - const CoordType2D &uv1, - const CoordType2D &uv2) -{ - #ifdef MIQ_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 t2(uv0,uv1,uv2); - return (!t2.IsCCW()); - #endif -} - -template -inline bool IsFlipped( FaceType &f) -{ - - typedef typename FaceType::ScalarType ScalarType; - typedef typename vcg::Point2 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 -void SelectFlippedFaces(MeshType &Tmesh) -{ - typedef typename MeshType::ScalarType ScalarType; - typedef typename vcg::Point2 CoordType2D; - vcg::tri::UpdateFlags::FaceClearS(Tmesh); - for (unsigned int i=0;i -typename FaceType::ScalarType LamdaDistortion(FaceType &f,typename FaceType::ScalarType h) -{ - typedef typename FaceType::CoordType CoordType3D; - typedef typename FaceType::ScalarType ScalarType; - typedef typename vcg::Point2 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((ScalarType)0.0, (ScalarType)(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; - if ((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 = LamdaDistortion(f, h); - ScalarType lapl=0; - for (int i=0;i<3;i++) - lapl += (mydist- LamdaDistortion(*f.FFp(i), h)); - return lapl; -} - -template< class MeshType> -void UpdateUVBox(MeshType &mesh) -{ - typedef typename MeshType::ScalarType ScalarType; - typename MeshType::template PerMeshAttributeHandle > Handle_UVBox; - - Handle_UVBox=vcg::tri::Allocator::template GetPerMeshAttribute >(mesh,std::string("UVBox")); - - Handle_UVBox().SetNull(); - for (unsigned int i=0;i -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[i],h); - if (dist>maxD)maxD=dist; - if (dist -#include -#include -#include - -#include "auxmath.h" -#include "sparsesystemdata.h" -#include "vertex_indexing.h" - -namespace vcg -{ -template -class PoissonSolver -{ - typedef typename MeshType::ScalarType ScalarType; - typedef typename MeshType::FaceType FaceType; - typedef typename MeshType::VertexType VertexType; - typedef typename MeshType::CoordType CoordType; - //typedef VertexIndexing VertexIndexingType; - //typedef typename VertexIndexingType::SeamInfo SeamInfo; - - typename MeshType::template PerFaceAttributeHandle HandleS_Index; - typename MeshType::template PerVertexAttributeHandle Handle_Singular; - typename MeshType::template PerFaceAttributeHandle Handle_Stiffness; - ///this handle for mesh - typename MeshType::template PerMeshAttributeHandle 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 Hard_constraints; - - ///vector of indexes to round - std::vector 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(*f,0); - //K2=CrossVector(*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;jIsD())continue; - //if (v->IsSingular()) - if (Handle_Singular[v]) - { - AddFixedVertex(v); - v->T().P()=vcg::Point2(0,0); - return; - } - } - ///if anything fixed fix the first - AddFixedVertex(&mesh.vert[0]); - mesh.vert[0].T().P()=vcg::Point2(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;iIsD()); - ///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;jIsD())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;jIsD()) - 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;jIsD())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 uv=vcg::Point2(U,V); - ///assing - //f->realUV[k]=uv; -// ScalarType factor=(ScalarType)SIZEQUADS/(ScalarType)SIZEPARA; -// uv*=factor; - ///assing - f->WT(k).P()=uv; - } - } - ///initialize the vector of integer variables to return their values -// Handle_SystemInfo().IntegerValues.resize(n_integer_vars*2); -// int baseIndex=(n_vert_vars)*2; -// int endIndex=baseIndex+n_integer_vars*2; -// int index=0; -// for (int i=baseIndex;i EdgeSeamInfo; - //VIndex.GetSeamInfo(EdgeSeamInfo); - for (unsigned int i=0;i((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::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::template GetPerFaceAttribute(mesh,"SystemIndex"); - Handle_Singular=vcg::tri::Allocator::template GetPerVertexAttribute(mesh,std::string("Singular")); - Handle_Stiffness=vcg::tri::Allocator::template GetPerFaceAttribute(mesh,std::string("Stiffness")); - Handle_SystemInfo=vcg::tri::Allocator::template GetPerMeshAttribute(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 diff --git a/wrap/miq/core/seams_initializer.h b/wrap/miq/core/seams_initializer.h deleted file mode 100644 index 9ea74638..00000000 --- a/wrap/miq/core/seams_initializer.h +++ /dev/null @@ -1,367 +0,0 @@ -#ifndef MIQ_SEAMS_INTIALIZER -#define MIQ_SEAMS_INTIALIZER -#include -#include -#include -#include - -template -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 Handle_MMatch; - ///per vertex singular or not - typename MeshType::template PerVertexAttributeHandle Handle_Singular; - ///per vertex degree of a singularity - typename MeshType::template PerVertexAttributeHandle Handle_SingularDegree; - ///seam per face - typename MeshType::template PerFaceAttributeHandle > Handle_Seams; - ///seam index per face - typename MeshType::template PerFaceAttributeHandle Handle_SeamsIndex; - - 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 VertexType &v,int &missmatch) - { - ///check that is on border.. - if (v.IsB()) return false; - - vcg::face::Pos pos(v.cVFp(), v.cVFi()); - std::vector > posVec; - vcg::face::VFOrderedStarFF(pos, posVec); - - missmatch=0; - for (unsigned int i=0;i::template GetPerFaceAttribute(*mesh,std::string("MissMatch")); - Handle_Singular=vcg::tri::Allocator::template GetPerVertexAttribute(*mesh,std::string("Singular")); - Handle_SingularDegree=vcg::tri::Allocator::template GetPerVertexAttribute(*mesh,std::string("SingularityDegree")); - Handle_Seams=vcg::tri::Allocator::template GetPerFaceAttribute >(*mesh,std::string("Seams")); - Handle_SeamsIndex=vcg::tri::Allocator::template GetPerFaceAttribute(*mesh,std::string("SeamsIndex")); - } - - - void FloodFill(FaceType* start) - { - std::deque 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()) ) - { - Handle_Seams[f][s]=false; - Handle_Seams[g][j]=false; - g->SetV(); - d.push_back(g); - } - } - } - } - - void Retract(){ - std::vector 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); - } - - void AddSeamsByMM() - { - for (unsigned int i=0;iface.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;ivert.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; - if (missmatch==3)missmatch=1; - else - if (missmatch==1)missmatch=3; - Handle_SingularDegree[i]=missmatch; - } - else - { - mesh->vert[i].ClearS(); - Handle_Singular[i]=false; - Handle_SingularDegree[i]=0; - } - } - } - - - int InitTopologycalCuts(){ - vcg::tri::UpdateFlags::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;iface.size();i++) - { - FaceType *curr_f=&mesh->face[i]; - for (int j=0;j<3;j++) - { - FaceType *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::MissMatchByCross(*curr_f,*opp_f); - } - } - } - - void InitSeamIndexes() - { - ///initialize seams indexes - for (unsigned int i=0;iface.size();i++) - { - FaceType *f=&mesh->face[i]; - for (int j=0;j<3;j++) - Handle_SeamsIndex[f][j]=-1; - } - - std::vector > > seamsVert; - seamsVert.resize(mesh->vert.size()); - for (unsigned int i=0;iface.size();i++) - { - FaceType *f=&mesh->face[i]; - for (int j=0;j<3;j++) - { - if (f->IsB(j))continue; - if (Handle_Seams[f][j]) - { - VertexType *v0=f->V0(j); - VertexType *v1=f->V1(j); - if (v0>v1)continue; - int index0=v0-&mesh->vert[0]; - int index1=v1-&mesh->vert[0]; - std::pair entry=std::pair(f,j); - seamsVert[index0].push_back(entry); - seamsVert[index1].push_back(entry); - } - } - } - - int curr_index=0; - for (unsigned int i=0;ivert.size();i++) - { - VertexType *v_seam=&mesh->vert[i]; - bool IsVertex=(Handle_Singular[i])||(seamsVert[i].size()>2)|| - (v_seam->IsB()); - - if(!IsVertex)continue; - - ///then follows the seam - for (unsigned int j=0;jFFp(edge); - int edge1=f->FFi(edge); - Handle_SeamsIndex[f1][edge1]=curr_index;///set current value - - assert((v0==f->V0(edge))|| - (v0==f->V1(edge))); - - if (f->V0(edge)==v0) - v1=f->V1(edge); - else - v1=f->V0(edge); - - assert(v0!=v1); - //int index0=v0-&mesh->vert[0]; - int index1=v1-&mesh->vert[0]; - - bool IsVertex=(Handle_Singular[v1])|| - (seamsVert[index1].size()>2)|| - (v1->IsB()); - if(IsVertex) - finished=true; - else - { - assert(seamsVert[index1].size()==2); - FaceType *f_new[2]; - int edge_new[2]; - - f_new[0]=seamsVert[index1][0].first; - edge_new[0]=seamsVert[index1][0].second; - f_new[1]=seamsVert[index1][1].first; - edge_new[1]=seamsVert[index1][1].second; - - assert((f_new[0]==f)||(f_new[1]==f)); - assert((edge_new[0]==edge)||(edge_new[1]==edge)); - - if ((f_new[0]==f)&&(edge_new[0]==edge)) - { - f=f_new[1]; - edge=edge_new[1]; - } - else - { - f=f_new[0]; - edge=edge_new[0]; - } - v0=v1; - } - } - while (!finished); - //return; - curr_index++; - } - } - - } - -public: - - - void Init(MeshType *_mesh, - bool orient_globally, - bool initMM, - bool initCuts) - { - - mesh=_mesh; - - vcg::tri::UpdateTopology::FaceFace(*_mesh); - vcg::tri::UpdateFlags::FaceBorderFromFF(*_mesh); - vcg::tri::UpdateFlags::VertexBorderFromFace(*_mesh); - - - AddAttributesIfNeeded(); - if (orient_globally) - vcg::tri::CrossField::MakeDirectionFaceCoherent(*mesh); - if (initMM) - InitMMatch(); - SelectSingularityByMM(); - if (initCuts) - { - InitTopologycalCuts(); - AddSeamsByMM(); - } - //InitSeamIndexes(); - } - - SeamsInitializer(){mesh=NULL;} -}; - -#endif diff --git a/wrap/miq/core/sparsesystemdata.h b/wrap/miq/core/sparsesystemdata.h deleted file mode 100644 index d9c87783..00000000 --- a/wrap/miq/core/sparsesystemdata.h +++ /dev/null @@ -1,256 +0,0 @@ -#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 diff --git a/wrap/miq/core/stiffening.h b/wrap/miq/core/stiffening.h deleted file mode 100644 index 180b914e..00000000 --- a/wrap/miq/core/stiffening.h +++ /dev/null @@ -1,170 +0,0 @@ -#ifndef MIQ_STIFFENING -#define MIQ_STIFFENING - - -template -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 StiffeningInitializer -{ - typedef typename MeshType::ScalarType ScalarType; - typedef typename MeshType::FaceType FaceType; - typedef typename MeshType::VertexType VertexType; - typedef typename MeshType::CoordType CoordType; - -public: - - static void colorByStiffening(MeshType & mesh, - typename MeshType::ScalarType MaxVal=16) - { - bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness")); - assert(hasStiffness); - typename MeshType::template PerFaceAttributeHandle Handle_Stiffness=vcg::tri::Allocator::template GetPerFaceAttribute(mesh,std::string("Stiffness")); - - for (unsigned int i=0;i Handle_Stiffness; - - bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness")); - if(!hasStiffness) - Handle_Stiffness=vcg::tri::Allocator::template AddPerFaceAttribute(mesh,std::string("Stiffness")); - else - Handle_Stiffness=vcg::tri::Allocator::template FindPerFaceAttribute(mesh,std::string("Stiffness")); - - bool hasSingular = vcg::tri::HasPerVertexAttribute(mesh,std::string("Singular")); - assert(hasSingular); - - typename MeshType::template PerVertexAttributeHandle Handle_Singular; - Handle_Singular=vcg::tri::Allocator:: template GetPerVertexAttribute(mesh,std::string("Singular")); - - std::vector to_stiff; - for(unsigned int i=0;iIsD())continue; - //if (!v->IsSingular())continue; - if (!Handle_Singular[v])continue; - to_stiff.push_back(v); - } - for(unsigned int i=0;iIsD())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()); - typename std::vector::iterator new_end; - 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 ring; - //mesh.GetConnectedFaces(v,r,ring); - VFExtendedStarVF(v,r,ring); - ///then set stiffening - for (unsigned int k=0;kstiffeningstiffening=stiffVal; - if (Handle_Stiffness[f] Handle_Stiffness; - - bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness")); - if(!hasStiffness) - Handle_Stiffness=vcg::tri::Allocator::template AddPerFaceAttribute(mesh,std::string("Stiffness")); - else - Handle_Stiffness=vcg::tri::Allocator::template FindPerFaceAttribute(mesh,std::string("Stiffness")); - - bool flipped = NumFlips(mesh)>0; - //if (h == 0.0) - // return flipped; - // - //assert(h != 0.0); - - if (!flipped) - return false; - ScalarType maxL=0; - 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) - { - ScalarType dist=LamdaDistortion(mesh.face[i],grad_size); - if (dist>maxD)maxD=dist; - ScalarType absLap=fabs(LaplaceDistortion(mesh.face[i], grad_size)); - if (absLap>maxL)maxL=absLap; - 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; - } - - static void InitDefaultStiffening(MeshType & mesh) - { - typename MeshType::template PerFaceAttributeHandle Handle_Stiffness; - - bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness")); - if(!hasStiffness) - Handle_Stiffness=vcg::tri::Allocator::template AddPerFaceAttribute(mesh,std::string("Stiffness")); - else - Handle_Stiffness=vcg::tri::Allocator::template FindPerFaceAttribute(mesh,std::string("Stiffness")); - - for(unsigned int i=0;istiffening=1; - Handle_Stiffness[f]=1; - } - } - -}; - -#endif diff --git a/wrap/miq/core/vertex_indexing.h b/wrap/miq/core/vertex_indexing.h deleted file mode 100644 index 6ffb7f1b..00000000 --- a/wrap/miq/core/vertex_indexing.h +++ /dev/null @@ -1,443 +0,0 @@ -#ifndef MIQ_VERTEX_INDEXING -#define MIQ_VERTEX_INDEXING -#include -#include -#include -#include - - - -struct SeamInfo -{ - int v0,v0p,v1,v1p; - int integerVar; - //unsigned char RotInt; - unsigned char MMatch; - - SeamInfo(int _v0, - int _v1, - int _v0p, - int _v1p, - int _MMatch, - int _integerVar) - { - v0=_v0; - v1=_v1; - v0p=_v0p; - v1p=_v1p; - integerVar=_integerVar; - MMatch=_MMatch; - } - - SeamInfo(const SeamInfo &S1) - { - v0=S1.v0; - v1=S1.v1; - v0p=S1.v0p; - v1p=S1.v1p; - integerVar=S1.integerVar; - 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 EdgeSeamInfo; - ///this are values of integer variables after optimization - std::vector IntegerValues; -}; - -template -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 IndexToVert; - ///this maps the integer for edge - typename MeshType::template PerFaceAttributeHandle Handle_Integer; - ///per face indexes of vertex in the solver - typename MeshType::template PerFaceAttributeHandle HandleS_Index; - ///per face per edge of mmatch in the solver - typename MeshType::template PerFaceAttributeHandle Handle_MMatch; - ///per vertex variable indexes - typename MeshType::template PerVertexAttributeHandle > HandleV_Integer; - ///per vertex singular or not - typename MeshType::template PerVertexAttributeHandle Handle_Singular; - ///per vertex degree of a singularity - typename MeshType::template PerVertexAttributeHandle Handle_SingularDegree; - ///seam per face - typename MeshType::template PerFaceAttributeHandle > Handle_Seams; - ///this handle for mesh - typename MeshType::template PerMeshAttributeHandle Handle_SystemInfo; - - ///this are used for drawing purposes - std::vector 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;icFFp(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(FaceType *f0,FaceType *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 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 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 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;ivert.size();i++) - if (!mesh->vert[i].IsD()) - InitMappingSeam(&mesh->vert[i]); - - for (unsigned int j=0;jvert.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 - Handle_SystemInfo=vcg::tri::Allocator::template GetPerMeshAttribute(*mesh,std::string("SystemInfo")); - - Handle_SystemInfo().num_scalar_variables=0; - Handle_SystemInfo().num_vert_variables=0; - Handle_SystemInfo().num_integer_cuts=0; - - duplicated.clear(); - - HandleS_Index = vcg::tri::Allocator::template GetPerFaceAttribute(*mesh,std::string("SystemIndex")); - - Handle_Integer= vcg::tri::Allocator::template GetPerFaceAttribute(*mesh,std::string("Integer")); - HandleV_Integer=vcg::tri::Allocator::template GetPerVertexAttribute >(*mesh,std::string("VertInteger")); - - - bool HasHandleMMatch=vcg::tri::HasPerFaceAttribute(*mesh,std::string("MissMatch")); - assert(HasHandleMMatch); - Handle_MMatch = vcg::tri::Allocator::template FindPerFaceAttribute(*mesh,std::string("MissMatch")); - - bool HasHandleSingularity=vcg::tri::HasPerVertexAttribute(*mesh,std::string("Singular")); - assert(HasHandleSingularity); - Handle_Singular=vcg::tri::Allocator::template FindPerVertexAttribute(*mesh,std::string("Singular")); - - bool HasHandleSingularityDegree=vcg::tri::HasPerVertexAttribute(*mesh,std::string("SingularityDegree")); - assert(HasHandleSingularityDegree); - Handle_SingularDegree=vcg::tri::Allocator::template FindPerVertexAttribute(*mesh,std::string("SingularityDegree")); - - bool HasHandleSeams=vcg::tri::HasPerFaceAttribute(*mesh,std::string("Seams")); - assert(HasHandleSeams); - Handle_Seams=vcg::tri::Allocator::template FindPerFaceAttribute >(*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 faces; - std::vector indexes; - - vcg::face::VFStarVF(v,faces,indexes); - for (unsigned int j=0;jIsD()); - 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;jface.size();j++) - { - FaceType *f=&mesh->face[j]; - if (f->IsD()) continue; - TestSeamMapping(f); - } - - ///TEST V-F MAPPING - for (unsigned int j=0;jvert.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(); - - 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;jface.size();j++) - { - FaceType *f=&mesh->face[j]; - if (f->IsD())continue; - for (int k=0;k<3;k++) - { - if (Handle_Seams[f][k]) - { - Handle_Integer[j][k]=Handle_SystemInfo().num_integer_cuts; - Handle_SystemInfo().num_integer_cuts++; - } - else - Handle_Integer[j][k]=-1; - } - } - } - - void InitSeamInfo() - { - Handle_SystemInfo().EdgeSeamInfo.clear(); - for (unsigned int j=0;jface.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 -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#define precisionQ 0.0000000001 - -namespace vcg { -namespace tri { -template -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; - - - struct InterpolationInfo - { - CoordType Pos3D; - vcg::Point2 PosUV; - ScalarType alpha; - bool to_split; - - InterpolationInfo() - { - Pos3D=CoordType(0,0,0); - PosUV=vcg::Point2(0,0); - to_split=false; - alpha=-1; - } - }; - - //the interpolation map that is saved once to be univoque per edge - typedef std::pair KeyEdgeType; - - std::map InterpMap; - - //ScalarType UVtolerance; - -private: - - bool ToSplit(const vcg::Point2 &uv0, - const vcg::Point2 &uv1, - int Dir, - ScalarType &alpha) - { - ScalarType val0=uv0.V(Dir); - ScalarType val1=uv1.V(Dir); - int IntegerLine0=floor(val0); - int IntegerLine1=floor(val1); - if (IntegerLine0==IntegerLine1) - return false;//no integer line pass throught the edge - - bool swapped=false; - if (IntegerLine0>IntegerLine1) - { - std::swap(IntegerLine0,IntegerLine1); - std::swap(val0,val1); - assert(val1>=val0); - swapped=true; - } - - //then get the first if extist that overcome the threshold - int IntegerSplit=IntegerLine0+1; - bool found=false; - ScalarType dist1,dist0; - for (int i=IntegerSplit;i<=IntegerLine1;i++) - { - dist1=fabs(val1-IntegerSplit); - dist0=fabs(val0-IntegerSplit); - -// if ((dist0>=UVtolerance)&& -// (dist1>=UVtolerance)) - if ((val0!=IntegerSplit)&& - (val1!=IntegerSplit)) - { - found=true; - break; - } - IntegerSplit++; - } - if (!found)return false; - - //have to check distance also in opposite direction - ScalarType lenght=val1-val0; - assert(lenght>=0); - //alpha=1.0-(dist/lenght); - alpha=(dist1/lenght); - if (swapped)alpha=1-alpha; - assert((alpha>0)&&(alpha<1)); - return true; - } - - void RoundInitial(TriMesh &to_split) - { - ScalarType minTolerance=precisionQ; - //first add all eddge - for (int i=0;i UV=f->WT(j).P(); - - int int0=floor(UV.X()+0.5); - int int1=floor(UV.Y()+0.5); - - ScalarType diff0=(fabs(UV.X()-(ScalarType)int0)); - ScalarType diff1=(fabs(UV.Y()-(ScalarType)int1)); - - if (diff0WT(j).P()=UV; - } - } - } - - void RoundSplits(TriMesh &to_split,int dir) - { - ScalarType minTolerance=precisionQ; - //first add all eddge - for (size_t i=0;iP0(j); - CoordType p1=f->P1(j); - KeyEdgeType k(p0,p1); - assert(InterpMap.count(k)==1); - if (!InterpMap[k].to_split)continue; - //then get the intepolated value - vcg::Point2 UV=InterpMap[k].PosUV; - - int int0=floor(UV.X()+0.5); - int int1=floor(UV.Y()+0.5); - - ScalarType diff0=(fabs(UV.X()-(ScalarType)int0)); - ScalarType diff1=(fabs(UV.Y()-(ScalarType)int1)); - - if (diff0P0(j); - CoordType p1=f->P1(j); - vcg::Point2 Uv0=f->V0(j)->T().P(); - vcg::Point2 Uv1=f->V1(j)->T().P(); - KeyEdgeType k(p0,p1); -// printf("p0 (%5.5f,%5.5f,%5.5f) p1(%5.5f,%5.5f,%5.5f) \n",p0.X(),p0.Y(),p0.Z(),p1.X(),p1.Y(),p1.Z()); -// printf("uv0 (%5.5f,%5.5f) uv1(%5.5f,%5.5f) \n",Uv0.X(),Uv0.Y(),Uv1.X(),Uv1.Y()); -// fflush(stdout); - assert(InterpMap.count(k)==0); - InterpMap[k]=InterpolationInfo(); - } - } - - //then set the ones to be splitted - for (size_t i=0;iP0(j); - CoordType p1=f->P1(j); - vcg::Point2 uv0=f->V0(j)->T().P(); - vcg::Point2 uv1=f->V1(j)->T().P(); - - ScalarType alpha; - if (!ToSplit(uv0,uv1,dir,alpha))continue; - - KeyEdgeType k(p0,p1); - assert(InterpMap.count(k)==1); - InterpMap[k].Pos3D=p0*alpha+p1*(1-alpha); - InterpMap[k].PosUV=uv0*alpha+uv1*(1-alpha); - InterpMap[k].to_split=true; - InterpMap[k].alpha=alpha; - } - } - - //then make them coherent - for (size_t i=0;iP0(j); - CoordType p1=f->P1(j); - vcg::Point2 uv0=f->V0(j)->T().P(); - vcg::Point2 uv1=f->V1(j)->T().P(); -// if (p0>p1)continue; //only one verse of coherence - - KeyEdgeType k0(p0,p1); - assert(InterpMap.count(k0)==1);//there should be already in the - //table and it should be coherent - - KeyEdgeType k1(p1,p0); - if(InterpMap.count(k1)==0)continue;//REAL border, no need for update - - bool to_split0=InterpMap[k0].to_split; - bool to_split1=InterpMap[k1].to_split; - - //the find all possible cases - if ((!to_split0)&&(!to_split1))continue; - - if ((to_split0)&&(to_split1)) - { - CoordType Pos3D=InterpMap[k1].Pos3D; - InterpMap[k0].Pos3D=Pos3D; - - //check if need to make coherent also the UV Position - //skip the fake border and do the rest - bool IsBorderFF=(f->FFp(j)==f); - - if (!IsBorderFF) //in this case they should have same UVs - InterpMap[k0].PosUV=InterpMap[k1].PosUV; - else - { - ScalarType alpha=InterpMap[k1].alpha; - assert((alpha>=0)&&(alpha<=1)); - alpha=1-alpha; - InterpMap[k0].PosUV=alpha*uv0+(1-alpha)*uv1; - InterpMap[k0].alpha=alpha; - } - - } - else - if ((!to_split0)&&(to_split1)) - { - CoordType Pos3D=InterpMap[k1].Pos3D; - InterpMap[k0].Pos3D=Pos3D; - - //check if need to make coherent also the UV Position - //skip the fake border and do the rest - bool IsBorderFF=(f->FFp(j)==f); - - InterpMap[k0].to_split=true; - - if (!IsBorderFF) //in this case they should have same UVs - InterpMap[k0].PosUV=InterpMap[k1].PosUV; - else //recalculate , it pass across a seam - { - ScalarType alpha=InterpMap[k1].alpha; - assert((alpha>=0)&&(alpha<=1)); - alpha=1-alpha; - InterpMap[k0].PosUV=alpha*uv0+(1-alpha)*uv1; - InterpMap[k0].alpha=alpha; - } - } - } - } - - RoundSplits(to_split,dir); - } - - // 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 - struct SplitMidPoint : public std::unary_function , typename MESH_TYPE::CoordType > - { - typedef typename MESH_TYPE::VertexType VertexType; - typedef typename MESH_TYPE::FaceType FaceType; - typedef typename MESH_TYPE::CoordType CoordType; - - std::map *MapEdge; - - void operator()(typename MESH_TYPE::VertexType &nv, - vcg::face::Pos ep) - { - VertexType* v0=ep.f->V0(ep.z); - VertexType* v1=ep.f->V1(ep.z); - assert(v0!=v1); - - CoordType p0=v0->P(); - CoordType p1=v1->P(); - assert(p0!=p1); - - KeyEdgeType k(p0,p1); - bool found=(MapEdge->count(k)==1); - assert(found); - bool to_split=(*MapEdge)[k].to_split; - assert(to_split); - - //get the value on which the edge must be splitted - nv.P()= (*MapEdge)[k].Pos3D; - //nv.N()= v0->N()*alpha+v1->N()*(1.0-alpha); - nv.T().P()=(*MapEdge)[k].PosUV; - } - - vcg::TexCoord2 WedgeInterp(vcg::TexCoord2 &t0, - vcg::TexCoord2 &t1) - { - return (vcg::TexCoord2(0,0)); - } - - SplitMidPoint(std::map *_MapEdge){MapEdge=_MapEdge;} - }; - - template - class EdgePredicate - { - typedef typename MESH_TYPE::VertexType VertexType; - typedef typename MESH_TYPE::FaceType FaceType; - typedef typename MESH_TYPE::ScalarType ScalarType; - - std::map *MapEdge; - - public: - - bool operator()(vcg::face::Pos ep) const - { - VertexType* v0=ep.f->V0(ep.z); - VertexType* v1=ep.f->V1(ep.z); - assert(v0!=v1); - - CoordType p0=v0->P(); - CoordType p1=v1->P(); - assert(p0!=p1); - - KeyEdgeType k(p0,p1); - bool found=(MapEdge->count(k)==1); - assert(found); - bool to_split=(*MapEdge)[k].to_split; - return(to_split); - } - - EdgePredicate(std::map *_MapEdge){MapEdge=_MapEdge;} - }; - - void SplitTrisDir(TriMesh &to_split, - int dir) - { - bool done=true; - //int step=0; - while (done) - { - printf("Number of Vertices %d \n",to_split.vn); - fflush(stdout); - - InitSplitMap(to_split,dir); - - SplitMidPoint splMd(&InterpMap); - EdgePredicate eP(&InterpMap); - - done=vcg::tri::RefineE,EdgePredicate >(to_split,splMd,eP); - - } - printf("Number of Vertices %d \n",to_split.vn); - fflush(stdout); - fflush(stdout); - } - - - bool IsOnIntegerLine(vcg::Point2 uv0, - vcg::Point2 uv1) - { - 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))>=UVtolerance)continue; -// if ((fabs(val1-(ScalarType)integer1))>=UVtolerance)continue; - if (val0!=(ScalarType)floor(val0))continue; - if (val1!=(ScalarType)floor(val1))continue; - return true; - } - return false; - } - - bool IsOnIntegerVertex(vcg::Point2 uv, - bool IsB) - { - int onIntegerL=0; - for (int dir=0;dir<2;dir++) - { - ScalarType val0=uv.V(dir); - int integer0=floor(val0+0.5); - //if ((fabs(val0-(ScalarType)integer0))0))return true; - return (onIntegerL==2); - } - - - void InitIntegerEdgesVert(TriMesh &Tmesh) - { - //IntegerEdges.clear(); - vcg::tri::UpdateFlags::FaceSetF(Tmesh); - vcg::tri::UpdateFlags::FaceClearS(Tmesh); - vcg::tri::UpdateFlags::VertexClearS(Tmesh); - - for (unsigned int i=0;iIsD())continue; - for (int j=0;j<3;j++) - { - bool IsBorder=f->IsB(j); - if (IsBorder) - f->ClearF(j); - else - { - vcg::Point2 uv0=f->WT(j).P(); - vcg::Point2 uv1=f->WT((j+1)%3).P(); - - if (IsOnIntegerLine(uv0,uv1)) - { - f->ClearF(j); - TriFaceType *f1=f->FFp(j); - int z=f->FFi(j); - assert(f1!=f); - f1->ClearF(z); - } - } - - bool BorderV=f->V(j)->IsB(); - - if (IsOnIntegerVertex(f->WT(j).P(),BorderV)) - f->V(j)->SetS(); - } - } - } - - short int AlignmentEdge(TriFaceType *f, - int edge_index) - { - vcg::Point2 uv0=f->WT(edge_index).P(); - vcg::Point2 uv1=f->WT((edge_index+1)%3).P(); - if (uv0.X()==uv1.X())return 0; - if (uv0.Y()==uv1.Y())return 1; - return -1; - } - - void FindPolygon(vcg::face::Pos &currPos, - std::vector &poly, - std::vector &UVpoly) - { - currPos.F()->SetV(); - currPos.F()->C()=vcg::Color4b(255,0,0,255); - poly.clear(); - assert(currPos.V()->IsS()); - TriVertexType *v_init=currPos.V(); - poly.push_back(currPos.V()); - - //retrieve UV - int indexV0=currPos.E(); - - short int Align=AlignmentEdge(currPos.F(),currPos.E()); - - std::vector TempUVpoly; - TempUVpoly.push_back(Align); - - do - { - currPos.NextNotFaux(); - currPos.F()->SetV(); - currPos.F()->C()=vcg::Color4b(255,0,0,255); - - if ((currPos.V()->IsS())&&(currPos.V()!=v_init)) - { - poly.push_back(currPos.V()); - - short int Align=AlignmentEdge(currPos.F(),currPos.E()); - - TempUVpoly.push_back(Align); - } - - }while (currPos.V()!=v_init); - - //then shift the order of UV by one - //to be consistent with edge ordering - int size=TempUVpoly.size(); - for (int i=0;i > &polygons, - std::vector > &UV) - { - vcg::tri::UpdateFlags::FaceClearV(Tmesh); - for (unsigned int i=0;iIsV())continue; - - for (int j=0;j<3;j++) - { - TriVertexType* v0=f->V0(j); - if (!v0->IsS())continue; - if (f->IsF(j))continue; - - vcg::face::Pos startPos(f,j); - - std::vector poly; - std::vector< short int> UVpoly; - - FindPolygon(startPos,poly,UVpoly); - - if (poly.size()>2) - { - assert(poly.size()==UVpoly.size()); - std::reverse(poly.begin(),poly.end()); -// std::reverse(UVpoly.begin(),UVpoly.end()); - - polygons.push_back(poly); - UV.push_back(UVpoly); - - } - //only one polygon per initial face - break; - } - } - - } - - //FUNCTIONS NEEDED BY "UV WEDGE TO VERTEX" FILTER - static void ExtractVertex(const TriMesh & srcMesh, - const TriFaceType & f, - int whichWedge, - const TriMesh & dstMesh, - TriVertexType & v) - { - (void)srcMesh; - (void)dstMesh; - // This is done to preserve every single perVertex property - // perVextex Texture Coordinate is instead obtained from perWedge one. - v.ImportData(*f.cV(whichWedge)); - v.T() = f.cWT(whichWedge); - } - - static bool CompareVertex(const TriMesh & m, - TriVertexType & vA, - TriVertexType & vB) - { - (void)m; - return (vA.cT() == vB.cT()); - } - - void ConvertWTtoVT(TriMesh &Tmesh) - { - int vn = Tmesh.vn; - vcg::tri::AttributeSeam::SplitVertex(Tmesh, ExtractVertex, CompareVertex); - vcg::tri::UpdateTopology::FaceFace(Tmesh); - // vcg::tri::UpdateFlags::FaceBorderFromFF(Tmesh); - } - - void ConvertVTtoWT(TriMesh &Tmesh) - { - vcg::tri::UpdateTexture::WedgeTexFromVertexTex(Tmesh); - vcg::tri::Clean::RemoveDuplicateVertex(Tmesh); - } - - void ReupdateMesh(TriMesh &Tmesh) - { - vcg::tri::UpdateNormal::PerFaceNormalized(Tmesh); // update Normals - vcg::tri::UpdateNormal::PerVertexNormalized(Tmesh);// update Normals - //compact the mesh - vcg::tri::Allocator::CompactVertexVector(Tmesh); - vcg::tri::Allocator::CompactFaceVector(Tmesh); - vcg::tri::UpdateTopology::FaceFace(Tmesh); // update Topology - vcg::tri::UpdateTopology::TestFaceFace(Tmesh); //and test it - //set flags - vcg::tri::UpdateFlags::VertexClearV(Tmesh); - vcg::tri::UpdateFlags::FaceBorderFromFF(Tmesh); - vcg::tri::UpdateFlags::VertexBorderFromFaceBorder(Tmesh); - } - -public: - - - void TestIsProper(TriMesh &Tmesh) - { - - - //test manifoldness - int test=vcg::tri::Clean::CountNonManifoldVertexFF(Tmesh); - //assert(test==0); - if (test != 0) - cerr << "Assertion failed: TestIsProper NonManifoldVertices!" << endl; - - test=vcg::tri::Clean::CountNonManifoldEdgeFF(Tmesh); - //assert(test==0); - if (test != 0) - cerr << "Assertion failed: TestIsProper NonManifoldEdges" << endl; - - for (unsigned int i=0;iIsD()); - 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, - std::vector< std::vector< short int> > &UV) - { - UV.clear(); - Pmesh.Clear(); - - TestIsProper(Tmesh); - - RoundInitial(Tmesh); - - //UVtolerance=tolerance; - - //split to per vert - ConvertWTtoVT(Tmesh); - - - vcg::tri::Allocator::CompactVertexVector(Tmesh); - vcg::tri::Allocator::CompactFaceVector(Tmesh); - vcg::tri::UpdateTopology::FaceFace(Tmesh); - (void)Pmesh; - //TestIsProper(Tmesh); - - //then split the tris along X - SplitTrisDir(Tmesh,0); - SplitTrisDir(Tmesh,1); - - //merge back the mesh and WT coords - ConvertVTtoWT(Tmesh); - - //CleanMesh(Pmesh); - - //update properties of the mesh - ReupdateMesh(Tmesh); - - //test manifoldness - TestIsProper(Tmesh); - - InitIntegerEdgesVert(Tmesh); - - for (int i=0;i > polygons; - FindPolygons(Tmesh,polygons,UV); - - //then add to the polygonal mesh - Pmesh.Clear(); - - int numV=vcg::tri::UpdateSelection::VertexCount(Tmesh); - - //first create vertices - vcg::tri::Allocator::AddVertices(Pmesh,numV); - - std::map VertMap; - int index=0; - for(unsigned int i=0;i UV=Tmesh.vert[i].T().P(); - Pmesh.vert[index].P()=typename PolyMesh::CoordType(pos.X(),pos.Y(),pos.Z()); - Pmesh.vert[index].N()=typename PolyMesh::CoordType(norm.X(),norm.Y(),norm.Z()); - Pmesh.vert[index].T().P()=UV; - VertMap[pos]=index; - index++; - } - - //then add polygonal mesh - vcg::tri::Allocator::AddFaces(Pmesh,polygons.size()); - for (unsigned int i=0;iP(); - assert(VertMap.count(pos)==1); - int index=VertMap[pos]; - Pmesh.face[i].V(j)=&(Pmesh.vert[index]); - } - } - - - } -}; -} -} -#endif