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] = ep[1] = ep[2] = false;
320  vp[0] = vp[1] = vp[2] = NULL;
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;
447  for(fi=m.face.begin();fi!=oldendf;++fi) if(!(*fi).IsD())
448  {
449  if(cb && (++step%PercStep)==0)
450  (*cb)(step/PercStep,"Refining...");
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] != NULL) ? 1 : 0) + ((vv[4] != NULL) ? 2 : 0) + ((vv[5] != NULL) ? 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  {
474  wtt[i]=(*fi).WT(i);
475  wtt[3+i]=mid.WedgeInterp((*fi).WT(i),(*fi).WT((i+1)%3));
476  }
477 
478  int orgflag = (*fi).Flags();
479  for (i=0; i<SplitTab[ind].TriNum; ++i)
480  for(j=0;j<3;++j)
481  {
482  (*nf[i]).V(j)=&*vv[SplitTab[ind].TV[i][j]];
483 
484  if(tri::HasPerWedgeTexCoord(m)) //analogo ai vertici...
485  (*nf[i]).WT(j) = wtt[SplitTab[ind].TV[i][j]];
486 
487  assert((*nf[i]).V(j)!=0);
488  if(SplitTab[ind].TE[i][j]!=3)
489  {
490  if(orgflag & (MESH_TYPE::FaceType::BORDER0<<(SplitTab[ind].TE[i][j])))
491  (*nf[i]).SetB(j);
492  else
493  (*nf[i]).ClearB(j);
494 
495  if(orgflag & (MESH_TYPE::FaceType::FACEEDGESEL0<<(SplitTab[ind].TE[i][j])))
496  (*nf[i]).SetFaceEdgeS(j);
497  else
498  (*nf[i]).ClearFaceEdgeS(j);
499  }
500  else
501  {
502  (*nf[i]).ClearB(j);
503  (*nf[i]).ClearFaceEdgeS(j);
504  }
505  }
506 
507  if(SplitTab[ind].TriNum==3 &&
508  SquaredDistance(vv[SplitTab[ind].swap[0][0]]->P(),vv[SplitTab[ind].swap[0][1]]->P()) <
509  SquaredDistance(vv[SplitTab[ind].swap[1][0]]->P(),vv[SplitTab[ind].swap[1][1]]->P()) )
510  { // swap the last two triangles
511  (*nf[2]).V(1)=(*nf[1]).V(0);
512  (*nf[1]).V(1)=(*nf[2]).V(0);
513  if(tri::HasPerWedgeTexCoord(m)){ //swap also textures coordinates
514  (*nf[2]).WT(1)=(*nf[1]).WT(0);
515  (*nf[1]).WT(1)=(*nf[2]).WT(0);
516  }
517 
518  if((*nf[1]).IsB(0)) (*nf[2]).SetB(1); else (*nf[2]).ClearB(1);
519  if((*nf[2]).IsB(0)) (*nf[1]).SetB(1); else (*nf[1]).ClearB(1);
520  (*nf[1]).ClearB(0);
521  (*nf[2]).ClearB(0);
522 
523  if((*nf[1]).IsFaceEdgeS(0)) (*nf[2]).SetFaceEdgeS(1); else (*nf[2]).ClearFaceEdgeS(1);
524  if((*nf[2]).IsFaceEdgeS(0)) (*nf[1]).SetFaceEdgeS(1); else (*nf[1]).ClearFaceEdgeS(1);
525  (*nf[1]).ClearFaceEdgeS(0);
526  (*nf[2]).ClearFaceEdgeS(0);
527  }
528  }
529 
530  assert(lastf==m.face.end()); // critical assert: we MUST have used all the faces that we forecasted we need and that we previously allocated.
531  assert(!m.vert.empty());
532  for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()){
533  assert((*fi).V(0)>=&*m.vert.begin() && (*fi).V(0)<=&m.vert.back() );
534  assert((*fi).V(1)>=&*m.vert.begin() && (*fi).V(1)<=&m.vert.back() );
535  assert((*fi).V(2)>=&*m.vert.begin() && (*fi).V(2)<=&m.vert.back() );
536  }
538 
539  tri::Allocator<MESH_TYPE> :: template DeletePerFaceAttribute<RefinedFaceData<VertexPointer> > (m,RD);
540 
541  return true;
542 }
543 
544 /*************************************************************************/
545 // simple wrapper of the base refine for lazy coder that do not need a edge predicate
546 
547 template<class MESH_TYPE,class MIDPOINT>
548 bool Refine(MESH_TYPE &m, MIDPOINT mid, typename MESH_TYPE::ScalarType thr=0,bool RefineSelected=false, CallBackPos *cb = 0)
549 {
550  EdgeLen <MESH_TYPE, typename MESH_TYPE::ScalarType> ep(thr);
551  return RefineE(m,mid,ep,RefineSelected,cb);
552 }
553 /*************************************************************************/
554 
555 /*
556 Modified Butterfly interpolation scheme,
557 as presented in
558 Zorin, Schroeder
559 Subdivision for modeling and animation
560 Siggraph 2000 Course Notes
561 */
562 
563 /*
564 
565  vul-------vu--------vur
566  \ / \ /
567  \ / \ /
568  \ / fu \ /
569  vl--------vr
570  / \ fd / \
571  / \ / \
572  / \ / \
573  vdl-------vd--------vdr
574 
575 */
576 
577 template<class MESH_TYPE>
578 struct MidPointButterfly : public std::unary_function<face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType>
579 {
580  MESH_TYPE &m;
581  MidPointButterfly(MESH_TYPE &_m):m(_m){}
582 
583  void operator()(typename MESH_TYPE::VertexType &nv, face::Pos<typename MESH_TYPE::FaceType> ep)
584  {
585  face::Pos<typename MESH_TYPE::FaceType> he(ep.f,ep.z,ep.f->V(ep.z));
586  typename MESH_TYPE::CoordType *vl,*vr;
587  typename MESH_TYPE::CoordType *vl0,*vr0;
588  typename MESH_TYPE::CoordType *vu,*vd,*vul,*vur,*vdl,*vdr;
589  vl=&he.v->P();
590  he.FlipV();
591  vr=&he.v->P();
592 
593  if( tri::HasPerVertexColor(m))
594  nv.C().lerp(ep.f->V(ep.z)->C(),ep.f->V1(ep.z)->C(),.5f);
595 
596  if(he.IsBorder())
597  {
598  he.NextB();
599  vr0=&he.v->P();
600  he.FlipV();
601  he.NextB();
602  assert(&he.v->P()==vl);
603  he.NextB();
604  vl0=&he.v->P();
605  nv.P()=((*vl)+(*vr))*(9.0/16.0)-((*vl0)+(*vr0))/16.0 ;
606  }
607  else
608  {
609  he.FlipE();he.FlipV();
610  vu=&he.v->P();
611  he.FlipF();he.FlipE();he.FlipV();
612  vur=&he.v->P();
613  he.FlipV();he.FlipE();he.FlipF(); assert(&he.v->P()==vu); // back to vu (on the right)
614  he.FlipE();
615  he.FlipF();he.FlipE();he.FlipV();
616  vul=&he.v->P();
617  he.FlipV();he.FlipE();he.FlipF(); assert(&he.v->P()==vu); // back to vu (on the left)
618  he.FlipV();he.FlipE();he.FlipF(); assert(&he.v->P()==vl);// again on vl (but under the edge)
619  he.FlipE();he.FlipV();
620  vd=&he.v->P();
621  he.FlipF();he.FlipE();he.FlipV();
622  vdl=&he.v->P();
623  he.FlipV();he.FlipE();he.FlipF(); assert(&he.v->P()==vd);// back to vd (on the right)
624  he.FlipE();
625  he.FlipF();he.FlipE();he.FlipV();
626  vdr=&he.v->P();
627 
628  nv.P()=((*vl)+(*vr))/2.0+((*vu)+(*vd))/8.0 - ((*vul)+(*vur)+(*vdl)+(*vdr))/16.0;
629  }
630  }
631 
633  Color4<typename MESH_TYPE::ScalarType> WedgeInterp(Color4<typename MESH_TYPE::ScalarType> &c0, Color4<typename MESH_TYPE::ScalarType> &c1)
634  {
635  Color4<typename MESH_TYPE::ScalarType> cc;
636  return cc.lerp(c0,c1,0.5f);
637  }
638 
639  template<class FL_TYPE>
640  TexCoord2<FL_TYPE,1> WedgeInterp(TexCoord2<FL_TYPE,1> &t0, TexCoord2<FL_TYPE,1> &t1)
641  {
642  TexCoord2<FL_TYPE,1> tmp;
643  assert(t0.n()== t1.n());
644  tmp.n()=t0.n();
645  tmp.t()=(t0.t()+t1.t())/2.0;
646  return tmp;
647  }
648 };
649 
650 
651 #if 0
652  int rule=0;
653  if(vr==vul) rule+=1;
654  if(vl==vur) rule+=2;
655  if(vl==vdr) rule+=4;
656  if(vr==vdl) rule+=8;
657  switch(rule){
658 /* */
659 /* */ case 0 : return ((*vl)+(*vr))/2.0+((*vu)+(*vd))/8.0 - ((*vul)+(*vur)+(*vdl)+(*vdr))/16.0;
660 /* ul */ case 1 : return (*vl*6 + *vr*10 + *vu + *vd*3 - *vur - *vdl -*vdr*2 )/16.0;
661 /* ur */ case 2 : return (*vr*6 + *vl*10 + *vu + *vd*3 - *vul - *vdr -*vdl*2 )/16.0;
662 /* dr */ case 4 : return (*vr*6 + *vl*10 + *vd + *vu*3 - *vdl - *vur -*vul*2 )/16.0;
663 /* dl */ case 8 : return (*vl*6 + *vr*10 + *vd + *vu*3 - *vdr - *vul -*vur*2 )/16.0;
664 /* ul,ur */ case 3 : return (*vl*4 + *vr*4 + *vd*2 + - *vdr - *vdl )/8.0;
665 /* dl,dr */ case 12 : return (*vl*4 + *vr*4 + *vu*2 + - *vur - *vul )/8.0;
666 
667 /* ul,dr */ case 5 :
668 /* ur,dl */ case 10 :
669  default:
670  return (*vl+ *vr)/2.0;
671  }
672 
673 
674 
675 #endif
676 /*
677  vul-------vu--------vur
678  \ / \ /
679  \ / \ /
680  \ / fu \ /
681  vl--------vr
682  / \ fd / \
683  / \ / \
684  / \ / \
685  vdl-------vd--------vdr
686 
687 */
688 
689 // Versione modificata per tenere di conto in meniara corretta dei vertici con valenza alta
690 
691 template<class MESH_TYPE>
692 struct MidPointButterfly2 : public std::unary_function<face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType>
693 {
694  typename MESH_TYPE::CoordType operator()(face::Pos<typename MESH_TYPE::FaceType> ep)
695  {
696 double Rules[11][10] =
697 {
698  {.0}, // valenza 0
699  {.0}, // valenza 1
700  {.0}, // valenza 2
701  { .4166666667, -.08333333333 , -.08333333333 }, // valenza 3
702  { .375 , .0 , -0.125 , .0 }, // valenza 4
703  { .35 , .03090169945 , -.08090169945 , -.08090169945, .03090169945 }, // valenza 5
704  { .5 , .125 , -0.0625 , .0 , -0.0625 , 0.125 }, // valenza 6
705  { .25 , .1088899050 , -.06042933822 , -.04846056675, -.04846056675, -.06042933822, .1088899050 }, // valenza 7
706  { .21875 , .1196383476 , -.03125 , -.05713834763, -.03125 , -.05713834763, -.03125 ,.1196383476 }, // valenza 8
707  { .1944444444, .1225409480 , -.00513312590 , -.05555555556, -.03407448880, -.03407448880, -.05555555556, -.00513312590, .1225409480 }, // valenza 9
708  { .175 , .1213525492 , .01545084973 , -.04635254918, -.04045084973, -.025 , -.04045084973, -.04635254918, .01545084973, .1213525492 } // valenza 10
709 };
710 
711 face::Pos<typename MESH_TYPE::FaceType> he(ep.f,ep.z,ep.f->V(ep.z));
712  typename MESH_TYPE::CoordType *vl,*vr;
713  vl=&he.v->P();
714  vr=&he.VFlip()->P();
715  if(he.IsBorder())
716  {he.FlipV();
717  typename MESH_TYPE::CoordType *vl0,*vr0;
718  he.NextB();
719  vr0=&he.v->P();
720  he.FlipV();
721  he.NextB();
722  assert(&he.v->P()==vl);
723  he.NextB();
724  vl0=&he.v->P();
725  return ((*vl)+(*vr))*(9.0/16.0)-((*vl0)+(*vr0))/16.0 ;
726  }
727 
728  int kl=0,kr=0; // valence of left and right vertices
729  bool bl=false,br=false; // if left and right vertices are of border
730  face::Pos<typename MESH_TYPE::FaceType> heStart=he;assert(he.v->P()==*vl);
731  do { // compute valence of left vertex
732  he.FlipE();he.FlipF();
733  if(he.IsBorder()) bl=true;
734  ++kl;
735  } while(he!=heStart);
736 
737  he.FlipV();heStart=he;assert(he.v->P()==*vr);
738  do { // compute valence of right vertex
739  he.FlipE();he.FlipF();
740  if(he.IsBorder()) br=true;
741  ++kr;
742  } while(he!=heStart);
743  if(br||bl) return MidPointButterfly<MESH_TYPE>()( ep );
744  if(kr==6 && kl==6) return MidPointButterfly<MESH_TYPE>()( ep );
745  // TRACE("odd vertex among valences of %i %i\n",kl,kr);
746  typename MESH_TYPE::CoordType newposl=*vl*.75, newposr=*vr*.75;
747  he.FlipV();heStart=he; assert(he.v->P()==*vl);
748  int i=0;
749  if(kl!=6)
750  do { // compute position of left vertex
751  newposl+= he.VFlip()->P() * Rules[kl][i];
752  he.FlipE();he.FlipF();
753  ++i;
754  } while(he!=heStart);
755  i=0;he.FlipV();heStart=he;assert(he.v->P()==*vr);
756  if(kr!=6)
757  do { // compute position of right vertex
758  newposr+=he.VFlip()->P()* Rules[kr][i];
759  he.FlipE();he.FlipF();
760  ++i;
761  } while(he!=heStart);
762  if(kr==6) return newposl;
763  if(kl==6) return newposr;
764  return newposl+newposr;
765  }
766 };
767 
768 /* 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.
769  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
770  initialized 0 and invoke the refine
771 
772  MyMesh A;
773  tri::UpdateQuality:MyMesh>::VertexFromPlane(plane);
774  QualityMidPointFunctor<MyMesh> slicingfunc(0.0);
775  QualityEdgePredicate<MyMesh> slicingpred(0.0);
776  tri::UpdateTopology<MyMesh>::FaceFace(A);
777  RefineE<MyMesh, QualityMidPointFunctor<MyMesh>, QualityEdgePredicate<MyMesh> > (A, slicingfunc, slicingpred, false);
778 
779  Note that they store in the vertex quality the plane distance.
780  */
781 
782 template<class MESH_TYPE>
783 class QualityMidPointFunctor : public std::unary_function<face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType>
784 {
785 public:
786  typedef Point3<typename MESH_TYPE::ScalarType> Point3x;
787  typedef typename MESH_TYPE::ScalarType ScalarType;
788 
789  ScalarType thr;
790 
791  QualityMidPointFunctor(ScalarType _thr):thr(_thr){}
792 
793 
794  void operator()(typename MESH_TYPE::VertexType &nv, const face::Pos<typename MESH_TYPE::FaceType> &ep){
795  Point3x p0=ep.f->V0(ep.z)->P();
796  Point3x p1=ep.f->V1(ep.z)->P();
797  ScalarType q0=ep.f->V0(ep.z)->Q()-thr;
798  ScalarType q1=ep.f->V1(ep.z)->Q()-thr;
799  double pp= q0/(q0-q1);
800  nv.P()=p1*pp + p0*(1.0-pp);
801  nv.Q()=thr;
802  }
803 
804  Color4<typename MESH_TYPE::ScalarType> WedgeInterp(Color4<typename MESH_TYPE::ScalarType> &c0, Color4<typename MESH_TYPE::ScalarType> &c1)
805  {
806  Color4<typename MESH_TYPE::ScalarType> cc;
807  return cc.lerp(c0,c1,0.5f);
808  }
809 
810  template<class FL_TYPE>
811  TexCoord2<FL_TYPE,1> WedgeInterp(TexCoord2<FL_TYPE,1> &t0, TexCoord2<FL_TYPE,1> &t1)
812  {
813  TexCoord2<FL_TYPE,1> tmp;
814  assert(t0.n()== t1.n());
815  tmp.n()=t0.n();
816  tmp.t()=(t0.t()+t1.t())/2.0;
817  return tmp;
818  }
819 };
820 
821 
822 template <class MESH_TYPE>
823 class QualityEdgePredicate
824 {
825  public:
826  typedef Point3<typename MESH_TYPE::ScalarType> Point3x;
827  typedef typename MESH_TYPE::ScalarType ScalarType;
828  ScalarType thr;
829  ScalarType tolerance;
830  QualityEdgePredicate(const ScalarType &thr,ScalarType _tolerance=0.02):thr(thr) {tolerance=_tolerance;}
831  bool operator()(face::Pos<typename MESH_TYPE::FaceType> ep)
832  {
833  ScalarType q0=ep.f->V0(ep.z)->Q()-thr;
834  ScalarType q1=ep.f->V1(ep.z)->Q()-thr;
835  if(q0>q1) std::swap(q0,q1);
836  if ( q0*q1 >= 0) return false;
837  // now a small check to be sure that we do not make too thin crossing.
838  double pp= q0/(q0-q1);
839  if ((fabs(pp)< tolerance)||(fabs(pp)> (1-tolerance))) return false;
840  return true;
841  }
842 };
843 
844 
845 template<class MESH_TYPE>
846 struct MidPointSphere : public std::unary_function<face::Pos<typename MESH_TYPE::FaceType> , typename MESH_TYPE::CoordType>
847 {
848  Sphere3<typename MESH_TYPE::ScalarType> sph;
849  typedef Point3<typename MESH_TYPE::ScalarType> Point3x;
850 
851  void operator()(typename MESH_TYPE::VertexType &nv, face::Pos<typename MESH_TYPE::FaceType> ep){
852  Point3x &p0=ep.f->V0(ep.z)->P();
853  Point3x &p1=ep.f->V1(ep.z)->P();
854  nv.P()= sph.c+((p0+p1)/2.0 - sph.c ).Normalize();
855  }
856 
857  Color4<typename MESH_TYPE::ScalarType> WedgeInterp(Color4<typename MESH_TYPE::ScalarType> &c0, Color4<typename MESH_TYPE::ScalarType> &c1)
858  {
859  Color4<typename MESH_TYPE::ScalarType> cc;
860  return cc.lerp(c0,c1,0.5f);
861  }
862 
863  template<class FL_TYPE>
864  TexCoord2<FL_TYPE,1> WedgeInterp(TexCoord2<FL_TYPE,1> &t0, TexCoord2<FL_TYPE,1> &t1)
865  {
866  TexCoord2<FL_TYPE,1> tmp;
867  assert(t0.n()== t1.n());
868  tmp.n()=t0.n();
869  tmp.t()=(t0.t()+t1.t())/2.0;
870  return tmp;
871  }
872 };
873 
874 
875 template <class FLT>
876 class EdgeSplSphere
877 {
878  public:
879  Sphere3<FLT> sph;
880  bool operator()(const Point3<FLT> &p0, const Point3<FLT> &p1) const
881  {
882  if(Distance(sph,p0)>0) {
883  if(Distance(sph,p1)>0) return false;
884  else return true;
885  }
886  else if(Distance(sph,p1)<=0) return false;
887  return true;
888  }
889 };
890 
891 template<class TRIMESH_TYPE>
892 struct CenterPointBarycenter : public std::unary_function<typename TRIMESH_TYPE::FacePointer, typename TRIMESH_TYPE::CoordType>
893 {
894  typename TRIMESH_TYPE::CoordType operator()(typename TRIMESH_TYPE::FacePointer f){
895  return vcg::Barycenter<typename TRIMESH_TYPE::FaceType>(*f);
896  }
897 };
898 
902 
903 template<class TRIMESH_TYPE, class CenterPoint=CenterPointBarycenter <TRIMESH_TYPE> >
904 class TriSplit
905 {
906 public:
907  typedef typename TRIMESH_TYPE::FaceType FaceType;
908  typedef typename TRIMESH_TYPE::VertexType VertexType;
909 
910  static void Apply(FaceType *f,
911  FaceType * f1,FaceType * f2,
912  VertexType * vB, CenterPoint Center)
913  {
914  vB->P() = Center(f);
915 
916  //i tre vertici della faccia da dividere
917  VertexType *V0,*V1,*V2;
918  V0 = f->V(0);
919  V1 = f->V(1);
920  V2 = f->V(2);
921 
922  //risistemo la faccia di partenza
923  (*f).V(2) = &(*vB);
924  //Faccia nuova #1
925  (*f1).V(0) = &(*vB);
926  (*f1).V(1) = V1;
927  (*f1).V(2) = V2;
928  //Faccia nuova #2
929  (*f2).V(0) = V0;
930  (*f2).V(1) = &(*vB);
931  (*f2).V(2) = V2;
932 
933  if(f->HasFFAdjacency())
934  {
935  //adiacenza delle facce adiacenti a quelle aggiunte
936  f->FFp(1)->FFp(f->FFi(1)) = f1;
937  f->FFp(2)->FFp(f->FFi(2)) = f2;
938 
939  //adiacenza ff
940  FaceType * FF0,*FF1,*FF2;
941  FF0 = f->FFp(0);
942  FF1 = f->FFp(1);
943  FF2 = f->FFp(2);
944 
945  //Indici di adiacenza ff
946  char FFi0,FFi1,FFi2;
947  FFi0 = f->FFi(0);
948  FFi1 = f->FFi(1);
949  FFi2 = f->FFi(2);
950 
951  //adiacenza della faccia di partenza
952  (*f).FFp(1) = &(*f1);
953  (*f).FFi(1) = 0;
954  (*f).FFp(2) = &(*f2);
955  (*f).FFi(2) = 0;
956 
957  //adiacenza della faccia #1
958  (*f1).FFp(0) = f;
959  (*f1).FFi(0) = 1;
960 
961  (*f1).FFp(1) = FF1;
962  (*f1).FFi(1) = FFi1;
963 
964  (*f1).FFp(2) = &(*f2);
965  (*f1).FFi(2) = 1;
966 
967  //adiacenza della faccia #2
968  (*f2).FFp(0) = f;
969  (*f2).FFi(0) = 2;
970 
971  (*f2).FFp(1) = &(*f1);
972  (*f2).FFi(1) = 2;
973 
974  (*f2).FFp(2) = FF2;
975  (*f2).FFi(2) = FFi2;
976  }
977  }
978 }; // end class TriSplit
979 
980 template <class MeshType>
981 void TrivialMidPointRefine(MeshType & m)
982 {
983  typedef typename MeshType::VertexIterator VertexIterator;
984  typedef typename MeshType::FaceIterator FaceIterator;
985  typedef typename MeshType::VertexPointer VertexPointer;
986  typedef typename MeshType::FacePointer FacePointer;
987 
989  int startFn = m.fn;
990  FaceIterator lastf = tri::Allocator<MeshType>::AddFaces(m,m.fn*3);
991  VertexIterator lastv = tri::Allocator<MeshType>::AddVertices(m,m.fn*3);
992 
993  /*
994  * v0
995  * / \
996  * / f0 \
997  * / \
998  * mp01----------mp02
999  * / \ f3 / \
1000  * / f1 \ / f2 \
1001  * / \ / \
1002  *v1 ---------- mp12------------v2
1003  *
1004  */
1005 
1006  for(int i=0;i<startFn;++i)
1007  {
1008  FacePointer f0= &m.face[i];
1009  FacePointer f1= &*lastf; ++lastf;
1010  FacePointer f2= &*lastf; ++lastf;
1011  FacePointer f3= &*lastf; ++lastf;
1012  VertexPointer v0 =m.face[i].V(0);
1013  VertexPointer v1 =m.face[i].V(1);
1014  VertexPointer v2 =m.face[i].V(2);
1015  VertexPointer mp01 = &*lastv; ++lastv;
1016  VertexPointer mp12 = &*lastv; ++lastv;
1017  VertexPointer mp02 = &*lastv; ++lastv;
1018 
1019  f0->V(0) = v0; f0->V(1) = mp01; f0->V(2) = mp02;
1020  f1->V(0) = v1; f1->V(1) = mp12; f1->V(2) = mp01;
1021  f2->V(0) = v2; f2->V(1) = mp02; f2->V(2) = mp12;
1022  f3->V(0) = mp12; f3->V(1) = mp02; f3->V(2) = mp01;
1023  mp01->P() = (v0>v1) ? (v0->P()+v1->P())/2.0 : (v1->P()+v0->P())/2.0;
1024  mp12->P() = (v1>v2) ? (v1->P()+v2->P())/2.0 : (v2->P()+v1->P())/2.0;
1025  mp02->P() = (v0>v2) ? (v0->P()+v2->P())/2.0 : (v2->P()+v0->P())/2.0;
1026  }
1027 
1029  printf("Vertex unification %i\n",vd);
1031  printf("Vertex unref %i\n",vu);
1032  Allocator<MeshType>::CompactEveryVector(m);
1033 }
1034 
1035 } // namespace tri
1036 } // namespace vcg
1037 
1038 #endif