manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

changed to version on igl linked on wrap

This commit is contained in:
Nico Pietroni 2014-11-12 15:42:43 +00:00
parent 836d1dd12d
commit c38743aaff
10 changed files with 0 additions and 3458 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

99
wrap/miq/MIQ.h
View File

@ -1,99 +0,0 @@
#ifndef __MIQ__
#define __MIQ__
#include <iostream>
#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 <vcg/complex/algorithms/clean.h>
#define USECOMISO
template <class MeshType>
class MIQ_parametrization{
public:
typename MeshType::template PerFaceAttributeHandle<float> 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<MeshType> PSolver(mesh);
if (mesh.fn==0)return;
StiffeningInitializer<MeshType>::InitDefaultStiffening(mesh);
if (stiffMode==GAUSSIAN)
{
StiffeningInitializer<MeshType>::AddGaussStiffening(mesh,Stiffness);
PSolver.SolvePoisson(GradientSize,1.f,DirectRound,localIter,DoRound);
}
else
if (stiffMode==ITERATIVE)
{
for (int i=0;i<iter;i++)
{
PSolver.SolvePoisson(GradientSize,1.f,DirectRound,localIter,DoRound);
int nflips=NumFlips(mesh);
bool folded=StiffeningInitializer<MeshType>::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<MeshType>::CountConnectedComponents(mesh);
int non_manifE=vcg::tri::Clean<MeshType>::CountNonManifoldEdgeFF(mesh);
int non_manifV=vcg::tri::Clean<MeshType>::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<MeshType> 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<MeshType> VInd;
VInd.Init(&mesh);
VInd.InitMapping();
VInd.InitFaceIntegerVal();
VInd.InitSeamInfo();
DoParameterize(mesh,stiffMode,Stiffness,GradientSize,DirectRound,iter,localIter , DoRound);
}
};
#endif

38
wrap/miq/core/auxmath.h
View File

@ -1,38 +0,0 @@
#ifndef AUXMATH_H
#define AUXMATH_H
#pragma once
#include <complex>
#include <cmath>
//#include "Claussen.h"
static const double pi = 3.1415926535897932384626433;
const double SQRT2_2 = 0.7071067811865; // sqrt(2)/2 = 1/sqrt(2)
const double EPSILON = 1.0e-8;
template <class T> T sqr( const T& a ) { return a*a; }
inline double l2s( const double angle ){
const double s = 2*fabs(sin(angle));
return s < 1e-40 ? -92.103403719761827360719658 : log(s);
}
//
//inline double Lob( const double angle ){
// return 0 <= angle && angle <= pi ? claussen(2*angle)/2 : -1e10;
//}
inline double PropLength( const double l, const double an, const double ad ){
return l*(sin(an)/sin(ad));
}
inline double cot( const double x ){
const double a = fabs(fmod(x+pi,2*pi)-pi)<1e-40 ? 1e-40 : x;
return 1./tan(a);
}
typedef std::complex<double> Cmplx;
#endif // AUXMATH_H

264
wrap/miq/core/glUtils.h
View File

@ -1,264 +0,0 @@
#ifndef MIQ_GL_UTILS
#define MIQ_GL_UTILS
//#include <wrap/gl/space.h>
#include "vertex_indexing.h"
class Miq_Gl_Utils
{
public:
///singular vertices should be selected
template <class MeshType>
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<bool> Handle_Singular;
typename MeshType::template PerVertexAttributeHandle<int> Handle_SingularDegree;
Handle_Singular=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<bool>(mesh,std::string("Singular"));
Handle_SingularDegree=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<int>(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;j<mesh.face.size();j++)
{
CFace *f=&mesh.face[j];
if (f->IsD())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 <class FaceType>
static void GLDrawFaceSeams(const FaceType &f,
vcg::Point3<bool> 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 <class MeshType>
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<vcg::Point3<bool> > SeamsHandleType;
typedef typename MeshType::template PerFaceAttributeHandle<vcg::Point3i > SeamsIndexHandleType;
typedef typename vcg::tri::Allocator<MeshType> SeamsAllocator;
SeamsHandleType Handle_Seam;
Handle_Seam=SeamsAllocator::template GetPerFaceAttribute<vcg::Point3<bool> >(mesh,std::string("Seams"));
SeamsIndexHandleType Handle_SeamIndex;
if (HasSeamIndex)
Handle_SeamIndex=SeamsAllocator::template GetPerFaceAttribute<vcg::Point3i >(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<mesh.face.size();i++)
{
if (mesh.face[i].IsD())continue;
vcg::Point3<bool> 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 <class VertexIndexingType>
static void GLDrawVertexIndexing(VertexIndexingType &VI)
{
typedef typename VertexIndexingType::ScalarType ScalarType;
glPushAttrib(GL_ALL_ATTRIB_BITS);
glEnable(GL_COLOR_MATERIAL);
glDisable(GL_LIGHTING);
glDepthRange(0,0.999);
typename VertexIndexingType::ScalarType size=5;
glPointSize(size);
vcg::glColor(vcg::Color4b(0,255,0,255));
glBegin(GL_POINTS);
for (unsigned int i=0;i<VI.duplicated.size();i++)
{
assert(!VI.duplicated[i]->IsD());
vcg::glVertex(VI.duplicated[i]->P());
}
glEnd();
glPopAttrib();
}
template <class QuadrangulatorType>
static void QuadGLDrawIntegerVertices(QuadrangulatorType &Quadr)
{
glPushAttrib(GL_ALL_ATTRIB_BITS);
glDisable(GL_LIGHTING);
glEnable(GL_COLOR_MATERIAL);
glDisable(GL_TEXTURE_2D);
glPointSize(8);
glDepthRange(0,0.997);
/*glColor3d(1,0,0);*/
glBegin(GL_POINTS);
for (int i=0;i<Quadr.IntegerVertex.size();i++)
{
typename QuadrangulatorType::TriVertexType* v=Quadr.IntegerVertex[i];
typename QuadrangulatorType::CoordType pos=v->P();
if (v->IsV())
glColor3d(1,0,0);
else
glColor3d(1,1,0);
glVertex(pos);
}
glEnd();
glPopAttrib();
}
template <class QuadrangulatorType>
static void GLDrawIntegerLines(QuadrangulatorType &Quadr)
{
glPushAttrib(GL_ALL_ATTRIB_BITS);
glDisable(GL_LIGHTING);
glEnable(GL_COLOR_MATERIAL);
glDisable(GL_TEXTURE_2D);
glLineWidth(2);
glColor3d(0,1,0);
glDepthRange(0,0.998);
for (int i=0;i<Quadr.IntegerLines.size();i++)
{
typename QuadrangulatorType::TriFaceType *f=Quadr.IntegerLines[i].first;
int edge=Quadr.IntegerLines[i].second;
typename QuadrangulatorType::TriVertexType* v0=f->V0(edge);
typename QuadrangulatorType::TriVertexType* v1=f->V1(edge);
glBegin(GL_LINES);
glVertex(v0->P());
glVertex(v1->P());
glEnd();
}
glPopAttrib();
}
template <class QuadrangulatorType>
static void GLDrawPolygons(QuadrangulatorType &Quadr)
{
glPushAttrib(GL_ALL_ATTRIB_BITS);
glEnable(GL_LIGHTING);
glEnable(GL_COLOR_MATERIAL);
glDisable(GL_TEXTURE_2D);
glColor3d(0.7,0.8,0.9);
//glFrontFace(GL_CW);
glDepthRange(0,0.998);
for (unsigned int i=0;i<Quadr.polygons.size();i++)
{
glBegin(GL_POLYGON);
for (unsigned int j=0;j<Quadr.polygons[i].size();j++)
{
typename QuadrangulatorType::TriVertexType* v=Quadr.polygons[i][j];
glNormal(v->N());
glVertex(v->P());
}
glEnd();
}
glDepthRange(0,0.997);
glDisable(GL_LIGHTING);
glEnable(GL_COLOR_MATERIAL);
glColor3d(0,0,0);
for (unsigned int i=0;i<Quadr.polygons.size();i++)
{
glBegin(GL_LINE_LOOP);
for (unsigned int j=0;j<Quadr.polygons[i].size();j++)
{
typename QuadrangulatorType::TriVertexType* v=Quadr.polygons[i][j];
glVertex(v->P());
}
glEnd();
}
glPopAttrib();
}
};
#endif

218
wrap/miq/core/param_stats.h
View File

@ -1,218 +0,0 @@
#ifndef PARAM_STATS_H
#define PARAM_STATS_H
#ifdef MIQ_USE_ROBUST
#include <predicates.h>
#endif
#include <vcg/space/triangle2.h>
template<typename CoordType2D>
inline bool IsFlipped(const CoordType2D &uv0,
const CoordType2D &uv1,
const CoordType2D &uv2)
{
#ifdef 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<typename CoordType2D::ScalarType> t2(uv0,uv1,uv2);
return (!t2.IsCCW());
#endif
}
template<typename FaceType>
inline bool IsFlipped( FaceType &f)
{
typedef typename FaceType::ScalarType ScalarType;
typedef typename vcg::Point2<ScalarType> CoordType2D;
CoordType2D uv0=f.WT(0).P();
CoordType2D uv1=f.WT(1).P();
CoordType2D uv2=f.WT(2).P();
return (IsFlipped(uv0,uv1,uv2));
}
template< class MeshType>
int NumFlips(MeshType &Tmesh)
{
int numFl=0;
for (unsigned int i=0;i<Tmesh.face.size();i++)
{
if (Tmesh.face[i].IsD())continue;
if (IsFlipped(Tmesh.face[i]))numFl++;
}
return numFl;
}
template< class MeshType>
void SelectFlippedFaces(MeshType &Tmesh)
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename vcg::Point2<ScalarType> CoordType2D;
vcg::tri::UpdateFlags<MeshType>::FaceClearS(Tmesh);
for (unsigned int i=0;i<Tmesh.face.size();i++)
{
CoordType2D uv0=Tmesh.face[i].WT(0).P();
CoordType2D uv1=Tmesh.face[i].WT(1).P();
CoordType2D uv2=Tmesh.face[i].WT(2).P();
if (IsFlipped(uv0,uv1,uv2) ) Tmesh.face[i].SetS();
}
}
// Compute the lambda (distortion) of a triangle in the current
// parameerization from the Mixed Integer paper.
// f facet on which to compute distortion
// h scaling factor applied to cross field
// return the triangle's distortion
template< class FaceType>
typename FaceType::ScalarType LamdaDistortion(FaceType &f,typename FaceType::ScalarType h)
{
typedef typename FaceType::CoordType CoordType3D;
typedef typename FaceType::ScalarType ScalarType;
typedef typename vcg::Point2<ScalarType> CoordType2D;
typedef typename FaceType::ScalarType ScalarType;
//Facet_const_handle self = f;
assert(h > 0);
//TriangleIndex ftri = m_flattenedTriangles[f.index()];
//TriangleIndex tri = m_triangles[f.index()];
CoordType2D uv0 = f.WT(0).P();
CoordType2D uv1 = f.WT(1).P();
CoordType2D uv2 = f.WT(2).P();
CoordType3D p0 = f.P(0);
CoordType3D p1 = f.P(1);
CoordType3D p2 = f.P(2);
CoordType3D norm = (p1 - p0) ^ (p2 - p0);
ScalarType area2 = (norm).Norm();
ScalarType area2_inv = 1.0 / area2;
norm *= area2_inv;
if (area2 > 0) {
// Singular values of the Jacobian
CoordType3D neg_t0 = norm ^ (p2 - p1);
CoordType3D neg_t1 = norm ^ ( p0 - p2);
CoordType3D neg_t2 = norm ^ ( p1 - p0);
/*CoordType3D diffu = area2_inv * (uv0.X() * neg_t0 +
uv1.X() * neg_t1 + uv2.X() * neg_t2);
CoordType3D diffv = area2_inv * (uv0.Y() * neg_t0 +
uv1.Y() * neg_t1 + uv2.Y() * neg_t2);*/
CoordType3D diffu = (neg_t0 * uv0.X() +neg_t1 *uv1.X() + neg_t2 * uv2.X() )*area2_inv;
CoordType3D diffv = (neg_t0 * uv0.Y() + neg_t1*uv1.Y() + neg_t2*uv2.Y() )*area2_inv;
// first fundamental form
ScalarType I00 = diffu*diffu; // guaranteed non-neg
ScalarType I01 = diffu*diffv; // I01 = I10
ScalarType I11 = diffv*diffv; // guaranteed non-neg
// eigenvalues of a 2x2 matrix
// [a00 a01]
// [a10 a11]
// 1/2 * [ (a00 + a11) +/- sqrt((a00 - a11)^2 + 4 a01 a10) ]
ScalarType trI = I00 + I11; // guaranteed non-neg
ScalarType diffDiag = I00 - I11; // guaranteed non-neg
ScalarType sqrtDet = sqrt(std::max((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<FaceType>(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<vcg::Box2<ScalarType> > Handle_UVBox;
Handle_UVBox=vcg::tri::Allocator<MeshType>::template GetPerMeshAttribute<vcg::Box2<ScalarType> >(mesh,std::string("UVBox"));
Handle_UVBox().SetNull();
for (unsigned int i=0;i<mesh.face.size();i++)
{
if (mesh.face[i].IsD())continue;
Handle_UVBox().Add(mesh.face[i].WT(0).P());
Handle_UVBox().Add(mesh.face[i].WT(1).P());
Handle_UVBox().Add(mesh.face[i].WT(2).P());
}
}
template< class MeshType>
void SetFaceQualityByDistortion(MeshType &Tmesh,
typename MeshType::ScalarType h)
{
typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::ScalarType ScalarType;
ScalarType minD=100;
ScalarType maxD=0;
///evaluate min and max
for (unsigned int i=0;i<Tmesh.face.size();i++)
{
ScalarType dist=LamdaDistortion<FaceType>(Tmesh.face[i],h);
if (dist>maxD)maxD=dist;
if (dist<minD)minD=dist;
Tmesh.face[i].Q()= dist;
}
for (unsigned int i=0;i<Tmesh.face.size();i++)
{
Tmesh.face[i].Q()= maxD-Tmesh.face[i].Q();
}
printf("Distortion MIN %4.4f \n",(float)minD);
printf("Distortion MAX %4.4f \n",(float)maxD);
}
#endif // PARAM_STATS_H

863
wrap/miq/core/poisson_solver.h
View File

@ -1,863 +0,0 @@
#ifndef MIQ_POISSON_SOLVER
#define MIQ_POISSON_SOLVER
#include <gmm/gmm.h>
#include <ConstrainedSolver.hh>
#include <MISolver.hh>
#include <GMM_Tools.hh>
#include "auxmath.h"
#include "sparsesystemdata.h"
#include "vertex_indexing.h"
namespace vcg
{
template <class MeshType>
class PoissonSolver
{
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::CoordType CoordType;
//typedef VertexIndexing<MeshType> VertexIndexingType;
//typedef typename VertexIndexingType::SeamInfo SeamInfo;
typename MeshType::template PerFaceAttributeHandle<vcg::Point3i> HandleS_Index;
typename MeshType::template PerVertexAttributeHandle<bool> Handle_Singular;
typename MeshType::template PerFaceAttributeHandle<float> Handle_Stiffness;
///this handle for mesh
typename MeshType::template PerMeshAttributeHandle<MeshSystemInfo> Handle_SystemInfo;
///range map data
MeshType &mesh;
///vertex indexing structure
//VertexIndexingType &VIndex;
///solver data
SparseSystemData S;
///vector of unknowns
std::vector< double > X;
////REAL PART
///number of fixed vertex
unsigned int n_fixed_vars;
///the number of REAL variables for vertices
unsigned int n_vert_vars;
///total number of variables of the system,
///do not consider constraints, but consider integer vars
unsigned int num_total_vars;
///the number of scalar variables
//unsigned int n_scalar_vars;
//////INTEGER PART
///the total number of integer variables
unsigned int n_integer_vars;
///CONSTRAINT PART
///number of cuts constraints
unsigned int num_cut_constraint;
///total number of constraints equations
unsigned int num_constraint_equations;
///total size of the system including constraints
unsigned int system_size;
///if you intend to make integer rotation
///and translations
bool integer_jumps_bary;
///vector of blocked vertices
std::vector<VertexType*> Hard_constraints;
///vector of indexes to round
std::vector<int> ids_to_round;
///boolean that is true if rounding to integer is needed
bool integer_rounding;
///START SYSTEM ACCESS METHODS
///add an entry to the LHS
void AddValA(int Xindex,
int Yindex,
ScalarType val)
{
int size=(int)S.nrows();
assert(0 <= Xindex && Xindex < size);
assert(0 <= Yindex && Yindex < size);
S.A().addEntryReal(Xindex,Yindex,val);
}
///add a complex entry to the LHS
void AddComplexA(int VarXindex,
int VarYindex,
Cmplx val)
{
int size=(int)S.nrows()/2;
assert(0 <= VarXindex && VarXindex < size);
assert(0 <= VarYindex && VarYindex < size);
S.A().addEntryCmplx(VarXindex,VarYindex,val);
}
///add a velue to the RHS
void AddValB(int Xindex,
ScalarType val)
{
int size=(int)S.nrows();
assert(0 <= Xindex && Xindex < size);
S.b()[Xindex] += val;
}
///add the area term, scalefactor is used to sum up
///and normalize on the overlap zones
void AddAreaTerm(int index[3][3][2],ScalarType ScaleFactor)
{
const ScalarType entry=0.5*ScaleFactor;
ScalarType val[3][3]= { {0, entry, -entry},
{-entry, 0, entry},
{entry, -entry, 0} };
for (int i=0;i<3;i++)
for (int j=0;j<3;j++)
{
///add for both u and v
int Xindex=index[i][j][0]*2;
int Yindex=index[i][j][1]*2;
AddValA(Xindex+1,Yindex,-val[i][j]);
AddValA(Xindex,Yindex+1,val[i][j]);
}
}
///set the diagonal of the matrix (which is zero at the beginning)
///such that the sum of a row or a colums is zero
void SetDiagonal(ScalarType val[3][3])
{
for (int i=0;i<3;i++)
{
ScalarType sum=0;
for (int j=0;j<3;j++)
sum+=val[i][j];
val[i][i]=-sum;
}
}
///given a vector of scalar values and
///a vector of indexes add such values
///as specified by the indexes
void AddRHS(ScalarType b[6],
int index[3])
{
for (int i=0;i<3;i++)
{
ScalarType valU=b[i*2];
ScalarType valV=b[(i*2)+1];
AddValB((index[i]*2),valU);
AddValB((index[i]*2)+1,valV);
}
}
///add a 3x3 block matrix to the system matrix...
///indexes are specified in the 3x3 matrix of x,y pairs
///indexes must be multiplied by 2 cause u and v
void Add33Block(ScalarType val[3][3],int index[3][3][2])
{
for (int i=0;i<3;i++)
for (int j=0;j<3;j++)
{
///add for both u and v
int Xindex=index[i][j][0]*2;
int Yindex=index[i][j][1]*2;
assert((unsigned)Xindex<(n_vert_vars*2));
assert((unsigned)Yindex<(n_vert_vars*2));
AddValA(Xindex,Yindex,val[i][j]);
AddValA(Xindex+1,Yindex+1,val[i][j]);
}
}
///add a 3x3 block matrix to the system matrix...
///indexes are specified in the 3x3 matrix of x,y pairs
///indexes must be multiplied by 2 cause u and v
void Add44Block(ScalarType val[4][4],int index[4][4][2])
{
for (int i=0;i<4;i++)
for (int j=0;j<4;j++)
{
///add for both u and v
int Xindex=index[i][j][0]*2;
int Yindex=index[i][j][1]*2;
assert((unsigned)Xindex<(n_vert_vars*2));
assert((unsigned)Yindex<(n_vert_vars*2));
AddValA(Xindex,Yindex,val[i][j]);
AddValA(Xindex+1,Yindex+1,val[i][j]);
}
}
///END SYSTEM ACCESS METHODS
///START COMMON MATH FUNCTIONS
///return the complex encoding the rotation
///for a given missmatch interval
Cmplx GetRotationComplex(int interval)
{
assert((interval>=0)&&(interval<4));
switch(interval)
{
case 0:return Cmplx(1,0);
case 1:return Cmplx(0,1);
case 2:return Cmplx(-1,0);
default:return Cmplx(0,-1);
}
}
vcg::Point2i Rotate(vcg::Point2i p,int interval)
{
assert(interval>=0);
assert(interval<4);
Cmplx rot=GetRotationComplex(interval);
/*
| real -imag| |p.x|
| |*| |
| imag real | |p.y|
*/
vcg::Point2i ret;
ret.X()=rot.real()*p.X()-rot.imag()*p.Y();
ret.Y()=rot.imag()*p.X()+rot.real()*p.Y();
return (ret);
}
///END COMMON MATH FUNCTIONS
///START ENERGY MINIMIZATION PART
///initialize the LHS for a given face
///for minimization of Dirichlet's energy
void perElementLHS(FaceType *f,
ScalarType val[3][3],
int index[3][3][2])
{
///initialize to zero
for (int x=0;x<3;x++)
for (int y=0;y<3;y++)
val[x][y]=0;
///get the vertices
VertexType *v[3];
v[0]=f->V(0);
v[1]=f->V(1);
v[2]=f->V(2);
///get the indexes of vertex instance (to consider cuts)
///for the current face
int Vindexes[3];
//Vindexes[0]=f->syst_index[0];
//Vindexes[1]=f->syst_index[1];
//Vindexes[2]=f->syst_index[2];
Vindexes[0]=HandleS_Index[f][0];
Vindexes[1]=HandleS_Index[f][1];
Vindexes[2]=HandleS_Index[f][2];
//VIndex.getIndexInfo(f,Vindexes);
///initialize the indexes for the block
for (int x=0;x<3;x++)
for (int y=0;y<3;y++)
{
index[x][y][0]=Vindexes[x];
index[x][y][1]=Vindexes[y];
}
///initialize edges
CoordType e[3];
for (int k=0;k<3;k++)
e[k]=v[(k+2)%3]->P()-v[(k+1)%3]->P();
///then consider area but also considering scale factor dur to overlaps
ScalarType areaT=((f->P(1)-f->P(0))^(f->P(2)-f->P(0))).Norm()/2.0;
//ScalarType ScaleFactor=f->area;
for (int x=0;x<3;x++)
for (int y=0;y<3;y++)
if (x!=y)
{
ScalarType num=(e[x]*e[y]);
val[x][y] =num/(4.0*areaT);//*ScaleFactor);//*(ScalarType)convex_size);
val[x][y]*=Handle_Stiffness[f];//f->stiffening;
//val[x][y]*=ScaleFactor;
}
///set the matrix as diagonal
SetDiagonal(val);
}
///initialize the RHS for a given face
///for minimization of Dirichlet's energy
void perElementRHS(FaceType *f,
ScalarType b[6],
ScalarType vector_field_scale=1)
{
/// then set the rhs
CoordType scaled_Kreal;
CoordType scaled_Kimag;
CoordType fNorm=f->N();
CoordType p[3];
p[0]=f->P0(0);
p[1]=f->P0(1);
p[2]=f->P0(2);
CoordType neg_t[3];
neg_t[0] = fNorm ^ (p[2] - p[1]);
neg_t[1] = fNorm ^ (p[0] - p[2]);
neg_t[2] = fNorm ^ (p[1] - p[0]);
CoordType K1,K2;
//K1=CrossVector<FaceType>(*f,0);
//K2=CrossVector<FaceType>(*f,1);
K1=f->PD1();
K2=f->PD2();
scaled_Kreal = K1*(vector_field_scale)/2;
scaled_Kimag = K2*(vector_field_scale)/2;
ScalarType stiff_val=(ScalarType)Handle_Stiffness[f];
b[0] = scaled_Kreal * neg_t[0]*stiff_val;
b[1] = scaled_Kimag * neg_t[0]*stiff_val;
b[2] = scaled_Kreal * neg_t[1]*stiff_val;
b[3] = scaled_Kimag * neg_t[1]*stiff_val;
b[4] = scaled_Kreal * neg_t[2]*stiff_val;
b[5] = scaled_Kimag * neg_t[2]*stiff_val;
}
///evaluate the LHS and RHS for a single face
///for minimization of Dirichlet's energy
void PerElementSystemReal(FaceType *f,
ScalarType val[3][3],
int index[3][3][2],
ScalarType b[6],
ScalarType vector_field_scale=1.0)
{
perElementLHS(f,val,index);
perElementRHS(f,b,vector_field_scale);
}
///END ENERGY MINIMIZATION PART
///START FIXING VERTICES
///set a given vertex as fixed
void AddFixedVertex(VertexType *v)
{
n_fixed_vars++;
Hard_constraints.push_back(v);
//v->blocked=true;
}
///find vertex to fix in case we're using
///a vector field NB: multiple components not handled
void FindFixedVertField()
{
Hard_constraints.clear();
n_fixed_vars=0;
///fix the first singularity
for (unsigned int j=0;j<mesh.vert.size();j++)
{
VertexType *v=&mesh.vert[j];
if (v->IsD())continue;
//if (v->IsSingular())
if (Handle_Singular[v])
{
AddFixedVertex(v);
v->T().P()=vcg::Point2<ScalarType>(0,0);
return;
}
}
///if anything fixed fix the first
AddFixedVertex(&mesh.vert[0]);
mesh.vert[0].T().P()=vcg::Point2<ScalarType>(0,0);
}
///find hard constraint depending if using or not
///a vector field
void FindFixedVert()
{
Hard_constraints.clear();
FindFixedVertField();
}
int GetFirstVertexIndex(VertexType *v)
{
FaceType *f=v->VFp();
int index=v->VFi();
//return (f->syst_index[index]);
return (HandleS_Index[f][index]);
//int indexV[3];
//VIndex.getIndexInfo(f,indexV);
//return indexV[index];
}
///fix the vertices which are flagged as fixed
void FixBlockedVertex()
{
int offset_row=n_vert_vars+num_cut_constraint;
unsigned int constr_num = 0;
for (unsigned int i=0;i<Hard_constraints.size();i++)
{
VertexType *v=Hard_constraints[i];
assert(!v->IsD());
///get first index of the vertex that must blocked
//int index=v->vertex_index[0];
int index=GetFirstVertexIndex(v);
///multiply times 2 because of uv
int indexvert=index*2;
///find the first free row to add the constraint
int indexRow=(offset_row+constr_num)*2;
int indexCol=indexRow;
///add fixing constraint LHS
AddValA(indexRow,indexvert,1);
AddValA(indexRow+1,indexvert+1,1);
///add fixing constraint RHS
//AddValB(indexCol,(ScalarType)v->range_index.X());
//AddValB(indexCol+1,(ScalarType)v->range_index.Y());
AddValB(indexCol,(ScalarType)v->T().P().X());
AddValB(indexCol+1,(ScalarType)v->T().P().Y());
//AddValB(indexCol,0);
//AddValB(indexCol+1,0);
constr_num++;
}
assert(constr_num==n_fixed_vars);
}
///END FIXING VERTICES
///HANDLING SINGULARITY
//set the singularity round to integer location
void AddSingularityRound()
{
for (unsigned int j=0;j<mesh.vert.size();j++)
{
VertexType *v=&mesh.vert[j];
if (v->IsD())continue;
//if (v->IsSingular())
if (Handle_Singular[v])
{
//assert(v->vertex_index.size()==1);
//ids_to_round.push_back(v->vertex_index[0]*2);
//ids_to_round.push_back((v->vertex_index[0]*2)+1);
int index0=GetFirstVertexIndex(v);
ids_to_round.push_back(index0*2);
ids_to_round.push_back((index0*2)+1);
}
}
}
///START GENERIC SYSTEM FUNCTIONS
//build the laplacian matrix cyclyng over all rangemaps
//and over all faces
void BuildLaplacianMatrix(double vfscale=1)
{
///then for each face
for (unsigned int j=0;j<mesh.face.size();j++)
{
FaceType *f=&mesh.face[j];
if (f->IsD())
continue;
int var_idx[3];//vertex variable indices
//for(int k = 0; k < 3; ++k)
// var_idx[k] = f->syst_index[k];
for(int k = 0; k < 3; ++k)
var_idx[k] = HandleS_Index[f][k];
//VIndex.getIndexInfo(f,var_idx);
///block of variables
ScalarType val[3][3];
///block of vertex indexes
int index[3][3][2];
///righe hand side
ScalarType b[6];
///compute the system for the given face
PerElementSystemReal(f, val,index, b, vfscale);
//Add the element to the matrix
Add33Block(val,index);
///add area term.. to test if needed
//if (!use_direction_field)
// AddAreaTerm(index,1.0);//f->area);
///add right hand side
//if (use_direction_field)
AddRHS(b,var_idx);
}
}
///find different sized of the system
void FindSizes()
{
///find the vertex that need to be fixed
FindFixedVert();
///REAL PART
n_vert_vars=Handle_SystemInfo().num_vert_variables;//VIndex.NumVertexVariables();
///INTEGER PART
///the total number of integer variables
n_integer_vars=Handle_SystemInfo().num_integer_cuts;//VIndex.NumInteger();
///CONSTRAINT PART
num_cut_constraint=Handle_SystemInfo().EdgeSeamInfo.size()*2;//VIndex.NumCutsConstraints();
num_constraint_equations=num_cut_constraint+n_fixed_vars;
///total variable of the system
num_total_vars=n_vert_vars+n_integer_vars;
///initialize matrix size
//get max constraint-integer size
int MaxAddSize=std::max(n_integer_vars,num_constraint_equations);
system_size = (n_vert_vars+MaxAddSize)*2;
printf("\n*** SYSTEM VARIABLES *** \n");
printf("* NUM REAL VERTEX VARIABLES %d \n",n_vert_vars);
printf("\n*** SINGULARITY *** \n ");
printf("* NUM SINGULARITY %d\n",(int)ids_to_round.size()/2);
printf("\n*** INTEGER VARIABLES *** \n");
printf("* NUM INTEGER VARIABLES %d \n",(int)n_integer_vars);
printf("\n*** CONSTRAINTS *** \n ");
printf("* NUM FIXED CONSTRAINTS %d\n",n_fixed_vars);
printf("* NUM CUTS CONSTRAINTS %d\n",num_cut_constraint);
printf("\n*** TOTAL SIZE *** \n");
printf("* TOTAL VARIABLE SIZE (WITH INTEGER TRASL) %d \n",num_total_vars);
printf("* TOTAL CONSTRAINTS %d \n",num_constraint_equations);
printf("* MATRIX SIZE %d \n",system_size);
}
void AllocateSystem()
{
S.initialize(system_size, system_size);
printf("\n INITIALIZED SPARSE MATRIX OF %d x %d \n",system_size, system_size);
// number of nonzero entries in the matrix
int num_facet=mesh.fn;
//unsigned int nentries_vert = 4*6*6*2*n_scalar_vars; // 6x6 for each facet
unsigned int nentries_vert = 6*6*num_facet*2;
unsigned int nentries_fixed = 2*n_fixed_vars; // 1 complex variable fixed in each constraint, i.e. 2 regular
unsigned int nentries_transition = 2*6*(n_integer_vars); // 3 complex variables involved in each constraint = 8 regular
unsigned int nentries_constraints = 4*4*6*(num_constraint_equations); //
int total_reserve=nentries_vert
+ nentries_fixed
+ nentries_transition
+nentries_constraints;
printf("\n*** SPACE ALLOCATION *** \n");
printf("* number of reserved vertices variables %d \n",nentries_vert);
printf("* number of reserved entries fixed %d \n",nentries_fixed);
printf("* number of reserved entries integer %d \n",nentries_transition);
printf("* number of reserved entries constraints %d \n",nentries_constraints);
printf("* total number of reserved entries %d \n",total_reserve);
S.A().reserve(2*total_reserve);
printf("\n*** ALLOCATED *** \n");
}
///intitialize the whole matrix
void InitMatrix()
{
///find singularities that must be rounded
//AddSingularityRound();
FindSizes();
AllocateSystem();
}
///map back coordinates after that
///the system has been solved
void MapCoords()
{
///map coords to faces
for (unsigned int j=0;j<mesh.face.size();j++)
{
FaceType *f=&mesh.face[j];
if (f->IsD())continue;
//int indexV[3];
//VIndex.getIndexInfo(f,indexV);
for (int k=0;k<3;k++)
{
//get the index of the variable in the system
//int indexUV=indexV[k];//f->syst_index[k];
int indexUV=HandleS_Index[f][k];
///then get U and V coords
double U=X[indexUV*2];
double V=X[indexUV*2+1];
vcg::Point2<ScalarType> uv=vcg::Point2<ScalarType>(U,V);
///assing
//f->realUV[k]=uv;
// ScalarType factor=(ScalarType)SIZEQUADS/(ScalarType)SIZEPARA;
// uv*=factor;
///assing
f->WT(k).P()=uv;
}
}
///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<endIndex;i++)
// {
// ///assert that the value is an integer value
// ScalarType value=X[i];
// ScalarType diff=value-(int)floor(value+0.5);
// assert(diff<0.00000001);
// Handle_SystemInfo().IntegerValues[index]=value;
// index++;
// }
}
///END GENERIC SYSTEM FUNCTIONS
///set the constraints for the inter-range cuts
void BuildSeamConstraintsExplicitTranslation()
{
///add constraint(s) for every seam edge (not halfedge)
int offset_row=n_vert_vars;
///current constraint row
int constr_row=offset_row;
///current constraint
unsigned int constr_num = 0;
///get the info for seams
//std::vector<SeamInfo> EdgeSeamInfo;
//VIndex.GetSeamInfo(EdgeSeamInfo);
for (unsigned int i=0;i<Handle_SystemInfo().EdgeSeamInfo.size();i++)
{
//if (i<25)continue;
unsigned char interval=Handle_SystemInfo().EdgeSeamInfo[i].MMatch;
if (interval==1)
interval=3;
else
if(interval==3)
interval=1;
//printf("%d\n",interval);
int p0 = Handle_SystemInfo().EdgeSeamInfo[i].v0;
int p1 = Handle_SystemInfo().EdgeSeamInfo[i].v1;
int p0p = Handle_SystemInfo().EdgeSeamInfo[i].v0p;
int p1p = Handle_SystemInfo().EdgeSeamInfo[i].v1p;
/*assert(p1!=p1p);
assert(p0!=p0p);*/
Cmplx rot=GetRotationComplex(interval);
///get the integer variable
int integerVar=offset_row+Handle_SystemInfo().EdgeSeamInfo[i].integerVar;
//Cmplx rotInt=GetRotationComplex(interval);
if (integer_rounding)
{
ids_to_round.push_back(integerVar*2);
ids_to_round.push_back(integerVar*2+1);
}
AddComplexA(constr_row, p0 , rot);
AddComplexA(constr_row, p0p, -1);
///then translation...considering the rotation
///due to substitution
AddComplexA(constr_row, integerVar, 1);
AddValB(2*constr_row,0);
AddValB(2*constr_row+1,0);
constr_row +=1;
constr_num++;
AddComplexA(constr_row, p1, rot);
AddComplexA(constr_row, p1p, -1);
///other translation
AddComplexA(constr_row, integerVar , 1);
AddValB(2*constr_row,0);
AddValB(2*constr_row+1,0);
constr_row +=1;
constr_num++;
//// r p1 - p1' + t = 0
}
//assert(constr_num==num_cut_constraint);
}
///call of the mixed integer solver
void MixedIntegerSolve(double cone_grid_res=1,
bool direct_round=true,
int localIter=0)
{
X=std::vector< double >((n_vert_vars+n_integer_vars)*2);
///variables part
int ScalarSize=(n_vert_vars)*2;
int SizeMatrix=(n_vert_vars+n_integer_vars)*2;
printf("\n ALLOCATED X \n");
///matrix A
gmm::col_matrix< gmm::wsvector< double > > A(SizeMatrix,SizeMatrix); // lhs matrix variables +
///constraints part
int CsizeX=num_constraint_equations*2;
int CsizeY=SizeMatrix+1;
gmm::row_matrix< gmm::wsvector< double > > C(CsizeX,CsizeY); // constraints
printf("\n ALLOCATED QMM STRUCTURES \n");
std::vector< double > rhs(SizeMatrix,0); // rhs
printf("\n ALLOCATED RHS STRUCTURES \n");
//// copy LHS
for(int i = 0; i < (int)S.A().nentries(); ++i)
{
int row = S.A().rowind()[i];
int col = S.A().colind()[i];
int size=(int)S.nrows();
assert(0 <= row && row < size);
assert(0 <= col && col < size);
// it's either part of the matrix
if (row < ScalarSize) {
A(row, col) += S.A().vals()[i];
}
// or it's a part of the constraint
else
{
//if (row<(n_vert_vars+num_constraint_equations)*2)
assert ((unsigned int)row<(n_vert_vars+num_constraint_equations)*2);
//{
int r = row - ScalarSize;
assert(r<CsizeX);
assert(col<CsizeY);
C(r , col ) += S.A().vals()[i];
//}
}
}
printf("\n SET %d INTEGER VALUES \n",n_integer_vars);
///add penalization term for integer variables
double penalization=0.000001;
int offline_index=ScalarSize;
for(unsigned int i = 0; i < (n_integer_vars)*2; ++i)
{
int index=offline_index+i;
A(index,index)=penalization;
}
printf("\n SET RHS \n");
// copy RHS
for(int i = 0; i < (int)ScalarSize; ++i)
{
rhs[i] = S.getRHSReal(i) * cone_grid_res;
}
// copy constraint RHS
printf("\n SET %d CONSTRAINTS \n",num_constraint_equations);
for(unsigned int i = 0; i < num_constraint_equations; ++i)
{
C(i, SizeMatrix) = -S.getRHSReal(ScalarSize + i) * cone_grid_res;
}
///copy values back into S
COMISO::ConstrainedSolver solver;
solver.misolver().set_local_iters(localIter);
solver.misolver().set_direct_rounding(direct_round);
std::sort(ids_to_round.begin(),ids_to_round.end());
std::vector<int>::iterator new_end=std::unique(ids_to_round.begin(),ids_to_round.end());
int dist=distance(ids_to_round.begin(),new_end);
ids_to_round.resize(dist);
solver.solve( C, A, X, rhs, ids_to_round, 0.0, true, true);
}
void GetAttributes()
{
// you can query if an attribute is present or not
bool hasSystIndex = vcg::tri::HasPerFaceAttribute(mesh,std::string("SystemIndex"));
bool hasSingular = vcg::tri::HasPerVertexAttribute(mesh,std::string("Singular"));
bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
bool HasSystemInfo=vcg::tri::HasPerMeshAttribute(mesh,std::string("SystemInfo"));
assert(hasSystIndex);
assert(hasSingular);
assert(hasStiffness);
assert(HasSystemInfo);
HandleS_Index = vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3i>(mesh,"SystemIndex");
Handle_Singular=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<bool>(mesh,std::string("Singular"));
Handle_Stiffness=vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<float>(mesh,std::string("Stiffness"));
Handle_SystemInfo=vcg::tri::Allocator<MeshType>::template GetPerMeshAttribute<MeshSystemInfo>(mesh,"SystemInfo");
}
public:
void SolvePoisson(ScalarType vector_field_scale=0.1f,
ScalarType grid_res=1.f,
bool direct_round=true,
int localIter=0,
bool _integer_rounding=true)
{
GetAttributes();
//initialization of flags and data structures
//ClearFlags();
integer_rounding=_integer_rounding;
ids_to_round.clear();
///Initializing Matrix
int t0=clock();
///initialize the matrix ALLOCATING SPACE
InitMatrix();
printf("\n ALLOCATED THE MATRIX \n");
///build the laplacian system
BuildLaplacianMatrix(vector_field_scale);
BuildSeamConstraintsExplicitTranslation();
////add the lagrange multiplier
FixBlockedVertex();
printf("\n BUILT THE MATRIX \n");
if (integer_rounding)
AddSingularityRound();
int t1=clock();
printf("\n time:%d \n",t1-t0);
printf("\n SOLVING \n");
MixedIntegerSolve(grid_res,direct_round,localIter);
int t2=clock();
printf("\n time:%d \n",t2-t1);
printf("\n ASSIGNING COORDS \n");
MapCoords();
int t3=clock();
printf("\n time:%d \n",t3-t2);
printf("\n FINISHED \n");
//TagFoldedFaces();
}
PoissonSolver(MeshType &_mesh):mesh(_mesh)
{}
};
};
#endif

367
wrap/miq/core/seams_initializer.h
View File

@ -1,367 +0,0 @@
#ifndef MIQ_SEAMS_INTIALIZER
#define MIQ_SEAMS_INTIALIZER
#include <vcg/complex/complex.h>
#include <vcg/simplex/face/pos.h>
#include <vcg/simplex/face/jumping_pos.h>
#include <wrap/io_trimesh/import_field.h>
template <class MeshType>
class SeamsInitializer
{
private:
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::FaceIterator FaceIterator;
MeshType *mesh;
///per face per edge of mmatch in the solver
typename MeshType::template PerFaceAttributeHandle<vcg::Point3i> Handle_MMatch;
///per vertex singular or not
typename MeshType::template PerVertexAttributeHandle<bool> Handle_Singular;
///per vertex degree of a singularity
typename MeshType::template PerVertexAttributeHandle<int> Handle_SingularDegree;
///seam per face
typename MeshType::template PerFaceAttributeHandle<vcg::Point3<bool> > Handle_Seams;
///seam index per face
typename MeshType::template PerFaceAttributeHandle<vcg::Point3i > 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<FaceType> pos(v.cVFp(), v.cVFi());
std::vector<vcg::face::Pos<FaceType> > posVec;
vcg::face::VFOrderedStarFF(pos, posVec);
missmatch=0;
for (unsigned int i=0;i<posVec.size();i++)
{
FaceType *curr_f=posVec[i].F();
int currMM=Handle_MMatch[curr_f][posVec[i].E()];
missmatch+=currMM;
}
missmatch=missmatch%4;
return(missmatch!=0);
}
///initialized mapping structures if are not already initialized
void AddAttributesIfNeeded()
{
Handle_MMatch = vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3i>(*mesh,std::string("MissMatch"));
Handle_Singular=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<bool>(*mesh,std::string("Singular"));
Handle_SingularDegree=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<int>(*mesh,std::string("SingularityDegree"));
Handle_Seams=vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3<bool> >(*mesh,std::string("Seams"));
Handle_SeamsIndex=vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3i >(*mesh,std::string("SeamsIndex"));
}
void FloodFill(FaceType* start)
{
std::deque<FaceType*> d;
///clean the visited flag
start->SetV();
d.push_back(start);
while (!d.empty()){
FaceType *f = d.at(0); d.pop_front();
for (int s = 0; s<3; s++)
{
FaceType *g = f->FFp(s);
int j = f->FFi(s);
if ((!(IsRotSeam(f,s))) && (!(IsRotSeam(g,j))) && (!g->IsV()) )
{
Handle_Seams[f][s]=false;
Handle_Seams[g][j]=false;
g->SetV();
d.push_back(g);
}
}
}
}
void Retract(){
std::vector<int> e(mesh->vert.size(),0); // number of edges per vert
VertexType *vb = &(mesh->vert[0]);
for (FaceIterator f = mesh->face.begin(); f!=mesh->face.end(); f++) if (!f->IsD()){
for (int s = 0; s<3; s++){
//if (f->seam[s])
if (Handle_Seams[f][s])
if (f->FFp(s)<=&*f) {
e[ f->V(s) - vb ] ++;
e[ f->V1(s) - vb ] ++;
}
}
}
bool over=true;
int guard = 0;
do {
over = true;
for (FaceIterator f = mesh->face.begin(); f!=mesh->face.end(); f++) if (!f->IsD()){
for (int s = 0; s<3; s++){
//if (f->seam[s])
if (Handle_Seams[f][s])
if (!(IsRotSeam(&(*f),s))) // never retract rot seams
//if (f->FFp(s)<=&*f)
{
if (e[ f->V(s) - vb ] == 1) {
// dissolve seam
//f->seam[s] = false;
Handle_Seams[f][s]=false;
//f->FFp(s)->seam[(int)f->FFi(s)] = false;
Handle_Seams[f->FFp(s)][(int)f->FFi(s)]=false;
e[ f->V(s) - vb ] --;
e[ f->V1(s) - vb ] --;
over = false;
}
}
}
}
if (guard++>10000) over = true;
} while (!over);
}
void AddSeamsByMM()
{
for (unsigned int i=0;i<mesh->face.size();i++)
{
FaceType *f=&mesh->face[i];
if (f->IsD())continue;
for (int j=0;j<3;j++)
{
if (IsRotSeam(f,j))
Handle_Seams[f][j]=true;
//f->SetSeam(j);
}
}
}
void SelectSingularityByMM()
{
for (unsigned int i=0;i<mesh->vert.size();i++)
{
if (mesh->vert[i].IsD())continue;
int missmatch;
bool isSing=IsSingularByMMatch(mesh->vert[i],missmatch);
if (isSing)
{
mesh->vert[i].SetS();
Handle_Singular[i]=true;
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<MeshType>::FaceClearV(*mesh);
for (FaceIterator f = mesh->face.begin(); f!=mesh->face.end(); f++)
if (!f->IsD())
{
Handle_Seams[f][0]=true;
Handle_Seams[f][1]=true;
Handle_Seams[f][2]=true;
}
int index=0;
for (FaceIterator f = mesh->face.begin(); f!=mesh->face.end(); f++)
if (!f->IsD())
{
if (!f->IsV())
{
index++;
FloodFill(&*f);
}
}
Retract();
return index;
}
void InitMMatch()
{
for (unsigned int i=0;i<mesh->face.size();i++)
{
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<MeshType>::MissMatchByCross(*curr_f,*opp_f);
}
}
}
void InitSeamIndexes()
{
///initialize seams indexes
for (unsigned int i=0;i<mesh->face.size();i++)
{
FaceType *f=&mesh->face[i];
for (int j=0;j<3;j++)
Handle_SeamsIndex[f][j]=-1;
}
std::vector<std::vector<std::pair<FaceType*,int> > > seamsVert;
seamsVert.resize(mesh->vert.size());
for (unsigned int i=0;i<mesh->face.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<FaceType*,int> entry=std::pair<FaceType*,int>(f,j);
seamsVert[index0].push_back(entry);
seamsVert[index1].push_back(entry);
}
}
}
int curr_index=0;
for (unsigned int i=0;i<mesh->vert.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;j<seamsVert[i].size();j++)
{
FaceType *f=seamsVert[i][j].first;
int edge=seamsVert[i][j].second;
VertexType *v0=v_seam;
VertexType *v1=NULL;
///the index is already initialized
if (Handle_SeamsIndex[f][edge]!=-1)continue;
bool finished=false;
do
{
///otherwise follow the index
Handle_SeamsIndex[f][edge]=curr_index;///set current value
FaceType *f1=f->FFp(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<MeshType>::FaceFace(*_mesh);
vcg::tri::UpdateFlags<MeshType>::FaceBorderFromFF(*_mesh);
vcg::tri::UpdateFlags<MeshType>::VertexBorderFromFace(*_mesh);
AddAttributesIfNeeded();
if (orient_globally)
vcg::tri::CrossField<MeshType>::MakeDirectionFaceCoherent(*mesh);
if (initMM)
InitMMatch();
SelectSingularityByMM();
if (initCuts)
{
InitTopologycalCuts();
AddSeamsByMM();
}
//InitSeamIndexes();
}
SeamsInitializer(){mesh=NULL;}
};
#endif

256
wrap/miq/core/sparsesystemdata.h
View File

@ -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

170
wrap/miq/core/stiffening.h
View File

@ -1,170 +0,0 @@
#ifndef MIQ_STIFFENING
#define MIQ_STIFFENING
template <class ScalarType>
ScalarType Gauss(ScalarType &value)
{
const ScalarType E_NEPER=2.71828;
const ScalarType den=sqrt(2.0*M_PI);
ScalarType exponent=-0.5*pow(value,2);
ScalarType res=((1.0/den)*pow(E_NEPER,exponent));
return res;
}
template <class MeshType>
class 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<float> Handle_Stiffness=vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<float>(mesh,std::string("Stiffness"));
for (unsigned int i=0;i<mesh.face.size();i++)
{
//MeshType::ScalarType val=MaxVal-mesh.face[i].stiffening+1;
ScalarType val=MaxVal-Handle_Stiffness[i]+1;
if (val<1)val=1;
mesh.face[i].C()=vcg::Color4b::ColorRamp(1.0,MaxVal,val);
}
}
static void AddGaussStiffening(MeshType & mesh,ScalarType C)
{
int radius=floor(C);
if (C<4)radius=4;
typename MeshType::template PerFaceAttributeHandle<float> Handle_Stiffness;
bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
if(!hasStiffness)
Handle_Stiffness=vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));
else
Handle_Stiffness=vcg::tri::Allocator<MeshType>::template FindPerFaceAttribute<float>(mesh,std::string("Stiffness"));
bool hasSingular = vcg::tri::HasPerVertexAttribute(mesh,std::string("Singular"));
assert(hasSingular);
typename MeshType::template PerVertexAttributeHandle<bool> Handle_Singular;
Handle_Singular=vcg::tri::Allocator<MeshType>:: template GetPerVertexAttribute<bool>(mesh,std::string("Singular"));
std::vector<VertexType*> to_stiff;
for(unsigned int i=0;i<mesh.vert.size();i++)
{
VertexType *v=&mesh.vert[i];
if (v->IsD())continue;
//if (!v->IsSingular())continue;
if (!Handle_Singular[v])continue;
to_stiff.push_back(v);
}
for(unsigned int i=0;i<mesh.face.size();i++)
{
FaceType *f=&(mesh.face[i]);
if (f->IsD())continue;
if (!f->IsV())continue;
to_stiff.push_back(f->V(0));
to_stiff.push_back(f->V(1));
to_stiff.push_back(f->V(2));
}
std::sort(to_stiff.begin(),to_stiff.end());
typename std::vector<VertexType*>::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<to_stiff.size();i++)
{
VertexType *v=to_stiff[i];
for (int r=0;r<radius;r++)
{
ScalarType stiffVal=((ScalarType)r)/(ScalarType)radius;//((ScalarType)(radius-r))/(ScalarType)radius;
stiffVal*=3.0;
stiffVal=Gauss(stiffVal)/0.4;
stiffVal=1+(stiffVal*C);
std::vector<FaceType*> ring;
//mesh.GetConnectedFaces(v,r,ring);
VFExtendedStarVF(v,r,ring);
///then set stiffening
for (unsigned int k=0;k<ring.size();k++)
{
FaceType* f=ring[k];
//if (f->stiffening<stiffVal)
// f->stiffening=stiffVal;
if (Handle_Stiffness[f]<stiffVal)
Handle_Stiffness[f]=stiffVal;
}
}
}
}
static bool updateStiffeningJacobianDistorsion(MeshType & mesh,ScalarType grad_size)
{
typename MeshType::template PerFaceAttributeHandle<float> Handle_Stiffness;
bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
if(!hasStiffness)
Handle_Stiffness=vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));
else
Handle_Stiffness=vcg::tri::Allocator<MeshType>::template FindPerFaceAttribute<float>(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<FaceType>(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<float> Handle_Stiffness;
bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
if(!hasStiffness)
Handle_Stiffness=vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));
else
Handle_Stiffness=vcg::tri::Allocator<MeshType>::template FindPerFaceAttribute<float>(mesh,std::string("Stiffness"));
for(unsigned int i=0;i<mesh.face.size();i++)
{
FaceType *f=&(mesh.face[i]);
//f->stiffening=1;
Handle_Stiffness[f]=1;
}
}
};
#endif

443
wrap/miq/core/vertex_indexing.h
View File

@ -1,443 +0,0 @@
#ifndef MIQ_VERTEX_INDEXING
#define MIQ_VERTEX_INDEXING
#include <vcg/complex/complex.h>
#include <vcg/simplex/face/pos.h>
#include <vcg/simplex/face/jumping_pos.h>
#include <vcg/simplex/face/topology.h>
struct SeamInfo
{
int v0,v0p,v1,v1p;
int 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<SeamInfo> EdgeSeamInfo;
///this are values of integer variables after optimization
std::vector<int> IntegerValues;
};
template <class MeshType>
class VertexIndexing
{
private:
typedef typename MeshType::ScalarType ScalarType;
typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::FaceIterator FaceIterator;
MeshType *mesh;
///this maps back index to vertices
std::vector<VertexType*> IndexToVert;
///this maps the integer for edge
typename MeshType::template PerFaceAttributeHandle<vcg::Point3i> Handle_Integer;
///per face indexes of vertex in the solver
typename MeshType::template PerFaceAttributeHandle<vcg::Point3i> HandleS_Index;
///per face per edge of mmatch in the solver
typename MeshType::template PerFaceAttributeHandle<vcg::Point3i> Handle_MMatch;
///per vertex variable indexes
typename MeshType::template PerVertexAttributeHandle<std::vector<int> > HandleV_Integer;
///per vertex singular or not
typename MeshType::template PerVertexAttributeHandle<bool> Handle_Singular;
///per vertex degree of a singularity
typename MeshType::template PerVertexAttributeHandle<int> Handle_SingularDegree;
///seam per face
typename MeshType::template PerFaceAttributeHandle<vcg::Point3<bool> > Handle_Seams;
///this handle for mesh
typename MeshType::template PerMeshAttributeHandle<MeshSystemInfo> Handle_SystemInfo;
///this are used for drawing purposes
std::vector<VertexType*> duplicated;
void FirstPos(const VertexType* v,FaceType *&f,int &edge)
{
f=v->cVFp();
edge=v->cVFi();
}
int AddNewIndex(VertexType* v0)
{
Handle_SystemInfo().num_scalar_variables++;
HandleV_Integer[v0].push_back(Handle_SystemInfo().num_scalar_variables);
IndexToVert.push_back(v0);
return Handle_SystemInfo().num_scalar_variables;
}
bool HasIndex(int indexVert,int indexVar)
{
for (unsigned int i=0;i<HandleV_Integer[indexVert].size();i++)
if (HandleV_Integer[indexVert][i]==indexVar)return true;
return false;
}
bool HasIndex(VertexType* v0,int indexVar)
{
for (unsigned int i=0;i<HandleV_Integer[v0].size();i++)
if (HandleV_Integer[v0][i]==indexVar)return true;
return false;
}
void GetSeamInfo(const FaceType *f0,
const FaceType *f1,
const int indexE,
int &v0,int &v1,
int &v0p,int &v1p,
unsigned char &_MMatch,
int &integerVar)
{
int edgef0=indexE;
v0=HandleS_Index[f0][edgef0];
v1=HandleS_Index[f0][(edgef0+1)%3];
////get the index on opposite side
assert(f0->cFFp(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<FaceType> VFI(f_init,edge_init);
bool vertexB=vert->IsB();
bool possible_split=false;
bool complete_turn=false;
do
{
FaceType *curr_f=VFI.F();
int curr_edge=VFI.E();
VFI.NextFE();
FaceType *next_f=VFI.F();
///test if I've just crossed a border
bool on_border=(curr_f->FFp(curr_edge)==curr_f);
//bool mismatch=false;
bool seam=false;
///or if I've just crossed a seam
///if I'm on a border I MUST start from the one next t othe border
if (!vertexB)
//seam=curr_f->IsSeam(next_f);
seam=IsSeam(curr_f,next_f);
if (vertexB)
assert(!Handle_Singular[vert]);
;
//assert(!vert->IsSingular());
possible_split=((on_border)||(seam));
complete_turn=(next_f==f_init);
}while ((!possible_split)&&(!complete_turn));
face=VFI.F();
edge=VFI.E();
///test that is not on a border
//assert(face->FFp(edge)!=face);
}
///intialize the mapping given an initial pos
///whih must be initialized with FindInitialPos
void MapIndexes(VertexType * vert,
int &edge_init,
FaceType *&f_init)
{
///check that is not on border..
///in such case maybe it's non manyfold
///insert an initial index
int curr_index=AddNewIndex(vert);
///and initialize the jumping pos
vcg::face::JumpingPos<FaceType> VFI(f_init,edge_init);
bool complete_turn=false;
do
{
FaceType *curr_f=VFI.F();
int curr_edge=VFI.E();
///assing the current index
HandleS_Index[curr_f][curr_edge]=curr_index;
VFI.NextFE();
FaceType *next_f=VFI.F();
///test if I've finiseh with the face exploration
complete_turn=(next_f==f_init);
///or if I've just crossed a mismatch
if (!complete_turn)
{
bool seam=false;
//seam=curr_f->IsSeam(next_f);
seam=IsSeam(curr_f,next_f);
if (seam)
{
///then add a new index
curr_index=AddNewIndex(vert);
}
}
}while (!complete_turn);
}
///intialize the mapping for a given vertex
void InitMappingSeam(VertexType *vert)
{
///first rotate until find the first pos after a mismatch
///or a border or return to the first position...
FaceType *f_init=vert->VFp();
int indexE=vert->VFi();
vcg::face::JumpingPos<FaceType> VFI(f_init,indexE);
int edge_init;
FaceType *face_init;
FindInitialPos(vert,edge_init,face_init);
MapIndexes(vert,edge_init,face_init);
}
///intialize the mapping for a given sampled mesh
void InitMappingSeam()
{
//num_scalar_variables=-1;
Handle_SystemInfo().num_scalar_variables=-1;
for (unsigned int i=0;i<mesh->vert.size();i++)
if (!mesh->vert[i].IsD())
InitMappingSeam(&mesh->vert[i]);
for (unsigned int j=0;j<mesh->vert.size();j++)
{
VertexType *v= &(mesh->vert[j]);
if (!v->IsD()){
assert(HandleV_Integer[j].size()>0);
if (HandleV_Integer[j].size()>1)
duplicated.push_back(v);
}
}
}
///initialized mapping structures if are not already initialized
void AddAttributesIfNeeded()
{
IndexToVert.clear();
///other per mesh attributes
///see if already exists otherwise it creates it
Handle_SystemInfo=vcg::tri::Allocator<MeshType>::template GetPerMeshAttribute<MeshSystemInfo>(*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<MeshType>::template GetPerFaceAttribute<vcg::Point3i>(*mesh,std::string("SystemIndex"));
Handle_Integer= vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<vcg::Point3i>(*mesh,std::string("Integer"));
HandleV_Integer=vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<std::vector<int> >(*mesh,std::string("VertInteger"));
bool HasHandleMMatch=vcg::tri::HasPerFaceAttribute(*mesh,std::string("MissMatch"));
assert(HasHandleMMatch);
Handle_MMatch = vcg::tri::Allocator<MeshType>::template FindPerFaceAttribute<vcg::Point3i>(*mesh,std::string("MissMatch"));
bool HasHandleSingularity=vcg::tri::HasPerVertexAttribute(*mesh,std::string("Singular"));
assert(HasHandleSingularity);
Handle_Singular=vcg::tri::Allocator<MeshType>::template FindPerVertexAttribute<bool>(*mesh,std::string("Singular"));
bool HasHandleSingularityDegree=vcg::tri::HasPerVertexAttribute(*mesh,std::string("SingularityDegree"));
assert(HasHandleSingularityDegree);
Handle_SingularDegree=vcg::tri::Allocator<MeshType>::template FindPerVertexAttribute<int>(*mesh,std::string("SingularityDegree"));
bool HasHandleSeams=vcg::tri::HasPerFaceAttribute(*mesh,std::string("Seams"));
assert(HasHandleSeams);
Handle_Seams=vcg::tri::Allocator<MeshType>::template FindPerFaceAttribute<vcg::Point3<bool> >(*mesh,std::string("Seams"));
}
///test consistency of face variables per vert mapping
void TestSeamMapping(FaceType *f)
{
for (int k=0;k<3;k++)
{
int indexV=HandleS_Index[f][k];
VertexType *v=f->V(k);
bool has_index=HasIndex(v,indexV);
assert(has_index);
}
}
///test consistency of face variables per vert mapping
void TestSeamMapping(int indexVert)
{
VertexType *v=&mesh->vert[indexVert];
for (unsigned int k=0;k<HandleV_Integer[indexVert].size();k++)
{
int indexV=HandleV_Integer[indexVert][k];
///get faces sharing vertex
std::vector<FaceType*> faces;
std::vector<int> indexes;
vcg::face::VFStarVF(v,faces,indexes);
for (unsigned int j=0;j<faces.size();j++)
{
FaceType *f=faces[j];
int index=indexes[j];
assert(!f->IsD());
assert(f->V(index)==v);
assert((index>=0)&&(index<3));
if (HandleS_Index[f][index]==indexV)
return;
}
}
assert(0);
}
///check consistency of variable mapping across seams
void TestSeamMapping()
{
printf("\n TESTING SEAM INDEXES \n");
///test F-V mapping
for (unsigned int j=0;j<mesh->face.size();j++)
{
FaceType *f=&mesh->face[j];
if (f->IsD()) continue;
TestSeamMapping(f);
}
///TEST V-F MAPPING
for (unsigned int j=0;j<mesh->vert.size();j++)
{
VertexType *v=&mesh->vert[j];
if (v->IsD()) continue;
TestSeamMapping(j);
}
}
public:
///vertex to variable mapping
void InitMapping()
{
//use_direction_field=_use_direction_field;
IndexToVert.clear();
duplicated.clear();
InitMappingSeam();
Handle_SystemInfo().num_vert_variables=Handle_SystemInfo().num_scalar_variables+1;
///end testing...
TestSeamMapping();
}
void InitFaceIntegerVal()
{
Handle_SystemInfo().num_integer_cuts=0;
for (unsigned int j=0;j<mesh->face.size();j++)
{
FaceType *f=&mesh->face[j];
if (f->IsD())continue;
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;j<mesh->face.size();j++)
{
FaceType *f0=&mesh->face[j];
if (f0->IsD())continue;
for (int k=0;k<3;k++)
{
FaceType *f1=f0->FFp(k);
if (f1==f0)continue;
bool seam=Handle_Seams[f0][k];//f0->IsSeam(k);
if ((seam) )//&&(f0<f1))
{
int v0,v0p,v1,v1p;
unsigned char MM;
int integerVar;
GetSeamInfo(f0,f1,k,v0,v1,v0p,v1p,MM,integerVar);
Handle_SystemInfo().EdgeSeamInfo.push_back(SeamInfo(v0,v1,v0p,v1p,MM,integerVar));
}
}
}
}
void Init(MeshType *_mesh){mesh=_mesh;AddAttributesIfNeeded();}
VertexIndexing(){mesh=NULL;}
};
#endif

740
wrap/miq/quadrangulator.h
View File

@ -1,740 +0,0 @@
#ifndef MIQ_QUADRANGULATOR_H
#define MIQ_QUADRANGULATOR_H
#include <vcg/complex/complex.h>
#include <vcg/simplex/face/pos.h>
#include <vcg/simplex/face/jumping_pos.h>
#include <vcg/complex/algorithms/attribute_seam.h>
#include <vcg/complex/algorithms/refine.h>
#include <vcg/complex/algorithms/smooth.h>
#include <vcg/complex/algorithms/clean.h>
#include <vcg/complex/algorithms/update/bounding.h>
#include <wrap/io_trimesh/export.h>
#include <vcg/complex/algorithms/update/texture.h>
#define precisionQ 0.0000000001
namespace vcg {
namespace tri {
template <class TriMesh,class PolyMesh>
class Quadrangulator
{
public:
typedef typename TriMesh::FaceType TriFaceType;
typedef typename TriMesh::VertexType TriVertexType;
typedef typename TriMesh::CoordType CoordType;
typedef typename TriMesh::ScalarType ScalarType;
typedef typename PolyMesh::FaceType PolyFaceType;
typedef typename PolyMesh::VertexType PolyVertexType;
typedef typename PolyMesh::CoordType PolyCoordType;
typedef typename PolyMesh::ScalarType PolyScalarType;
struct InterpolationInfo
{
CoordType Pos3D;
vcg::Point2<ScalarType> PosUV;
ScalarType alpha;
bool to_split;
InterpolationInfo()
{
Pos3D=CoordType(0,0,0);
PosUV=vcg::Point2<ScalarType>(0,0);
to_split=false;
alpha=-1;
}
};
//the interpolation map that is saved once to be univoque per edge
typedef std::pair<CoordType,CoordType > KeyEdgeType;
std::map<KeyEdgeType,InterpolationInfo> InterpMap;
//ScalarType UVtolerance;
private:
bool ToSplit(const vcg::Point2<ScalarType> &uv0,
const vcg::Point2<ScalarType> &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<to_split.face.size();i++)
{
TriFaceType *f=&to_split.face[i];
for (int j =0;j<3;j++)
{
vcg::Point2<ScalarType> 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 (diff0<minTolerance)
UV.X()=(ScalarType)int0;
if (diff1<minTolerance)
UV.Y()=(ScalarType)int1;
f->WT(j).P()=UV;
}
}
}
void RoundSplits(TriMesh &to_split,int dir)
{
ScalarType minTolerance=precisionQ;
//first add all eddge
for (size_t i=0;i<to_split.face.size();i++)
{
TriFaceType *f=&to_split.face[i];
for (int j =0;j<3;j++)
{
CoordType p0=f->P0(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<ScalarType> 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 (diff0<minTolerance)
UV.X()=(ScalarType)int0;
if (diff1<minTolerance)
UV.Y()=(ScalarType)int1;
InterpMap[k].PosUV=UV;
}
}
}
void InitSplitMap(TriMesh &to_split,
int dir)
{
assert((dir==0)||(dir==1));
InterpMap.clear();
//printf("direction %d\n",dir );
//first add all eddge
for (int i=0;i<to_split.face.size();i++)
{
TriFaceType *f=&to_split.face[i];
for (int j =0;j<3;j++)
{
CoordType p0=f->P0(j);
CoordType p1=f->P1(j);
vcg::Point2<ScalarType> Uv0=f->V0(j)->T().P();
vcg::Point2<ScalarType> 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;i<to_split.face.size();i++)
{
TriFaceType *f=&to_split.face[i];
for (int j =0;j<3;j++)
{
CoordType p0=f->P0(j);
CoordType p1=f->P1(j);
vcg::Point2<ScalarType> uv0=f->V0(j)->T().P();
vcg::Point2<ScalarType> 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;i<to_split.face.size();i++)
{
TriFaceType *f=&to_split.face[i];
for (int j =0;j<3;j++)
{
CoordType p0=f->P0(j);
CoordType p1=f->P1(j);
vcg::Point2<ScalarType> uv0=f->V0(j)->T().P();
vcg::Point2<ScalarType> 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<class MESH_TYPE>
struct SplitMidPoint : public std::unary_function<vcg::face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType >
{
typedef typename MESH_TYPE::VertexType VertexType;
typedef typename MESH_TYPE::FaceType FaceType;
typedef typename MESH_TYPE::CoordType CoordType;
std::map<KeyEdgeType,InterpolationInfo> *MapEdge;
void operator()(typename MESH_TYPE::VertexType &nv,
vcg::face::Pos<typename MESH_TYPE::FaceType> 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<ScalarType> WedgeInterp(vcg::TexCoord2<ScalarType> &t0,
vcg::TexCoord2<ScalarType> &t1)
{
return (vcg::TexCoord2<ScalarType>(0,0));
}
SplitMidPoint(std::map<KeyEdgeType,InterpolationInfo> *_MapEdge){MapEdge=_MapEdge;}
};
template <class MESH_TYPE>
class EdgePredicate
{
typedef typename MESH_TYPE::VertexType VertexType;
typedef typename MESH_TYPE::FaceType FaceType;
typedef typename MESH_TYPE::ScalarType ScalarType;
std::map<KeyEdgeType,InterpolationInfo> *MapEdge;
public:
bool operator()(vcg::face::Pos<typename MESH_TYPE::FaceType> 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<KeyEdgeType,InterpolationInfo> *_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<TriMesh> splMd(&InterpMap);
EdgePredicate<TriMesh> eP(&InterpMap);
done=vcg::tri::RefineE<TriMesh,SplitMidPoint<TriMesh>,EdgePredicate<TriMesh> >(to_split,splMd,eP);
}
printf("Number of Vertices %d \n",to_split.vn);
fflush(stdout);
fflush(stdout);
}
bool IsOnIntegerLine(vcg::Point2<ScalarType> uv0,
vcg::Point2<ScalarType> 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<ScalarType> 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))<UVtolerance)onIntegerL++;
if (val0==(ScalarType)floor(val0))onIntegerL++;
}
if ((IsB)&&(onIntegerL>0))return true;
return (onIntegerL==2);
}
void InitIntegerEdgesVert(TriMesh &Tmesh)
{
//IntegerEdges.clear();
vcg::tri::UpdateFlags<TriMesh>::FaceSetF(Tmesh);
vcg::tri::UpdateFlags<TriMesh>::FaceClearS(Tmesh);
vcg::tri::UpdateFlags<TriMesh>::VertexClearS(Tmesh);
for (unsigned int i=0;i<Tmesh.face.size();i++)
{
TriFaceType *f=&Tmesh.face[i];
if (f->IsD())continue;
for (int j=0;j<3;j++)
{
bool IsBorder=f->IsB(j);
if (IsBorder)
f->ClearF(j);
else
{
vcg::Point2<ScalarType> uv0=f->WT(j).P();
vcg::Point2<ScalarType> 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<ScalarType> uv0=f->WT(edge_index).P();
vcg::Point2<ScalarType> 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<TriFaceType> &currPos,
std::vector<TriVertexType *> &poly,
std::vector<short int> &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<short int> 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<size;i++)
UVpoly.push_back(TempUVpoly[(i+1)%size]);
}
void FindPolygons(TriMesh &Tmesh,
std::vector<std::vector<TriVertexType *> > &polygons,
std::vector<std::vector<short int> > &UV)
{
vcg::tri::UpdateFlags<TriMesh>::FaceClearV(Tmesh);
for (unsigned int i=0;i<Tmesh.face.size();i++)
{
TriFaceType * f=&Tmesh.face[i];
if (f->IsV())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<TriFaceType> startPos(f,j);
std::vector<TriVertexType *> 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<TriMesh>::FaceFace(Tmesh);
// vcg::tri::UpdateFlags<TriMesh>::FaceBorderFromFF(Tmesh);
}
void ConvertVTtoWT(TriMesh &Tmesh)
{
vcg::tri::UpdateTexture<TriMesh>::WedgeTexFromVertexTex(Tmesh);
vcg::tri::Clean<TriMesh>::RemoveDuplicateVertex(Tmesh);
}
void ReupdateMesh(TriMesh &Tmesh)
{
vcg::tri::UpdateNormal<TriMesh>::PerFaceNormalized(Tmesh); // update Normals
vcg::tri::UpdateNormal<TriMesh>::PerVertexNormalized(Tmesh);// update Normals
//compact the mesh
vcg::tri::Allocator<TriMesh>::CompactVertexVector(Tmesh);
vcg::tri::Allocator<TriMesh>::CompactFaceVector(Tmesh);
vcg::tri::UpdateTopology<TriMesh>::FaceFace(Tmesh); // update Topology
vcg::tri::UpdateTopology<TriMesh>::TestFaceFace(Tmesh); //and test it
//set flags
vcg::tri::UpdateFlags<TriMesh>::VertexClearV(Tmesh);
vcg::tri::UpdateFlags<TriMesh>::FaceBorderFromFF(Tmesh);
vcg::tri::UpdateFlags<TriMesh>::VertexBorderFromFaceBorder(Tmesh);
}
public:
void TestIsProper(TriMesh &Tmesh)
{
//test manifoldness
int test=vcg::tri::Clean<TriMesh>::CountNonManifoldVertexFF(Tmesh);
//assert(test==0);
if (test != 0)
cerr << "Assertion failed: TestIsProper NonManifoldVertices!" << endl;
test=vcg::tri::Clean<TriMesh>::CountNonManifoldEdgeFF(Tmesh);
//assert(test==0);
if (test != 0)
cerr << "Assertion failed: TestIsProper NonManifoldEdges" << endl;
for (unsigned int i=0;i<Tmesh.face.size();i++)
{
TriFaceType *f=&Tmesh.face[i];
assert (!f->IsD());
for (int z=0;z<3;z++)
{
//int indexOpp=f->FFi(z);
TriFaceType *Fopp=f->FFp(z);
if (Fopp==f) continue;
//assert( f->FFp(z)->FFp(f->FFi(z))==f );
if (f->FFp(z)->FFp(f->FFi(z))!=f)
cerr << "Assertion failed: TestIsProper f->FFp(z)->FFp(f->FFi(z))!=f " << endl;
}
}
}
void Quadrangulate(TriMesh &Tmesh,
PolyMesh &Pmesh,
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<TriMesh>::CompactVertexVector(Tmesh);
vcg::tri::Allocator<TriMesh>::CompactFaceVector(Tmesh);
vcg::tri::UpdateTopology<TriMesh>::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<Tmesh.face.size();i++)
Tmesh.face[i].C()=vcg::Color4b(255,255,255,255);
std::vector<std::vector<TriVertexType *> > polygons;
FindPolygons(Tmesh,polygons,UV);
//then add to the polygonal mesh
Pmesh.Clear();
int numV=vcg::tri::UpdateSelection<TriMesh>::VertexCount(Tmesh);
//first create vertices
vcg::tri::Allocator<PolyMesh>::AddVertices(Pmesh,numV);
std::map<CoordType,int> VertMap;
int index=0;
for(unsigned int i=0;i<Tmesh.vert.size();i++)
{
if (!Tmesh.vert[i].IsS())continue;
CoordType pos=Tmesh.vert[i].P();
CoordType norm=Tmesh.vert[i].N();
vcg::Point2<ScalarType> 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<PolyMesh>::AddFaces(Pmesh,polygons.size());
for (unsigned int i=0;i<polygons.size();i++)
{
int size=polygons[i].size();
Pmesh.face[i].Alloc(size);
for (int j=0;j<size;j++)
{
CoordType pos=(polygons[i][j])->P();
assert(VertMap.count(pos)==1);
int index=VertMap[pos];
Pmesh.face[i].V(j)=&(Pmesh.vert[index]);
}
}
}
};
}
}
#endif