diff --git a/apps/sample/trimesh_topological_cut/trimesh_topological_cut.cpp b/apps/sample/trimesh_topological_cut/trimesh_topological_cut.cpp index 67777387..6f9c63b3 100644 --- a/apps/sample/trimesh_topological_cut/trimesh_topological_cut.cpp +++ b/apps/sample/trimesh_topological_cut/trimesh_topological_cut.cpp @@ -26,6 +26,7 @@ #include +#include #include #include @@ -44,47 +45,62 @@ class MyMesh : public tri::TriMesh< std::vector, std::vector, int main(int argc,char ** argv ) { - MyMesh base,basecopy, poly; + MyMesh base, basecopy, poly; int ret0 = tri::io::Importer::Open(base,argv[1]); - tri::UpdateBounding::Box(base); - tri::Append::MeshCopy(basecopy,base); - printf( "Mesh %s has %i vert and %i faces\n", argv[1], basecopy.VN(), basecopy.FN() ); if(ret0 != 0 ) { printf("Failed Loading\n"); exit(-1); - } + } + + tri::UpdateBounding::Box(base); + printf( "Mesh %s has %i vert and %i faces\n", argv[1], base.VN(), base.FN() ); + srand(time(0)); + tri::CutTree ct(base); + ct.BuildVisitTree(poly,rand()%base.fn); + tri::io::ExporterPLY::Save(poly,"tree.ply",tri::io::Mask::IOM_EDGEINDEX); + tri::CoM cc(base); cc.Init(); - cc.BuildVisitTree(poly); - tri::UpdateBounding::Box(poly); - tri::io::ExporterPLY::Save(poly,"tree.ply",tri::io::Mask::IOM_EDGEINDEX); - while(cc.OptimizeTree(poly)); - tri::io::ExporterPLY::Save(poly,"tree1.ply",tri::io::Mask::IOM_EDGEINDEX); - - cc.MarkFauxEdgeWithPolyLine(basecopy,poly); + cc.MarkFauxEdgeWithPolyLine(poly); + tri::Append::MeshCopy(basecopy,base); tri::UpdateTopology::FaceFace(basecopy); tri::CutMeshAlongNonFauxEdges(basecopy); - tri::io::ExporterPLY::Save(basecopy,"basecut.ply"); + tri::io::ExporterPLY::Save(basecopy,"base_cut_with_tree.ply"); + + cc.par.surfDistThr = base.bbox.Diag()/100.0; - cc.par.maxSimpEdgeLen = base.bbox.Diag()/60.0; - cc.SmoothProject(poly,5,0.8, 0.2); + cc.par.maxSimpEdgeLen = base.bbox.Diag()/40.0; + cc.par.minRefEdgeLen = base.bbox.Diag()/80.0; + cc.SmoothProject(poly,40,0.5, .7); + Distribution dist; cc.EvaluateHausdorffDistance(poly, dist ); - tri::io::ExporterPLY::Save(poly,"poly_adapted.ply",tri::io::Mask::IOM_EDGEINDEX+tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); - cc.par.surfDistThr = base.bbox.Diag()/1000.0; - cc.par.maxSimpEdgeLen = base.bbox.Diag()/500.0; - cc.SmoothProject(poly,5,0.3, 0.7); + +// tri::io::ExporterPLY::Save(poly,"poly_adapted.ply",tri::io::Mask::IOM_EDGEINDEX+tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); +// cc.par.surfDistThr = base.bbox.Diag()/2000.0; +// cc.par.maxSimpEdgeLen = base.bbox.Diag()/100.0; +// cc.par.minRefEdgeLen = base.bbox.Diag()/200.0; +// cc.SmoothProject(poly,10,0.3, 0.7); +// cc.SmoothProject(poly,1,0.01, 0.99); tri::io::ExporterPLY::Save(poly,"poly_adapted2.ply",tri::io::Mask::IOM_EDGEINDEX+tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); - cc.par.minRefEdgeLen = base.bbox.Diag()/200.0; - cc.Refine(poly,true); - cc.Refine(poly,true); + std::vector newVertVec; cc.SnapPolyline(poly); - tri::io::ExporterPLY::Save(poly,"poly_adapted2_snap.ply",tri::io::Mask::IOM_EDGEINDEX+tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); - cc.par.maxSimpEdgeLen = base.bbox.Diag()/200.0; - cc.SmoothProject(poly,10 ,0.8, 0.2); - cc.RefineBaseMesh(poly); - tri::io::ExporterPLY::Save(poly,"poly_adapted2_snapSmooth.ply",tri::io::Mask::IOM_EDGEINDEX+tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); + tri::io::ExporterPLY::Save(poly,"poly_snap.ply",tri::io::Mask::IOM_EDGEINDEX+tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); + cc.RefineCurveByBaseMesh(poly); + tri::io::ExporterPLY::Save(poly,"poly_refined.ply",tri::io::Mask::IOM_EDGEINDEX+tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); + + cc.SplitMeshWithPolyline(poly); + cc.RefineCurveByBaseMesh(poly); + tri::io::ExporterPLY::Save(base,"base_refined.ply",tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); + cc.SplitMeshWithPolyline(poly); + cc.RefineCurveByBaseMesh(poly); + cc.SplitMeshWithPolyline(poly); + cc.RefineCurveByBaseMesh(poly); + tri::io::ExporterPLY::Save(base,"base_refined2.ply",tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY+tri::io::Mask::IOM_FACEFLAGS); + cc.MarkFauxEdgeWithPolyLine(poly); + CutMeshAlongNonFauxEdges(base); + tri::io::ExporterPLY::Save(base,"base_refined2_cut.ply",tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); return 0; } diff --git a/vcg/complex/algorithms/curve_on_manifold.h b/vcg/complex/algorithms/curve_on_manifold.h index 66016ae3..e9189fe5 100644 --- a/vcg/complex/algorithms/curve_on_manifold.h +++ b/vcg/complex/algorithms/curve_on_manifold.h @@ -36,7 +36,6 @@ #include #include #include -#include #include #include @@ -67,7 +66,7 @@ public: typedef typename vcg::GridStaticPtr MeshGrid; typedef typename vcg::GridStaticPtr EdgeGrid; typedef typename face::Pos PosType; - + typedef typename tri::UpdateTopology::PEdge PEdge; class Param { public: @@ -79,6 +78,7 @@ public: ScalarType maxSmoothDelta; // The maximum movement that is admitted during smoothing. ScalarType maxSnapThr; // The maximum distance allowed when snapping a vertex of the polyline onto a mesh vertex ScalarType gridBailout; // The maximum distance bailout used in grid sampling + ScalarType barycentricSnapThr; // The maximum distance bailout used in grid sampling Param(MeshType &m) { Default(m);} @@ -87,10 +87,11 @@ public: surfDistThr = m.bbox.Diag()/1000.0; polyDistThr = m.bbox.Diag()/5000.0; minRefEdgeLen = m.bbox.Diag()/16000.0; - maxSimpEdgeLen = m.bbox.Diag()/1000.0; + maxSimpEdgeLen = m.bbox.Diag()/100.0; maxSmoothDelta = m.bbox.Diag()/100.0; - maxSnapThr = m.bbox.Diag()/10000.0; + maxSnapThr = m.bbox.Diag()/1000.0; gridBailout = m.bbox.Diag()/20.0; + barycentricSnapThr = 0.05; } void Dump() const { @@ -112,336 +113,99 @@ public: Param par; CoM(MeshType &_m) :base(_m),par(_m){} - - //******** CUT TREE GENERATION - - - // Given two points return true if on the base mesh there exist an edge with that two coords - // if return true the pos indicate the found edge. - bool ExistEdge(KdTree &kdtree, CoordType &p0, CoordType &p1, PosType &fpos) + FaceType *GetClosestFace(const CoordType &p) { - ScalarType locEps = SquaredDistance(p0,p1)/100000.0; - - VertexType *v0=0,*v1=0; - unsigned int veInd; - ScalarType sqdist; - kdtree.doQueryClosest(p0,veInd,sqdist); - if(sqdist vfi(v0); - while(!vfi.End()) - { - if(vfi.V1()==v1) - { - fpos = PosType(vfi.F(),vfi.I(), v0); - return true; - } - if(vfi.V2()==v1) - { - fpos = PosType(vfi.F(),(vfi.I()+1)%3, v1); - return true; - } - ++vfi; - } + float closestDist; + CoordType closestP,closestN; + return vcg::tri::GetClosestFaceBase(base,uniformGrid,p, this->par.gridBailout, closestDist, closestP,closestN,ip); } + + FaceType *GetClosestFacePoint(const CoordType &p, CoordType &closestP) + { + float closestDist; + return vcg::tri::GetClosestFaceBase(base,uniformGrid,p, this->par.gridBailout, closestDist, closestP); + } + + bool IsSnappedEdge(CoordType &ip, int &ei) + { + for(int i=0;i<3;++i) + if(ip[i]>0.0 && ip[(i+1)%3]>0.0 && ip[(i+2)%3]==0.0 ) { + ei=i; + return true; + } + ei=-1; return false; } - - bool OptimizeTree(MeshType &t) - { - tri::Allocator::CompactEveryVector(t); - tri::UpdateTopology::VertexEdge(t); - tri::UpdateTopology::VertexFace(base); - VertexConstDataWrapper vdw(base); - KdTree kdtree(vdw); - - // First simple loop that search for 2->1 moves. - for(VertexIterator vi=t.vert.begin();vi!=t.vert.end();++vi) - { - std::vector starVec; - edge::VVStarVE(&*vi,starVec); - if(starVec.size()==2) - { - PosType pos; - if(ExistEdge(kdtree,starVec[0]->P(),starVec[1]->P(),pos)) - edge::VEEdgeCollapse(t,&*vi); - } - } - return (t.en < t.edge.size()); - } - - int findNonVisitedEdgesDuringRetract(VertexType * vp, EdgeType * &ep) - { - std::vector starVec; - edge::VEStarVE(&*vp,starVec); - int cnt =0; - for(size_t i=0;iIsV()) { - cnt++; - ep = starVec[i]; - } - } - return cnt; - } - - void Retract(MeshType &t) - { - tri::Clean::RemoveDuplicateVertex(t); - printf("Retracting a tree of %i edges and %i vertices\n",t.en,t.vn); - tri::UpdateTopology::VertexEdge(t); - - std::stack vertStack; - - // Put on the stack all the vertex with just a single incident edge. - for(VertexIterator vi=t.vert.begin();vi!=t.vert.end();++vi) - { - std::vector starVec; - edge::VEStarVE(&*vi,starVec); - if(starVec.size()==1) - vertStack.push(&*vi); - } - - tri::UpdateFlags::EdgeClearV(t); - tri::UpdateFlags::VertexClearV(t); - - int unvisitedEdgeNum = t.en; - while((!vertStack.empty()) && (unvisitedEdgeNum > 2) ) - { - VertexType *vp = vertStack.top(); - vertStack.pop(); - vp->C()=Color4b::Blue; - EdgeType *ep=0; - int eCnt = findNonVisitedEdgesDuringRetract(vp,ep); - if(eCnt==1) // We have only one non visited edge over vp - { - assert(!ep->IsV()); - ep->SetV(); - --unvisitedEdgeNum; - VertexType *otherVertP; - if(ep->V(0)==vp) otherVertP = ep->V(1); - else otherVertP = ep->V(0); - vertStack.push(otherVertP); - } - } - assert(unvisitedEdgeNum >0); - - for(size_t i =0; i::DeleteEdge(t,t.edge[i]); - assert(t.en >0); - tri::Clean::RemoveUnreferencedVertex(t); - tri::Allocator::CompactEveryVector(t); - } - - void CleanSpuriousDanglingEdges(MeshType &poly) - { - EdgeGrid edgeGrid; - edgeGrid.Set(poly.edge.begin(), poly.edge.end()); - std::vector< std::pair > mergeVec; - UpdateFlags::FaceClearV(base); - UpdateTopology::FaceFace(base); - for(int i=0;iIsV()) - { - for(int j=0;j<3;++j) - if(face::IsBorder(*fp,j)) - { - face::Pos startPos(fp,int(j)); - assert(startPos.IsBorder()); - face::Pos curPos=startPos; - face::Pos prevPos=startPos; - int edgeCnt=0; - do - { - prevPos=curPos; - curPos.F()->SetV(); - curPos.NextB(); - edgeCnt++; - if(Distance(curPos.V()->P(),prevPos.VFlip()->P())<0.000000001f) - { - Point3f endp = curPos.VFlip()->P(); - float minDist=par.gridBailout; - Point3f closestP; - EdgeType *cep = vcg::tri::GetClosestEdgeBase(poly,edgeGrid,endp,par.gridBailout,minDist,closestP); - if(minDist > par.polyDistThr){ - mergeVec.push_back(std::make_pair(curPos.V(),prevPos.VFlip())); - printf("Vertex %i and %i should be merged\n",tri::Index(base,curPos.V()),tri::Index(base,prevPos.VFlip())); - } - } - - } while(curPos!=startPos); - printf("walked along a border of %i edges\n",edgeCnt); - break; - } - } - } - printf("Found %i vertex pairs to be merged\n",mergeVec.size()); - for(int k=0;kV(j)==vdel) fp->V(j)=vmerge; - } - Allocator::DeleteVertex(base,*vdel); - } - Allocator::CompactEveryVector(base); - - - for(int i=0;iV0(j)==fp->V1(j)) - { - Allocator::DeleteFace(base,*fp); - printf("Deleted face %i\n",tri::Index(base,fp)); - } - - } - Allocator::CompactEveryVector(base); - } - - - // \brief This function build a cut tree. - // - // First we make a bread first FF face visit. - // Each time that we encounter a visited face we cut add to the tree the edge - // that brings to the already visited face. - // this structure build a dense graph and we retract this graph retracting each - // leaf until we remains with just the loops that cuts the object. - - void BuildVisitTree(MeshType &dualMesh) - { - tri::UpdateTopology::FaceFace(base); - tri::UpdateFlags::FaceClearV(base); - - std::vector > visitStack; // the stack contain the pos on the 'starting' face. - - - MeshType treeMesh; // this mesh will contain the set of the non traversed edge - - base.face[0].SetV(); - for(int i=0;i<3;++i) - visitStack.push_back(PosType(&(base.face[0]),i,base.face[0].V(i))); - - int cnt=1; - - while(!visitStack.empty()) - { - std::swap(visitStack.back(),visitStack[rand()%visitStack.size()]); - PosType c = visitStack.back(); - visitStack.pop_back(); - assert(c.F()->IsV()); - c.F()->C() = Color4b::ColorRamp(0,base.fn,cnt); - c.FlipF(); - if(!c.F()->IsV()) - { - tri::Allocator::AddEdge(treeMesh,Barycenter(*(c.FFlip())),Barycenter(*(c.F()))); - ++cnt; - c.F()->SetV(); - c.FlipE();c.FlipV(); - visitStack.push_back(c); - c.FlipE();c.FlipV(); - visitStack.push_back(c); - } - else - { - tri::Allocator::AddEdge(dualMesh,c.V()->P(),c.VFlip()->P()); - } - } - assert(cnt==base.fn); - - Retract(dualMesh); -} - - /* - * Given a mesh and a polyline whose edges are a subset of edges of m - * This function marks the edges of m as non-faux when they coincide with the polyline ones. - * - * Use this function toghether with the CutMeshAlongCrease function to actually cut the mesh. - */ - - void MarkFauxEdgeWithPolyLine(MeshType &m, MeshType &e) - { - tri::UpdateFlags::FaceSetF(m); - tri::UpdateTopology::VertexEdge(e); - VertexConstDataWrapper vdw(e); - KdTree edgeTree(vdw); - - for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) - { - for(int i=0;i<3;++i) - { - ScalarType locEps = SquaredDistance(fi->P0(i), fi->P1(i))/100000.0; - unsigned int veInd; - ScalarType sqdist; - edgeTree.doQueryClosest(fi->P(i),veInd,sqdist); - if(sqdist starVecVp; - edge::VVStarVE(&(e.vert[veInd]),starVecVp); - for(size_t j=0;jP(),fi->P1(i))< locEps) - fi->ClearF(i); - } - } - } - } - - } - - - - - void SnapPolyline(MeshType &t) - { - printf("Selected vert num %i\n",tri::UpdateSelection::VertexCount(t)); - - tri::UpdateFlags::VertexClearS(base); - tri::UpdateFlags::VertexClearS(t); - tri::Allocator::CompactEveryVector(t); - VertexConstDataWrapper vdw(base); - KdTree kdtree(vdw); - for(VertexIterator vi=t.vert.begin();vi!=t.vert.end();++vi) - { - unsigned int veInd; - ScalarType sqdist; - kdtree.doQueryClosest(vi->P(),veInd,sqdist); - VertexPointer vp = &base.vert[veInd]; - if(!vp->IsS()) - { - ScalarType dist = sqrt(sqdist); - std::vector starVecVp; - face::VVStarVF(vp,starVecVp); - ScalarType minEdgeLen = std::numeric_limits::max(); - for(int i=0;iP(),starVecVp[i]->P()),minEdgeLen); - if(dist < minEdgeLen/3.0) - { - vi->P() = vp->P(); - vi->SetS(); - vp->SetS(); - } + // Given a baricentric coordinate finds that we assume that snaps onto an edge, it finds the vertex on which it is snapping + bool IsSnappedVertex(CoordType &ip, int &vi) + { + for(int i=0;i<3;++i) + if(ip[i]==1.0 && ip[(i+1)%3]==0.0 && ip[(i+2)%3]==0.0 ) { + vi=i; + return true; } - } - printf("Selected vert num %i\n",tri::UpdateSelection::VertexCount(t)); + vi=-1; + return false; + } + + // Given a baricentric coordinate finds that we assume that snaps onto an edge, it finds the vertex on which it is snapping + VertexPointer FindVertexSnap(FacePointer fp, CoordType &ip) + { + for(int i=0;i<3;++i) + if(ip[i]==1.0 && ip[(i+1)%3]==0.0 && ip[(i+2)%3]==0.0 ) return fp->V(i); + return 0; } - + /** + * @brief MarkFauxEdgeWithPolyLine marks the edges of basemesh as non-faux when they coincide with the polyline ones * + * @param poly + * @return true if all the edges of the polyline are snapped onto the mesh. + * + * Use this function toghether with the CutMeshAlongCrease function to actually cut the mesh with a snapped polyline. + * + */ + + bool MarkFauxEdgeWithPolyLine( MeshType &poly) + { + tri::UpdateFlags::FaceSetF(base); + tri::UpdateTopology::VertexFace(base); + tri::UpdateTopology::FaceFace(base); + + for(EdgeIterator ei=poly.edge.begin(); ei!=poly.edge.end();++ei) + { + Point3f ip0,ip1; + FaceType *f0 = GetClosestFaceIP(ei->cP(0),ip0); + FaceType *f1 = GetClosestFaceIP(ei->cP(1),ip1); + + if(BarycentricSnap(ip0) && BarycentricSnap(ip1)) + { + VertexPointer v0 = FindVertexSnap(f0,ip0); + VertexPointer v1 = FindVertexSnap(f1,ip1); + assert(v1>=0 && v0>=0 && v0!=v1); + FacePointer ff0,ff1; + int e0,e1; + bool ret=face::FindSharedFaces(v0,v1,ff0,ff1,e0,e1); + assert(ret); + assert(ff0->V(e0)==v0 || ff0->V(e0)==v1); + ff0->ClearF(e0); + ff1->ClearF(e1); + } + else { + assert(0); + } + } + } + float MinDistOnEdge(Point3f samplePnt, EdgeGrid &edgeGrid, MeshType &poly, Point3f &closestPoint) { float polyDist; @@ -511,411 +275,54 @@ public: return v0->P()*w0 + v1->P()*w1; } - class QualitySign + + /** + * @brief SnapPolyline snaps the vertexes of a polyline onto the base mesh + * @param poly + * @param newVertVec the vector of the indexes of the snapped vertices + * + * Polyline vertices can be snapped either on vertexes or on edges. + * + */ + + void SnapPolyline(MeshType &poly) { - public: - EdgeGrid &edgeGrid; - MeshType &poly; - CoM &com; - QualitySign(EdgeGrid &_e,MeshType &_poly, CoM &_com):edgeGrid(_e),poly(_poly),com(_com) {}; - bool operator()(face::Pos ep) const - { - VertexType *v0 = ep.V(); - VertexType *v1 = ep.VFlip(); - if(v0->Q() * v1->Q() < 0) - { - Point3f pp = QLerp(v0,v1); - Point3f closestP; - if(com.MinDistOnEdge(pp,edgeGrid,poly,closestP) , Point3f> - { - EdgeGrid &edgeGrid; - MeshType &poly; - CoM &com; - vector &newVertVec; - - QualitySignSplit(EdgeGrid &_e,MeshType &_p, CoM &_com, vector &_vec):edgeGrid(_e),poly(_p),com(_com),newVertVec(_vec) {}; - void operator()(VertexType &nv, face::Pos ep) - { - VertexType *v0 = ep.V(); - VertexType *v1 = ep.VFlip(); - Point3f pp = QLerp(v0,v1); - Point3f closestP; - com.MinDistOnEdge(pp,edgeGrid,poly,closestP); - -// float minDist = MinDistOnEdge(v0,v1,edgeGrid,poly,closestP); - nv.P()=closestP; - nv.Q()=0; - newVertVec.push_back(tri::Index(com.base,&nv)); -// nv.P() = CoS::QLerp(v0,v1); - } - Color4b WedgeInterp(Color4b &c0, Color4b &c1) - { - Color4b cc; - cc.lerp(c0,c1,0.5f); - return Color4b::Red; - } - TexCoord2f WedgeInterp(TexCoord2f &t0, TexCoord2f &t1) - { - TexCoord2f tmp; - assert(t0.n()== t1.n()); - tmp.n()=t0.n(); - tmp.t()=(t0.t()+t1.t())/2.0; - return tmp; - } - }; - - class EdgePointPred - { - public: - CoM &com; - KdTree &kdtree; - - EdgePointPred(CoM &_com, KdTree &_kdtree):com(_com),kdtree(_kdtree) {}; - bool operator()(face::Pos ep) const - { - CoordType p0 = ep.V()->P(); - CoordType p1 = ep.VFlip()->P(); - float stepNum=100.0; - float locEps = Distance(p0,p1)/stepNum; - for(float j=0;j , Point3f> - { - CoM &com; - KdTree &kdtree; - MeshType &poly; - - EdgePointSplit(CoM &_com, KdTree &_kdtree, MeshType &_poly):com(_com),kdtree(_kdtree),poly(_poly) {}; - void operator()(VertexType &nv, face::Pos ep) - { - CoordType p0 = ep.V()->P(); - CoordType p1 = ep.VFlip()->P(); - float stepNum=100.0; - float locEps = Distance(p0,p1)/stepNum; - for(float j=0;j &planeVec, int i =0) - { - MeshType full; - for(int i=0;i::Mesh(full,t); - } - char buf[100]; - sprintf(buf,"planes%03i.ply",i); - tri::io::ExporterPLY::Save(full,buf); - } - - Plane3f ComputeEdgePlane(VertexType *v0, VertexType *v1) - { - Plane3f pl; - Point3f mid = (v0->P()+v1->P())/2.0; - Point3f delta = (v1->P()-v0->P()).Normalize(); - Point3f n0 = (v0->N() ^ delta).Normalize(); - Point3f n1 = (v1->N() ^ delta).Normalize(); - Point3f n = (n0+n1)/2.0; - pl.Init(mid,n); - return pl; - } - - - void ComputePlaneField(MeshType &poly, EdgeGrid &edgeGrid, int ind) - { - // First Compute per-edge planes - std::vector planeVec(poly.en); - for(int i=0;i p = vi->P(); - - float minDist=par.gridBailout; - Point3f closestP; - EdgeType *cep = vcg::tri::GetClosestEdgeBase(poly,edgeGrid,p,par.gridBailout,minDist,closestP); - if(cep) - { - int ind = tri::Index(poly,cep); - vi->Q() = SignedDistancePlanePoint(planeVec[ind],p); - Point3f proj = planeVec[ind].Projection(p); - -// if(Distance(proj,closestP)>edge::Length(*cep)) -// { -// vi->Q() =1; -// vi->SetS(); -// } - } - else { - vi->Q() =1; - } - } - } - - - void CutAlongPolyLineUsingField(MeshType &poly,EdgeGrid &edgeGrid,std::vector &newVertVec) -{ - QualitySign qsPred(edgeGrid,poly,*this); - QualitySignSplit qsSplit(edgeGrid,poly,*this,newVertVec); - tri::UpdateTopology::FaceFace(base); - tri::RefineE(base,qsSplit,qsPred); - tri::UpdateTopology::FaceFace(base); - - - for(int i=0;iIsD()) - { - for(int j=0;j<3;++j) - { - if(Distance(fp->P0(j),fp->P1(j)) < par.polyDistThr) - { - if(face::FFLinkCondition(*fp,j)) - { -// if(fp->V0(j)->Q()==0) fp->V1(j)->Q()=0; -// face::FFEdgeCollapse(base,*fp,j); - break; - } - } - } - } - } - tri::Allocator::CompactEveryVector(base); - - for(int i=0;iV(0)->Q()==0) && - (fp->V(1)->Q()==0) && - (fp->V(2)->Q()==0) ) - { - ScalarType maxDist = 0; - int maxInd = -1; - for(int j=0;j<3;++j) - { - Point3f closestPt; - ScalarType d = MinDistOnEdge(fp->P(j),edgeGrid,poly,closestPt); - if(d>maxDist) - { - maxDist= d; - maxInd=j; - } - } -// assert(maxInd!=-1); -// if(maxInd>=0 && maxDist > par.surfDistThr) -// fp->V(maxInd)->Q() = maxDist; - } - } - - for(int i=0;iV(0)->Q()>=0) && - (fp->V(1)->Q()>=0) && - (fp->V(2)->Q()>=0) ) - fp->C() = Color4b::Blue; - if( (fp->V(0)->Q()<=0) && - (fp->V(1)->Q()<=0) && - (fp->V(2)->Q()<=0) ) - fp->C() = Color4b::Red; - if( (fp->V(0)->Q()==0) && - (fp->V(1)->Q()==0) && - (fp->V(2)->Q()==0) ) - fp->C() = Color4b::Green; - - if( (fp->V(0)->Q()>0) && - (fp->V(1)->Q()>0) && - (fp->V(2)->Q()>0) ) - fp->C() = Color4b::White; - if( (fp->V(0)->Q()<0) && - (fp->V(1)->Q()<0) && - (fp->V(2)->Q()<0) ) - fp->C() = Color4b::White; - } - tri::AttributeSeam::SplitVertex(base, ExtractVertex, CompareVertex); - - } - - void WalkAlongPolyLine(MeshType &poly, std::vector &ptVec) - { - // Search a starting vertex - VertexType *startVert; - for(int i=0;iP()) < par.polyDistThr) - { - startVert = &base.vert[i]; - break; - } - } - tri::UpdateTopology::VertexFace(base); - tri::UpdateTopology::FaceFace(base); - - - - } - - - - -// } -// } - -/** - * - * - */ - - void CutWithPolyLine(MeshType &poly) - { - std::vector newVertVec; - SnapPolyline(poly, &newVertVec); - tri::io::ExporterPLY::Save(poly,"poly_snapped.ply",tri::io::Mask::IOM_EDGEINDEX+tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); - - DecomposeNonManifoldPolyline(poly); - tri::io::ExporterPLY::Save(poly,"poly_manif.ply",tri::io::Mask::IOM_EDGEINDEX+tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); - std::vector< std::vector< int> > ccVec; - BuildConnectedComponentVectors(poly,ccVec); - printf("PolyLine of %i edges decomposed into %i manifold components\n",poly.en,ccVec.size()); - Reorient(poly,ccVec); - char buf[1024]; - for(int i=0;i ptVec; - FindTerminalPoints(subPoly,ptVec); - printf("Component %i (%i edges) has %i terminal points\n",i,subPoly.en, ptVec.size());fflush(stdout); - SplitMeshWithPoints(base,ptVec,newVertVec); -// sprintf(buf,"CuttingPoly%02i.ply",i); -// tri::io::ExporterPLY::Save(subPoly, buf,tri::io::Mask::IOM_EDGEINDEX+tri::io::Mask::IOM_VERTCOLOR+tri::io::Mask::IOM_VERTQUALITY); - EdgeGrid edgeGrid; - ComputePlaneField(subPoly, edgeGrid,i); - sprintf(buf,"PlaneField%02i.ply",i); - tri::io::ExporterPLY::Save(base,buf,tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY ); - CutAlongPolyLineUsingField(subPoly,edgeGrid,newVertVec); - sprintf(buf,"PlaneCut%02i.ply",i); - tri::io::ExporterPLY::Save(base,buf,tri::io::Mask::IOM_FACECOLOR + tri::io::Mask::IOM_VERTQUALITY ); - } - -// printf("Added %i vertices\n",newVertVec.size()); -// for(int i=0;i::Save(base,"base_cut.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY ); - CleanSpuriousDanglingEdges(poly); - tri::io::ExporterPLY::Save(base,"base_cut_clean.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY ); - - } - - - - - void SnapPolyline(MeshType &poly, std::vector *newVertVec) - { - const float maxDist = base.bbox.Diag()/100.0; - const ScalarType interpEps = 0.0001; + tri::Allocator::CompactEveryVector(poly); + tri::UpdateTopology::VertexEdge(poly); int vertSnapCnt=0; int edgeSnapCnt=0; + int borderCnt=0,midCnt=0,nonmanifCnt=0; for(VertexIterator vi=poly.vert.begin(); vi!=poly.vert.end();++vi) { - float closestDist; - Point3f closestP,closestN,ip; - FaceType *f = vcg::tri::GetClosestFaceBase(base,uniformGrid,vi->P(),maxDist, closestDist, closestP, closestN,ip); - assert(f); - VertexType *closestVp=0; - int indIp = -1; - ScalarType minDist = std::numeric_limits::max(); - ScalarType minIp = minDist; - for(int i=0;i<3;++i) + Point3f ip; + FaceType *f = GetClosestFaceIP(vi->cP(),ip); + if(BarycentricSnap(ip)) { - if(Distance(vi->P(),f->P(i))P(),f->P(i)); - closestVp = f->V(i); - } - if(minIp > ip[i]) - { - indIp = i; - minIp=ip[i]; - } + if(ip[0]>0 && ip[1]>0) { vi->P() = f->P(0)*ip[0]+f->P(1)*ip[1]; edgeSnapCnt++; assert(ip[2]==0); vi->C()=Color4b::White;} + if(ip[0]>0 && ip[2]>0) { vi->P() = f->P(0)*ip[0]+f->P(2)*ip[2]; edgeSnapCnt++; assert(ip[1]==0); vi->C()=Color4b::White;} + if(ip[1]>0 && ip[2]>0) { vi->P() = f->P(1)*ip[1]+f->P(2)*ip[2]; edgeSnapCnt++; assert(ip[0]==0); vi->C()=Color4b::White;} + + if(ip[0]==1.0) { vi->P() = f->P(0); vertSnapCnt++; assert(ip[1]==0 && ip[2]==0); vi->C()=Color4b::Black; } + if(ip[1]==1.0) { vi->P() = f->P(1); vertSnapCnt++; assert(ip[0]==0 && ip[2]==0); vi->C()=Color4b::Black;} + if(ip[2]==1.0) { vi->P() = f->P(2); vertSnapCnt++; assert(ip[0]==0 && ip[1]==0); vi->C()=Color4b::Black;} } - assert(closestVp && (indIp!=-1)); - - - if(minDist < par.maxSnapThr) { // First Case: Snap to vertex; - vi->P() = closestVp->P(); - vertSnapCnt++; - if(newVertVec) - newVertVec->push_back(tri::Index(base,closestVp)); - } else { - if(minIp < interpEps) { // Second Case: Snap to edge; - ScalarType T1 = ip[(indIp+1)%3]; - ScalarType T2 = ip[(indIp+2)%3]; - vi->P() = (f->V1(indIp)->P() * T1 + f->V2(indIp)->P() * T2)/(T1+T2); - edgeSnapCnt++; - } + else + { + int deg = edge::VEDegree(&*vi); + if (deg > 2) { nonmanifCnt++; vi->C()=Color4b::Magenta; } + if (deg < 2) { borderCnt++; vi->C()=Color4b::Green;} + if (deg== 2) { midCnt++; vi->C()=Color4b::Blue;} } } - printf("Snapped %i onto vert and %i onto edges\n",vertSnapCnt, edgeSnapCnt); + printf("SnapPolyline %i vertices: snapped %i onto vert and %i onto edges %i nonmanif, %i border, %i mid\n", + poly.vn, vertSnapCnt, edgeSnapCnt, nonmanifCnt,borderCnt,midCnt); fflush(stdout); + int dupCnt=tri::Clean::RemoveDuplicateVertex(poly); + tri::Allocator::CompactEveryVector(poly); + if(dupCnt) printf("SnapPolyline: Removed %i Duplicated vertices\n",dupCnt); } - - + + /* * Make an edge mesh 1-manifold by splitting all the * vertexes that have more than two incident edges @@ -975,278 +382,190 @@ public: assert(firstVi == poly.vert.end()); } - /* - * Given a manifold edgemesh it returns the boundary/terminal endpoints. + + + + /** + * @brief SplitMeshWithPolyline + * @param poly + * + * First it splits the base mesh with all the non snapped points doing a standard 1 to 3 split; + * */ - static void FindTerminalPoints(MeshType &poly, std::vector &vec) - { - tri::UpdateTopology::VertexEdge(poly); + void SplitMeshWithPolyline(MeshType &poly) + { + std::vector< std::pair > toSplitVec; // the index of the face to be split and the poly vertex to be used + for(VertexIterator vi=poly.vert.begin(); vi!=poly.vert.end();++vi) { - if(edge::VEDegree(&*vi)==1) - vec.push_back(&*vi); - } - } - - // This function will decompose the input edge mesh into a set of - // connected components. - // the vector will contain, for each connected component, a vector with all the edge indexes. - void BuildConnectedComponentVectors(MeshType &poly, std::vector< std::vector< int> > &ccVec) - { - UpdateTopology::VertexEdge(poly); - for(size_t i=0;i(&(poly.vert[i])) <=2); + Point3f ip; + FaceType *f = GetClosestFaceIP(vi->cP(),ip); + if(!BarycentricSnap(ip)) + toSplitVec.push_back(std::make_pair(tri::Index(base,f),&*vi)); } - tri::UpdateTopology::EdgeEdge(poly); - tri::UpdateFlags::EdgeClearV(poly); + printf("SplitMeshWithPolyline found %i non snapped points\n",toSplitVec.size());fflush(stdout); - int visitedEdgeNum=0 ; - int ccCnt=0; - EdgeIterator eIt = poly.edge.begin(); + FaceIterator newFi = tri::Allocator::AddFaces(base,toSplitVec.size()*2); + VertexIterator newVi = tri::Allocator::AddVertices(base,toSplitVec.size()); + tri::UpdateColor::PerVertexConstant(base,Color4b::White); - while(visitedEdgeNum < poly.en) + for(size_t i =0; iIsV()) ++eIt; -// printf("Starting component from edge %i\n",tri::Index(poly,&*eIt)); - assert(eIt != poly.edge.end()); - edge::Pos startPos(&*eIt,0); - edge::Pos curPos(&*eIt,0); - do - { -// printf("(%i %i %i)-",tri::Index(poly,curPos.VFlip()), tri::Index(poly,curPos.E()) ,tri::Index(poly,curPos.V())); - curPos.NextE(); - } - while(curPos!=startPos && !curPos.IsBorder()) ; - - curPos.FlipV(); - assert(!curPos.IsBorder()); - do - { -// printf("<%i %i %i>-",tri::Index(poly,curPos.VFlip()), tri::Index(poly,curPos.E()) ,tri::Index(poly,curPos.V())); - curPos.E()->SetV(); - visitedEdgeNum++; - ccVec[ccCnt].push_back(tri::Index(poly,curPos.E())); - curPos.NextE(); - } while(!curPos.E()->IsV()); - printf("Completed component %i of %i edges\n",ccCnt, ccVec[ccCnt].size()); - ccCnt++; - } - } - - // This function will decompose the input edge mesh into a set of - // connected components. - // the vector will contain, for each connected component, a vector with all the edge indexes. - void BuildConnectedComponentVectorsOld(MeshType &poly, std::vector< std::vector< int> > &ccVec) - { - tri::UpdateTopology::EdgeEdge(poly); - tri::UpdateTopology::VertexEdge(poly); - tri::UpdateFlags::EdgeClearV(poly); - - int visitedEdgeNum=0 ; - int ccCnt=0; - - EdgeIterator eIt = poly.edge.begin(); - while(visitedEdgeNum < poly.en) - { - ccVec.resize(ccVec.size()+1); - while((eIt != poly.edge.end()) && eIt->IsV()) ++eIt; - EdgeType *startE = &*eIt; - - EdgeType *curEp = &*eIt; - int curEi = 0; - printf("Starting Visit of connected Component %i from edge %i\n",ccCnt,tri::Index(poly,*eIt)); - while( (curEp->EEp(curEi) != startE) && - (curEp->EEp(curEi) != curEp) ) - { - EdgeType *nextEp = curEp->EEp(curEi); - int nextEi = (curEp->EEi(curEi)+1)%2; - curEp = nextEp; - curEi = nextEi; - } - - curEp->SetV(); - curEi = (curEi +1)%2; // Flip the visit direction! - visitedEdgeNum++; - ccVec[ccCnt].push_back(tri::Index(poly,curEp)); - while(! curEp->EEp(curEi)->IsV()) - { - EdgeType *nextEp = curEp->EEp(curEi); - int nextEi = (curEp->EEi(curEi)+1)%2; - curEp->SetV(); - curEp = nextEp; - curEi = nextEi; - curEp->V(curEi)->C() = Color4b::Scatter(30,ccCnt%30); - if(!curEp->IsV()) { - ccVec[ccCnt].push_back(tri::Index(poly,curEp)); - visitedEdgeNum++; - } - } - printf("Completed visit of component of size %i\n",ccVec[ccCnt].size()); - ccCnt++; + newVi->P() = toSplitVec[i].second->P(); + newVi->C()=Color4b::Green; + face::TriSplit(&base.face[toSplitVec[i].first],&*(newFi++),&*(newFi++),&*(newVi++)); } - printf("en %i - VisitedEdgeNum = %i\n",poly.en, visitedEdgeNum); + Init(); // need to reset everthing + SnapPolyline(poly); - } - - void ExtractSubMesh(MeshType &poly, std::vector &ind, MeshType &subPoly) - { - subPoly.Clear(); - std::vector remap(poly.vert.size(),-1); - for(size_t i=0;i, VertexPointer> edgeToPolyVertMap; + for(VertexIterator vi=poly.vert.begin(); vi!=poly.vert.end();++vi) { - int v0 = tri::Index(poly,poly.edge[ind[i]].V(0)); - int v1 = tri::Index(poly,poly.edge[ind[i]].V(1)); - if(remap[v0]==-1) + Point3f ip; + FaceType *f = GetClosestFaceIP(vi->cP(),ip); + if(!BarycentricSnap(ip)) { assert(0); } + for(int i=0;i<3;++i) { - remap[v0]=subPoly.vn; - tri::Allocator::AddVertex(subPoly,poly.vert[v0].P(),poly.vert[v0].N()); + if(ip[i]>0 && ip[(i+1)%3]>0 && ip[(i+2)%3]==0 ) + { + CoordType p0=f->P0(i); + CoordType p1=f->P1(i); + if (p0>p1) std::swap(p0,p1); + if(edgeToPolyVertMap[make_pair(p0,p1)]) printf("Found an already used Edge %i - %i %i!!!\n", tri::Index(base,f->V0(i)),tri::Index(base,f->V1(i)),tri::Index(poly,&*vi)); + edgeToPolyVertMap[make_pair(p0,p1)]=&*vi; + } } - if(remap[v1]==-1) - { - remap[v1]=subPoly.vn; - tri::Allocator::AddVertex(subPoly,poly.vert[v1].P(),poly.vert[v1].N()); - } - tri::Allocator::AddEdge(subPoly, &subPoly.vert[remap[v0]], &subPoly.vert[remap[v1]]); - } - } - - // It takes a vector of vector of connected components and cohorently reorient each one of them. - // it usese the EE adjacency and requires that the input edgemesh is 1manifold. - - void Reorient(MeshType &poly, std::vector< std::vector< int> > &ccVec) - { - UpdateTopology::VertexEdge(poly); - for(size_t i=0;i(&(poly.vert[i])) <=2); } - UpdateTopology::EdgeEdge(poly); + printf("SplitMeshWithPolyline: Created a map of %i edges to be split\n",edgeToPolyVertMap.size()); + EdgePointPred ePred(edgeToPolyVertMap); + EdgePointSplit eSplit(edgeToPolyVertMap); + tri::UpdateTopology::FaceFace(base); + tri::RefineE(base,eSplit,ePred); + Init(); // need to reset everthing + } - for(size_t i=0;i toFlipVec(ccVec[i].size(),false); - - for(int j=0;jEEp(0) == cur) - prev = cur; // boundary - else - prev = & poly.edge[ccVec[i].back()]; // cc is a loop - } - else prev = & poly.edge[ccVec[i][j-1]]; - if(cur->EEp(0) != prev) - { - toFlipVec[j] = true; - assert(cur->EEp(1) == prev || j==0); - } - } - for(int j=0;j &vec, std::vector &newVertVec) - { - int faceToAdd=0; - int vertToAdd=0; - - // For each splitting point we save the index of the face to be splitten and the "kind" of split to do: - // 3 -> means classical 1 to 3 face split - // 2 -> means edge split. - // 0 -> means no need of a split (e.g. the point is coincident with a vertex) - - std::vector< std::pair > toSplitVec(vec.size(), std::make_pair(0,0)); - MeshGrid uniformGrid; - uniformGrid.Set(m.face.begin(), m.face.end()); - - for(size_t i =0; iP(); - float closestDist; - Point3f closestP; - FaceType *f = vcg::tri::GetClosestFaceBase(m,uniformGrid,newP,par.gridBailout, closestDist, closestP); - assert(f); - VertexType *closestVp=0; - ScalarType minDist = std::numeric_limits::max(); - for(int i=0;i<3;++i) { - if(Distance(newP,f->P(i))P(i)); - closestVp = f->V(i); - } - } - assert(closestVp); - if(minDist < par.maxSnapThr) { - vec[i]->P() = closestVp->P(); - } - else - { - toSplitVec[i].first = tri::Index(m,f); - toSplitVec[i].second = 3; - faceToAdd += 2; - vertToAdd += 1; - } - } -// printf("Splitting with %i points: adding %i faces and %i vertices\n",vec.size(), faceToAdd,vertToAdd); - FaceIterator newFi = tri::Allocator::AddFaces(m,faceToAdd); - VertexIterator newVi = tri::Allocator::AddVertices(m,vertToAdd); - - tri::UpdateColor::PerFaceConstant(m,Color4b::White); - - for(size_t i =0; iP() = vec[i]->P(); - VertexType *vp0 = fp0->V(0); - VertexType *vp1 = fp0->V(1); - VertexType *vp2 = fp0->V(2); - - fp0->V(0) = vp0; fp0->V(1) = vp1; fp0->V(2) = vp; - fp1->V(0) = vp1; fp1->V(1) = vp2; fp1->V(2) = vp; - fp2->V(0) = vp2; fp2->V(1) = vp0; fp2->V(2) = vp; - - fp0->C() = Color4b::Green; - fp1->C() = Color4b::Green; - fp2->C() = Color4b::Green; - } - } - } - void Init() { // Construction of the uniform grid - uniformGrid.Set(base.face.begin(), base.face.end()); UpdateNormal::PerFaceNormalized(base); UpdateTopology::FaceFace(base); uniformGrid.Set(base.face.begin(), base.face.end()); } + void SimplifyMidEdge(MeshType &poly) + { + int startVn; + int midEdgeCollapseCnt=0; + tri::Allocator::CompactEveryVector(poly); + do + { + startVn = poly.vn; + for(int ei =0; eiP(),ip0); + FaceType *f1=GetClosestFaceIP(v1->P(),ip1); + + bool snap0=BarycentricSnap(ip0); + bool snap1=BarycentricSnap(ip1); + int e0i,e1i; + bool e0 = IsSnappedEdge(ip0,e0i); + bool e1 = IsSnappedEdge(ip1,e1i); + if(e0 && e1) + if( ( f0 == f1 && e0i == e1i) || + ( f0 == f1->FFp(e1i) && e0i == f1->FFi(e1i)) || + (f0->FFp(e0i) == f1 && f0->FFi(e0i) == e1i) || + (f0->FFp(e0i) == f1->FFp(e1i) && f0->FFi(e0i) == f1->FFi(e1i)) ) + { + CoordType newp = (v0->P()+v1->P())/2.0; + v0->P()=newp; + v1->P()=newp; + midEdgeCollapseCnt++; + } + } + tri::Clean::RemoveDuplicateVertex(poly); + tri::Allocator::CompactEveryVector(poly); + printf("SimplifyMidEdge %5i -> %5i %i mid %i ve \n",startVn,poly.vn,midEdgeCollapseCnt); + } while(startVn>poly.vn); + } + + /** + * @brief SimplifyMidFace remove all the vertices that in the mid of a face and between two of the points snapped onto the edges of the same face + * @param poly + * + * It assumes that the mesh has been snapped and refined by the BaseMesh + */ + void SimplifyMidFace(MeshType &poly) + { + int startVn; + do + { + tri::Allocator::CompactEveryVector(poly); + startVn = poly.vn; + UpdateTopology::VertexEdge(poly); + int midFaceCollapseCnt=0; + int vertexEdgeCollapseCnt=0; + for(int i =0; i starVecVp; + edge::VVStarVE(&(poly.vert[i]),starVecVp); + if( (starVecVp.size()==2) ) + { + CoordType ipP, ipN, ipI; + FacePointer fpP = GetClosestFaceIP(starVecVp[0]->P(),ipP); + FacePointer fpN = GetClosestFaceIP(starVecVp[1]->P(),ipN); + FacePointer fpI = GetClosestFaceIP(poly.vert[i].P(), ipI); + + bool snapP = (BarycentricSnap(ipP)); + bool snapN = (BarycentricSnap(ipN)); + bool snapI = (BarycentricSnap(ipI)); + VertexPointer vertexSnapP = 0; + VertexPointer vertexSnapN = 0; + VertexPointer vertexSnapI = 0; + for(int j=0;j<3;++j) + { + if(ipP[j]==1.0) vertexSnapP=fpP->V(j); + if(ipN[j]==1.0) vertexSnapN=fpN->V(j); + if(ipI[j]==1.0) vertexSnapI=fpI->V(j); + } + + bool collapseFlag=false; + + if((!snapI && snapP && snapN) || // First case a vertex that is not snapped between two snapped vertexes + (!snapI && !snapP && fpI==fpP) || // Or a two vertex not snapped but on the same face + (!snapI && !snapN && fpI==fpN) ) + { + collapseFlag=true; + midFaceCollapseCnt++; + } + + else // case 2) a vertex snap and edge snap we have to check that the edge do not share the same vertex of the vertex snap + if(snapI && snapP && snapN && vertexSnapI==0 && (vertexSnapP!=0 || vertexSnapN!=0) ) + { + for(int j=0;j<3;++j) { + if(ipI[j]!=0 && (fpI->V(j)==vertexSnapP || fpI->V(j)==vertexSnapN)) { + collapseFlag=true; + vertexEdgeCollapseCnt++; + } + } + } + + if(collapseFlag) + edge::VEEdgeCollapse(poly,&(poly.vert[i])); + } + } + printf("SimplifyMidFace %5i -> %5i %i mid %i ve \n",startVn,poly.vn,midFaceCollapseCnt,vertexEdgeCollapseCnt); + } while(startVn>poly.vn); + } + void Simplify( MeshType &poly) { int startEn = poly.en; @@ -1270,10 +589,9 @@ public: Point3f fp,fn; float maxSurfDist = MaxSegDist(starVecVp[0], starVecVp[1],fp,fn); - if(maxSurfDist < par.surfDistThr && (newSegLen < par.maxSimpEdgeLen) ) + if((maxSurfDist < par.surfDistThr) && (newSegLen < par.maxSimpEdgeLen) ) { - edge::VEEdgeCollapse(poly,&(poly.vert[i])); - + edge::VEEdgeCollapse(poly,&(poly.vert[i])); } } } @@ -1303,9 +621,280 @@ public: tri::UpdateColor::PerVertexQualityRamp(poly,0,dist.Max()); } + + /** + * @brief BarycentricSnap + * @param ip the baricentric coords to be snapped + * @return true if they have been snapped. + * + * This is the VERY important function that is used everywhere. + * Given a barycentric coord of a point inside a triangle it decides if it should be "snapped" either onto an edge or on a vertex. * + * + */ + bool BarycentricSnap(CoordType &ip) + { + for(int i=0;i<3;++i) + { + if(ip[i] <= par.barycentricSnapThr) ip[i]=0; + if(ip[i] >= 1.0-par.barycentricSnapThr) ip[i]=1; + } + ScalarType sum = ip[0]+ip[1]+ip[2]; + + for(int i=0;i<3;++i) + if(ip[i]!=1) ip[i]/=sum; + + if(ip[0]==0 || ip[1]==0 || ip[2]==0) return true; + return false; + } + + + /** + * @brief TestSplitSegWithMesh Given a poly segment decide if it should be split along elements of base mesh. + * @param v0 + * @param v1 + * @param splitPt + * @return true if it should be split + * + * We make a few samples onto the edge and if some of them snaps onto a an edge we use it. + * In case there are more than one candidate we choose the sample closeset to its snapping point. + * We explicitly avoid snapping twice on the same edge by checking the starting and ending edges. + * + * Two cases: + * - poly edge pass near a vertex of the mesh + * - poly edge cross one or more edges + * + * Note that we have to check the case where + */ + bool TestSplitSegWithMesh(VertexType *v0, VertexType *v1, Point3f &splitPt) + { + Segment3f segPoly(v0->P(),v1->P()); + const float sampleNum = 40; + CoordType ip0,ip1; + + FaceType *f0=GetClosestFaceIP(v0->P(),ip0); + FaceType *f1=GetClosestFaceIP(v1->P(),ip1); + if(f0==f1) return false; + + bool snap0=false,snap1=false; // true if the segment start/end on a edge/vert + + Segment3f seg0; // The two segments to be avoided + Segment3f seg1; // from which the current poly segment can start + VertexPointer vertexSnap0 = 0; + VertexPointer vertexSnap1 = 0; + if(BarycentricSnap(ip0)) { + snap0=true; + for(int i=0;i<3;++i) { + if(ip0[i]==1.0) vertexSnap0=f0->V(i); + if(ip0[i]==0.0) seg0=Segment3f(f0->P1(i),f0->P2(i)); + } + } + if(BarycentricSnap(ip1)) { + snap1=true; + for(int i=0;i<3;++i){ + if(ip1[i]==1.0) vertexSnap1=f1->V(i); + if(ip1[i]==0.0) seg1=Segment3f(f1->P1(i),f1->P2(i)); + } + } + + CoordType bestSplitPt(0,0,0); + ScalarType bestDist = std::numeric_limits::max(); + for(float k = 1;kV(i); + CoordType closestPt = f->P(0)*ip[0]+f->P(1)*ip[1]+f->P(2)*ip[2]; + if(Distance(samplePnt,closestPt) < bestDist ) + { + ScalarType dist0=std::numeric_limits::max(); + ScalarType dist1=std::numeric_limits::max(); + CoordType closestSegPt; + if(snap0) SegmentPointDistance(seg0,closestPt,closestSegPt,dist0); + if(snap1) SegmentPointDistance(seg1,closestPt,closestSegPt,dist1); + if( (!vertexSnapI && (dist0 > par.surfDistThr/1000 && dist1>par.surfDistThr/1000) ) || + ( vertexSnapI!=vertexSnap0 && vertexSnapI!=vertexSnap1) ) + { + bestDist = Distance(samplePnt,closestPt); + bestSplitPt = closestPt; + } + } + } + } + if(bestDist < par.surfDistThr*100) + { + splitPt = bestSplitPt; + return true; + } + + return false; + } + /** + * @brief SnappedOnSameFace Return true if the two points are snapped to a common face; + * @param f0 + * @param i0 + * @param f1 + * @param i0 + * @return + * + * Require FFAdj. se assume that both SNAPPED. Three cases: + * - Edge Edge - true iff the two edges belongs to a common face. + * - Vert Edge - true iff there is one of the two snapped edge faces has the vert as non-edge face; + * - Vert Vert + * + */ + bool SnappedOnSameFace(FacePointer f0, CoordType i0, FacePointer f1, CoordType i1) + { + if(f0==f1) return true; + int e0,e1; + int v0,v1; + bool e0Snap = IsSnappedEdge(i0,e0); + bool e1Snap = IsSnappedEdge(i1,e1); + bool v0Snap = IsSnappedVertex(i0,v0); + bool v1Snap = IsSnappedVertex(i1,v1); + FacePointer f0p=0; int e0p=-1; // When Edge snap the other face and the index of the snapped edge on the other face + FacePointer f1p=0; int e1p=-1; + assert((e0Snap != v0Snap) && (e1Snap != v1Snap)); + // For EdgeSnap compute the 'other' face stuff + if(e0Snap){ + f0p = f0->FFp(e0); e0p=f0->FFi(e0); assert(f0p->FFp(e0p)==f0); + } + if(e1Snap){ + f1p = f1->FFp(e1); e1p=f1->FFi(e1); assert(f1p->FFp(e1p)==f1); + } + + if(e0Snap && e1Snap) { + if(f0==f1p || f0p==f1p || f0p==f1 || f0==f1) return true; + } + + if(e0Snap && v1Snap) { + assert(v1>=0 && v1<3 && v0==-1 && e1==-1); + if(f0->V2(e0) ==f1->V(v1)) return true; + if(f0p->V2(e0p)==f1->V(v1)) return true; + } + + if(e1Snap && v0Snap) { + assert(v0>=0 && v0<3 && v1==-1 && e0==-1); + if(f1->V2(e1) ==f0->V(v0)) return true; + if(f1p->V2(e1p)==f0->V(v0)) return true; + } + + if(v1Snap && v0Snap) { + PosType startPos(f0,f0->V(v0)); + PosType curPos=startPos; + do + { + if(curPos.VFlip()==f1->V(v1)) return true; + curPos.FlipE(); + curPos.FlipF(); + } + while(curPos!=startPos); + } + return false; + } + + /** + * @brief TestSplitSegWithMesh Given a poly segment decide if it should be split along elements of base mesh. + * @param v0 + * @param v1 + * @param splitPt + * @return true if it should be split + * + * We make a few samples onto the edge and if some of them snaps onto a an edge we use it. + * In case there are more than one candidate we choose the sample closeset to its snapping point. + * We explicitly avoid snapping twice on the same edge by checking the starting and ending edges. + * + * Two cases: + * - poly edge pass near a vertex of the mesh + * - poly edge cross one or more edges + * + * Note that we have to check the case where + */ + bool TestSplitSegWithMeshAdapt(VertexType *v0, VertexType *v1, Point3f &splitPt) + { + splitPt=(v0->P()+v1->P())/2.0; + + CoordType ip0,ip1,ipm; + FaceType *f0=GetClosestFaceIP(v0->P(),ip0); + FaceType *f1=GetClosestFaceIP(v1->P(),ip1); + FaceType *fm=GetClosestFaceIP(splitPt,ipm); + + if(f0==f1) return false; + + bool snap0=BarycentricSnap(ip0); + bool snap1=BarycentricSnap(ip1); + bool snapm=BarycentricSnap(ipm); + + if(!snap0 && !snap1) { + assert(f0!=f1); + return true; + } + if(snap0 && snap1) + { + if(SnappedOnSameFace(f0,ip0,f1,ip1)) + return false; + + } + + if(snap0) { + int e0,v0; + if (IsSnappedEdge(ip0,e0)) { + if(f0->FFp(e0) == f1) return false; + } + if(IsSnappedVertex(ip0,v0)) { + for(int i=0;i<3;++i) + if(f1->V(i)==f0->V(v0)) return false; + } + } + return true; + } + + + bool TestSplitSegWithMeshAdaptOld(VertexType *v0, VertexType *v1, Point3f &splitPt) + { + Segment3f segPoly(v0->P(),v1->P()); + const float sampleNum = 40; + CoordType ip0,ip1; + FaceType *f0=GetClosestFaceIP(v0->P(),ip0); + FaceType *f1=GetClosestFaceIP(v1->P(),ip1); + if(f0==f1) return false; + + bool snap0=BarycentricSnap(ip0); + bool snap1=BarycentricSnap(ip1); + + if(!snap0 && !snap1) { + assert(f0!=f1); + splitPt=(v0->P()+v1->P())/2.0; + return true; + } + if(snap0 && snap1) + { + if(SnappedOnSameFace(f0,ip0,f1,ip1)) + return false; + } + + if(snap0) { + int e0,v0; + if (IsSnappedEdge(ip0,e0)) { + if(f0->FFp(e0) == f1) return false; + } + if(IsSnappedVertex(ip0,v0)) { + for(int i=0;i<3;++i) + if(f1->V(i)==f0->V(v0)) return false; + } + } + splitPt=(v0->P()+v1->P())/2.0; + return true; + } // Given a segment find the maximum distance from it to the original surface. + // It is used to evaluate the Haustdorff distance of a Segment from the mesh. float MaxSegDist(VertexType *v0, VertexType *v1, Point3f &farthestPointOnSurf, Point3f &farthestN, Distribution *dist=0) { float maxSurfDist = 0; @@ -1330,7 +919,19 @@ public: return maxSurfDist; } - void Refine(MeshType &poly, bool uniformFlag = false) + + /** + * @brief RefineCurve + * @param poly the curve to be refined + * @param uniformFlag + * + * Make one pass of refinement for all the edges of the curve that are distant from the basemesh + * uses two parameters: + * - par.minRefEdgeLen + * - par.surfDistThr + */ + + void RefineCurveByDistance(MeshType &poly) { tri::Allocator::CompactEveryVector(poly); int startEdgeSize = poly.en; @@ -1345,70 +946,36 @@ public: { edge::VEEdgeSplit(poly, &ei, farthestP, farthestN); } - else if(uniformFlag) - { - edge::VEEdgeSplit(poly,&ei,(ei.P(0)+ei.P(1))/2.0,(ei.V(0)->N()+ei.V(1)->N())/2.0); - } } } // tri::Allocator::CompactEveryVector(poly); printf("Refine %i -> %i\n",startEdgeSize,poly.en);fflush(stdout); } - // Take a poly and find the position of edge -edge crossing. - void RefineBaseMesh(MeshType &poly) + /** + * @brief RefineCurveByBaseMesh + * @param poly + */ + + void RefineCurveByBaseMesh(MeshType &poly) { - std::vector newPtVec; - for(int i=0;i::CompactEveryVector(poly); + int lastEn; + do { + lastEn = poly.en; + for(int i =0; iFFp(ei) == lastFace) - { - Point3f ps0,ps1; - float segDist; - bool par; - Segment3f faceSeg(fp->P0(ei),fp->P1(ei)); - float angDeg = fabs(math::ToDeg(Angle(faceSeg.Direction(),seg.Direction()))); - - SegmentSegmentDistance(seg,faceSeg,segDist,par,ps0,ps1); - if(!par && (angDeg > 5 && angDeg<175)) newPtVec.push_back(ps1); - } - } - } - lastFace=fp; - lastqp = qp; - } - } - } - MeshType meshPt; tri::BuildMeshFromCoordVector(meshPt,newPtVec); - tri::io::ExporterPLY::Save(meshPt,"newPoints.ply"); - - VertexConstDataWrapper vdw(meshPt); - KdTree kdtree(vdw); - - EdgePointPred epPred(*this,kdtree); - EdgePointSplit epSplit(*this,kdtree,meshPt); - tri::UpdateTopology::FaceFace(base); - tri::RefineE(base,epSplit,epPred); - tri::UpdateTopology::FaceFace(base); - tri::io::ExporterPLY::Save(base,"newbase.ply"); + CoordType splitPt; + if(TestSplitSegWithMeshAdapt(poly.edge[i].V(0),poly.edge[i].V(1),splitPt)) + edge::VEEdgeSplit(poly, &poly.edge[i], splitPt); + } + printf("RefineCurveByBaseMesh %i en -> %i en\n",lastEn,poly.en); fflush(stdout); + tri::Allocator::CompactEveryVector(poly); + } while(lastEn < poly.en); + SimplifyMidFace(poly); + SimplifyMidEdge(poly); + SnapPolyline(poly); } @@ -1434,6 +1001,8 @@ public: assert(poly.en>0 && base.fn>0); for(int k=0;k posVec(poly.vn,Point3f(0,0,0)); std::vector cntVec(poly.vn,0); @@ -1447,7 +1016,7 @@ public: } } - const float maxDist = base.bbox.Diag()/10.0; + const float maxDist = base.bbox.Diag()/10.0; for(int i=0; iN(); } -// Refine(poly); + // Refine(poly); tri::UpdateTopology::TestVertexEdge(poly); - Refine(poly); + RefineCurveByDistance(poly); tri::UpdateTopology::TestVertexEdge(poly); - Simplify(poly); + Simplify(poly); tri::UpdateTopology::TestVertexEdge(poly); int dupVertNum = Clean::RemoveDuplicateVertex(poly); if(dupVertNum) { @@ -1485,7 +1054,56 @@ public: } +class EdgePointPred +{ +public: + std::map, VertexPointer> &edgeToPolyVertMap; + + EdgePointPred(std::map, VertexPointer> &_edgeToPolyVertMap):edgeToPolyVertMap(_edgeToPolyVertMap){}; + bool operator()(face::Pos ep) const + { + CoordType p0 = ep.V()->P(); + CoordType p1 = ep.VFlip()->P(); + if (p0>p1) std::swap(p0,p1); + VertexPointer vp=edgeToPolyVertMap[make_pair(p0,p1)]; + return vp!=0; + } }; + +struct EdgePointSplit : public std::unary_function , Point3f> +{ +public: + std::map, VertexPointer> &edgeToPolyVertMap; + + EdgePointSplit(std::map, VertexPointer> &_edgeToPolyVertMap):edgeToPolyVertMap(_edgeToPolyVertMap){}; + void operator()(VertexType &nv, face::Pos ep) + { + CoordType p0 = ep.V()->P(); + CoordType p1 = ep.VFlip()->P(); + if (p0>p1) std::swap(p0,p1); + VertexPointer vp=edgeToPolyVertMap[make_pair(p0,p1)]; + assert(vp); + nv.P()=vp->P(); + return; + } + Color4b WedgeInterp(Color4b &c0, Color4b &c1) + { + Color4b cc; + cc.lerp(c0,c1,0.5f); + return Color4b::Red; + } + TexCoord2f WedgeInterp(TexCoord2f &t0, TexCoord2f &t1) + { + TexCoord2f tmp; + assert(t0.n()== t1.n()); + tmp.n()=t0.n(); + tmp.t()=(t0.t()+t1.t())/2.0; + return tmp; + } +}; + +}; + } // end namespace tri } // end namespace vcg diff --git a/vcg/complex/algorithms/cut_tree.h b/vcg/complex/algorithms/cut_tree.h new file mode 100644 index 00000000..258ed4b3 --- /dev/null +++ b/vcg/complex/algorithms/cut_tree.h @@ -0,0 +1,238 @@ +/**************************************************************************** +* VCGLib o o * +* Visual and Computer Graphics Library o o * +* _ O _ * +* Copyright(C) 2004-2016 \/)\/ * +* Visual Computing Lab /\/| * +* ISTI - Italian National Research Council | * +* \ * +* All rights reserved. * +* * +* This program is free software; you can redistribute it and/or modify * +* it under the terms of the GNU General Public License as published by * +* the Free Software Foundation; either version 2 of the License, or * +* (at your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, * +* but WITHOUT ANY WARRANTY; without even the implied warranty of * +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * +* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * +* for more details. * +* * +****************************************************************************/ +#ifndef CUT_TREE_H +#define CUT_TREE_H + +namespace vcg { +namespace tri { + +template +class CutTree +{ +public: + typedef typename MeshType::ScalarType ScalarType; + typedef typename MeshType::CoordType CoordType; + typedef typename MeshType::VertexType VertexType; + typedef typename MeshType::VertexPointer VertexPointer; + typedef typename MeshType::VertexIterator VertexIterator; + typedef typename MeshType::EdgeIterator EdgeIterator; + typedef typename MeshType::EdgeType EdgeType; + typedef typename MeshType::FaceType FaceType; + typedef typename MeshType::FacePointer FacePointer; + typedef typename MeshType::FaceIterator FaceIterator; + typedef Box3 Box3Type; + typedef typename vcg::GridStaticPtr MeshGrid; + typedef typename vcg::GridStaticPtr EdgeGrid; + typedef typename face::Pos PosType; + typedef typename tri::UpdateTopology::PEdge PEdge; + + MeshType &base; +// MeshGrid uniformGrid; + +// Param par; + CutTree(MeshType &_m) :base(_m){} + + +void OptimizeTree(MeshType &t) +{ + tri::Allocator::CompactEveryVector(t); + int lastEn=t.en; + do + { + lastEn=t.en; + tri::UpdateTopology::VertexEdge(t); + tri::UpdateTopology::VertexFace(base); + VertexConstDataWrapper vdw(base); + KdTree kdtree(vdw); + + // First simple loop that search for 2->1 moves. + for(VertexIterator vi=t.vert.begin();vi!=t.vert.end();++vi) + { + std::vector starVec; + edge::VVStarVE(&*vi,starVec); + if(starVec.size()==2) + { + PosType pos; + if(ExistEdge(kdtree,starVec[0]->P(),starVec[1]->P(),pos)) + edge::VEEdgeCollapse(t,&*vi); + } + } + tri::Allocator::CompactEveryVector(t); + } + while(t.en &kdtree, CoordType &p0, CoordType &p1, PosType &fpos) +{ + ScalarType locEps = SquaredDistance(p0,p1)/100000.0; + + VertexType *v0=0,*v1=0; + unsigned int veInd; + ScalarType sqdist; + kdtree.doQueryClosest(p0,veInd,sqdist); + if(sqdist vfi(v0); + while(!vfi.End()) + { + if(vfi.V1()==v1) + { + fpos = PosType(vfi.F(),vfi.I(), v0); + return true; + } + if(vfi.V2()==v1) + { + fpos = PosType(vfi.F(),(vfi.I()+1)%3, v1); + return true; + } + ++vfi; + } + } + return false; +} + + +int findNonVisitedEdgesDuringRetract(VertexType * vp, EdgeType * &ep) +{ + std::vector starVec; + edge::VEStarVE(&*vp,starVec); + int cnt =0; + for(size_t i=0;iIsV()) { + cnt++; + ep = starVec[i]; + } + } + return cnt; +} + +void Retract(MeshType &t) +{ + printf("Retracting a tree of %i edges and %i vertices\n",t.en,t.vn); + tri::UpdateTopology::VertexEdge(t); + + std::stack vertStack; + + // Put on the stack all the vertex with just a single incident edge. + for(VertexIterator vi=t.vert.begin();vi!=t.vert.end();++vi) + { + std::vector starVec; + edge::VEStarVE(&*vi,starVec); + if(starVec.size()==1) + vertStack.push(&*vi); + } + + tri::UpdateFlags::EdgeClearV(t); + tri::UpdateFlags::VertexClearV(t); + + int unvisitedEdgeNum = t.en; + while((!vertStack.empty()) && (unvisitedEdgeNum > 2) ) + { + VertexType *vp = vertStack.top(); + vertStack.pop(); + vp->C()=Color4b::Blue; + EdgeType *ep=0; + int eCnt = findNonVisitedEdgesDuringRetract(vp,ep); + if(eCnt==1) // We have only one non visited edge over vp + { + assert(!ep->IsV()); + ep->SetV(); + --unvisitedEdgeNum; + VertexType *otherVertP; + if(ep->V(0)==vp) otherVertP = ep->V(1); + else otherVertP = ep->V(0); + vertStack.push(otherVertP); + } + } + assert(unvisitedEdgeNum >0); + + for(size_t i =0; i::DeleteEdge(t,t.edge[i]); + assert(t.en >0); + tri::Clean::RemoveUnreferencedVertex(t); + tri::Allocator::CompactEveryVector(t); +} + + +// \brief This function build a cut tree. +// +// First we make a bread first FF face visit. +// Each time that we encounter a visited face we cut add to the tree the edge +// that brings to the already visited face. +// this structure build a dense graph and we retract this graph retracting each +// leaf until we remains with just the loops that cuts the object. + +void BuildVisitTree(MeshType &dualMesh, int startingFaceInd=0) +{ + tri::UpdateTopology::FaceFace(base); + tri::UpdateFlags::FaceClearV(base); + std::vector > visitStack; // the stack contain the pos on the 'starting' face. + + base.face[startingFaceInd].SetV(); + for(int i=0;i<3;++i) + visitStack.push_back(PosType(&(base.face[startingFaceInd]),i,base.face[startingFaceInd].V(i))); + + int cnt=1; + + while(!visitStack.empty()) + { + std::swap(visitStack.back(),visitStack[rand()%visitStack.size()]); + PosType c = visitStack.back(); + visitStack.pop_back(); + assert(c.F()->IsV()); + c.F()->C() = Color4b::ColorRamp(0,base.fn,cnt); + c.FlipF(); + if(!c.F()->IsV()) + { + ++cnt; + c.F()->SetV(); + c.FlipE();c.FlipV(); + visitStack.push_back(c); + c.FlipE();c.FlipV(); + visitStack.push_back(c); + } + else + { + tri::Allocator::AddEdge(dualMesh,c.V()->P(),c.VFlip()->P()); + } + } + assert(cnt==base.fn); + + tri::Clean::RemoveDuplicateVertex(dualMesh); + Retract(dualMesh); + OptimizeTree(dualMesh); + tri::UpdateBounding::Box(dualMesh); +} + +}; +} // end namespace tri +} // end namespace vcg + +#endif // CUT_TREE_H