manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
vcglib/wrap/miq/core/poisson_solver.h

861 lines
25 KiB
C++
Raw Blame History

#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"
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
Reference in New Issue View Git Blame Copy Permalink