VCG Library
refine.h
1 /***********F*****************************************************************
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 
24 #ifndef __VCGLIB_REFINE
25 #define __VCGLIB_REFINE
26 
27 #include <functional>
28 #include <map>
29 #include <vcg/space/sphere3.h>
30 #include <vcg/space/plane3.h>
31 #include <vcg/complex/algorithms/clean.h>
32 #include <vcg/space/texcoord2.h>
33 #include <vcg/space/triangle3.h>
34 
35 namespace vcg{
36 namespace tri{
37 
38 /* A very short intro about the generic refinement framework,
39  the main fuction is the
40 
41  template<class MESH_TYPE,class MIDPOINT, class EDGEPRED>
42  bool RefineE(MESH_TYPE &m, MIDPOINT mid, EDGEPRED ep,bool RefineSelected=false, CallBackPos *cb = 0)
43 
44  You have to provide two functor objects to this, one for deciding what edge has to be spltted and one to decide position and new values for the attributes of the new point.
45 
46  for example the minimal EDGEPRED is
47 
48  template <class MESH_TYPE, class FLT> class EdgeLen
49  {
50  public:
51  FLT thr2;
52  bool operator()(face::Pos<typename MESH_TYPE::FaceType> ep) const
53  {
54  return SquaredDistance(ep.f->V(ep.z)->P(), ep.f->V1(ep.z)->P())>thr2;
55  }
56  };
57 
58  With a bit of patience you can customize to make also slicing operation.
59 
60 */
61 
62 
63 /* The table which encodes how to subdivide a triangle depending
64  on the splitted edges is organized as such:
65 
66  TriNum (the first number): encodes the number of triangles
67  TV (the following 4 triples): encodes the resulting triangles where
68  0, 1, 2 are the original vertices of the triangles and 3, 4, 5
69  (mp01, mp12, mp20) are the midpoints of the three edges.
70 
71  In the case two edges are splitted the triangle has 2 possible splittings:
72 we need to choose a diagonal of the resulting trapezoid.
73 'swap' encodes the two diagonals to test: if diag1 < diag2 we swap the diagonal
74 like this (140, 504 -> 150, 514) (the second vertex of each triangles is replaced
75  by the first vertex of the other one).
76  2
77  / \
78  5---4
79  / \
80  0-------1
81 
82 */
83 
84 class Split {
85 public:
86  int TriNum; // number of triangles
87  int TV[4][3]; // The triangles coded as the following convention
88  // 0..2 vertici originali del triangolo
89  // 3..5 mp01, mp12, mp20 midpoints of the three edges
90  int swap[2][2]; // the two diagonals to test for swapping
91  int TE[4][3]; // the edge-edge correspondence between refined triangles and the old one
92  // (3) means the edge of the new triangle is internal;
93 };
94 
95 const Split SplitTab[8]={
96 /* m20 m12 m01 */
97 /* 0 0 0 */ {1, {{0,1,2},{0,0,0},{0,0,0},{0,0,0}}, {{0,0},{0,0}}, {{0,1,2},{0,0,0},{0,0,0},{0,0,0}} },
98 /* 0 0 1 */ {2, {{0,3,2},{3,1,2},{0,0,0},{0,0,0}}, {{0,0},{0,0}}, {{0,3,2},{0,1,3},{0,0,0},{0,0,0}} },
99 /* 0 1 0 */ {2, {{0,1,4},{0,4,2},{0,0,0},{0,0,0}}, {{0,0},{0,0}}, {{0,1,3},{3,1,2},{0,0,0},{0,0,0}} },
100 /* 0 1 1 */ {3, {{3,1,4},{0,3,2},{4,2,3},{0,0,0}}, {{0,4},{3,2}}, {{0,1,3},{0,3,2},{1,3,3},{0,0,0}} },
101 /* 1 0 0 */ {2, {{0,1,5},{5,1,2},{0,0,0},{0,0,0}}, {{0,0},{0,0}}, {{0,3,2},{3,1,2},{0,0,0},{0,0,0}} },
102 /* 1 0 1 */ {3, {{0,3,5},{3,1,5},{2,5,1},{0,0,0}}, {{3,2},{5,1}}, {{0,3,2},{0,3,3},{2,3,1},{0,0,0}} },
103 /* 1 1 0 */ {3, {{2,5,4},{0,1,5},{4,5,1},{0,0,0}}, {{0,4},{5,1}}, {{2,3,1},{0,3,2},{3,3,1},{0,0,0}} },
104 /* 1 1 1 */ //{4, {{0,3,5},{3,1,4},{5,4,2},{3,4,5}}, {{0,0},{0,0}}, {{0,3,2},{0,1,3},{3,1,2},{3,3,3}} },
105 /* 1 1 1 */ {4, {{3,4,5},{0,3,5},{3,1,4},{5,4,2}}, {{0,0},{0,0}}, {{3,3,3},{0,3,2},{0,1,3},{3,1,2}} },
106 };
107 
108 
109 template <class MeshType>
110 struct BaseInterpolator
111 {
112  typedef typename face::Pos<typename MeshType::FaceType> PosType;
113  typedef typename MeshType::VertexType VertexType;
114  void operator()(VertexType &, PosType ){}
115 };
116 
117 // Basic subdivision class
118 // This class must provide methods for finding the position of the newly created vertices
119 // In this implemenation we simply put the new vertex in the MidPoint position.
120 // Color and TexCoords are interpolated accordingly.
121 // This subdivision class allow also the correct interpolation of userdefined data by
122 // providing, in the constructor, an interpolator functor that will be called for each new vertex to be created.
123 
124 template<class MESH_TYPE, class InterpolatorFunctorType = BaseInterpolator< MESH_TYPE> >
125 struct MidPoint : public std::unary_function<face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType >
126 {
127  typedef typename face::Pos<typename MESH_TYPE::FaceType> PosType;
128  typedef typename MESH_TYPE::VertexType VertexType;
129 
130  MidPoint(MESH_TYPE *_mp,
131  InterpolatorFunctorType *_intFunc=0) {
132  mp=_mp;
133  intFunc =_intFunc;
134  }
135 
136  MESH_TYPE *mp;
137  InterpolatorFunctorType *intFunc;
138 
139  void operator()(VertexType &nv, PosType ep){
140  assert(mp);
141  VertexType *V0 = ep.V() ;
142  VertexType *V1 = ep.VFlip() ;
143  if(V0 > V1) std::swap(V1,V0);
144 
145  nv.P()= (V0->P()+V1->P())/2.0;
146 
147  if( tri::HasPerVertexNormal(*mp))
148  nv.N()= (V0->N()+V1->N()).normalized();
149 
150  if( tri::HasPerVertexColor(*mp))
151  nv.C().lerp(V0->C(),V1->C(),.5f);
152 
153  if( tri::HasPerVertexQuality(*mp))
154  nv.Q() = (V0->Q()+V1->Q()) / 2.0;
155 
156  if( tri::HasPerVertexTexCoord(*mp))
157  nv.T().P() = (V0->T().P()+V1->T().P()) / 2.0;
158  if(intFunc)
159  (*intFunc)(nv,ep);
160  }
161 
162  Color4<typename MESH_TYPE::ScalarType> WedgeInterp(Color4<typename MESH_TYPE::ScalarType> &c0, Color4<typename MESH_TYPE::ScalarType> &c1)
163  {
164  Color4<typename MESH_TYPE::ScalarType> cc;
165  return cc.lerp(c0,c1,0.5f);
166  }
167 
168  template<class FL_TYPE>
169  TexCoord2<FL_TYPE,1> WedgeInterp(TexCoord2<FL_TYPE,1> &t0, TexCoord2<FL_TYPE,1> &t1)
170  {
171  TexCoord2<FL_TYPE,1> tmp;
172  assert(t0.n()== t1.n());
173  tmp.n()=t0.n();
174  tmp.t()=(t0.t()+t1.t())/2.0;
175  return tmp;
176  }
177 };
178 
179 
180 
181 template<class MESH_TYPE>
182 struct MidPointArc : public std::unary_function<face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType>
183 {
184  void operator()(typename MESH_TYPE::VertexType &nv, face::Pos<typename MESH_TYPE::FaceType> ep)
185  {
186  const typename MESH_TYPE::ScalarType EPS =1e-10;
187  typename MESH_TYPE::CoordType vp = (ep.f->V(ep.z)->P()+ep.f->V1(ep.z)->P())/2.0;
188  typename MESH_TYPE::CoordType n = (ep.f->V(ep.z)->N()+ep.f->V1(ep.z)->N())/2.0;
189  typename MESH_TYPE::ScalarType w =n.Norm();
190  if(w<EPS) { nv.P()=(ep.f->V(ep.z)->P()+ep.f->V1(ep.z)->P())/2.0; return;}
191  n/=w;
192  typename MESH_TYPE::CoordType d0 = ep.f->V(ep.z)->P() - vp;
193  typename MESH_TYPE::CoordType d1 = ep.f->V1(ep.z)->P()- vp;
194  typename MESH_TYPE::ScalarType d=Distance(ep.f->V(ep.z)->P(),ep.f->V1(ep.z)->P())/2.0;
195 
196  typename MESH_TYPE::CoordType nn = ep.f->V1(ep.z)->N() ^ ep.f->V(ep.z)->N();
197  typename MESH_TYPE::CoordType np = n ^ d0; //vector perpendicular to the edge plane, normal is interpolated
198  np.Normalize();
199  double sign=1;
200  if(np*nn<0) sign=-1; // se le normali non divergono sposta il punto nella direzione opposta
201 
202  typename MESH_TYPE::CoordType n0=ep.f->V(ep.z)->N() -np*(ep.f->V(ep.z)->N()*np);
203  n0.Normalize();
204  typename MESH_TYPE::CoordType n1=ep.f->V1(ep.z)->N()-np*(ep.f->V1(ep.z)->N()*np);
205  assert(n1.Norm()>EPS);
206  n1.Normalize();
207  typename MESH_TYPE::ScalarType cosa0=n0*n;
208  typename MESH_TYPE::ScalarType cosa1=n1*n;
209  if(2-cosa0-cosa1<EPS) {nv.P()=(ep.f->V(ep.z)->P()+ep.f->V1(ep.z)->P())/2.0;return;}
210  typename MESH_TYPE::ScalarType cosb0=(d0*n)/d;
211  typename MESH_TYPE::ScalarType cosb1=(d1*n)/d;
212  assert(1+cosa0>EPS);
213  assert(1+cosa1>EPS);
214  typename MESH_TYPE::ScalarType delta0=d*(cosb0 +sqrt( ((1-cosb0*cosb0)*(1-cosa0))/(1+cosa0)) );
215  typename MESH_TYPE::ScalarType delta1=d*(cosb1 +sqrt( ((1-cosb1*cosb1)*(1-cosa1))/(1+cosa1)) );
216  assert(delta0+delta1<2*d);
217  nv.P()=vp+n*sign*(delta0+delta1)/2.0;
218  return ;
219  }
220 
221  // Aggiunte in modo grezzo le due wedgeinterp
222  Color4<typename MESH_TYPE::ScalarType> WedgeInterp(Color4<typename MESH_TYPE::ScalarType> &c0, Color4<typename MESH_TYPE::ScalarType> &c1)
223  {
224  Color4<typename MESH_TYPE::ScalarType> cc;
225  return cc.lerp(c0,c1,0.5f);
226  }
227 
228  template<class FL_TYPE>
229  TexCoord2<FL_TYPE,1> WedgeInterp(TexCoord2<FL_TYPE,1> &t0, TexCoord2<FL_TYPE,1> &t1)
230  {
231  TexCoord2<FL_TYPE,1> tmp;
232  assert(t0.n()== t1.n());
233  tmp.n()=t0.n();
234  tmp.t()=(t0.t()+t1.t())/2.0;
235  return tmp;
236  }
237 
238 };
239 
240 /*
241 Versione Della Midpoint basata sul paper:
242 S. Karbacher, S. Seeger, G. Hausler
243 A non linear subdivision scheme for triangle meshes
244 
245  Non funziona!
246  Almeno due problemi:
247  1) il verso delle normali influenza il risultato (e.g. si funziona solo se le normali divergono)
248  Risolvibile controllando se le normali divergono
249  2) gli archi vanno calcolati sul piano definito dalla normale interpolata e l'edge.
250  funziona molto meglio nelle zone di sella e non semplici.
251 
252 */
253 template<class MESH_TYPE>
254 struct MidPointArcNaive : public std::unary_function< face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType>
255 {
256  typename MESH_TYPE::CoordType operator()(face::Pos<typename MESH_TYPE::FaceType> ep)
257  {
258 
259  typename MESH_TYPE::CoordType vp = (ep.f->V(ep.z)->P()+ep.f->V1(ep.z)->P())/2.0;
260  typename MESH_TYPE::CoordType n = (ep.f->V(ep.z)->N()+ep.f->V1(ep.z)->N())/2.0;
261  n.Normalize();
262  typename MESH_TYPE::CoordType d0 = ep.f->V(ep.z)->P() - vp;
263  typename MESH_TYPE::CoordType d1 = ep.f->V1(ep.z)->P()- vp;
264  typename MESH_TYPE::ScalarType d=Distance(ep.f->V(ep.z)->P(),ep.f->V1(ep.z)->P())/2.0;
265 
266  typename MESH_TYPE::ScalarType cosa0=ep.f->V(ep.z)->N()*n;
267  typename MESH_TYPE::ScalarType cosa1=ep.f->V1(ep.z)->N()*n;
268  typename MESH_TYPE::ScalarType cosb0=(d0*n)/d;
269  typename MESH_TYPE::ScalarType cosb1=(d1*n)/d;
270 
271  typename MESH_TYPE::ScalarType delta0=d*(cosb0 +sqrt( ((1-cosb0*cosb0)*(1-cosa0))/(1+cosa0)) );
272  typename MESH_TYPE::ScalarType delta1=d*(cosb1 +sqrt( ((1-cosb1*cosb1)*(1-cosa1))/(1+cosa1)) );
273 
274  return vp+n*(delta0+delta1)/2.0;
275  }
276 };
277 
278 
279 // Basic Predicate that tells if a given edge must be splitted.
280 // the constructure requires the threshold.
281 // VERY IMPORTANT REQUIREMENT: this function must be symmetric
282 // e.g. it must return the same value if the Pos is VFlipped.
283 // If this function is not symmetric the Refine can crash.
284 
285 template <class MESH_TYPE, class FLT>
286 class EdgeLen
287 {
288  FLT squaredThr;
289 public:
290  EdgeLen(){};
291  EdgeLen(FLT threshold) {setThr(threshold);}
292  void setThr(FLT threshold) {squaredThr = threshold*threshold; }
293  bool operator()(face::Pos<typename MESH_TYPE::FaceType> ep) const
294  {
295  return SquaredDistance(ep.V()->P(), ep.VFlip()->P())>squaredThr;
296  }
297 };
298 
299 /*********************************************************/
300 /*********************************************************
301 
302 Given a mesh the following function refines it according to two functor objects:
303 
304 - a predicate that tells if a given edge must be splitted
305 
306 - a functor that gives you the new poistion of the created vertices (starting from an edge)
307 
308 If RefineSelected is true only selected faces are taken into account for being splitted.
309 
310 Requirement: FF Adjacency and Manifoldness
311 
312 **********************************************************/
313 /*********************************************************/
314 template <class VertexPointer>
315 class RefinedFaceData
316  {
317  public:
318  RefinedFaceData(){
319  ep[0]=0;ep[1]=0;ep[2]=0;
320  vp[0]=0;vp[1]=0;vp[2]=0;
321  }
322  bool ep[3];
323  VertexPointer vp[3];
324  };
325 
326 template<class MESH_TYPE,class MIDPOINT, class EDGEPRED>
327 bool RefineE(MESH_TYPE &m, MIDPOINT &mid, EDGEPRED &ep,bool RefineSelected=false, CallBackPos *cb = 0)
328 {
329  // common typenames
330  typedef typename MESH_TYPE::VertexIterator VertexIterator;
331  typedef typename MESH_TYPE::FaceIterator FaceIterator;
332  typedef typename MESH_TYPE::VertexPointer VertexPointer;
333  typedef typename MESH_TYPE::FacePointer FacePointer;
334  typedef typename MESH_TYPE::FaceType FaceType;
335  typedef typename MESH_TYPE::FaceType::TexCoordType TexCoordType;
336  assert(tri::HasFFAdjacency(m));
338  typedef face::Pos<FaceType> PosType;
339 
340  int j,NewVertNum=0,NewFaceNum=0;
341 
342  typedef RefinedFaceData<VertexPointer> RFD;
343  typedef typename MESH_TYPE :: template PerFaceAttributeHandle<RFD> HandleType;
344  HandleType RD = tri::Allocator<MESH_TYPE>:: template AddPerFaceAttribute<RFD> (m,std::string("RefineData"));
345 
346  // Callback stuff
347  int step=0;
348  int PercStep=std::max(1,m.fn/33);
349 
350  // First Loop: We analyze the mesh to compute the number of the new faces and new vertices
351  FaceIterator fi;
352  for(fi=m.face.begin(),j=0;fi!=m.face.end();++fi) if(!(*fi).IsD())
353  {
354  if(cb && (++step%PercStep)==0) (*cb)(step/PercStep,"Refining...");
355  // skip unselected faces if necessary
356  if(RefineSelected && !(*fi).IsS()) continue;
357 
358  for(j=0;j<3;j++)
359  {
360  if(RD[fi].ep[j]) continue;
361 
362  PosType edgeCur(&*fi,j);
363  if(RefineSelected && ! edgeCur.FFlip()->IsS()) continue;
364  if(!ep(edgeCur)) continue;
365 
366  RD[edgeCur.F()].ep[edgeCur.E()]=true;
367  ++NewFaceNum;
368  ++NewVertNum;
369  assert(edgeCur.IsManifold());
370  if(!edgeCur.IsBorder())
371  {
372  edgeCur.FlipF();
373  edgeCur.F()->SetV();
374  RD[edgeCur.F()].ep[edgeCur.E()]=true;
375  ++NewFaceNum;
376  }
377  }
378 
379  } // end face loop
380 
381  if(NewVertNum ==0 )
382  {
383  tri::Allocator<MESH_TYPE> :: template DeletePerFaceAttribute<RefinedFaceData<VertexPointer> > (m,RD);
384  return false;
385  }
386  VertexIterator lastv = tri::Allocator<MESH_TYPE>::AddVertices(m,NewVertNum);
387 
388  // Secondo loop: We initialize a edge->vertex map
389 
390  for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
391  {
392  if(cb && (++step%PercStep)==0)(*cb)(step/PercStep,"Refining...");
393  for(j=0;j<3;j++)
394  {
395  // skip unselected faces if necessary
396  if(RefineSelected && !(*fi).IsS()) continue;
397  for(j=0;j<3;j++)
398  {
399  PosType edgeCur(&*fi,j);
400  if(RefineSelected && ! edgeCur.FFlip()->IsS()) continue;
401 
402  if( RD[edgeCur.F()].ep[edgeCur.E()] && RD[edgeCur.F()].vp[edgeCur.E()] ==0 )
403  {
404  RD[edgeCur.F()].vp[edgeCur.E()] = &*lastv;
405  mid(*lastv,edgeCur);
406  if(!edgeCur.IsBorder())
407  {
408  edgeCur.FlipF();
409  assert(RD[edgeCur.F()].ep[edgeCur.E()]);
410  RD[edgeCur.F()].vp[edgeCur.E()] = &*lastv;
411  }
412  ++lastv;
413  }
414  }
415  }
416  }
417 
418  assert(lastv==m.vert.end()); // critical assert: we MUST have used all the vertex that we forecasted we need
419 
420  FaceIterator lastf = tri::Allocator<MESH_TYPE>::AddFaces(m,NewFaceNum);
421  FaceIterator oldendf = lastf;
422 
423 /*
424  * v0
425  *
426  *
427  * f0
428  *
429  * mp01 f3 mp02
430  *
431  *
432  * f1 f2
433  *
434  *v1 mp12 v2
435  *
436 */
437 
438  VertexPointer vv[6]; // The six vertices that arise in the single triangle splitting
439  // 0..2 Original triangle vertices
440  // 3..5 mp01, mp12, mp20 midpoints of the three edges
441  FacePointer nf[4]; // The (up to) four faces that are created.
442 
443  TexCoordType wtt[6]; // per ogni faccia sono al piu' tre i nuovi valori
444  // di texture per wedge (uno per ogni edge)
445 
446  int fca=0,fcn =0;
447  for(fi=m.face.begin();fi!=oldendf;++fi) if(!(*fi).IsD())
448  {
449  if(cb && (++step%PercStep)==0)(*cb)(step/PercStep,"Refining...");
450  fcn++;
451  vv[0]=(*fi).V(0);
452  vv[1]=(*fi).V(1);
453  vv[2]=(*fi).V(2);
454  vv[3] = RD[fi].vp[0];
455  vv[4] = RD[fi].vp[1];
456  vv[5] = RD[fi].vp[2];
457 
458  int ind=((&*vv[3])?1:0)+((&*vv[4])?2:0)+((&*vv[5])?4:0);
459 
460  nf[0]=&*fi;
461  int i;
462  for(i=1;i<SplitTab[ind].TriNum;++i){
463  nf[i]=&*lastf; ++lastf; fca++;
464  if(RefineSelected || (*fi).IsS()) (*nf[i]).SetS();
465  nf[i]->ImportData(*fi);
466 // if(tri::HasPerFaceColor(m))
467 // nf[i]->C()=(*fi).cC();
468  }
469 
470 
471  if(tri::HasPerWedgeTexCoord(m))
472  for(i=0;i<3;++i) {
473  wtt[i]=(*fi).WT(i);
474  wtt[3+i]=mid.WedgeInterp((*fi).WT(i),(*fi).WT((i+1)%3));
475  }
476 
477  int orgflag= (*fi).Flags();
478  for(i=0;i<SplitTab[ind].TriNum;++i)
479  for(j=0;j<3;++j){
480  (*nf[i]).V(j)=&*vv[SplitTab[ind].TV[i][j]];
481 
482  if(tri::HasPerWedgeTexCoord(m)) //analogo ai vertici...
483  (*nf[i]).WT(j)=wtt[SplitTab[ind].TV[i][j]];
484 
485  assert((*nf[i]).V(j)!=0);
486  if(SplitTab[ind].TE[i][j]!=3){
487  if(orgflag & (MESH_TYPE::FaceType::BORDER0<<(SplitTab[ind].TE[i][j])))
488  (*nf[i]).SetB(j);
489  else
490  (*nf[i]).ClearB(j);
491  }
492  else (*nf[i]).ClearB(j);
493  }
494 
495  if(SplitTab[ind].TriNum==3 &&
496  SquaredDistance(vv[SplitTab[ind].swap[0][0]]->P(),vv[SplitTab[ind].swap[0][1]]->P()) <
497  SquaredDistance(vv[SplitTab[ind].swap[1][0]]->P(),vv[SplitTab[ind].swap[1][1]]->P()) )
498  { // swap the last two triangles
499  (*nf[2]).V(1)=(*nf[1]).V(0);
500  (*nf[1]).V(1)=(*nf[2]).V(0);
501  if(tri::HasPerWedgeTexCoord(m)){ //swap also textures coordinates
502  (*nf[2]).WT(1)=(*nf[1]).WT(0);
503  (*nf[1]).WT(1)=(*nf[2]).WT(0);
504  }
505 
506  if((*nf[1]).IsB(0)) (*nf[2]).SetB(1); else (*nf[2]).ClearB(1);
507  if((*nf[2]).IsB(0)) (*nf[1]).SetB(1); else (*nf[1]).ClearB(1);
508  (*nf[1]).ClearB(0);
509  (*nf[2]).ClearB(0);
510  }
511  }
512 
513  assert(lastf==m.face.end()); // critical assert: we MUST have used all the faces that we forecasted we need and that we previously allocated.
514  assert(!m.vert.empty());
515  for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()){
516  assert((*fi).V(0)>=&*m.vert.begin() && (*fi).V(0)<=&m.vert.back() );
517  assert((*fi).V(1)>=&*m.vert.begin() && (*fi).V(1)<=&m.vert.back() );
518  assert((*fi).V(2)>=&*m.vert.begin() && (*fi).V(2)<=&m.vert.back() );
519  }
521 
522  tri::Allocator<MESH_TYPE> :: template DeletePerFaceAttribute<RefinedFaceData<VertexPointer> > (m,RD);
523 
524  return true;
525 }
526 
527 /*************************************************************************/
528 // simple wrapper of the base refine for lazy coder that do not need a edge predicate
529 
530 template<class MESH_TYPE,class MIDPOINT>
531 bool Refine(MESH_TYPE &m, MIDPOINT mid, typename MESH_TYPE::ScalarType thr=0,bool RefineSelected=false, CallBackPos *cb = 0)
532 {
533  EdgeLen <MESH_TYPE, typename MESH_TYPE::ScalarType> ep(thr);
534  return RefineE(m,mid,ep,RefineSelected,cb);
535 }
536 /*************************************************************************/
537 
538 /*
539 Modified Butterfly interpolation scheme,
540 as presented in
541 Zorin, Schroeder
542 Subdivision for modeling and animation
543 Siggraph 2000 Course Notes
544 */
545 
546 /*
547 
548  vul-------vu--------vur
549  \ / \ /
550  \ / \ /
551  \ / fu \ /
552  vl--------vr
553  / \ fd / \
554  / \ / \
555  / \ / \
556  vdl-------vd--------vdr
557 
558 */
559 
560 template<class MESH_TYPE>
561 struct MidPointButterfly : public std::unary_function<face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType>
562 {
563  MESH_TYPE &m;
564  MidPointButterfly(MESH_TYPE &_m):m(_m){}
565 
566  void operator()(typename MESH_TYPE::VertexType &nv, face::Pos<typename MESH_TYPE::FaceType> ep)
567  {
568  face::Pos<typename MESH_TYPE::FaceType> he(ep.f,ep.z,ep.f->V(ep.z));
569  typename MESH_TYPE::CoordType *vl,*vr;
570  typename MESH_TYPE::CoordType *vl0,*vr0;
571  typename MESH_TYPE::CoordType *vu,*vd,*vul,*vur,*vdl,*vdr;
572  vl=&he.v->P();
573  he.FlipV();
574  vr=&he.v->P();
575 
576  if( tri::HasPerVertexColor(m))
577  nv.C().lerp(ep.f->V(ep.z)->C(),ep.f->V1(ep.z)->C(),.5f);
578 
579  if(he.IsBorder())
580  {
581  he.NextB();
582  vr0=&he.v->P();
583  he.FlipV();
584  he.NextB();
585  assert(&he.v->P()==vl);
586  he.NextB();
587  vl0=&he.v->P();
588  nv.P()=((*vl)+(*vr))*(9.0/16.0)-((*vl0)+(*vr0))/16.0 ;
589  }
590  else
591  {
592  he.FlipE();he.FlipV();
593  vu=&he.v->P();
594  he.FlipF();he.FlipE();he.FlipV();
595  vur=&he.v->P();
596  he.FlipV();he.FlipE();he.FlipF(); assert(&he.v->P()==vu); // back to vu (on the right)
597  he.FlipE();
598  he.FlipF();he.FlipE();he.FlipV();
599  vul=&he.v->P();
600  he.FlipV();he.FlipE();he.FlipF(); assert(&he.v->P()==vu); // back to vu (on the left)
601  he.FlipV();he.FlipE();he.FlipF(); assert(&he.v->P()==vl);// again on vl (but under the edge)
602  he.FlipE();he.FlipV();
603  vd=&he.v->P();
604  he.FlipF();he.FlipE();he.FlipV();
605  vdl=&he.v->P();
606  he.FlipV();he.FlipE();he.FlipF(); assert(&he.v->P()==vd);// back to vd (on the right)
607  he.FlipE();
608  he.FlipF();he.FlipE();he.FlipV();
609  vdr=&he.v->P();
610 
611  nv.P()=((*vl)+(*vr))/2.0+((*vu)+(*vd))/8.0 - ((*vul)+(*vur)+(*vdl)+(*vdr))/16.0;
612  }
613  }
614 
616  Color4<typename MESH_TYPE::ScalarType> WedgeInterp(Color4<typename MESH_TYPE::ScalarType> &c0, Color4<typename MESH_TYPE::ScalarType> &c1)
617  {
618  Color4<typename MESH_TYPE::ScalarType> cc;
619  return cc.lerp(c0,c1,0.5f);
620  }
621 
622  template<class FL_TYPE>
623  TexCoord2<FL_TYPE,1> WedgeInterp(TexCoord2<FL_TYPE,1> &t0, TexCoord2<FL_TYPE,1> &t1)
624  {
625  TexCoord2<FL_TYPE,1> tmp;
626  assert(t0.n()== t1.n());
627  tmp.n()=t0.n();
628  tmp.t()=(t0.t()+t1.t())/2.0;
629  return tmp;
630  }
631 };
632 
633 
634 #if 0
635  int rule=0;
636  if(vr==vul) rule+=1;
637  if(vl==vur) rule+=2;
638  if(vl==vdr) rule+=4;
639  if(vr==vdl) rule+=8;
640  switch(rule){
641 /* */
642 /* */ case 0 : return ((*vl)+(*vr))/2.0+((*vu)+(*vd))/8.0 - ((*vul)+(*vur)+(*vdl)+(*vdr))/16.0;
643 /* ul */ case 1 : return (*vl*6 + *vr*10 + *vu + *vd*3 - *vur - *vdl -*vdr*2 )/16.0;
644 /* ur */ case 2 : return (*vr*6 + *vl*10 + *vu + *vd*3 - *vul - *vdr -*vdl*2 )/16.0;
645 /* dr */ case 4 : return (*vr*6 + *vl*10 + *vd + *vu*3 - *vdl - *vur -*vul*2 )/16.0;
646 /* dl */ case 8 : return (*vl*6 + *vr*10 + *vd + *vu*3 - *vdr - *vul -*vur*2 )/16.0;
647 /* ul,ur */ case 3 : return (*vl*4 + *vr*4 + *vd*2 + - *vdr - *vdl )/8.0;
648 /* dl,dr */ case 12 : return (*vl*4 + *vr*4 + *vu*2 + - *vur - *vul )/8.0;
649 
650 /* ul,dr */ case 5 :
651 /* ur,dl */ case 10 :
652  default:
653  return (*vl+ *vr)/2.0;
654  }
655 
656 
657 
658 #endif
659 /*
660  vul-------vu--------vur
661  \ / \ /
662  \ / \ /
663  \ / fu \ /
664  vl--------vr
665  / \ fd / \
666  / \ / \
667  / \ / \
668  vdl-------vd--------vdr
669 
670 */
671 
672 // Versione modificata per tenere di conto in meniara corretta dei vertici con valenza alta
673 
674 template<class MESH_TYPE>
675 struct MidPointButterfly2 : public std::unary_function<face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType>
676 {
677  typename MESH_TYPE::CoordType operator()(face::Pos<typename MESH_TYPE::FaceType> ep)
678  {
679 double Rules[11][10] =
680 {
681  {.0}, // valenza 0
682  {.0}, // valenza 1
683  {.0}, // valenza 2
684  { .4166666667, -.08333333333 , -.08333333333 }, // valenza 3
685  { .375 , .0 , -0.125 , .0 }, // valenza 4
686  { .35 , .03090169945 , -.08090169945 , -.08090169945, .03090169945 }, // valenza 5
687  { .5 , .125 , -0.0625 , .0 , -0.0625 , 0.125 }, // valenza 6
688  { .25 , .1088899050 , -.06042933822 , -.04846056675, -.04846056675, -.06042933822, .1088899050 }, // valenza 7
689  { .21875 , .1196383476 , -.03125 , -.05713834763, -.03125 , -.05713834763, -.03125 ,.1196383476 }, // valenza 8
690  { .1944444444, .1225409480 , -.00513312590 , -.05555555556, -.03407448880, -.03407448880, -.05555555556, -.00513312590, .1225409480 }, // valenza 9
691  { .175 , .1213525492 , .01545084973 , -.04635254918, -.04045084973, -.025 , -.04045084973, -.04635254918, .01545084973, .1213525492 } // valenza 10
692 };
693 
694 face::Pos<typename MESH_TYPE::FaceType> he(ep.f,ep.z,ep.f->V(ep.z));
695  typename MESH_TYPE::CoordType *vl,*vr;
696  vl=&he.v->P();
697  vr=&he.VFlip()->P();
698  if(he.IsBorder())
699  {he.FlipV();
700  typename MESH_TYPE::CoordType *vl0,*vr0;
701  he.NextB();
702  vr0=&he.v->P();
703  he.FlipV();
704  he.NextB();
705  assert(&he.v->P()==vl);
706  he.NextB();
707  vl0=&he.v->P();
708  return ((*vl)+(*vr))*(9.0/16.0)-((*vl0)+(*vr0))/16.0 ;
709  }
710 
711  int kl=0,kr=0; // valence of left and right vertices
712  bool bl=false,br=false; // if left and right vertices are of border
713  face::Pos<typename MESH_TYPE::FaceType> heStart=he;assert(he.v->P()==*vl);
714  do { // compute valence of left vertex
715  he.FlipE();he.FlipF();
716  if(he.IsBorder()) bl=true;
717  ++kl;
718  } while(he!=heStart);
719 
720  he.FlipV();heStart=he;assert(he.v->P()==*vr);
721  do { // compute valence of right vertex
722  he.FlipE();he.FlipF();
723  if(he.IsBorder()) br=true;
724  ++kr;
725  } while(he!=heStart);
726  if(br||bl) return MidPointButterfly<MESH_TYPE>()( ep );
727  if(kr==6 && kl==6) return MidPointButterfly<MESH_TYPE>()( ep );
728  // TRACE("odd vertex among valences of %i %i\n",kl,kr);
729  typename MESH_TYPE::CoordType newposl=*vl*.75, newposr=*vr*.75;
730  he.FlipV();heStart=he; assert(he.v->P()==*vl);
731  int i=0;
732  if(kl!=6)
733  do { // compute position of left vertex
734  newposl+= he.VFlip()->P() * Rules[kl][i];
735  he.FlipE();he.FlipF();
736  ++i;
737  } while(he!=heStart);
738  i=0;he.FlipV();heStart=he;assert(he.v->P()==*vr);
739  if(kr!=6)
740  do { // compute position of right vertex
741  newposr+=he.VFlip()->P()* Rules[kr][i];
742  he.FlipE();he.FlipF();
743  ++i;
744  } while(he!=heStart);
745  if(kr==6) return newposl;
746  if(kl==6) return newposr;
747  return newposl+newposr;
748  }
749 };
750 
751 /* The two following classes are the functor and the predicate that you need for using the refine framework to cut a mesh along a linear interpolation of the quality.
752  This can be used for example to slice a mesh with a plane. Just set the quality value as distance from plane and then functor and predicate
753  initialized 0 and invoke the refine
754 
755  MyMesh A;
756  tri::UpdateQuality:MyMesh>::VertexFromPlane(plane);
757  QualityMidPointFunctor<MyMesh> slicingfunc(0.0);
758  QualityEdgePredicate<MyMesh> slicingpred(0.0);
759  tri::UpdateTopology<MyMesh>::FaceFace(A);
760  RefineE<MyMesh, QualityMidPointFunctor<MyMesh>, QualityEdgePredicate<MyMesh> > (A, slicingfunc, slicingpred, false);
761 
762  Note that they store in the vertex quality the plane distance.
763  */
764 
765 template<class MESH_TYPE>
766 class QualityMidPointFunctor : public std::unary_function<face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType>
767 {
768 public:
769  typedef Point3<typename MESH_TYPE::ScalarType> Point3x;
770  typedef typename MESH_TYPE::ScalarType ScalarType;
771 
772  ScalarType thr;
773 
774  QualityMidPointFunctor(ScalarType _thr):thr(_thr){}
775 
776 
777  void operator()(typename MESH_TYPE::VertexType &nv, const face::Pos<typename MESH_TYPE::FaceType> &ep){
778  Point3x p0=ep.f->V0(ep.z)->P();
779  Point3x p1=ep.f->V1(ep.z)->P();
780  ScalarType q0=ep.f->V0(ep.z)->Q()-thr;
781  ScalarType q1=ep.f->V1(ep.z)->Q()-thr;
782  double pp= q0/(q0-q1);
783  nv.P()=p1*pp + p0*(1.0-pp);
784  nv.Q()=thr;
785  }
786 
787  Color4<typename MESH_TYPE::ScalarType> WedgeInterp(Color4<typename MESH_TYPE::ScalarType> &c0, Color4<typename MESH_TYPE::ScalarType> &c1)
788  {
789  Color4<typename MESH_TYPE::ScalarType> cc;
790  return cc.lerp(c0,c1,0.5f);
791  }
792 
793  template<class FL_TYPE>
794  TexCoord2<FL_TYPE,1> WedgeInterp(TexCoord2<FL_TYPE,1> &t0, TexCoord2<FL_TYPE,1> &t1)
795  {
796  TexCoord2<FL_TYPE,1> tmp;
797  assert(t0.n()== t1.n());
798  tmp.n()=t0.n();
799  tmp.t()=(t0.t()+t1.t())/2.0;
800  return tmp;
801  }
802 };
803 
804 
805 template <class MESH_TYPE>
806 class QualityEdgePredicate
807 {
808  public:
809  typedef Point3<typename MESH_TYPE::ScalarType> Point3x;
810  typedef typename MESH_TYPE::ScalarType ScalarType;
811  ScalarType thr;
812  ScalarType tolerance;
813  QualityEdgePredicate(const ScalarType &thr,ScalarType _tolerance=0.02):thr(thr) {tolerance=_tolerance;}
814  bool operator()(face::Pos<typename MESH_TYPE::FaceType> ep)
815  {
816  ScalarType q0=ep.f->V0(ep.z)->Q()-thr;
817  ScalarType q1=ep.f->V1(ep.z)->Q()-thr;
818  if(q0>q1) std::swap(q0,q1);
819  if ( q0*q1 >= 0) return false;
820  // now a small check to be sure that we do not make too thin crossing.
821  double pp= q0/(q0-q1);
822  if ((fabs(pp)< tolerance)||(fabs(pp)> (1-tolerance))) return false;
823  return true;
824  }
825 };
826 
827 
828 template<class MESH_TYPE>
829 struct MidPointSphere : public std::unary_function<face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType>
830 {
831  Sphere3<typename MESH_TYPE::ScalarType> sph;
832  typedef Point3<typename MESH_TYPE::ScalarType> Point3x;
833 
834  void operator()(typename MESH_TYPE::VertexType &nv, face::Pos<typename MESH_TYPE::FaceType> ep){
835  Point3x &p0=ep.f->V0(ep.z)->P();
836  Point3x &p1=ep.f->V1(ep.z)->P();
837  nv.P()= sph.c+((p0+p1)/2.0 - sph.c ).Normalize();
838  }
839 
840  Color4<typename MESH_TYPE::ScalarType> WedgeInterp(Color4<typename MESH_TYPE::ScalarType> &c0, Color4<typename MESH_TYPE::ScalarType> &c1)
841  {
842  Color4<typename MESH_TYPE::ScalarType> cc;
843  return cc.lerp(c0,c1,0.5f);
844  }
845 
846  template<class FL_TYPE>
847  TexCoord2<FL_TYPE,1> WedgeInterp(TexCoord2<FL_TYPE,1> &t0, TexCoord2<FL_TYPE,1> &t1)
848  {
849  TexCoord2<FL_TYPE,1> tmp;
850  assert(t0.n()== t1.n());
851  tmp.n()=t0.n();
852  tmp.t()=(t0.t()+t1.t())/2.0;
853  return tmp;
854  }
855 };
856 
857 
858 template <class FLT>
859 class EdgeSplSphere
860 {
861  public:
862  Sphere3<FLT> sph;
863  bool operator()(const Point3<FLT> &p0, const Point3<FLT> &p1) const
864  {
865  if(Distance(sph,p0)>0) {
866  if(Distance(sph,p1)>0) return false;
867  else return true;
868  }
869  else if(Distance(sph,p1)<=0) return false;
870  return true;
871  }
872 };
873 
874 template<class TRIMESH_TYPE>
875 struct CenterPointBarycenter : public std::unary_function<typename TRIMESH_TYPE::FacePointer, typename TRIMESH_TYPE::CoordType>
876 {
877  typename TRIMESH_TYPE::CoordType operator()(typename TRIMESH_TYPE::FacePointer f){
878  return vcg::Barycenter<typename TRIMESH_TYPE::FaceType>(*f);
879  }
880 };
881 
885 
886 template<class TRIMESH_TYPE, class CenterPoint=CenterPointBarycenter <TRIMESH_TYPE> >
887 class TriSplit
888 {
889 public:
890  typedef typename TRIMESH_TYPE::FaceType FaceType;
891  typedef typename TRIMESH_TYPE::VertexType VertexType;
892 
893  static void Apply(FaceType *f,
894  FaceType * f1,FaceType * f2,
895  VertexType * vB, CenterPoint Center)
896  {
897  vB->P() = Center(f);
898 
899  //i tre vertici della faccia da dividere
900  VertexType *V0,*V1,*V2;
901  V0 = f->V(0);
902  V1 = f->V(1);
903  V2 = f->V(2);
904 
905  //risistemo la faccia di partenza
906  (*f).V(2) = &(*vB);
907  //Faccia nuova #1
908  (*f1).V(0) = &(*vB);
909  (*f1).V(1) = V1;
910  (*f1).V(2) = V2;
911  //Faccia nuova #2
912  (*f2).V(0) = V0;
913  (*f2).V(1) = &(*vB);
914  (*f2).V(2) = V2;
915 
916  if(f->HasFFAdjacency())
917  {
918  //adiacenza delle facce adiacenti a quelle aggiunte
919  f->FFp(1)->FFp(f->FFi(1)) = f1;
920  f->FFp(2)->FFp(f->FFi(2)) = f2;
921 
922  //adiacenza ff
923  FaceType * FF0,*FF1,*FF2;
924  FF0 = f->FFp(0);
925  FF1 = f->FFp(1);
926  FF2 = f->FFp(2);
927 
928  //Indici di adiacenza ff
929  char FFi0,FFi1,FFi2;
930  FFi0 = f->FFi(0);
931  FFi1 = f->FFi(1);
932  FFi2 = f->FFi(2);
933 
934  //adiacenza della faccia di partenza
935  (*f).FFp(1) = &(*f1);
936  (*f).FFi(1) = 0;
937  (*f).FFp(2) = &(*f2);
938  (*f).FFi(2) = 0;
939 
940  //adiacenza della faccia #1
941  (*f1).FFp(0) = f;
942  (*f1).FFi(0) = 1;
943 
944  (*f1).FFp(1) = FF1;
945  (*f1).FFi(1) = FFi1;
946 
947  (*f1).FFp(2) = &(*f2);
948  (*f1).FFi(2) = 1;
949 
950  //adiacenza della faccia #2
951  (*f2).FFp(0) = f;
952  (*f2).FFi(0) = 2;
953 
954  (*f2).FFp(1) = &(*f1);
955  (*f2).FFi(1) = 2;
956 
957  (*f2).FFp(2) = FF2;
958  (*f2).FFi(2) = FFi2;
959  }
960  }
961 }; // end class TriSplit
962 
963 template <class MeshType>
964 void TrivialMidPointRefine(MeshType & m)
965 {
966  typedef typename MeshType::VertexIterator VertexIterator;
967  typedef typename MeshType::FaceIterator FaceIterator;
968  typedef typename MeshType::VertexPointer VertexPointer;
969  typedef typename MeshType::FacePointer FacePointer;
970 
972  int startFn = m.fn;
973  FaceIterator lastf = tri::Allocator<MeshType>::AddFaces(m,m.fn*3);
974  VertexIterator lastv = tri::Allocator<MeshType>::AddVertices(m,m.fn*3);
975 
976  /*
977  * v0
978  * / \
979  * / f0 \
980  * / \
981  * mp01----------mp02
982  * / \ f3 / \
983  * / f1 \ / f2 \
984  * / \ / \
985  *v1 ---------- mp12------------v2
986  *
987  */
988 
989  for(int i=0;i<startFn;++i)
990  {
991  FacePointer f0= &m.face[i];
992  FacePointer f1= &*lastf; ++lastf;
993  FacePointer f2= &*lastf; ++lastf;
994  FacePointer f3= &*lastf; ++lastf;
995  VertexPointer v0 =m.face[i].V(0);
996  VertexPointer v1 =m.face[i].V(1);
997  VertexPointer v2 =m.face[i].V(2);
998  VertexPointer mp01 = &*lastv; ++lastv;
999  VertexPointer mp12 = &*lastv; ++lastv;
1000  VertexPointer mp02 = &*lastv; ++lastv;
1001 
1002  f0->V(0) = v0; f0->V(1) = mp01; f0->V(2) = mp02;
1003  f1->V(0) = v1; f1->V(1) = mp12; f1->V(2) = mp01;
1004  f2->V(0) = v2; f2->V(1) = mp02; f2->V(2) = mp12;
1005  f3->V(0) = mp12; f3->V(1) = mp02; f3->V(2) = mp01;
1006  mp01->P() = (v0>v1) ? (v0->P()+v1->P())/2.0 : (v1->P()+v0->P())/2.0;
1007  mp12->P() = (v1>v2) ? (v1->P()+v2->P())/2.0 : (v2->P()+v1->P())/2.0;
1008  mp02->P() = (v0>v2) ? (v0->P()+v2->P())/2.0 : (v2->P()+v0->P())/2.0;
1009  }
1010 
1012  printf("Vertex unification %i\n",vd);
1014  printf("Vertex unref %i\n",vu);
1015  Allocator<MeshType>::CompactEveryVector(m);
1016 }
1017 
1018 } // namespace tri
1019 } // namespace vcg
1020 
1021 #endif