/**************************************************************************** * MeshLab o o * * A versatile mesh processing toolbox o o * * _ O _ * * Copyright(C) 2005 \/)\/ * * 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 VORONOI_PROCESSING_H #define VORONOI_PROCESSING_H #include #include #include #include #include #include namespace vcg { namespace tri { struct VoronoiProcessingParameter { enum { None=0, DistanceFromSeed=1, DistanceFromBorder=2, RegionArea=3 }; VoronoiProcessingParameter() { colorStrategy = DistanceFromSeed; areaThresholdPerc=0; deleteUnreachedRegionFlag=false; constrainSelectedSeed=false; preserveFixedSeed=false; collapseShortEdge=false; collapseShortEdgePerc = 0.01f; triangulateRegion=false; unbiasedSeedFlag = true; geodesicRelaxFlag = true; relaxOnlyConstrainedFlag=false; refinementRatio = 5.0f; seedPerturbationProbability=0; seedPerturbationAmount = 0.001f; } int colorStrategy; float areaThresholdPerc; bool deleteUnreachedRegionFlag; bool unbiasedSeedFlag; bool constrainSelectedSeed; /// If true the selected vertexes define a constraining domain: /// During relaxation all selected seeds are constrained to move /// only on other selected vertices. /// In this way you can constrain some seed to move only on certain /// domains, for example moving only along some linear features /// like border of creases. bool relaxOnlyConstrainedFlag; bool preserveFixedSeed; /// If true the 'fixed' seeds are not moved during relaxation. /// \see FixVertexVector function to see how to fix a set of seeds. float refinementRatio; /// It defines how much the input mesh has to be refined in order to have a supporting /// triangulation that is dense enough to well approximate the voronoi diagram. /// reasonable values are in the range 4..10. It is used by PreprocessForVoronoi and this value /// says how many triangles you should expect in a voronoi region of a given radius. float seedPerturbationProbability; /// if true at each iteration step each seed has the given probability to be perturbed a little. float seedPerturbationAmount; /// As a bbox diag fraction (e.g. in the 0..1 range). // Convertion to Voronoi Diagram Parameters bool triangulateRegion; /// If true when building the voronoi diagram mesh each region is a /// triangulated polygon. Otherwise it each voronoi region is a star /// triangulation with the original seed in the center. bool collapseShortEdge; float collapseShortEdgePerc; bool geodesicRelaxFlag; }; template > class VoronoiProcessing { typedef typename MeshType::CoordType CoordType; typedef typename MeshType::ScalarType ScalarType; typedef typename MeshType::VertexType VertexType; typedef typename MeshType::VertexPointer VertexPointer; typedef typename MeshType::VertexIterator VertexIterator; typedef typename MeshType::FacePointer FacePointer; typedef typename MeshType::FaceIterator FaceIterator; typedef typename MeshType::FaceType FaceType; typedef typename MeshType::FaceContainer FaceContainer; typedef typename tri::Geodesic::VertDist VertDist; static math::MarsenneTwisterRNG &RandomGenerator() { static math::MarsenneTwisterRNG rnd; return rnd; } public: typedef typename MeshType::template PerVertexAttributeHandle PerVertexPointerHandle; typedef typename MeshType::template PerVertexAttributeHandle PerVertexBoolHandle; typedef typename MeshType::template PerVertexAttributeHandle PerVertexFloatHandle; typedef typename MeshType::template PerFaceAttributeHandle PerFacePointerHandle; // Given a vector of point3f it finds the closest vertices on the mesh. static void SeedToVertexConversion(MeshType &m,std::vector &seedPVec,std::vector &seedVVec, bool compactFlag = true) { typedef typename vcg::SpatialHashTable HashVertexGrid; seedVVec.clear(); HashVertexGrid HG; HG.Set(m.vert.begin(),m.vert.end()); const float dist_upper_bound=m.bbox.Diag()/10.0; typename std::vector::iterator pi; for(pi=seedPVec.begin();pi!=seedPVec.end();++pi) { ScalarType dist; VertexPointer vp; vp=tri::GetClosestVertex(m, HG, *pi, dist_upper_bound, dist); if(vp) { seedVVec.push_back(vp); } } if(compactFlag) { std::sort(seedVVec.begin(),seedVVec.end()); typename std::vector::iterator vi = std::unique(seedVVec.begin(),seedVVec.end()); seedVVec.resize( std::distance(seedVVec.begin(),vi) ); } } static void ComputePerVertexSources(MeshType &m, std::vector &seedVec, DistanceFunctor &df) { tri::Allocator::DeletePerVertexAttribute(m,"sources"); // delete any conflicting handle regardless of the type... PerVertexPointerHandle vertexSources = tri::Allocator:: template AddPerVertexAttribute (m,"sources"); tri::Allocator::DeletePerFaceAttribute(m,"sources"); // delete any conflicting handle regardless of the type... tri::Allocator::template AddPerFaceAttribute (m,"sources"); assert(tri::Allocator::IsValidHandle(m,vertexSources)); tri::Geodesic::Compute(m,seedVec,df,std::numeric_limits::max(),0,&vertexSources); } static void VoronoiColoring(MeshType &m, bool frontierFlag=true) { PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); assert(tri::Allocator::IsValidHandle(m,sources)); if(frontierFlag) { //static_cast(NULL) has been introduced just to avoid an error in the MSVS2010's compiler confusing pointer with int. You could use nullptr to avoid it, but it's not supported by all compilers. //The error should have been removed from MSVS2012 std::pair zz(0.0f,static_cast(NULL)); std::vector< std::pair > regionArea(m.vert.size(),zz); std::vector frontierVec; GetAreaAndFrontier(m, sources, regionArea, frontierVec); tri::Geodesic::Compute(m,frontierVec); } tri::UpdateColor::PerVertexQualityRamp(m); } static void VoronoiAreaColoring(MeshType &m,std::vector &seedVec, std::vector< std::pair > ®ionArea) { PerVertexPointerHandle vertexSources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); float meshArea = tri::Stat::ComputeMeshArea(m); float expectedArea = meshArea/float(seedVec.size()); for(size_t i=0;i:: template GetPerFaceAttribute (m,"sources"); PerVertexPointerHandle vertexSources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) { faceSources[fi]=0; std::vector vp(3); for(int i=0;i<3;++i) vp[i]=vertexSources[fi->V(i)]; for(int i=0;i<3;++i) // First try to associate to the most reached vertex { if(vp[0]==vp[1] && vp[0]==vp[2]) faceSources[fi] = vp[0]; else { if(vp[0]==vp[1] && vp[0]->Q()< vp[2]->Q()) faceSources[fi] = vp[0]; if(vp[0]==vp[2] && vp[0]->Q()< vp[1]->Q()) faceSources[fi] = vp[0]; if(vp[1]==vp[2] && vp[1]->Q()< vp[0]->Q()) faceSources[fi] = vp[1]; } } } tri::UpdateTopology::FaceFace(m); int unassCnt=0; do { unassCnt=0; for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) { if(faceSources[fi]==0) { std::vector vp(3); for(int i=0;i<3;++i) vp[i]=faceSources[fi->FFp(i)]; if(vp[0]!=0 && (vp[0]==vp[1] || vp[0]==vp[2])) faceSources[fi] = vp[0]; else if(vp[1]!=0 && (vp[1]==vp[2])) faceSources[fi] = vp[1]; else faceSources[fi] = std::max(vp[0],std::max(vp[1],vp[2])); if(faceSources[fi]==0) unassCnt++; } } } while(unassCnt>0); } // Select all the faces with a given source vertex // It reads the PerFace attribute 'sources' static int FaceSelectAssociateRegion(MeshType &m, VertexPointer vp) { PerFacePointerHandle sources = tri::Allocator:: template FindPerFaceAttribute (m,"sources"); assert(tri::Allocator::IsValidHandle(m,sources)); tri::UpdateSelection::Clear(m); int selCnt=0; for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) { if(sources[fi]==vp) { fi->SetS(); ++selCnt; } } return selCnt; } // Given a seed , it selects all the faces that have the minimal distance vertex sourced by the given . // can be null (it search for unreached faces...) // returns the number of selected faces; // // It reads the PerVertex attribute 'sources' static int FaceSelectRegion(MeshType &m, VertexPointer vp) { PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); assert(tri::Allocator::IsValidHandle(m,sources)); tri::UpdateSelection::Clear(m); int selCnt=0; for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) { int minInd = 0; float minVal=std::numeric_limits::max(); for(int i=0;i<3;++i) { if((*fi).V(i)->Q()Q(); } } if( sources[(*fi).V(minInd)] == vp) { fi->SetS(); selCnt++; } } return selCnt; } /// Given a mesh with for each vertex the link to the closest seed /// (e.g. for all vertexes we know what is the corresponding voronoi region) /// we compute: /// area of all the voronoi regions /// the vector of the frontier vertexes (e.g. vert of faces shared by at least two regions) /// /// Area is computed only for triangles that fully belong to a given source. static void GetAreaAndFrontier(MeshType &m, PerVertexPointerHandle &sources, std::vector< std::pair > ®ionArea, // for each seed we store area std::vector &frontierVec) { tri::UpdateFlags::VertexClearV(m); frontierVec.clear(); for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) { VertexPointer s0 = sources[(*fi).V(0)]; VertexPointer s1 = sources[(*fi).V(1)]; VertexPointer s2 = sources[(*fi).V(2)]; assert(s0 && s1 && s2); if((s0 != s1) || (s0 != s2) ) { for(int i=0;i<3;++i) if(!fi->V(i)->IsV()) { frontierVec.push_back(fi->V(i)); fi->V(i)->SetV(); } } else // the face belongs to a single region; accumulate area; { if(s0 != 0) { int seedIndex = tri::Index(m,s0); regionArea[seedIndex].first+=DoubleArea(*fi)*0.5f; regionArea[seedIndex].second=s0; } } } } /// Given a mesh with for each vertex the link to the closest seed /// we compute: /// the vector of the corner faces (ie the faces shared exactly by three regions) /// the vector of the frontier faces that are on the boundary. static void GetFaceCornerVec(MeshType &m, PerVertexPointerHandle &sources, std::vector &cornerVec, std::vector &borderCornerVec) { tri::UpdateFlags::VertexClearV(m); cornerVec.clear(); for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) { VertexPointer s0 = sources[(*fi).V(0)]; VertexPointer s1 = sources[(*fi).V(1)]; VertexPointer s2 = sources[(*fi).V(2)]; assert(s0 && s1 && s2); if(s1!=s2 && s0!=s1 && s0!=s2) { cornerVec.push_back(&*fi); } else { if(isBorderCorner(&*fi,sources)) borderCornerVec.push_back(&*fi); } } } static bool isBorderCorner(FaceType *f, typename MeshType::template PerVertexAttributeHandle &sources) { for(int i=0;i<3;++i) { if(sources[(*f).V0(i)] != sources[(*f).V1(i)] && f->IsB(i)) return true; } return false; } // Given two supposedly adjacent border corner faces it finds the source common to them; static VertexPointer CommonSourceBetweenBorderCorner(FacePointer f0, FacePointer f1, typename MeshType::template PerVertexAttributeHandle &sources) { assert(isBorderCorner(f0,sources)); assert(isBorderCorner(f1,sources)); int b0 =-1,b1=-1; for(int i=0;i<3;++i) { if(face::IsBorder(*f0,i)) b0=i; if(face::IsBorder(*f1,i)) b1=i; } assert(b0!=-1 && b1!=-1); if( (sources[f0->V0(b0)] == sources[f1->V0(b1)]) || (sources[f0->V0(b0)] == sources[f1->V1(b1)]) ) return sources[f0->V0(b0)]; if( (sources[f0->V1(b0)] == sources[f1->V0(b1)]) || (sources[f0->V1(b0)] == sources[f1->V1(b1)]) ) return sources[f0->V1(b0)]; assert(0); return 0; } static void ConvertVoronoiDiagramToMesh(MeshType &m, MeshType &outMesh, MeshType &outPoly, std::vector &seedVec, VoronoiProcessingParameter &vpp ) { tri::RequirePerVertexAttribute(m,"sources"); PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); outMesh.Clear(); outPoly.Clear(); tri::UpdateTopology::FaceFace(m); tri::UpdateFlags::FaceBorderFromFF(m); std::vector innerCornerVec, // Faces adjacent to three different regions borderCornerVec; // Faces that are on the border and adjacent to at least two regions. GetFaceCornerVec(m, sources, innerCornerVec, borderCornerVec); // For each seed collect all the vertices and build for(size_t i=0;i::AddVertex(outMesh,seedVec[i]->P(),Color4b::DarkGray); for(size_t i=0;i pt; for(size_t j=0;jV(qq)] == curSeed) { pt.push_back(Barycenter(*innerCornerVec[j])); break; } for(size_t j=0;jV(qq)] == curSeed) { CoordType edgeCenter; for(int jj=0;jj<3;++jj) if(face::IsBorder(*(borderCornerVec[j]),jj)) edgeCenter=(borderCornerVec[j]->P0(jj)+borderCornerVec[j]->P1(jj))/2.0f; pt.push_back(edgeCenter); break; } Plane3 pl; pt.push_back(curSeed->P()); FitPlaneToPointSet(pt,pl); pt.pop_back(); CoordType nZ = pl.Direction(); CoordType nX = (pt[0]-curSeed->P()).Normalize(); CoordType nY = (nX^nZ).Normalize(); vector > angleVec(pt.size()); for(size_t j=0;jP()).Normalize(); float angle = 180.0f+math::ToDeg(atan2(p*nY,p*nX)); angleVec[j] = make_pair(angle,j); } std::sort(angleVec.begin(),angleVec.end()); // Now build another piece of mesh. int curRegionStart=outMesh.vert.size(); for(size_t j=0;j::AddVertex(outMesh,pt[angleVec[j].second],Color4b::LightGray); for(size_t j=0;j::AddFace(outMesh, &outMesh.vert[i ], &outMesh.vert[curRegionStart + j ], &outMesh.vert[curRegionStart + ((j+1)%pt.size())]); outMesh.face.back().SetF(0); outMesh.face.back().SetF(2); } } // end for each seed. tri::Clean::RemoveDuplicateVertex(outMesh); tri::UpdateTopology::FaceFace(outMesh); bool oriented,orientable; tri::Clean::OrientCoherentlyMesh(outMesh,oriented,orientable); tri::UpdateTopology::FaceFace(outMesh); // last loop to remove faux edges bit that are now on the boundary. for(FaceIterator fi=outMesh.face.begin();fi!=outMesh.face.end();++fi) for(int i=0;i<3;++i) if(face::IsBorder(*fi,i) && fi->IsF(i)) fi->ClearF(i); std::vector< typename tri::UpdateTopology::PEdge> EdgeVec; // ******************* star to tri conversion ********* // If requested the voronoi regions are converted from a star arragned polygon // with vertex on the seed to a simple triangulated polygon by mean of a simple edge collapse if(vpp.triangulateRegion) { tri::UpdateFlags::FaceBorderFromFF(outMesh); tri::UpdateFlags::VertexBorderFromFaceBorder(outMesh); for(FaceIterator fi=outMesh.face.begin();fi!=outMesh.face.end();++fi) if(!fi->IsD()) { for(int i=0;i<3;++i) { bool b0 = fi->V0(i)->IsB(); bool b1 = fi->V1(i)->IsB(); if( ((b0 && b1) || (fi->IsF(i) && !b0) ) && tri::Index(outMesh,fi->V0(i))V0(i))]->IsS()) if(face::FFLinkCondition(*fi, i)) { face::FFEdgeCollapse(outMesh, *fi,i); // we delete vertex fi->V0(i) break; } } } } } // Now a plain conversion of the non faux edges into a polygonal mesh tri::UpdateTopology::FillUniqueEdgeVector(outMesh,EdgeVec,false); tri::UpdateTopology::AllocateEdge(outMesh); for(size_t i=0;i::AddVertex(outPoly,outMesh.vert[i].P()); for(size_t i=0;i::AddEdge(outPoly,&(outPoly.vert[e0]),&(outPoly.vert[e1])); } } /// \brief Build a mesh of voronoi diagram from the given seeds /// /// This function assumes that you have just run a geodesic like algorithm over your mesh using /// a seed set as starting points and that there is an PerVertex Attribute called 'sources' /// with pointers to the seed source. Usually you can initialize it with something like /// /// DistanceFunctor &df, /// tri::Geodesic::Compute(m, seedVec, df, std::numeric_limits::max(),0,&sources); /// static void ConvertVoronoiDiagramToMeshOld(MeshType &m, MeshType &outMesh, MeshType &outPoly, std::vector &seedVec, VoronoiProcessingParameter &vpp ) { tri::RequirePerVertexAttribute(m,"sources"); PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); outMesh.Clear(); outPoly.Clear(); tri::UpdateTopology::FaceFace(m); tri::UpdateFlags::FaceBorderFromFF(m); std::map seedMap; // It says if a given vertex of m is a seed (and what position it has in the seed vector) for(size_t i=0;i=0) && (ind innerCornerVec, // Faces adjacent to three different regions borderCornerVec; // Faces that are on the border and adjacent to at least two regions. GetFaceCornerVec(m, sources, innerCornerVec, borderCornerVec); std::map vertexIndCornerMap; // Given a cornerFace (border or inner) what is the corresponding vertex? for(size_t i=0;i::AddVertex(outMesh, seedVec[i]->P(),Color4b::White); for(size_t i=0;i::AddVertex(outMesh, vcg::Barycenter(*(innerCornerVec[i])),Color4b::Gray); vertexIndCornerMap[innerCornerVec[i]] = outMesh.vn-1; } for(size_t i=0;iP0(j)+borderCornerVec[i]->P1(j))/2.0f; tri::Allocator::AddVertex(outMesh, edgeCenter,Color4b::Gray); vertexIndCornerMap[borderCornerVec[i]] = outMesh.vn-1; } tri::Append::MeshCopy(outPoly,outMesh); // There is a voronoi edge if there are two corner face that share two sources. // In such a case we add a pair of triangles with an edge connecting these two corner faces // and with the two involved sources // For each pair of adjacent seed we store the first of the two corner that we encounter. std::map, FacePointer > VoronoiEdge; // 1) Build internal triangles // Loop build all the triangles connecting seeds with internal corners // we loop over the all the voronoi corner (triangles with three different sources) // we build for(size_t i=0;iV0(j)]; VertexPointer v1 = sources[innerCornerVec[i]->V1(j)]; assert(seedMap[v0]>=0);assert(seedMap[v1]>=0); if(v1::AddFace(outMesh,&(outMesh.vert[seedMap[v0]]), corner0, corner1); tri::Allocator::AddFace(outMesh,&(outMesh.vert[seedMap[v1]]), corner1, corner0); } } } // 2) build the boundary facets: // We loop over border corners and build triangles with seed vertex // we do **only** triangles with a bordercorner and a internal 'corner' for(size_t i=0;iV(0)]; // All bordercorner faces have only two different regions VertexPointer s1 = sources[borderCornerVec[i]->V(1)]; if(s1==s0) s1 = sources[borderCornerVec[i]->V(2)]; if(s1::AddFace(outMesh,&(outMesh.vert[seedMap[s0]]), corner0, corner1); tri::Allocator::AddFace(outMesh,&(outMesh.vert[seedMap[s1]]), corner0, corner1); } } // Final pass tri::UpdateFlags::FaceClearV(m); bool AllFaceVisited = false; while(!AllFaceVisited) { // search for a unvisited boundary face face::Pos pos,startPos; AllFaceVisited=true; for(size_t i=0; (AllFaceVisited) && (iIsV()) { for(int j=0;j<3;++j) if(face::IsBorder(*(borderCornerVec[i]),j)) { pos.Set(borderCornerVec[i],j,borderCornerVec[i]->V(j)); AllFaceVisited =false; } } if(AllFaceVisited) break; assert(pos.IsBorder()); startPos=pos; bool foundBorderSeed=false; FacePointer curBorderCorner = pos.F(); do { pos.F()->SetV(); pos.NextB(); if(sources[pos.V()]==pos.V()) foundBorderSeed=true; assert(isBorderCorner(curBorderCorner,sources)); if(isBorderCorner(pos.F(),sources)) if(pos.F() != curBorderCorner) { VertexPointer curReg = CommonSourceBetweenBorderCorner(curBorderCorner, pos.F(),sources); VertexPointer curSeed = &(outMesh.vert[seedMap[curReg]]); int otherCorner0 = vertexIndCornerMap[pos.F() ]; int otherCorner1 = vertexIndCornerMap[curBorderCorner]; VertexPointer corner0 = &(outMesh.vert[otherCorner0]); VertexPointer corner1 = &(outMesh.vert[otherCorner1]); if(!foundBorderSeed) tri::Allocator::AddFace(outMesh,curSeed,corner0,corner1); foundBorderSeed=false; curBorderCorner=pos.F(); } } while(pos!=startPos); } //**************** CLEANING *************** // 1) reorient bool oriented,orientable; tri::UpdateTopology::FaceFace(outMesh); tri::Clean::OrientCoherentlyMesh(outMesh,oriented,orientable); // assert(orientable); // check that the normal of the input mesh are consistent with the result tri::UpdateNormal::PerVertexNormalizedPerFaceNormalized(outMesh); tri::UpdateNormal::PerVertexNormalizedPerFaceNormalized(m); if(seedVec[0]->N() * outMesh.vert[0].N() < 0 ) tri::Clean::FlipMesh(outMesh); tri::UpdateTopology::FaceFace(outMesh); tri::UpdateFlags::FaceBorderFromFF(outMesh); // 2) Remove Flips tri::UpdateNormal::PerFaceNormalized(outMesh); tri::UpdateFlags::FaceClearV(outMesh); for(FaceIterator fi=outMesh.face.begin();fi!=outMesh.face.end();++fi) { int badDiedralCnt=0; for(int i=0;i<3;++i) if(fi->N() * fi->FFp(i)->N() <0 ) badDiedralCnt++; if(badDiedralCnt == 2) fi->SetV(); } for(FaceIterator fi=outMesh.face.begin();fi!=outMesh.face.end();++fi) if(fi->IsV()) Allocator::DeleteFace(outMesh,*fi); tri::Allocator::CompactEveryVector(outMesh); tri::UpdateTopology::FaceFace(outMesh); tri::UpdateFlags::FaceBorderFromFF(outMesh); tri::UpdateFlags::VertexBorderFromFaceBorder(outMesh); // 3) set up faux bits for(FaceIterator fi=outMesh.face.begin();fi!=outMesh.face.end();++fi) for(int i=0;i<3;++i) { size_t v0 = tri::Index(outMesh,fi->V0(i) ); size_t v1 = tri::Index(outMesh,fi->V1(i) ); if (v0 < seedVec.size() && !(seedVec[v0]->IsB() && fi->IsB(i))) fi->SetF(i); if (v1 < seedVec.size() && !(seedVec[v1]->IsB() && fi->IsB(i))) fi->SetF(i); } if(vpp.collapseShortEdge) { float distThr = m.bbox.Diag() * vpp.collapseShortEdgePerc; for(FaceIterator fi=outMesh.face.begin();fi!=outMesh.face.end();++fi) if(!fi->IsD()) { for(int i=0;i<3;++i) if((Distance(fi->P0(i),fi->P1(i))IsF(i)) { // printf("Collapsing face %i:%i e%i \n",tri::Index(outMesh,*fi),tri::Index(outMesh,fi->FFp(i)),i); if ((!fi->V(i)->IsB())&&(face::FFLinkCondition(*fi, i))) face::FFEdgeCollapse(outMesh, *fi,i); break; } } } //******************** END OF CLEANING **************** // ******************* star to tri conversion ********* // If requested the voronoi regions are converted from a star arragned polygon // with vertex on the seed to a simple triangulated polygon by mean of a simple edge collapse if(vpp.triangulateRegion) { for(FaceIterator fi=outMesh.face.begin();fi!=outMesh.face.end();++fi) if(!fi->IsD()) { for(int i=0;i<3;++i) { bool b0 = fi->V0(i)->IsB(); bool b1 = fi->V1(i)->IsB(); if( ((b0 && b1) || (fi->IsF(i) && !b0 && !b1) ) && tri::Index(outMesh,fi->V(i))V(i))]->IsS()) if(face::FFLinkCondition(*fi, i)) { face::FFEdgeCollapse(outMesh, *fi,i); break; } } } } } // Now a plain conversion of the non faux edges into a polygonal mesh std::vector< typename tri::UpdateTopology::PEdge> EdgeVec; tri::UpdateTopology::FillUniqueEdgeVector(outMesh,EdgeVec,false); tri::UpdateTopology::AllocateEdge(outMesh); for(size_t i=0;i::AddEdge(outPoly,&(outPoly.vert[e0]),&(outPoly.vert[e1])); } } class VoronoiEdge { public: VertexPointer r0,r1; FacePointer f0,f1; bool operator == (const VoronoiEdge &ve) const {return ve.r0==r0 && ve.r1==r1; } bool operator < (const VoronoiEdge &ve) const { return (ve.r0==r0)?ve.r1 &edgeVec) { PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); edgeVec.clear(); std::vector cornerVec; std::vector borderCornerVec; GetFaceCornerVec(m,sources,cornerVec,borderCornerVec); // Now find all the voronoi edges: each edge (a *face pair) is identified by two voronoi regions typedef std::map< std::pair, std::pair > EdgeMapType; EdgeMapType EdgeMap; printf("cornerVec.size() %i\n",(int)cornerVec.size()); for(size_t i=0;iV0(j)]; VertexPointer v1 = sources[cornerVec[i]->V1(j)]; assert(v0!=v1); if(v0>v1) std::swap(v1,v0); std::pair adjRegion = std::make_pair(v0,v1); if(EdgeMap[adjRegion].first==0) EdgeMap[adjRegion].first = cornerVec[i]; else EdgeMap[adjRegion].second = cornerVec[i]; } } for(size_t i=0;iV(0)]; VertexPointer v1 = sources[borderCornerVec[i]->V(1)]; if(v0==v1) v1 = sources[borderCornerVec[i]->V(2)]; assert(v0!=v1); if(v0>v1) std::swap(v1,v0); std::pair adjRegion = std::make_pair(v0,v1); if(EdgeMap[adjRegion].first==0) EdgeMap[adjRegion].first = borderCornerVec[i]; else EdgeMap[adjRegion].second = borderCornerVec[i]; } typename EdgeMapType::iterator mi; for(mi=EdgeMap.begin();mi!=EdgeMap.end();++mi) { if((*mi).second.first && (*mi).second.second) { assert((*mi).first.first && (*mi).first.second); edgeVec.push_back(VoronoiEdge()); edgeVec.back().r0 = (*mi).first.first; edgeVec.back().r1 = (*mi).first.second; edgeVec.back().f0 = (*mi).second.first; edgeVec.back().f1 = (*mi).second.second; } } } static void BuildBiasedSeedVec(MeshType &m, DistanceFunctor &df, std::vector &seedVec, std::vector &frontierVec, std::vector &biasedFrontierVec, VoronoiProcessingParameter &vpp) { (void)df; biasedFrontierVec.clear(); if(vpp.unbiasedSeedFlag) { for(size_t i=0;i edgeVec; BuildVoronoiEdgeVec(m,edgeVec); printf("Found %lu edges on a diagram of %lu seeds\n",edgeVec.size(),seedVec.size()); std::map > SeedToEdgeVecMap; std::map< std::pair, VoronoiEdge *> SeedPairToEdgeMap; float totalLen=0; for(size_t i=0;i regionPerymeter; for(size_t i=0;iLen(); } printf("perimeter of region %i is %f\n",(int)i,regionPerymeter[seedVec[i]]); } PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); // The real bias for each edge is (perim)/(edge) // each source can belong to two edges max. so the weight is std::map weight; std::map cnt; float biasSum = totalLen/5.0f; for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) { for(int i=0;i<3;++i) { VertexPointer s0 = sources[(*fi).V0(i)]; VertexPointer s1 = sources[(*fi).V1(i)]; if(s0!=s1) { if(s0>s1) std::swap(s0,s1); VoronoiEdge *ve = SeedPairToEdgeMap[std::make_pair(s0,s1)]; if(!ve) printf("v %i %i \n",(int)tri::Index(m,s0),(int)tri::Index(m,s1)); assert(ve); float el = ve->Len(); weight[(*fi).V0(i)] += (regionPerymeter[s0]+biasSum)/(el+biasSum) ; weight[(*fi).V1(i)] += (regionPerymeter[s1]+biasSum)/(el+biasSum) ; cnt[(*fi).V0(i)]++; cnt[(*fi).V1(i)]++; } } } for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) { if(cnt[&*vi]>0) { // float bias = weight[&*vi]/float(cnt[&*vi]); float bias = weight[&*vi]/float(cnt[&*vi]) + totalLen; biasedFrontierVec.push_back(VertDist(&*vi, bias)); } } printf("Collected %i frontier vertexes\n",(int)biasedFrontierVec.size()); } static void DeleteUnreachedRegions(MeshType &m, PerVertexPointerHandle &sources) { tri::UpdateFlags::VertexClearV(m); for(size_t i=0;iV(0)->IsV() || fi->V(1)->IsV() || fi->V(2)->IsV() ) { face::VFDetach(*fi); tri::Allocator::DeleteFace(m,*fi); } // qDebug("Deleted faces not reached: %i -> %i",int(m.face.size()),m.fn); tri::Clean::RemoveUnreferencedVertex(m); tri::Allocator::CompactEveryVector(m); } /// Let f_p(q) be the squared distance of q from p /// f_p(q) = (p_x-q_x)^2 + (p_y-q_y)^2 + (p_z-q_z)^2 /// f_p(q) = p_x^2 -2p_xq_x +q_x^2 + ... + p_z^2 -2p_zq_z +q_z^2 /// struct QuadricSumDistance { ScalarType a; ScalarType c; CoordType b; QuadricSumDistance() {a=0; c=0; b[0]=0; b[1]=0; b[2]=0;} void AddPoint(CoordType p) { a+=1; assert(c>=0); c+=p*p; b[0]+= -2.0f*p[0]; b[1]+= -2.0f*p[1]; b[2]+= -2.0f*p[2]; } ScalarType Eval(CoordType p) const { ScalarType d = a*(p*p) + b*p + c; assert(d>=0); return d; } CoordType Min() const { return b * -0.5f; } }; /// \brief Relax the seeds of a Voronoi diagram according to the quadric distance rule. /// /// For each region it search the vertex that minimize the sum of the squared distance /// from all the points of the region. /// /// It uses a vector of QuadricSumDistances; /// for simplicity it is sized as the vertex vector even if only the ones of the quadric /// corresponding to seeds are actually used. /// /// It return true if at least one seed changed position. /// static bool QuadricRelax(MeshType &m, std::vector &seedVec, std::vector &frontierVec, std::vector &newSeeds, DistanceFunctor &df, VoronoiProcessingParameter &vpp) { (void)seedVec; (void)frontierVec; (void)df; newSeeds.clear(); PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); PerVertexBoolHandle fixed = tri::Allocator:: template GetPerVertexAttribute (m,"fixed"); QuadricSumDistance dz; std::vector dVec(m.vert.size(),dz); for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) { assert(sources[vi]!=0); int seedIndex = tri::Index(m,sources[vi]); // When constraining seeds movement we move selected seeds only onto other selected vertices if(vpp.constrainSelectedSeed) { // So we sum only the contribs of the selected vertices if( (sources[vi]->IsS() && vi->IsS()) || (!sources[vi]->IsS())) dVec[seedIndex].AddPoint(vi->P()); } else dVec[seedIndex].AddPoint(vi->P()); } // Search the local maxima for each region and use them as new seeds std::pair zz(std::numeric_limits::max(), static_cast(0)); std::vector< std::pair > seedMaximaVec(m.vert.size(),zz); for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) { assert(sources[vi]!=0); int seedIndex = tri::Index(m,sources[vi]); ScalarType val = dVec[seedIndex].Eval(vi->P()); vi->Q()=val; // if constrainSelectedSeed we search only among selected vertices if(!vpp.constrainSelectedSeed || !sources[vi]->IsS() || vi->IsS()) { if(seedMaximaVec[seedIndex].first > val) { seedMaximaVec[seedIndex].first = val; seedMaximaVec[seedIndex].second = &*vi; } } } if(vpp.colorStrategy==VoronoiProcessingParameter::DistanceFromBorder) tri::UpdateColor::PerVertexQualityRamp(m); // tri::io::ExporterPLY::Save(m,"last.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY ); bool seedChanged=false; // update the seedvector with the new maxima (For the vertex not fixed) for(size_t i=0;i &seedVec, std::vector &frontierVec, std::vector &newSeeds, DistanceFunctor &df, VoronoiProcessingParameter &vpp) { newSeeds.clear(); typename MeshType::template PerVertexAttributeHandle sources; sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); typename MeshType::template PerVertexAttributeHandle fixed; fixed = tri::Allocator:: template GetPerVertexAttribute (m,"fixed"); std::vector::VertDist> biasedFrontierVec; BuildBiasedSeedVec(m,df,seedVec,frontierVec,biasedFrontierVec,vpp); tri::Geodesic::Visit(m,biasedFrontierVec,df); if(vpp.colorStrategy == VoronoiProcessingParameter::DistanceFromSeed) tri::UpdateColor::PerVertexQualityRamp(m); // tri::io::ExporterPLY::Save(m,"last.ply",tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY ); if(vpp.colorStrategy == VoronoiProcessingParameter::DistanceFromBorder) tri::UpdateColor::PerVertexQualityRamp(m); // Search the local maxima for each region and use them as new seeds std::pair zz(0.0f,static_cast(NULL)); std::vector< std::pair > seedMaximaVec(m.vert.size(),zz); for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) { assert(sources[vi]!=0); int seedIndex = tri::Index(m,sources[vi]); if(!vpp.constrainSelectedSeed || !sources[vi]->IsS() || vi->IsS()) { if(seedMaximaVec[seedIndex].first < (*vi).Q()) { seedMaximaVec[seedIndex].first=(*vi).Q(); seedMaximaVec[seedIndex].second=&*vi; } } } bool seedChanged=false; // update the seedvector with the new maxima (For the vertex not selected) for(size_t i=0;i &seedVec, std::vector< std::pair > ®ionArea, VoronoiProcessingParameter &vpp) { // Smaller area region are discarded Distribution H; for(size_t i=0;i newSeedVec; // update the seedvector with the new maxima (For the vertex not selected) for(size_t i=0;i= areaThreshold) newSeedVec.push_back(seedVec[i]); } swap(seedVec,newSeedVec); } /// \brief Mark a vector of seeds to be fixed. /// /// Vertex pointers must belong to the mesh. /// The framework use a boolean attribute called "fixed" to store this info. /// static void FixVertexVector(MeshType &m, std::vector &vertToFixVec) { typename MeshType::template PerVertexAttributeHandle fixed; fixed = tri::Allocator:: template GetPerVertexAttribute (m,"fixed"); for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) fixed[vi]=false; for(size_t i=0;i &seedPosVec, std::vector &fixedVec, int relaxStep, VoronoiProcessingParameter &vpp, vcg::CallBackPos *cb=0) { PerVertexFloatHandle area = tri::Allocator:: template GetPerVertexAttribute (m,"area"); for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) area[vi]=0; for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi) { ScalarType a3 = DoubleArea(*fi)/6.0; for(int i=0;i<3;++i) area[fi->V(i)]+=a3; } assert(m.vn > seedPosVec.size()*20); int i; ScalarType perturb = m.bbox.Diag()*vpp.seedPerturbationAmount; for(i=0;i > vdw(seedPosVec); KdTree seedTree(vdw); std::vector > sumVec(seedPosVec.size(),std::make_pair(0,CoordType(0,0,0))); for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) { unsigned int seedInd; ScalarType sqdist; seedTree.doQueryClosest(vi->P(),seedInd,sqdist); vi->Q()=sqrt(sqdist); sumVec[seedInd].first+=area[vi]; sumVec[seedInd].second+=vi->cP()*area[vi]; } vector newseedVec; vector newfixedVec; for(size_t i=0;i RandomGenerator().generate01()) newseedVec.back()+=math::GeneratePointInUnitBallUniform( RandomGenerator())*perturb; newfixedVec.push_back(false); } } } std::swap(seedPosVec,newseedVec); std::swap(fixedVec,newfixedVec); tri::UpdateColor::PerVertexQualityRamp(m); } return relaxStep; } /// \brief Perform a Lloyd relaxation cycle over a mesh /// It uses two conventions: /// 1) a few vertexes can remain fixed, you have to set a per vertex bool attribute named 'fixed' /// 2) /// static int VoronoiRelaxing(MeshType &m, std::vector &seedVec, int relaxIter, DistanceFunctor &df, VoronoiProcessingParameter &vpp, vcg::CallBackPos *cb=0) { tri::RequireVFAdjacency(m); tri::RequireCompactness(m); for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi) assert(vi->VFp() && "Require mesh without unreferenced vertexes\n"); std::vector selectedVec; if(vpp.relaxOnlyConstrainedFlag) { for(size_t i=0;iIsS()) selectedVec.push_back(seedVec[i]); std::swap(seedVec,selectedVec); } tri::UpdateFlags::FaceBorderFromVF(m); tri::UpdateFlags::VertexBorderFromFaceBorder(m); PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); PerVertexBoolHandle fixed = tri::Allocator:: template GetPerVertexAttribute (m,"fixed"); int iter; for(iter=0;iter::Compute(m, seedVec, df,std::numeric_limits::max(),0,&sources); if(vpp.colorStrategy == VoronoiProcessingParameter::DistanceFromSeed) tri::UpdateColor::PerVertexQualityRamp(m); // Delete all the (hopefully) small regions that have not been reached by the seeds; if(vpp.deleteUnreachedRegionFlag) DeleteUnreachedRegions(m,sources); std::pair zz(0.0f,static_cast(NULL)); std::vector< std::pair > regionArea(m.vert.size(),zz); std::vector frontierVec; GetAreaAndFrontier(m, sources, regionArea, frontierVec); assert(frontierVec.size()>0); if(vpp.colorStrategy == VoronoiProcessingParameter::RegionArea) VoronoiAreaColoring(m, seedVec, regionArea); // qDebug("We have found %i regions range (%f %f), avg area is %f, Variance is %f 10perc is %f",(int)seedVec.size(),H.Min(),H.Max(),H.Avg(),H.StandardDeviation(),areaThreshold); if(cb) cb(iter*100/relaxIter,"Voronoi Lloyd Relaxation: Searching New Seeds"); std::vector newSeedVec; bool changed; if(vpp.geodesicRelaxFlag) changed = GeodesicRelax(m,seedVec, frontierVec, newSeedVec, df,vpp); else changed = QuadricRelax(m,seedVec,frontierVec, newSeedVec, df,vpp); assert(newSeedVec.size() == seedVec.size()); PruneSeedByRegionArea(newSeedVec,regionArea,vpp); for(size_t i=0;iC() = Color4b::Gray; for(size_t i=0;iC() = Color4b::Black; for(size_t i=0;iC() = Color4b::White; swap(newSeedVec,seedVec); if(!changed) break; } // Last run: Needed if we have changed the seed set to leave the sources handle correct. if(iter==relaxIter) tri::Geodesic::Compute(m, seedVec, df,std::numeric_limits::max(),0,&sources); if(vpp.relaxOnlyConstrainedFlag) { std::swap(seedVec,selectedVec); size_t i,j; for(i=0,j=0;iIsS()) { seedVec[i]=selectedVec[j]; fixed[seedVec[i]]=true; ++j; } } } return iter; } // Base vertex voronoi coloring algorithm. // it assumes VF adjacency. No attempt of computing real geodesic distnace is done. Just a BFS visit starting from the seeds. static void TopologicalVertexColoring(MeshType &m, std::vector &seedVec) { std::queue VQ; tri::UpdateQuality::VertexConstant(m,0); for(size_t i=0;iQ()=i+1; } while(!VQ.empty()) { VertexPointer vp = VQ.front(); VQ.pop(); std::vector vertStar; vcg::face::VVStarVF(vp,vertStar); for(typename std::vector::iterator vv = vertStar.begin();vv!=vertStar.end();++vv) { if((*vv)->Q()==0) { (*vv)->Q()=vp->Q(); VQ.push(*vv); } } } // end while(!VQ.empty()) } template static std::pair ordered_pair(const genericType &a, const genericType &b) { if(a, VertexPointer > &midMap) { PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); for(FaceIterator fi = m.face.begin(); fi!=m.face.end(); ++fi) { VertexPointer vp[3],sp[3]; vp[0] = (*fi).V(0); vp[1] = (*fi).V(1); vp[2] = (*fi).V(2); sp[0] = sources[vp[0]]; sp[1] = sources[vp[1]]; sp[2] = sources[vp[2]]; if((sp[0] == sp[1]) && (sp[0] == sp[2])) continue; // skip internal faces // if((sp[0] != sp[1]) && (sp[0] != sp[2]) && (sp[1] != sp[2])) continue; // skip corner faces for(int i=0;i<3;++i) // for each edge of a frontier face { int i0 = i; int i1 = (i+1)%3; // if((sp[i0]->IsS() && sp[i1]->IsS()) && !( vp[i0]->IsS() || vp[i1]->IsS() ) ) continue; VertexPointer closestVert = vp[i0]; if( vp[i1]->Q() < closestVert->Q()) closestVert = vp[i1]; if(sp[i0]->IsS() && sp[i1]->IsS()) { if ( (vp[i0]->IsS()) && !(vp[i1]->IsS()) ) closestVert = vp[i0]; if (!(vp[i0]->IsS()) && (vp[i1]->IsS()) ) closestVert = vp[i1]; if ( (vp[i0]->IsS()) && (vp[i1]->IsS()) ) closestVert = (vp[i0]->Q() < vp[i1]->Q()) ? vp[i0]:vp[i1]; } if(midMap[ordered_pair(sp[i0],sp[i1])] == 0 ) { midMap[ordered_pair(sp[i0],sp[i1])] = closestVert; } else { if(sp[i0]->IsS() && sp[i1]->IsS()) // constrained edge { if(!(midMap[ordered_pair(sp[i0],sp[i1])]->IsS()) && closestVert->IsS()) midMap[ordered_pair(sp[i0],sp[i1])] = closestVert; if( midMap[ordered_pair(sp[i0],sp[i1])]->IsS() && closestVert->IsS() && closestVert->Q() < midMap[ordered_pair(sp[i0],sp[i1])]->Q()) { midMap[ordered_pair(sp[i0],sp[i1])] = closestVert; } } else // UNCOSTRAINED EDGE { if(closestVert->Q() < midMap[ordered_pair(sp[i0],sp[i1])]->Q()) midMap[ordered_pair(sp[i0],sp[i1])] = closestVert; } } } } } /// \brief Check the topological correcteness of the induced Voronoi diagram /// /// This function assumes that you have just run a geodesic like algorithm over your mesh using /// a seed set as starting points and that there is an PerVertex Attribute called 'sources' /// with pointers to the seed source. Usually you can initialize it with something like /// /// DistanceFunctor &df, /// tri::Geodesic::Compute(m, seedVec, df, std::numeric_limits::max(),0,&sources); static bool CheckVoronoiTopology(MeshType& m,std::vector &seedVec) { tri::RequirePerVertexAttribute(m,"sources"); tri::RequireCompactness(m); typename MeshType::template PerVertexAttributeHandle sources; sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); std::map seedMap; // It says if a given vertex of m is a seed (and its index in seedVec) BuildSeedMap(m,seedVec,seedMap); // Very basic check: each vertex must have a source that is a seed. for(int i=0;i regionVec(seedVec.size(),0); for(size_t i=0; i< seedVec.size();i++) regionVec[i] = new MeshType; for(int i=0;i=0 && vi1>=0 && vi2>=0); tri::Allocator::AddFace(*regionVec[vi0], m.face[i].cP(0),m.face[i].cP(1),m.face[i].cP(2)); if(vi1 != vi0) tri::Allocator::AddFace(*regionVec[vi1], m.face[i].cP(0),m.face[i].cP(1),m.face[i].cP(2)); if((vi2 != vi0) && (vi2 != vi1) ) tri::Allocator::AddFace(*regionVec[vi2], m.face[i].cP(0),m.face[i].cP(1),m.face[i].cP(2)); } bool AllDiskRegion=true; for(size_t i=0; i< seedVec.size();i++) { MeshType &rm = *(regionVec[i]); tri::Clean::RemoveDuplicateVertex(rm); tri::Allocator::CompactEveryVector(rm); tri::UpdateTopology::FaceFace(rm); // char buf[100]; sprintf(buf,"disk%04i.ply",i); tri::io::ExporterPLY::Save(rm,buf,tri::io::Mask::IOM_VERTCOLOR + tri::io::Mask::IOM_VERTQUALITY ); int NNmanifoldE=tri::Clean::CountNonManifoldEdgeFF(rm); if (NNmanifoldE!=0) AllDiskRegion= false; int G=tri::Clean::MeshGenus(rm); int numholes=tri::Clean::CountHoles(rm); if (numholes!=1) AllDiskRegion= false; if(G!=0) AllDiskRegion= false; delete regionVec[i]; } if(!AllDiskRegion) return false; // **** Final step build a rough delaunay tri and check that it is manifold MeshType delaMesh; std::vector innerCornerVec, // Faces adjacent to three different regions borderCornerVec; // Faces that are on the border and adjacent to at least two regions. GetFaceCornerVec(m, sources, innerCornerVec, borderCornerVec); // First add all the needed vertices: seeds and corners for(size_t i=0;i::AddVertex(delaMesh, seedVec[i]->P()); // Now just add one face for each inner corner for(size_t i=0;iV(0)]]]; VertexPointer v1 = & delaMesh.vert[seedMap[sources[innerCornerVec[i]->V(1)]]]; VertexPointer v2 = & delaMesh.vert[seedMap[sources[innerCornerVec[i]->V(2)]]]; tri::Allocator::AddFace(delaMesh,v0,v1,v2); } Clean::RemoveUnreferencedVertex(delaMesh); tri::Allocator::CompactVertexVector(delaMesh); tri::UpdateTopology::FaceFace(delaMesh); int nonManif = tri::Clean::CountNonManifoldEdgeFF(delaMesh); if(nonManif>0) return false; return true; } static void BuildSeedMap(MeshType &m, std::vector &seedVec, std::map &seedMap) { seedMap.clear(); for(size_t i=0;i=0 && tri::Index(m,seedVec[i])::Compute(m, seedVec, df, std::numeric_limits::max(),0,&sources); /// /// The function can also static void ConvertDelaunayTriangulationToMesh(MeshType &m, MeshType &outMesh, std::vector &seedVec, bool refineFlag=true) { tri::RequirePerVertexAttribute(m,"sources"); tri::RequireCompactness(m); tri::RequireVFAdjacency(m); PerVertexPointerHandle sources = tri::Allocator:: template GetPerVertexAttribute (m,"sources"); outMesh.Clear(); tri::UpdateTopology::FaceFace(m); tri::UpdateFlags::FaceBorderFromFF(m); std::map seedMap; // It says if a given vertex of m is a seed (and its index in seedVec) BuildSeedMap(m,seedVec,seedMap); std::vector innerCornerVec, // Faces adjacent to three different regions borderCornerVec; // Faces that are on the border and adjacent to at least two regions. GetFaceCornerVec(m, sources, innerCornerVec, borderCornerVec); // First add all the needed vertices: seeds and corners for(size_t i=0;i::AddVertex(outMesh, seedVec[i]->P(),Color4b::White); map, int > midMapInd; // Given a pair of sources gives the index of the mid vertex map, VertexPointer > midMapPt; if(refineFlag) { GenerateMidPointMap(m, midMapPt); typename std::map, VertexPointer >::iterator mi; for(mi=midMapPt.begin(); mi!=midMapPt.end(); ++mi) { midMapInd[ordered_pair(mi->first.first, mi->first.second)]=outMesh.vert.size(); tri::Allocator::AddVertex(outMesh, mi->second->cP(), Color4b::LightBlue); } } // Now just add one (or four) face for each inner corner for(size_t i=0;iV(0)]; VertexPointer s1 = sources[innerCornerVec[i]->V(1)]; VertexPointer s2 = sources[innerCornerVec[i]->V(2)]; assert ( (s0!=s1) && (s0!=s2) && (s1!=s2) ); VertexPointer v0 = & outMesh.vert[seedMap[s0]]; VertexPointer v1 = & outMesh.vert[seedMap[s1]]; VertexPointer v2 = & outMesh.vert[seedMap[s2]]; if(refineFlag) { VertexPointer mp01 = & outMesh.vert[ midMapInd[ordered_pair(s0,s1)]]; VertexPointer mp02 = & outMesh.vert[ midMapInd[ordered_pair(s0,s2)]]; VertexPointer mp12 = & outMesh.vert[ midMapInd[ordered_pair(s1,s2)]]; assert ( (mp01!=mp02) && (mp01!=mp12) && (mp02!=mp12) ); tri::Allocator::AddFace(outMesh,v0,mp01,mp02); tri::Allocator::AddFace(outMesh,v1,mp12,mp01); tri::Allocator::AddFace(outMesh,v2,mp02,mp12); tri::Allocator::AddFace(outMesh,mp01,mp12,mp02); } else tri::Allocator::AddFace(outMesh,v0,v1,v2); } Clean::RemoveUnreferencedVertex(outMesh); tri::Allocator::CompactVertexVector(outMesh); } template static void PreprocessForVoronoi(MeshType &m, ScalarType radius, MidPointType mid, VoronoiProcessingParameter &vpp) { const int maxSubDiv = 10; tri::RequireFFAdjacency(m); tri::UpdateTopology::FaceFace(m); tri::Clean::RemoveUnreferencedVertex(m); ScalarType edgeLen = tri::Stat::ComputeFaceEdgeLengthAverage(m); for(int i=0;i(m,mid,min(edgeLen*2.0f,radius/vpp.refinementRatio)); if(!ret) break; } tri::Allocator::CompactEveryVector(m); tri::UpdateTopology::VertexFace(m); } static void PreprocessForVoronoi(MeshType &m, float radius, VoronoiProcessingParameter &vpp) { tri::MidPoint mid(&m); PreprocessForVoronoi >(m, radius,mid,vpp); } static void RelaxRefineTriangulationSpring(MeshType &m, MeshType &delaMesh, int refineStep=3, int relaxStep=10 ) { tri::RequireCompactness(m); tri::RequireCompactness(delaMesh); tri::RequireVFAdjacency(delaMesh); tri::RequireFFAdjacency(delaMesh); tri::RequirePerFaceMark(delaMesh); const float convergenceThr = 0.001f; const float eulerStep = 0.1f; tri::UpdateNormal::PerVertexNormalizedPerFaceNormalized(m); typedef GridStaticPtr TriMeshGrid; TriMeshGrid ug; ug.Set(m.face.begin(),m.face.end()); typedef typename vcg::SpatialHashTable HashVertexGrid; HashVertexGrid HG; HG.Set(m.vert.begin(),m.vert.end()); PerVertexBoolHandle fixed = tri::Allocator:: template GetPerVertexAttribute (m,"fixed"); const ScalarType maxDist = m.bbox.Diag()/4.f; for(int kk=0;kk::FaceFace(delaMesh); if(kk!=0) // first step do not refine; { int nonManif = tri::Clean::CountNonManifoldEdgeFF(delaMesh); if(nonManif) return; tri::Refine >(delaMesh,tri::MidPoint(&delaMesh)); } tri::UpdateTopology::VertexFace(delaMesh); const float dist_upper_bound=m.bbox.Diag()/10.0; float dist; for(int k=0;k avgForce(delaMesh.vn); std::vector avgLenVec(delaMesh.vn,0); for(int i=0;i starVec; face::VVStarVF(&delaMesh.vert[i],starVec); for(size_t j=0;jcP()); avgLenVec[i] /= float(starVec.size()); avgForce[i] =Point3f(0,0,0); for(size_t j=0;jcP(); float len = force.Norm(); force.Normalize(); avgForce[i] += force * (avgLenVec[i]-len); } } bool changed=false; for(int i=0;i(m, HG, delaMesh.vert[i].P(), dist_upper_bound, dist); if(!fixed[vp] && !(vp->IsS())) // update only non fixed vertices { delaMesh.vert[i].P() += (avgForce[i]*eulerStep); CoordType closest; float dist; tri::GetClosestFaceBase(m,ug,delaMesh.vert[i].cP(), maxDist,dist,closest); assert(dist!=maxDist); if(Distance(closest,delaMesh.vert[i].P()) > avgLenVec[i]*convergenceThr) changed = true; delaMesh.vert[i].P()=closest; } } if(!changed) k=relaxStep; } // end for k } } static void RelaxRefineTriangulationLaplacian(MeshType &m, MeshType &delaMesh, int refineStep=3, int relaxStep=10 ) { tri::RequireCompactness(m); tri::RequireCompactness(delaMesh); tri::RequireFFAdjacency(delaMesh); tri::RequirePerFaceMark(delaMesh); tri::UpdateTopology::FaceFace(delaMesh); typedef GridStaticPtr TriMeshGrid; TriMeshGrid ug; ug.Set(m.face.begin(),m.face.end()); const ScalarType maxDist = m.bbox.Diag()/4.f; int origVertNum = delaMesh.vn; for(int k=0;k::VertexClear(delaMesh); tri::Refine >(delaMesh,tri::MidPoint(&delaMesh)); for(int j=0;j::VertexCoordLaplacian(delaMesh,1,true); for(int i=origVertNum;i::VertexCoordLaplacianBlend(delaMesh,1,0.2f,true); } } for(int i=origVertNum;i