1 /****************************************************************************
2 * VCGLib o o *
3 * Visual and Computer Graphics Library o o *
4 * _ O _ *
5 * Copyright(C) 2004-2016 \/)\/ *
6 * Visual Computing Lab /\/| *
7 * ISTI - Italian National Research Council | *
8 * \ *
9 * All rights reserved. *
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 * This program is distributed in the hope that it will be useful, *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
20 * for more details. *
21 * *
22 ****************************************************************************/
23 #ifndef CUT_TREE_H
24 #define CUT_TREE_H
25
26 #include<vcg/complex/complex.h>
27 #include <vcg/space/index/kdtree/kdtree.h>
28 #include<vcg/complex/algorithms/update/quality.h>
29 #include<vcg/complex/algorithms/update/color.h>
30
31 namespace vcg {
32 namespace tri {
33
34 template <class MeshType>
35 class CutTree
36 {
37 public:
38 typedef typename MeshType::ScalarType ScalarType;
39 typedef typename MeshType::CoordType CoordType;
40 typedef typename MeshType::VertexType VertexType;
41 typedef typename MeshType::VertexPointer VertexPointer;
42 typedef typename MeshType::VertexIterator VertexIterator;
43 typedef typename MeshType::EdgeIterator EdgeIterator;
44 typedef typename MeshType::EdgeType EdgeType;
45 typedef typename MeshType::FaceType FaceType;
46 typedef typename MeshType::FacePointer FacePointer;
47 typedef typename MeshType::FaceIterator FaceIterator;
48 typedef Box3<ScalarType> Box3Type;
49 typedef typename face::Pos<FaceType> PosType;
50 typedef typename tri::UpdateTopology<MeshType>::PEdge PEdge;
51
52 MeshType &base;
53
CutTree(MeshType & _m)54 CutTree(MeshType &_m) :base(_m){}
55
56
57 // Perform a simple optimization of the three applying simple shortcuts:
58 // if the endpoints of two consecutive edges are connected by an edge existing on base mesh just use that edges
59
OptimizeTree(KdTree<ScalarType> & kdtree,MeshType & t)60 void OptimizeTree(KdTree<ScalarType> &kdtree, MeshType &t)
61 {
62 tri::Allocator<MeshType>::CompactEveryVector(t);
63 int lastEn=t.en;
64 do
65 {
66 lastEn=t.en;
67 tri::UpdateTopology<MeshType>::VertexEdge(t);
68
69 // First simple loop that search for 2->1 moves.
70 for(VertexIterator vi=t.vert.begin();vi!=t.vert.end();++vi)
71 {
72 std::vector<VertexType *> starVec;
73 edge::VVStarVE(&*vi,starVec);
74 if(starVec.size()==2) // middle vertex has to be 1-manifold
75 {
76 PosType pos;
77 if(ExistEdge(kdtree,starVec[0]->P(),starVec[1]->P(),pos))
78 edge::VEEdgeCollapse(t,&*vi);
79 }
80 }
81 tri::Allocator<MeshType>::CompactEveryVector(t);
82 }
83 while(t.en<lastEn);
84 }
85
86 // Given two points return true if on the base mesh there exist an edge with that two coords
87 // if return true the pos indicate the found edge.
ExistEdge(KdTree<ScalarType> & kdtree,CoordType & p0,CoordType & p1,PosType & fpos)88 bool ExistEdge(KdTree<ScalarType> &kdtree, CoordType &p0, CoordType &p1, PosType &fpos)
89 {
90 ScalarType locEps = SquaredDistance(p0,p1)/100000.0;
91
92 VertexType *v0=0,*v1=0;
93 unsigned int veInd;
94 ScalarType sqdist;
95 kdtree.doQueryClosest(p0,veInd,sqdist);
96 if(sqdist<locEps)
97 v0 = &base.vert[veInd];
98 kdtree.doQueryClosest(p1,veInd,sqdist);
99 if(sqdist<locEps)
100 v1 = &base.vert[veInd];
101 if(v0 && v1)
102 {
103 fpos =PosType(v0->VFp(),v0);
104 assert(fpos.V()==v0);
105 PosType startPos=fpos;
106 do
107 {
108 fpos.FlipE(); fpos.FlipF();
109 if(fpos.VFlip()== v1) return true;
110 } while(startPos!=fpos);
111 }
112 return false;
113 }
114
115
findNonVisitedEdgesDuringRetract(VertexType * vp,EdgeType * & ep)116 int findNonVisitedEdgesDuringRetract(VertexType * vp, EdgeType * &ep)
117 {
118 std::vector<EdgeType *> starVec;
119 edge::VEStarVE(&*vp,starVec);
120 int cnt =0;
121 for(size_t i=0;i<starVec.size();++i) {
122 if(!starVec[i]->IsV()) {
123 cnt++;
124 ep = starVec[i];
125 }
126 }
127 return cnt;
128 }
129
IsBoundaryVertexOnBase(KdTree<ScalarType> & kdtree,const CoordType & p)130 bool IsBoundaryVertexOnBase(KdTree<ScalarType> &kdtree, const CoordType &p)
131 {
132 VertexType *v0=0;
133 unsigned int veInd;
134 ScalarType sqdist;
135 kdtree.doQueryClosest(p,veInd,sqdist);
136 if(sqdist>0) { assert(0); }
137 v0 = &base.vert[veInd];
138 return v0->IsB();
139 }
140
141 /**
142 * @brief Retract
143 * @param t the edgemesh containing the visit tree.
144 *
145 */
Retract(KdTree<ScalarType> & kdtree,MeshType & t)146 void Retract(KdTree<ScalarType> &kdtree, MeshType &t)
147 {
148 printf("Retracting a tree of %i edges and %i vertices\n",t.en,t.vn);
149 tri::UpdateTopology<MeshType>::VertexEdge(t);
150 tri::Allocator<MeshType>::CompactEveryVector(t);
151 std::stack<VertexType *> vertStack;
152
153 // Put on the stack all the vertex with just a single incident edge.
154 ForEachVertex(t, [&](VertexType &v){
155 if(edge::VEDegree<EdgeType>(&v) ==1)
156 vertStack.push(&v);
157 });
158
159 tri::UpdateFlags<MeshType>::EdgeClearV(t);
160 tri::UpdateFlags<MeshType>::VertexClearV(t);
161
162 int unvisitedEdgeNum = t.en;
163 while((!vertStack.empty()) && (unvisitedEdgeNum > 2) )
164 {
165 VertexType *vp = vertStack.top();
166 vertStack.pop();
167 vp->C()=Color4b::Blue;
168 EdgeType *ep=0;
169 int eCnt = findNonVisitedEdgesDuringRetract(vp,ep);
170 if(eCnt==1) // We have only one non visited edge over vp
171 {
172 assert(!ep->IsV());
173 ep->SetV();
174 --unvisitedEdgeNum;
175 VertexType *otherVertP;
176 if(ep->V(0)==vp) otherVertP = ep->V(1);
177 else otherVertP = ep->V(0);
178 vertStack.push(otherVertP);
179 }
180 }
181 assert(unvisitedEdgeNum >0);
182 for(size_t i =0; i<t.edge.size();++i){
183 PosType fpos;
184 if( ExistEdge(kdtree, t.edge[i].P(0), t.edge[i].P(1), fpos)){
185 if(fpos.IsBorder()) {
186 t.edge[i].SetV();
187 }
188 }
189 else assert(0);
190 }
191
192 // All the boundary edges are in the initial tree so the clean boundary loops chains remains as irreducible loops
193 // We delete them (leaving dangling edges with a vertex on the boundary)
194 for(size_t i =0; i<t.edge.size();++i){
195 if (t.edge[i].IsV())
196 tri::Allocator<MeshType>::DeleteEdge(t,t.edge[i]) ;
197 }
198 assert(t.en >0);
199 tri::Clean<MeshType>::RemoveUnreferencedVertex(t);
200 tri::Allocator<MeshType>::CompactEveryVector(t);
201 }
202
203 /** \brief Main function
204 *
205 * It builds a cut tree that open the mesh into a topological disk
206 *
207 *
208 */
209 void Build(MeshType &dualMesh, int startingFaceInd=0)
210 {
211 tri::UpdateTopology<MeshType>::FaceFace(base);
212 tri::UpdateTopology<MeshType>::VertexFace(base);
213
214 BuildVisitTree(dualMesh,startingFaceInd);
215 // BuildDijkstraVisitTree(dualMesh,startingFaceInd);
216
217 VertexConstDataWrapper<MeshType > vdw(base);
218 KdTree<ScalarType> kdtree(vdw);
219 Retract(kdtree,dualMesh);
220 OptimizeTree(kdtree, dualMesh);
221 tri::UpdateBounding<MeshType>::Box(dualMesh);
222 }
223
224 /* Auxiliary class for keeping the heap of vertices to visit and their estimated distance */
225 struct FaceDist{
FaceDistFaceDist226 FaceDist(FacePointer _f):f(_f),dist(_f->Q()){}
227 FacePointer f;
228 ScalarType dist;
229 bool operator < (const FaceDist &o) const
230 {
231 if( dist != o.dist)
232 return dist > o.dist;
233 return f<o.f;
234 }
235 };
236
237
238 void BuildDijkstraVisitTree(MeshType &dualMesh, int startingFaceInd=0, ScalarType maxDistanceThr=std::numeric_limits<ScalarType>::max())
239 {
240 tri::RequireFFAdjacency(base);
241 tri::RequirePerFaceMark(base);
242 tri::RequirePerFaceQuality(base);
243 typename MeshType::template PerFaceAttributeHandle<FacePointer> parentHandle
244 = tri::Allocator<MeshType>::template GetPerFaceAttribute<FacePointer>(base, "parent");
245
246 std::vector<FacePointer> seedVec;
247 seedVec.push_back(&base.face[startingFaceInd]);
248
249 std::vector<FaceDist> Heap;
250 tri::UnMarkAll(base);
251 tri::UpdateQuality<MeshType>::FaceConstant(base,0);
252 ForEachVertex(base, [&](VertexType &v){
253 tri::Allocator<MeshType>::AddVertex(dualMesh,v.cP());
254 });
255
256 // Initialize the face heap;
257 // All faces in the heap are already marked; Q() store the distance from the source faces;
258 for(size_t i=0;i<seedVec.size();++i)
259 {
260 seedVec[i]->Q()=0;
261 Heap.push_back(FaceDist(seedVec[i]));
262 }
263 // Main Loop
264 int boundary=0;
265 std::make_heap(Heap.begin(),Heap.end());
266
267 int vCnt=0;
268 int eCnt=0;
269 int fCnt=0;
270
271 // The main idea is that in the heap we maintain all the faces to be visited.
272 int nonDiskCnt=0;
273 while(!Heap.empty() && nonDiskCnt<10)
274 {
275 int eulerChi= vCnt-eCnt+fCnt;
276 if(eulerChi==1) nonDiskCnt=0;
277 else ++nonDiskCnt;
278 // printf("HeapSize %i: %i - %i + %i = %i\n",Heap.size(), vCnt,eCnt,fCnt,eulerChi);
279 pop_heap(Heap.begin(),Heap.end());
280 FacePointer currFp = (Heap.back()).f;
281 if(tri::IsMarked(base,currFp))
282 {
283 // printf("Found an already visited face %f %f \n",Heap.back().dist, Heap.back().f->Q());
284 //assert(Heap.back().dist != currFp->Q());
285
286 Heap.pop_back();
287 continue;
288 }
289 Heap.pop_back();
290 ++fCnt;
291 eCnt+=3;
292 tri::Mark(base,currFp);
293
294 // printf("pop face %i \n", tri::Index(base,currFp));
295 for(int i=0;i<3;++i)
296 {
297 if(!currFp->V(i)->IsV()) {++vCnt; currFp->V(i)->SetV();}
298
299 FacePointer nextFp = currFp->FFp(i);
300 if( tri::IsMarked(base,nextFp) )
301 {
302 eCnt-=1;
303 printf("is marked\n");
304 if(nextFp != parentHandle[currFp] )
305 {
306 if(currFp>nextFp){
307 tri::Allocator<MeshType>::AddEdge(dualMesh,tri::Index(base,currFp->V0(i)), tri::Index(base,currFp->V1(i)));
308 }
309 }
310 }
311 else // add it to the heap;
312 {
313 // printf("is NOT marked\n");
314 parentHandle[nextFp] = currFp;
315 ScalarType nextDist = currFp->Q() + Distance(Barycenter(*currFp),Barycenter(*nextFp));
316 int adjMarkedNum=0;
317 for(int k=0;k<3;++k) if(tri::IsMarked(base,nextFp->FFp(k))) ++adjMarkedNum;
318 if(nextDist < maxDistanceThr || adjMarkedNum>1)
319 {
320 nextFp->Q() = nextDist;
321 Heap.push_back(FaceDist(nextFp));
322 push_heap(Heap.begin(),Heap.end());
323 }
324 else {
325 // printf("boundary %i\n",++boundary);
326 tri::Allocator<MeshType>::AddEdge(dualMesh,tri::Index(base,currFp->V0(i)), tri::Index(base,currFp->V1(i)));
327 }
328 }
329 }
330 } // End while
331 printf("fulltree %i vn %i en \n",dualMesh.vn, dualMesh.en);
332 int dupVert=tri::Clean<MeshType>::RemoveDuplicateVertex(dualMesh,false); printf("Removed %i dup vert\n",dupVert);
333 int dupEdge=tri::Clean<MeshType>::RemoveDuplicateEdge(dualMesh); printf("Removed %i dup edges %i\n",dupEdge,dualMesh.EN());
334 tri::Clean<MeshType>::RemoveUnreferencedVertex(dualMesh);
335
336 tri::io::ExporterPLY<MeshType>::Save(dualMesh,"fulltree.ply",tri::io::Mask::IOM_EDGEINDEX);
337 tri::UpdateColor<MeshType>::PerFaceQualityRamp(base);
338 tri::io::ExporterPLY<MeshType>::Save(base,"colored_Bydistance.ply",tri::io::Mask::IOM_FACECOLOR);
339 }
340
341 // \brief This function build a cut tree.
342 //
343 // First we make a bread first FF face visit.
344 // Each time that we encounter a visited face we add to the tree the edge
345 // that brings to the already visited face.
346 // this structure build a dense graph and we retract this graph retracting each
347 // leaf until we remains with just the loops that cuts the object.
348
349 void BuildVisitTree(MeshType &dualMesh, int startingFaceInd=0)
350 {
351 tri::UpdateFlags<MeshType>::FaceClearV(base);
352 tri::UpdateFlags<MeshType>::VertexBorderFromFaceAdj(base);
353 std::vector<face::Pos<FaceType> > visitStack; // the stack contain the pos on the 'starting' face.
354
355 base.face[startingFaceInd].SetV();
356 for(int i=0;i<3;++i)
357 visitStack.push_back(PosType(&(base.face[startingFaceInd]),i,base.face[startingFaceInd].V(i)));
358
359 int cnt=1;
360
361 while(!visitStack.empty())
362 {
363 std::swap(visitStack.back(),visitStack[rand()%visitStack.size()]);
364 PosType c = visitStack.back();
365 visitStack.pop_back();
366 assert(c.F()->IsV());
367 c.F()->C() = Color4b::ColorRamp(0,base.fn,cnt);
368 c.FlipF();
369 if(!c.F()->IsV())
370 {
371 ++cnt;
372 c.F()->SetV();
373 c.FlipE();c.FlipV();
374 visitStack.push_back(c);
375 c.FlipE();c.FlipV();
376 visitStack.push_back(c);
377 }
378 else
379 {
380 tri::Allocator<MeshType>::AddEdge(dualMesh,c.V()->P(),c.VFlip()->P());
381 }
382 }
383 assert(cnt==base.fn);
384
385 tri::Clean<MeshType>::RemoveDuplicateVertex(dualMesh);
386 tri::io::ExporterPLY<MeshType>::Save(dualMesh,"fulltree.ply",tri::io::Mask::IOM_EDGEINDEX);
387 }
388
389 };
390 } // end namespace tri
391 } // end namespace vcg
392
393 #endif // CUT_TREE_H
394