diff --git a/vcg/complex/algorithms/parametrization/poisson_solver.h b/vcg/complex/algorithms/parametrization/poisson_solver.h index 19560d98..a34a7d67 100644 --- a/vcg/complex/algorithms/parametrization/poisson_solver.h +++ b/vcg/complex/algorithms/parametrization/poisson_solver.h @@ -26,6 +26,8 @@ #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET #include +//#include +//#include #include #include @@ -34,785 +36,800 @@ #include #include #include +#include namespace vcg { - namespace tri{ +namespace tri{ template class PoissonSolver { - typedef typename MeshType::ScalarType ScalarType; - typedef typename MeshType::FaceType FaceType; - typedef typename MeshType::VertexType VertexType; - typedef typename MeshType::CoordType CoordType; - typedef typename MeshType:: template PerFaceAttributeHandle PerFaceCoordHandle; + typedef typename MeshType::ScalarType ScalarType; + typedef typename MeshType::FaceType FaceType; + typedef typename MeshType::VertexType VertexType; + typedef typename MeshType::CoordType CoordType; + typedef typename MeshType:: template PerFaceAttributeHandle PerFaceCoordHandle; - ///the mesh itself - MeshType &mesh; + ///the mesh itself + MeshType &mesh; - ///solver data - - std::map VertexToInd; - std::map IndToVertex; - - ///vertices to fix - std::vector to_fix; - - ///unknown vector - - Eigen:: DynamicSparseMatrix A; // A - Eigen::VectorXd b,x;// x and b - - //number of variables - unsigned int n_vert_vars; - ///total system size - unsigned int total_size; - ///number of fixed variables - unsigned int n_fixed_vars; + ///solver data - ///if you intend to follow the cross field - bool use_direction_field,fix_selected,correct_fixed; - ///size of the scalar field - ScalarType fieldScale; - ///handle per direction field - PerFaceCoordHandle Fh0,Fh1; + std::map VertexToInd; + std::map IndToVertex; - int VertexIndex(VertexType* v) - { - typename std::map::iterator iteMap=VertexToInd.find(v); - assert(iteMap!=VertexToInd.end()); - return ((*iteMap).second); - } - - VertexType* IndexVertex(int index) - { - typename std::map::iterator iteMap=IndToVertex.find(index); - assert(iteMap!=IndToVertex.end()); - return ((*iteMap).second); - } + ///vertices to fix + std::vector to_fix; - void AddVertexIndex(VertexType* v,int index) - { - VertexToInd.insert(std::pair(v,index)); - IndToVertex.insert(std::pair(index,v)); - } - ///set the value of A of the system Ax=b - void SetValA(int Xindex,int Yindex,ScalarType val) - { - //int size=(int)S.nrows(); - assert(0 <= Xindex && Xindex < int(total_size)); - assert(0 <= Yindex && Yindex < int(total_size)); - //S.A().addEntryReal(Xindex,Yindex,val); - //if (Xindex>=Yindex) - A.coeffRef(Xindex,Yindex) +=val; - - } - - /*void FindFarestVert(VertexType* &v0,VertexType* &v1) - { - UpdateBounding::Box(mesh); - ScalarType d0=mesh.bbox.Diag(); - ScalarType d1=d0; - v0=NULL; - v1=NULL; - for (unsigned int j=0;jIsD()) - { - ScalarType d_test=(v->P()-mesh.bbox.min).Norm(); - if (d_testP()-mesh.bbox.max).Norm(); - if (d_test::Box(mesh); + ///unknown vector - tri::UpdateTopology::FaceFace(mesh); - tri::UpdateFlags::FaceBorderFromFF(mesh); - tri::UpdateFlags::VertexBorderFromFace(mesh); + Eigen::DynamicSparseMatrix A; // A + Eigen::VectorXd b,x;// x and b - ScalarType dmax=0; - v0=NULL; - v1=NULL; - for (unsigned int i=0;iIsD())continue; - if (vt1->IsD())continue; - if (!vt0->IsB())continue; - if (!vt1->IsB())continue; - ScalarType d_test=(vt0->P()-vt1->P()).Norm(); - if (d_test>dmax) - { - dmax=d_test; - v0=vt0; - v1=vt1; - } - } - assert(v0!=NULL); - assert(v1!=NULL); - } + //number of variables + unsigned int n_vert_vars; + ///total system size + unsigned int total_size; + ///number of fixed variables + unsigned int n_fixed_vars; - ///set the value of b of the system Ax=b - void SetValB(int Xindex, - ScalarType val) - { - /*S.b()[Xindex] += val;*/ - b[Xindex] += val; - } + ///if you intend to follow the cross field + bool use_direction_field,fix_selected,correct_fixed; + ///size of the scalar field + ScalarType fieldScale; + ///handle per direction field + PerFaceCoordHandle Fh0,Fh1; - ///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} }; + int VertexIndex(VertexType* v) + { + typename std::map::iterator iteMap=VertexToInd.find(v); + assert(iteMap!=VertexToInd.end()); + return ((*iteMap).second); + } - 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; + VertexType* IndexVertex(int index) + { + typename std::map::iterator iteMap=IndToVertex.find(index); + assert(iteMap!=IndToVertex.end()); + return ((*iteMap).second); + } - SetValA(Xindex+1,Yindex,-val[i][j]); - SetValA(Xindex,Yindex+1,val[i][j]); - - } - } + void AddVertexIndex(VertexType* v,int index) + { + VertexToInd.insert(std::pair(v,index)); + IndToVertex.insert(std::pair(index,v)); + } + ///set the value of A of the system Ax=b + void SetValA(int Xindex,int Yindex,ScalarType val) + { + //int size=(int)S.nrows(); + assert(0 <= Xindex && Xindex < int(total_size)); + assert(0 <= Yindex && Yindex < int(total_size)); + //S.A().addEntryReal(Xindex,Yindex,val); + //if (Xindex>=Yindex) + A.coeffRef(Xindex,Yindex) +=val; - ///set the diagonal of the matrix (which is zero at the beginning) - ///as the sum of the other element inverted by sign - 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; - } - } + } - ///add this values to the right hand side - 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]; - SetValB((index[i]*2),valU); - SetValB((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(Xindex::Box(mesh); + ScalarType d0=mesh.bbox.Diag(); + ScalarType d1=d0; + v0=NULL; + v1=NULL; + for (unsigned int j=0;jIsD()) + { + ScalarType d_test=(v->P()-mesh.bbox.min).Norm(); + if (d_testP()-mesh.bbox.max).Norm(); + if (d_test::Box(mesh); - ///get the vertices - VertexType *v[3]; - v[0]=f->V(0); - v[1]=f->V(1); - v[2]=f->V(2); + tri::UpdateTopology::FaceFace(mesh); + tri::UpdateFlags::FaceBorderFromFF(mesh); + tri::UpdateFlags::VertexBorderFromFace(mesh); - ///get the indexes of vertex instance (to consider cuts) - ///for the current face - int Vindexes[3]; - Vindexes[0]=VertexIndex(f->V(0)); - Vindexes[1]=VertexIndex(f->V(1)); - Vindexes[2]=VertexIndex(f->V(2)); + ScalarType dmax=0; + v0=NULL; + v1=NULL; + for (unsigned int i=0;iIsD())continue; + if (vt1->IsD())continue; + if (!vt0->IsB())continue; + if (!vt1->IsB())continue; +// ScalarType d_test=(vt0->P()-vt1->P()).Norm(); + ScalarType Dx=fabs(vt0->P().X()-vt1->P().X()); + ScalarType Dy=fabs(vt0->P().Y()-vt1->P().Y()); + ScalarType Dz=fabs(vt0->P().Z()-vt1->P().Z()); - ///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]; - } +// ScalarType d_test=std::max(Dx,std::max(Dy,Dz)); + ScalarType d_test=std::max(fabs(Dx-Dy),std::max(fabs(Dx-Dz),fabs(Dy-Dz))); + if (d_test>dmax) + { + dmax=d_test; + v0=vt0; + v1=vt1; + } + } + assert(v0!=NULL); + assert(v1!=NULL); + } - ///initialize edges - CoordType e[3]; - for (int k=0;k<3;k++) - e[k]=v[(k+2)%3]->P()-v[(k+1)%3]->P(); + ///set the value of b of the system Ax=b + void SetValB(int Xindex, + ScalarType val) + { + /*S.b()[Xindex] += val;*/ + b[Xindex] += val; + } - ///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; - 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); - } + ///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} }; - ///set the matrix as diagonal - SetDiagonal(val); - } + 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; - ///return the RHS for a given face - void perElementRHS(FaceType *f, - ScalarType b[6], - ScalarType vector_field_scale=1) - { + SetValA(Xindex+1,Yindex,-val[i][j]); + SetValA(Xindex,Yindex+1,val[i][j]); - /// then set the rhs - CoordType scaled_Kreal; - CoordType scaled_Kimag; - CoordType fNorm=f->N(); - fNorm.Normalize(); - 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; - /*MyMesh::PerFaceCoordHandle Fh = tri::Allocator::AddPerVertexAttribute (m,std::string("Irradiance")); - bool CrossDir0 = tri::HasPerVertexAttribute(mesh,"CrossDir0"); - bool CrossDir1 = tri::HasPerVertexAttribute(mesh,"CrossDir1"); - assert(CrossDir0); - assert(CrossDir1);*/ - - //K1=f->Q3(); - K1=Fh0[f]; - K1.Normalize(); - //K2=fNorm^K1; - K2=Fh1[f]; - K2.Normalize(); - - scaled_Kreal = K1*(vector_field_scale);///2); - scaled_Kimag = K2*(vector_field_scale);///2); - - b[0] = scaled_Kreal * neg_t[0]; - b[1] = scaled_Kimag * neg_t[0]; - b[2] = scaled_Kreal * neg_t[1]; - b[3] = scaled_Kimag * neg_t[1]; - b[4] = scaled_Kreal * neg_t[2]; - b[5] = scaled_Kimag * neg_t[2]; - ////fine codice mio - } + } + } - ///return the LHS and RHS for a given face - 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); + ///set the diagonal of the matrix (which is zero at the beginning) + ///as the sum of the other element inverted by sign + 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; + } + } - if (use_direction_field) - perElementRHS(f,b,vector_field_scale); - } + ///add this values to the right hand side + 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]; + SetValB((index[i]*2),valU); + SetValB((index[i]*2)+1,valV); + } + } - void FixPointLSquares() - { - ScalarType penalization=1000; - int offset_row=n_vert_vars; - assert(to_fix.size()>0); - for (size_t i=0;iIsD()); - int index=VertexIndex(v); - //v->vertex_index[0]; - int indexvert=index*2; - int indexRow=(offset_row+i)*2; + ///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(XindexT().U()*penalization; - ScalarType V=v->T().V()*penalization; - SetValB(indexRow,U); - SetValB(indexRow+1,V); + } - /*///set upper right part - SetValA(indexvert,indexCol,penalization); - SetValA(indexvert+1,indexCol+1,penalization);*/ + ///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(Xindex<(n_vert_vars*2)); + assert(Yindex<(n_vert_vars*2)); + SetValA(Xindex,Yindex,val[i][j]); + SetValA(Xindex+1,Yindex+1,val[i][j]); + } - SetValA(indexvert,indexvert,penalization); - SetValA(indexvert+1,indexvert+1,penalization); - SetValA(indexRow,indexRow,penalization); - SetValA(indexRow+1,indexRow+1,penalization); - SetValA(indexvert,indexRow,-penalization); - SetValA(indexvert+1,indexRow+1,-penalization); - SetValA(indexRow,indexvert,-penalization); - SetValA(indexRow+1,indexvert+1,-penalization); - //SetValA(indexCol+1,indexCol+1,-1); - } - } + } - //build the laplacian matrix cyclyng over all rangemaps - //and over all faces - void BuildLaplacianMatrix(double vfscale=1) - { + ///return the LHS for a given face + 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; - ///then for each face - for (unsigned int j=0;jV(0); + v[1]=f->V(1); + v[2]=f->V(2); - FaceType *f=&mesh.face[j]; - if (f->IsD()) - continue; + ///get the indexes of vertex instance (to consider cuts) + ///for the current face + int Vindexes[3]; + Vindexes[0]=VertexIndex(f->V(0)); + Vindexes[1]=VertexIndex(f->V(1)); + Vindexes[2]=VertexIndex(f->V(2)); - int var_idx[3];//vertex variable indices - for(int k = 0; k < 3; ++k) - { - VertexType *v=f->V(k); - var_idx[k] = VertexIndex(v); - } - ScalarType val[3][3]; - int index[3][3][2]; - ScalarType b[6]; - PerElementSystemReal(f, val,index, b, vfscale); + ///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]; + } - //Add the element to the matrix - Add33Block(val,index); + ///initialize edges + CoordType e[3]; + for (int k=0;k<3;k++) + e[k]=v[(k+2)%3]->P()-v[(k+1)%3]->P(); - /////add area term.. to test if needed - /*if (!use_direction_field) - AddAreaTerm(index,1.0);//f->area);*/ - /*ScalarType area=((f->P(1)-f->P(0))^(f->P(2)-f->P(0))).Norm(); - if (!use_direction_field) - AddAreaTerm(index,area);*/ + ///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; + 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); + } - //ScalarType area=((f->P(1)-f->P(0))^(f->P(2)-f->P(0))).Norm(); - if (!use_direction_field) - AddAreaTerm(index,1); + ///set the matrix as diagonal + SetDiagonal(val); + } - ///add right hand side - if (use_direction_field) - AddRHS(b,var_idx); - } - } + ///return the RHS for a given face + 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(); + fNorm.Normalize(); + 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; + /*MyMesh::PerFaceCoordHandle Fh = tri::Allocator::AddPerVertexAttribute (m,std::string("Irradiance")); + bool CrossDir0 = tri::HasPerVertexAttribute(mesh,"CrossDir0"); + bool CrossDir1 = tri::HasPerVertexAttribute(mesh,"CrossDir1"); + assert(CrossDir0); + assert(CrossDir1);*/ + + //K1=f->Q3(); + K1=Fh0[f]; + K1.Normalize(); + //K2=fNorm^K1; + K2=Fh1[f]; + K2.Normalize(); + + scaled_Kreal = K1*(vector_field_scale);///2); + scaled_Kimag = K2*(vector_field_scale);///2); + + b[0] = scaled_Kreal * neg_t[0]; + b[1] = scaled_Kimag * neg_t[0]; + b[2] = scaled_Kreal * neg_t[1]; + b[3] = scaled_Kimag * neg_t[1]; + b[4] = scaled_Kreal * neg_t[2]; + b[5] = scaled_Kimag * neg_t[2]; + ////fine codice mio + } + + ///return the LHS and RHS for a given face + 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); + + if (use_direction_field) + perElementRHS(f,b,vector_field_scale); + } + + void FixPointLSquares() + { + ScalarType penalization=1000; + int offset_row=n_vert_vars; + assert(to_fix.size()>0); + for (size_t i=0;iIsD()); + int index=VertexIndex(v); + //v->vertex_index[0]; + int indexvert=index*2; + int indexRow=(offset_row+i)*2; + + SetValA(indexRow,indexRow,penalization); + SetValA(indexRow+1,indexRow+1,penalization); + + ///add values to the B vector + ScalarType U=v->T().U()*penalization; + ScalarType V=v->T().V()*penalization; + SetValB(indexRow,U); + SetValB(indexRow+1,V); + + /*///set upper right part + SetValA(indexvert,indexCol,penalization); + SetValA(indexvert+1,indexCol+1,penalization);*/ + + SetValA(indexvert,indexvert,penalization); + SetValA(indexvert+1,indexvert+1,penalization); + SetValA(indexRow,indexRow,penalization); + SetValA(indexRow+1,indexRow+1,penalization); + SetValA(indexvert,indexRow,-penalization); + SetValA(indexvert+1,indexRow+1,-penalization); + SetValA(indexRow,indexvert,-penalization); + SetValA(indexRow+1,indexvert+1,-penalization); + //SetValA(indexCol+1,indexCol+1,-1); + } + } + + //build the laplacian matrix cyclyng over all rangemaps + //and over all faces + void BuildLaplacianMatrix(double vfscale=1) + { + + ///then for each face + for (unsigned int j=0;jIsD()) + continue; + + int var_idx[3];//vertex variable indices + for(int k = 0; k < 3; ++k) + { + VertexType *v=f->V(k); + var_idx[k] = VertexIndex(v); + } + ScalarType val[3][3]; + int index[3][3][2]; + ScalarType b[6]; + 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);*/ + /*ScalarType area=((f->P(1)-f->P(0))^(f->P(2)-f->P(0))).Norm(); + if (!use_direction_field) + AddAreaTerm(index,area);*/ + + //ScalarType area=((f->P(1)-f->P(0))^(f->P(2)-f->P(0))).Norm(); + if (!use_direction_field) + AddAreaTerm(index,1); + + ///add right hand side + if (use_direction_field) + AddRHS(b,var_idx); + } + } - void FindSizes() - { - // tag vertices and compute numbers of equations to determine the number of rows in the matrix - //TagVertices_Constrained(n_vert_vars, n_transition_eqs, n_align_sharp_eqs); - n_vert_vars=mesh.vn; + void FindSizes() + { + // tag vertices and compute numbers of equations to determine the number of rows in the matrix + //TagVertices_Constrained(n_vert_vars, n_transition_eqs, n_align_sharp_eqs); + n_vert_vars=mesh.vn; - ///initialize matrix size - total_size = (n_fixed_vars + n_vert_vars)*2;///must be multiplied by 2 becasue of u and v - - } + ///initialize matrix size + total_size = (n_fixed_vars + n_vert_vars)*2;///must be multiplied by 2 becasue of u and v - void AllocateSystem() - { - //--- Allocates the data for Ax=b - A=Eigen::DynamicSparseMatrix(total_size, total_size); // A - b = Eigen::VectorXd::Zero(total_size); // x and b - } + } - + void AllocateSystem() + { + //--- Allocates the data for Ax=b + A=Eigen::DynamicSparseMatrix(total_size, total_size); // A + b = Eigen::VectorXd::Zero(total_size); // x and b + } - ///intitialize the whole matrix - void InitMatrix() - { - FindSizes(); - AllocateSystem(); - } - - bool Solve() - { - //A.finalize(); - //Eigen::SparseMatrix As=Eigen::SparseMatrix(A); - //As.finalize(); - // - //SparseLDLT,Taucs> ldlt_of_A(As); - //if(!ldlt_of_A.succeeded()) - //{ - // printf("Decomposition Failed \n"); - // return false; - //} - ///*printf("\n");*/ - - //ldlt_of_A.solveInPlace(b); - - //return true; - A.finalize(); - Eigen::SparseMatrix As=Eigen::SparseMatrix(A); - As.finalize(); - Eigen::SimplicialCholesky > solver(As); - x = solver.solve(b); - return (solver.info()==Eigen::Success); - } - - - void InitIndex() - { - for (size_t i=0;i::VertexClearV(mesh); - //set fixed to V - for (size_t i=0;iSetV(); + ///intitialize the whole matrix + void InitMatrix() + { + FindSizes(); + AllocateSystem(); + } - Box2 bbox; - if (normalize) - { - for (size_t i=0;i(U,V)); - } - } + bool Solve() + { + //A.finalize(); + //Eigen::SparseMatrix As=Eigen::SparseMatrix(A); + //As.finalize(); + // + //SparseLDLT,Taucs> ldlt_of_A(As); + //if(!ldlt_of_A.succeeded()) + //{ + // printf("Decomposition Failed \n"); + // return false; + //} + ///*printf("\n");*/ - //for each vertex - for (size_t i=0;i p; - if (!v->IsV()) - p=Point2(U,V); - else - p=v->T().P(); - //p/=fieldScale; - if (normalize) - { - p-=bbox.min; - p*=1/bbox.Diag(); - } - - v->T().P()=p; - } + //ldlt_of_A.solveInPlace(b); - ///then copy to faces - for (size_t i=0;iV(j); - Point2 p=v->T().P(); - f->WT(j).P()=p; - } - } - } + //return true; + A.finalize(); + Eigen::SparseMatrix As=Eigen::SparseMatrix(A); + As.finalize(); + + Eigen::SimplicialCholesky > solver(As); + x = solver.solve(b); + return (solver.info()==Eigen::Success); + } + + + void InitIndex() + { + for (size_t i=0;i::VertexClearV(mesh); + //set fixed to V + for (size_t i=0;iSetV(); + + Box2 bbox; + if (normalize) + { + for (size_t i=0;i(U,V)); + } + } + + //for each vertex + for (size_t i=0;i p; + if (!v->IsV()) + p=Point2(U,V); + else + p=v->T().P(); + //p/=fieldScale; + if (normalize) + { + p-=bbox.min; + p*=1/bbox.Diag(); + } + + v->T().P()=p; + } + + ///then copy to faces + for (size_t i=0;iV(j); + Point2 p=v->T().P(); + f->WT(j).P()=p; + } + } + } public: - ///return true if is possible to - bool IsFeaseable() - { - tri::UpdateTopology::FaceFace(mesh); - int NNmanifoldE=tri::Clean::CountNonManifoldEdgeFF(mesh); - if (NNmanifoldE!=0)return false; - /*int NNmanifoldV=tri::Clean::CountNonManifoldVertexFF(mesh); - if (NNmanifoldV!=0)return false;*/ - int G=tri::Clean::MeshGenus(mesh); - int numholes=tri::Clean::CountHoles(mesh); - if (numholes==0) return false; - return (G==0); - } + ///return true if is possible to + bool IsFeaseable() + { + tri::UpdateTopology::FaceFace(mesh); + int NNmanifoldE=tri::Clean::CountNonManifoldEdgeFF(mesh); + if (NNmanifoldE!=0) + { + printf("Non Manifold"); + return false; + } + /*int NNmanifoldV=tri::Clean::CountNonManifoldVertexFF(mesh); + if (NNmanifoldV!=0)return false;*/ + int G=tri::Clean::MeshGenus(mesh); + int numholes=tri::Clean::CountHoles(mesh); + if (numholes==0) + { + printf("Non omeomorph to a disc"); + return false; + } + return (G==0); + } - ///set the border as fixed - void SetBorderAsFixed() - { - for (size_t i=0;iIsD())continue; - if(v->IsB())to_fix.push_back(v); - } - std::sort(to_fix.begin(),to_fix.end()); - typename std::vector::iterator new_end=std::unique(to_fix.begin(),to_fix.end()); - int dist=distance(to_fix.begin(),new_end); - to_fix.resize(dist); - } - - ///set selected vertices as fixed - void SetSelectedAsFixed() - { - for (int i=0;iIsD())continue; - if(v->IsS())to_fix.push_back(v); - } - std::sort(to_fix.begin(),to_fix.end()); - typename std::vector::iterator new_end=std::unique(to_fix.begin(),to_fix.end()); - int dist=distance(to_fix.begin(),new_end); - to_fix.resize(dist); - } - - /*///fix default vertices no need if already border on other vertices are fixed - ///you need at least 2 fixed for solving without field , - ///while only 1 if you conforms to a given cross field - void FixDefaultVertices() - { - ///in this case there are already vertices fixed, so no need to fix by default - assert(to_fix.size()==0); - ///then fix only one vertex - if (use_direction_field) - { - for (size_t i=0;i(0,0); - to_fix.push_back(&mesh.vert[i]); - return; - } - } - ///then fix 2 vertices - else - { - VertexType *v0; - VertexType *v1; - FindFarestVert(v0,v1); - if (v0==v1) - { - tri::io::ExporterPLY::Save(mesh,"./parametrized.ply"); - assert(0); - } - v0->T().P()=Point2(0,0); - v1->T().P()=Point2(1,0); - to_fix.push_back(v0); - to_fix.push_back(v1); - return; - } - }*/ - - ///fix default vertices no need if already border on other vertices are fixed - ///you need at least 2 fixed for solving without field , - ///while only 1 if you conforms to a given cross field - void FixDefaultVertices() - { - ///in this case there are already vertices fixed, so no need to fix by default - assert(to_fix.size()==0); - ///then fix only one vertex - if (use_direction_field) - { - for (size_t i=0;i(0,0); - to_fix.push_back(&mesh.vert[i]); - return; - } - } - ///then fix 2 vertices - else - { - VertexType *v0; - VertexType *v1; - FindFarthestVert(v0,v1); - if (v0==v1) - { -// tri::io::ExporterPLY::Save(mesh,"./parametrized.ply"); - assert(0); - } - v0->T().P()=Point2(0,0); - v1->T().P()=Point2(1,0); - to_fix.push_back(v0); - to_fix.push_back(v1); - return; - } - } - ///intialize parameters and setup fixed vertices vector - void Init(bool _use_direction_field=false, - bool _correct_fixed=true, - ScalarType _fieldScale=1.0) - { - use_direction_field=_use_direction_field; - //query if an attribute is present or not - if (use_direction_field) - { - bool CrossDir0 = tri::HasPerFaceAttribute(mesh,"CrossDir0"); - bool CrossDir1 = tri::HasPerFaceAttribute(mesh,"CrossDir1"); - assert(CrossDir0); - assert(CrossDir1); + ///set the border as fixed + void SetBorderAsFixed() + { + for (size_t i=0;iIsD())continue; + if(v->IsB())to_fix.push_back(v); + } + std::sort(to_fix.begin(),to_fix.end()); + typename std::vector::iterator new_end=std::unique(to_fix.begin(),to_fix.end()); + int dist=distance(to_fix.begin(),new_end); + to_fix.resize(dist); + } + + ///set selected vertices as fixed + void SetSelectedAsFixed() + { + for (int i=0;iIsD())continue; + if(v->IsS())to_fix.push_back(v); + } + std::sort(to_fix.begin(),to_fix.end()); + typename std::vector::iterator new_end=std::unique(to_fix.begin(),to_fix.end()); + int dist=distance(to_fix.begin(),new_end); + to_fix.resize(dist); + } + + /*///fix default vertices no need if already border on other vertices are fixed + ///you need at least 2 fixed for solving without field , + ///while only 1 if you conforms to a given cross field + void FixDefaultVertices() + { + ///in this case there are already vertices fixed, so no need to fix by default + assert(to_fix.size()==0); + ///then fix only one vertex + if (use_direction_field) + { + for (size_t i=0;i(0,0); + to_fix.push_back(&mesh.vert[i]); + return; + } + } + ///then fix 2 vertices + else + { + VertexType *v0; + VertexType *v1; + FindFarestVert(v0,v1); + if (v0==v1) + { + tri::io::ExporterPLY::Save(mesh,"./parametrized.ply"); + assert(0); + } + v0->T().P()=Point2(0,0); + v1->T().P()=Point2(1,0); + to_fix.push_back(v0); + to_fix.push_back(v1); + return; + } + }*/ + + ///fix default vertices no need if already border on other vertices are fixed + ///you need at least 2 fixed for solving without field , + ///while only 1 if you conforms to a given cross field + void FixDefaultVertices() + { + ///in this case there are already vertices fixed, so no need to fix by default + assert(to_fix.size()==0); + ///then fix only one vertex + if (use_direction_field) + { + for (size_t i=0;i(0,0); + to_fix.push_back(&mesh.vert[i]); + return; + } + } + ///then fix 2 vertices + else + { + VertexType *v0; + VertexType *v1; + FindFarthestVert(v0,v1); + if (v0==v1) + { + // tri::io::ExporterPLY::Save(mesh,"./parametrized.ply"); + assert(0); + } + v0->T().P()=Point2(0,0); + v1->T().P()=Point2(1,0); + to_fix.push_back(v0); + to_fix.push_back(v1); + return; + } + } + ///intialize parameters and setup fixed vertices vector + void Init(bool _use_direction_field=false, + bool _correct_fixed=true, + ScalarType _fieldScale=1.0) + { + use_direction_field=_use_direction_field; + //query if an attribute is present or not + if (use_direction_field) + { + bool CrossDir0 = tri::HasPerFaceAttribute(mesh,"CrossDir0"); + bool CrossDir1 = tri::HasPerFaceAttribute(mesh,"CrossDir1"); + assert(CrossDir0); + assert(CrossDir1); Fh0= tri::Allocator :: template GetPerFaceAttribute(mesh,std::string("CrossDir0")); Fh1= tri::Allocator :: template GetPerFaceAttribute(mesh,std::string("CrossDir1")); - } - correct_fixed=_correct_fixed; - fieldScale=_fieldScale; - to_fix.clear(); - } + } + correct_fixed=_correct_fixed; + fieldScale=_fieldScale; + to_fix.clear(); + } - ///solve the system, it return false if the matrix is singular - bool SolvePoisson(bool _write_messages=false, - ScalarType fieldScale=1.0, - bool solve_global_fold=true) - { - int t0,t1,t2,t3; + ///solve the system, it return false if the matrix is singular + bool SolvePoisson(bool _write_messages=false, + ScalarType fieldScale=1.0, + bool solve_global_fold=true) + { + int t0,t1,t2,t3; - ///Initializing Matrix - if (_write_messages) - { - printf("\n INITIALIZING THE MATRIX \n"); - t0=clock(); - } - - ///set vertex indexes - InitIndex(); + ///Initializing Matrix + if (_write_messages) + { + printf("\n INITIALIZING THE MATRIX \n"); + t0=clock(); + } - /*///find vertex to fix - std::vector to_fix; - FindFixedVertices(to_fix); - n_fixed_vars=to_fix.size();*/ - if (use_direction_field) - { - assert(to_fix.size()>0); - } - else - { - assert(to_fix.size()>1); - } + ///set vertex indexes + InitIndex(); - n_fixed_vars=to_fix.size(); - ///initialize the matrix ALLOCATING SPACE - InitMatrix(); - - if (use_direction_field) - { - bool CrossDir0 = tri::HasPerFaceAttribute(mesh,"CrossDir0"); - bool CrossDir1 = tri::HasPerFaceAttribute(mesh,"CrossDir1"); - assert(CrossDir0); - assert(CrossDir1); - } + /*///find vertex to fix + std::vector to_fix; + FindFixedVertices(to_fix); + n_fixed_vars=to_fix.size();*/ + if (use_direction_field) + { + assert(to_fix.size()>0); + } + else + { + assert(to_fix.size()>1); + } - ///build the laplacian system - BuildLaplacianMatrix(fieldScale); - - ////add the lagrange multiplier - FixPointLSquares(); + n_fixed_vars=to_fix.size(); + ///initialize the matrix ALLOCATING SPACE + InitMatrix(); - if (_write_messages) - { - t1=clock(); - printf("\n time:%d \n",t1-t0); - printf("\n SOLVING \n"); - } - - //int n_vars=(n_vert_vars)*2; - //int integer_constr_size=(n_transition_vars+n_fixed_vars+n_bary_transition_vars)*2; - //X=std::vector< double >(n_vars+n_fixed_vars*2); - bool done=Solve(); - if (!done) - return false; - if (_write_messages) - { - t2=clock(); - printf("\n time:%d \n",t2-t1); - printf("\n ASSIGNING COORDS \n"); - } + if (use_direction_field) + { + bool CrossDir0 = tri::HasPerFaceAttribute(mesh,"CrossDir0"); + bool CrossDir1 = tri::HasPerFaceAttribute(mesh,"CrossDir1"); + assert(CrossDir0); + assert(CrossDir1); + } - MapCoords(false,fieldScale); - if (_write_messages) - { - t3=clock(); - printf("\n time:%d \n",t3-t2); - } + ///build the laplacian system + BuildLaplacianMatrix(fieldScale); - ///then check if majority of faces are folded - if (!solve_global_fold) return true; - if (tri::Distortion::GloballyUnFolded(mesh)) - { - tri::UV_Utils::GloballyMirrorX(mesh); - bool isUnfolded = tri::Distortion::GloballyUnFolded(mesh); - assert( ! isUnfolded); - } - return true; - } + ////add the lagrange multiplier + FixPointLSquares(); - PoissonSolver(MeshType &_mesh):mesh(_mesh) - { - assert(mesh.vert.size()>3); - assert(mesh.face.size()>1); - } - - - }; // end class - } //End Namespace Tri + if (_write_messages) + { + t1=clock(); + printf("\n time:%d \n",t1-t0); + printf("\n SOLVING \n"); + } + + //int n_vars=(n_vert_vars)*2; + //int integer_constr_size=(n_transition_vars+n_fixed_vars+n_bary_transition_vars)*2; + //X=std::vector< double >(n_vars+n_fixed_vars*2); + bool done=Solve(); + if (!done) + return false; + if (_write_messages) + { + t2=clock(); + printf("\n time:%d \n",t2-t1); + printf("\n ASSIGNING COORDS \n"); + } + + MapCoords(false,fieldScale); + if (_write_messages) + { + t3=clock(); + printf("\n time:%d \n",t3-t2); + } + + ///then check if majority of faces are folded + if (!solve_global_fold) return true; + if (tri::Distortion::GloballyUnFolded(mesh)) + { + tri::UV_Utils::GloballyMirrorX(mesh); + bool isUnfolded = tri::Distortion::GloballyUnFolded(mesh); + assert( ! isUnfolded); + } + return true; + } + + PoissonSolver(MeshType &_mesh):mesh(_mesh) + { + assert(mesh.vert.size()>3); + assert(mesh.face.size()>1); + } + + +}; // end class +} //End Namespace Tri } // End Namespace vcg #endif