diff --git a/wrap/miq/MIQ.h b/wrap/miq/MIQ.h new file mode 100644 index 00000000..c9638806 --- /dev/null +++ b/wrap/miq/MIQ.h @@ -0,0 +1,351 @@ +#ifndef __MIQ__ +#define __MIQ__ + + +#define SIZEQUADS 512 +#define V_SIZE 1 +#define SIZEPARA 1024 + +#include +#include +#include "mesh_type.h" +#include "quadrangulator.h" +#include "poisson_solver.h" +#include "param_stats.h" +#include "seams_initializer.h" +#include "vertex_indexing.h" + +#include +#include + +#define USECOMISO + +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 MIQ{ + +public: + typename MeshType::template PerFaceAttributeHandle 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::Info objInf; + int mask; + vcg::tri::io::ImporterOBJ::LoadMask(PathMesh.c_str(),mask); + err=vcg::tri::io::ImporterOBJ::Open(mesh,PathMesh.c_str(),mask); + assert(position!=-1); + if (err!=0)return false; + } + else + { + err=vcg::tri::io::ImporterPLY::Open(mesh,PathMesh.c_str()); + if (err!=ply::E_NOERROR)return false; + } + ///UPDATE MESH STUFF + vcg::tri::UpdateBounding::Box(mesh); + vcg::tri::UpdateNormal::PerVertexNormalizedPerFace(mesh); + vcg::tri::UpdateNormal::PerFaceNormalized(mesh); + vcg::tri::UpdateTopology::FaceFace(mesh); + vcg::tri::UpdateTopology::VertexFace(mesh); + vcg::tri::UpdateFlags::FaceBorderFromFF(mesh); + vcg::tri::UpdateFlags::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 >& V, std::vector< std::vector >& F); + +// void exportQuadMesh(std::string out); + + //bool LoadData(std::string PathMesh, + // std::string PathFField); + // Bounding box of the param domain + vcg::Box2 UVBox; + + void ScaleGLtoInt() + { + double factor=(double)SIZEPARA/(double)SIZEQUADS; + for (unsigned int i=0;i 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::AddPerFaceAttribute(mesh,std::string("Stiffness")); + + bool hasSingular = vcg::tri::HasPerVertexAttribute(mesh,std::string("Singular")); + assert(hasSingular); + CMesh::PerVertexAttributeHandle Handle_Singular; + Handle_Singular=vcg::tri::Allocator::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()); + std::vector::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 ring; + //mesh.GetConnectedFaces(v,r,ring); + VFExtendedStarVF(v,r,ring); + ///then set stiffening + for (unsigned int k=0;kstiffeningstiffening=stiffVal; + if (Handle_Stiffness[f]::AddPerFaceAttribute(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::AddPerFaceAttribute(mesh,std::string("Stiffness")); + + for(unsigned int i=0;istiffening=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 SInit; + + // Vertex indexing class used for the solver + VertexIndexing VInd; + + // Poisson solver + PoissonSolver PSolver; + +}; + +#endif diff --git a/wrap/miq/auxmath.h b/wrap/miq/auxmath.h new file mode 100644 index 00000000..cbbc032a --- /dev/null +++ b/wrap/miq/auxmath.h @@ -0,0 +1,38 @@ +#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/glUtils.h b/wrap/miq/glUtils.h new file mode 100644 index 00000000..ac555a15 --- /dev/null +++ b/wrap/miq/glUtils.h @@ -0,0 +1,150 @@ +#ifndef MIQ_GL_UTILS +#define MIQ_GL_UTILS +//#include +#include "vertex_indexing.h" +#include "quadrangulator.h" + +class Miq_Gl_Utils +{ +public: + + 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 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 + 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/param_stats.h b/wrap/miq/param_stats.h new file mode 100644 index 00000000..e0775fab --- /dev/null +++ b/wrap/miq/param_stats.h @@ -0,0 +1,199 @@ +#ifndef PARAM_STATS_H +#define PARAM_STATS_H +#ifdef MIQ_USE_ROBUST +#include "predicates.h" +#endif +#include + +template +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 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 Distortion(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(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[i],h); + if (dist>maxD)maxD=dist; + if (dist +#include +#include +#include + +#include "auxmath.h" +#include "sparsesystemdata.h" +#include "mesh_type.h" +#include "vertex_indexing.h" + +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; + } + } + } + + ///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 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/quadrangulator.h b/wrap/miq/quadrangulator.h new file mode 100644 index 00000000..7169ae03 --- /dev/null +++ b/wrap/miq/quadrangulator.h @@ -0,0 +1,517 @@ +#ifndef MIQ_QUADRANGULATOR_H +#define MIQ_QUADRANGULATOR_H + + +#include +#include +#include +#include +#include +#include + +///return the list of faces that share a specified vertex +///together with indexes those faces are sorted in ccw +template +void SortedStar(const vcg::face::Pos &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 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 +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 +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 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 > IntegerEdges; + + ///the set of all integer vertices and the other vertices on integer lines which is connectes to + std::map > IntegerLineAdj; + ///the set of integer vertices + std::set IntegerVertices; + ///temporary polygons + std::vector > polygons; + + ///drawing debug structures + std::vector > IntegerLines; + std::vector IntegerVertex; + + static bool ToSplit(const vcg::Point2 &uv0, + const vcg::Point2 &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 &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 uv0=v0->T().P()*factor; + vcg::Point2 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 + 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; + + ScalarType factor; + ScalarType tolerance; + ScalarType alpha; + + void operator()(typename MESH_TYPE::VertexType &nv, + vcg::face::Pos 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 WedgeInterp(vcg::TexCoord2 &t0, + vcg::TexCoord2 &t1) + { + vcg::TexCoord2 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 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 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 splMd; + EdgePredicate eP; + + splMd.tolerance=tolerance; + splMd.factor=factor; + eP.tolerance=tolerance; + eP.factor=factor; + + while (done) + done=vcg::tri::RefineE,EdgePredicate >(to_split,splMd,eP); + + for (unsigned int i=0;iT().P(); + } + + + bool IsOnIntegerLine(vcg::Point2 uv0, + vcg::Point2 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 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 >(IntegerEdges.begin(),IntegerEdges.end()); + IntegerVertex=std::vector (IntegerVertices.begin(),IntegerVertices.end()); + } + + void EraseIntegerEdge(const vcg::face::Pos &ep) + { + std::pair edge(ep.F(),ep.E()); + assert(IntegerEdges.count(edge)!=0); + IntegerEdges.erase(edge); + } + + void EraseIntegerEdge(const std::vector > &to_erase) + { + for (unsigned int i=0;i 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 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;iIsD())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 uv0=f->WT(j).P()*factor; + vcg::Point2 uv1=f->WT((j+1)%3).P()*factor; + if (IsOnIntegerLine(uv0,uv1,tolerance)||IsBorder) + { + //IntegerEdges.insert(std::pair(v0,v1)); + IntegerEdges.insert(std::pair(f,j)); + if (!IsBorder) + IntegerEdges.insert(std::pair(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 &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 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 &ep, + std::vector &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 > to_erase; + to_erase.push_back(std::pair(ep.F(),ep.E())); + do + { + bool done=OneIntegerStep(ep); + if (!done) + { + EraseIntegerEdge(to_erase); + return false; + } + to_erase.push_back(std::pair(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 > &polygons) + { + //int limit=2; + for (unsigned int i=0;iV0(j); + //TriVertexType* v1=f->V1(j); + + //std::pair edge(v0,v1);*/ + std::pair 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 ep(f,j); + std::vector 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::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, + ScalarType factor=1.0, + ScalarType tolerance=0.000001) + { + TestIsProper(Tmesh); + vcg::tri::AttributeSeam::SplitVertex(Tmesh, ExtractVertex, CompareVertex); + vcg::tri::Allocator::CompactVertexVector(Tmesh); + vcg::tri::Allocator::CompactFaceVector(Tmesh); + vcg::tri::UpdateTopology::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::MergeCloseVertex(Tmesh,EPS); + + 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::VertexBorderFromFace(Tmesh); + ///test manifoldness + TestIsProper(Tmesh); + + vcg::tri::UpdateFlags::VertexClearV(Tmesh); + + InitIntegerEdgesVert(Tmesh,factor,tolerance); + InitVertexQuadMesh(Tmesh); + } +}; + +#endif diff --git a/wrap/miq/seams_initializer.h b/wrap/miq/seams_initializer.h new file mode 100644 index 00000000..a3dc115f --- /dev/null +++ b/wrap/miq/seams_initializer.h @@ -0,0 +1,325 @@ +#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; + + 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 faces; + std::vector edges; + //SortedFaces(v,faces); + vcg::face::VFOrderedStarVF_FF(v,faces,edges); + + missmatch=0; + for (unsigned int i=0;i::template AddPerFaceAttribute(*mesh,std::string("MissMatch")); + else + Handle_MMatch = vcg::tri::Allocator::template GetPerFaceAttribute(*mesh,std::string("MissMatch")); + + bool HasHandleSingularity=vcg::tri::HasPerVertexAttribute(*mesh,std::string("Singular")); + if (!HasHandleSingularity) + Handle_Singular=vcg::tri::Allocator::template AddPerVertexAttribute(*mesh,std::string("Singular")); + else + Handle_Singular=vcg::tri::Allocator::template GetPerVertexAttribute(*mesh,std::string("Singular")); + + bool HasHandleSingularityDegree=vcg::tri::HasPerVertexAttribute(*mesh,std::string("SingularityDegree")); + if (!HasHandleSingularityDegree) + Handle_SingularDegree=vcg::tri::Allocator::template AddPerVertexAttribute(*mesh,std::string("SingularityDegree")); + else + Handle_SingularDegree=vcg::tri::Allocator::template GetPerVertexAttribute(*mesh,std::string("SingularityDegree")); + + bool HasHandleSeams=vcg::tri::HasPerFaceAttribute(*mesh,std::string("Seams")); + if (!HasHandleSeams) + Handle_Seams=vcg::tri::Allocator::template AddPerFaceAttribute >(*mesh,std::string("Seams")); + else + Handle_Seams=vcg::tri::Allocator::template GetPerFaceAttribute >(*mesh,std::string("Seams")); + } + + + 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()) ) + { + //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 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;iface.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;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; + 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++) + { + 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::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::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::LoadFFIELD(mesh,PathFField.c_str()); + if (!field_loaded)return false; + vcg::tri::CrossField::MakeDirectionFaceCoherent(*mesh); + InitMMatch(); + SelectSingularityByMM(); + InitTopologycalCuts(); + AddSeamsByMM(); + return true; + } + + void Init(MeshType *_mesh){mesh=_mesh;}//AllocateMappingStructures();} + + SeamsInitializer(){mesh=NULL;} +}; + +#endif diff --git a/wrap/miq/sparsesystemdata.h b/wrap/miq/sparsesystemdata.h new file mode 100644 index 00000000..d9c87783 --- /dev/null +++ b/wrap/miq/sparsesystemdata.h @@ -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 diff --git a/wrap/miq/vertex_indexing.h b/wrap/miq/vertex_indexing.h new file mode 100644 index 00000000..155339b7 --- /dev/null +++ b/wrap/miq/vertex_indexing.h @@ -0,0 +1,462 @@ +#ifndef MIQ_VERTEX_INDEXING +#define MIQ_VERTEX_INDEXING +#include +#include +#include +#include + +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 EdgeSeamInfo; +}; + +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;iFFp(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 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 + bool HasSystemInfo=vcg::tri::HasPerMeshAttribute(*mesh,std::string("SystemInfo")); + if (!HasSystemInfo) + Handle_SystemInfo=vcg::tri::Allocator::template AddPerMeshAttribute(*mesh,std::string("SystemInfo")); + else + Handle_SystemInfo=vcg::tri::Allocator::template GetPerMeshAttribute(*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::template AddPerFaceAttribute(*mesh,std::string("SystemIndex")); + else + HandleS_Index = vcg::tri::Allocator::template GetPerFaceAttribute(*mesh,std::string("SystemIndex")); + + bool HasHandleInteger=vcg::tri::HasPerFaceAttribute(*mesh,std::string("Integer")); + if (!HasHandleInteger) + Handle_Integer= vcg::tri::Allocator::template AddPerFaceAttribute(*mesh,std::string("Integer")); + else + Handle_Integer= vcg::tri::Allocator::template GetPerFaceAttribute(*mesh,std::string("Integer")); + + bool HasHandleVertexInteger=vcg::tri::HasPerVertexAttribute(*mesh,std::string("VertInteger")); + if (!HasHandleVertexInteger) + HandleV_Integer=vcg::tri::Allocator::template AddPerVertexAttribute >(*mesh,std::string("VertInteger")); + else + 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 GetPerFaceAttribute(*mesh,std::string("MissMatch")); + + bool HasHandleSingularity=vcg::tri::HasPerVertexAttribute(*mesh,std::string("Singular")); + assert(HasHandleSingularity); + Handle_Singular=vcg::tri::Allocator::template GetPerVertexAttribute(*mesh,std::string("Singular")); + + bool HasHandleSingularityDegree=vcg::tri::HasPerVertexAttribute(*mesh,std::string("SingularityDegree")); + assert(HasHandleSingularityDegree); + Handle_SingularDegree=vcg::tri::Allocator::template GetPerVertexAttribute(*mesh,std::string("SingularityDegree")); + + bool HasHandleSeams=vcg::tri::HasPerFaceAttribute(*mesh,std::string("Seams")); + assert(HasHandleSeams); + Handle_Seams=vcg::tri::Allocator::template GetPerFaceAttribute >(*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(); + + //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;jface.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;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