From d231b9d0215f7217d4ce7bacc879773e346c081d Mon Sep 17 00:00:00 2001
From: mtarini <marco.tarini@isti.cnr.it>
Date: Fri, 28 Aug 2009 15:17:23 +0000
Subject: [PATCH] Encapsulated everything in a static class. Also, templated
 with Interpolator "single-method static class" functor to make custom vertex
 interpolations during collapses.

---
 vcg/complex/trimesh/bitquad_creation.h     | 181 ++++---
 vcg/complex/trimesh/bitquad_optimization.h | 160 +++---
 vcg/complex/trimesh/bitquad_support.h      | 548 +++++++++++----------
 3 files changed, 430 insertions(+), 459 deletions(-)

diff --git a/vcg/complex/trimesh/bitquad_creation.h b/vcg/complex/trimesh/bitquad_creation.h
index 80882172..d012c0a7 100644
--- a/vcg/complex/trimesh/bitquad_creation.h
+++ b/vcg/complex/trimesh/bitquad_creation.h
@@ -3,7 +3,7 @@
 
 /** BIT-QUAD creation support:
     a collection of methods that,
-    starting from a triangular mesh, will create your quad-only or quad-domainant mesh.
+    starting from a triangular mesh, will create your quad-pure or quad-domainant mesh.
     
     They all require:
       - per face Q, and FF connectivity, 2-manyfold meshes, 
@@ -12,37 +12,37 @@
 
 [ list of available methods: ]
   
-void MakeBitQuadOnlyByRefine(Mesh &m) 
+void MakePureByRefine(Mesh &m) 
    - adds a vertex for each tri or quad present
    - thus, miminal complexity increase is the mesh is quad-dominant already
    - old non-border edges are made faux
    - never fails
 
-void MakeBitQuadOnlyByCatmullClark(Mesh &m) 
+void MakePureByCatmullClark(MeshType &m) 
    - adds a vertex in each (non-faux) edge. 
    - twice complexity increase w.r.t. "ByRefine" method.
    - preserves edges: old edges are still edges 
    - never fails
       
-bool MakeBitQuadOnlyByFlip(Mesh &m [, int maxdist] )  
+bool MakePureByFlip(MeshType &m [, int maxdist] )  
    - does not increase # vertices, just flips edges
    - call in a loop until it returns true  (temporary hack)
    - fails if number of triangle is odd (only happens in open meshes)
    - add "StepByStep" to method name if you want it to make a single step (debugging purposes)
 
-bool MakeTriEvenBySplit(Mesh& m)
-bool MakeTriEvenByDelete(Mesh& m)
+bool MakeTriEvenBySplit(MeshType& m)
+bool MakeTriEvenByDelete(MeshType& m)
    - two simple variants that either delete or split *at most one* border face
      so that the number of tris will be made even. Return true if it did it.
    - useful to use the previous method, when mesh is still all triangle
 
-void MakeBitQuadDominant(Mesh &m, int level)
+void MakeDominant(MeshType &m, int level)
    - just merges traingle pairs into quads, trying its best
    - various heuristic available, see descr. for parameter "level"
    - provides good starting point for make-Quad-Only methods
    - uses an ad-hoc measure for "quad quality" (which is hard-wired, for now)
 
-void MakeBitTriOnly(Mesh &m)
+void MakeBitTriOnly(MeshType &m)
    - inverse process: returns to tri-only mesh
 
    
@@ -52,14 +52,28 @@ void MakeBitTriOnly(Mesh &m)
 
 namespace vcg{namespace tri{
 
+template <class _MeshType>
+class BitQuadCreation{
+  
+public:
+  
+typedef _MeshType MeshType;
+typedef typename MeshType::ScalarType ScalarType;
+typedef typename MeshType::CoordType CoordType;
+typedef typename MeshType::FaceType FaceType;
+typedef typename MeshType::FaceType* FaceTypeP;
+typedef typename MeshType::VertexType VertexType;
+typedef typename MeshType::FaceIterator FaceIterator;
+typedef typename MeshType::VertexIterator VertexIterator;
 
+typedef BitQuad<MeshType> BQ; // static class to make basic quad operations
+
+  
 // helper function:
 // given a triangle, merge it with its best neightboord to form a quad
-template <class Face, bool override>
-static void selectBestQuadDiag(Face *fi){
+template <bool override>
+static void selectBestDiag(FaceType *fi){
   
-  typedef typename Face::ScalarType ScalarType;
-  typedef typename Face::VertexType VertexType;
   if (!override) {
     if (fi->IsAnyF()) return;
   }
@@ -79,7 +93,7 @@ static void selectBestQuadDiag(Face *fi){
     }
     if (fi->FFp(k)==fi) continue; // never make a border faux
     
-    ScalarType score = quadQuality( &*fi, k );
+    ScalarType score = BQ::quadQuality( &*fi, k );
     if (override) {
       // don't override anyway iff other face has a better match
       if (score < fi->FFp(k)->Q()) continue; 
@@ -129,33 +143,27 @@ static void selectBestQuadDiag(Face *fi){
 
 // helper funcion:
 // a pass though all triangles to merge triangle pairs into quads
-template <class Mesh, bool override> // override previous decisions?
-static void MakeQuadDominantPass(Mesh &m){
-  typedef typename Mesh::FaceType Face;
-  typedef typename Mesh::FaceIterator FaceIterator;
+template <bool override> // override previous decisions?
+static void MakeDominantPass(MeshType &m){
   
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) {
-    selectBestQuadDiag<Face,override>(&(*fi));
+    selectBestDiag<override>(&(*fi));
   }
   
 }
 
 // make tri count even by splitting a single triangle...
-template <class Mesh>
-bool MakeTriEvenBySplit(Mesh& m){
+static bool MakeTriEvenBySplit(MeshType& m){
   if (m.fn%2==0) return false; // it's already Even
   assert(0); // todo!
 }
 
 // make tri count even by delete...
-template <class Mesh>
-bool MakeTriEvenByDelete(Mesh& m)
+static bool MakeTriEvenByDelete(MeshType& m)
 {
   
   if (m.fn%2==0) return false; // it's already Even
   
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType Face;
   
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) {
     for (int k=0; k<3; k++) {
@@ -165,7 +173,7 @@ bool MakeTriEvenByDelete(Mesh& m)
         for (int h=1; h<3; h++) {
           int kh=(k+h)%3;
           int j = fi->FFi( kh );
-          Face *f = fi->FFp(kh);
+          FaceType *f = fi->FFp(kh);
           if (f != &* fi) {
             f->FFp( j ) = f;
             f->FFi( j ) = j;
@@ -174,7 +182,7 @@ bool MakeTriEvenByDelete(Mesh& m)
         }
         
         // delete found face
-        Allocator<Mesh>::DeleteFace(m,*fi);
+        Allocator<MeshType>::DeleteFace(m,*fi);
         return true;
       }
     }
@@ -187,9 +195,7 @@ bool MakeTriEvenByDelete(Mesh& m)
 /** 
   Given a mesh, makes it bit trianglular (makes all edges NOT faux)
 */
-template <class Mesh>
-void MakeBitTriOnly(Mesh &m){
-  typedef typename Mesh::FaceIterator FaceIterator;
+static void MakeBitTriOnly(MeshType &m){
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) {
     fi->ClearAllF();
   }
@@ -201,17 +207,13 @@ void MakeBitTriOnly(Mesh &m){
  * Updates: per wedge attributes, if any
  * Other connectivity structures, and per edge and per wedge flags are ignored
  */
-template <class Mesh>
-bool MakeBitTriQuadConventional(Mesh &m){
+static bool MakeBitTriQuadConventional(MeshType &m){
   assert(0); // todo
 }
 
 /* returns true if mesh is a "conventional" quad mesh.
    I.e. if it is all quads, with third edge faux fora all triangles*/
-template <class Mesh>
-bool IsBitTriQuadConventional(Mesh &m){
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
+static bool IsBitTriQuadConventional(MeshType &m){
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) {
     if (fi->IsAnyF())
     if ( fi->Flags() & ( FaceType::FAUX012 ) != FaceType::FAUX2 ) {
@@ -226,16 +228,10 @@ bool IsBitTriQuadConventional(Mesh &m){
  previous diags
  requires that the mesh is made only of quads and tris.
 */
-template <class Mesh> 
-void MakeBitQuadOnlyByRefine(Mesh &m){
+static void MakePureByRefine(MeshType &m){
   
   // todo: update VF connectivity if present
 
-
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::VertexIterator VertexIterator;
-  typedef typename Mesh::FaceType Face;
-  typedef typename Mesh::VertexType Vert;
   
   int ev = 0; // EXTRA vertices (times 2)
   int ef = 0; // EXTRA faces  
@@ -269,8 +265,8 @@ void MakeBitQuadOnlyByRefine(Mesh &m){
   ev/=2; // I was counting each of them twice
 
   int originalFaceNum = m.fn;
-  FaceIterator nfi = tri::Allocator<Mesh>::AddFaces(m,ef); 
-  VertexIterator nvi = tri::Allocator<Mesh>::AddVertices(m,ev); 
+  FaceIterator nfi = tri::Allocator<MeshType>::AddFaces(m,ef); 
+  VertexIterator nvi = tri::Allocator<MeshType>::AddVertices(m,ev); 
 
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD())  fi->ClearV();
 
@@ -291,12 +287,12 @@ void MakeBitQuadOnlyByRefine(Mesh &m){
       if (k==0) // add a vertex in the center of the face, splitting it in 3
       {
         assert(nvi!=m.vert.end());
-        Vert *nv = &*nvi; nvi++;
+        VertexType *nv = &*nvi; nvi++;
         *nv = *fi->V0( 0 ); // lazy: copy everything from the old vertex
         nv->P() = ( fi->V(0)->P() + fi->V(1)->P() + fi->V(2)->P() )  /3.0;
-        Face *fa = &*fi;
-        Face *fb = &*nfi; nfi++;
-        Face *fc = &*nfi; nfi++;
+        FaceType *fa = &*fi;
+        FaceType *fb = &*nfi; nfi++;
+        FaceType *fc = &*nfi; nfi++;
         *fb = *fc = *fa;  // lazy: copy everything from the old faces
         fa->V(0) = nv;
         fb->V(1) = nv;
@@ -342,9 +338,9 @@ void MakeBitQuadOnlyByRefine(Mesh &m){
     }
     else {      
       // assuming is a part of quad (not a penta, etc), i.e. only one faux
-      Face *fa = &*fi;
-      int ea2 = FauxIndex(fa); // index of the only faux edge
-      Face *fb = fa->FFp(ea2);
+      FaceType *fa = &*fi;
+      int ea2 = BQ::FauxIndex(fa); // index of the only faux edge
+      FaceType *fb = fa->FFp(ea2);
       int eb2 = fa->FFi(ea2);
       assert(fb->FFp(eb2)==fa) ;
       assert(fa->IsF(ea2));
@@ -357,15 +353,15 @@ void MakeBitQuadOnlyByRefine(Mesh &m){
 
       // create new vert in center of faux edge
       assert(nvi!=m.vert.end());
-      Vert *nv = &*nvi; nvi++;
+      VertexType *nv = &*nvi; nvi++;
       *nv = * fa->V0( ea2 );
       nv->P() = ( fa->V(ea2)->P() + fa->V(ea0)->P() ) /2.0;
           
       // split faces: add 2 faces (one per side)
       assert(nfi!=m.face.end());
-      Face *fc = &*nfi; nfi++;
+      FaceType *fc = &*nfi; nfi++;
       assert(nfi!=m.face.end());
-      Face *fd = &*nfi; nfi++;
+      FaceType *fd = &*nfi; nfi++;
       *fc = *fa;
       *fd = *fb;
           
@@ -443,10 +439,10 @@ void MakeBitQuadOnlyByRefine(Mesh &m){
   
   // last pass: add vertex on faux border faces... (if any)
   if (nsplit>0) { 
-    FaceIterator nfi = tri::Allocator<Mesh>::AddFaces(m,nsplit); 
-    VertexIterator nvi = tri::Allocator<Mesh>::AddVertices(m,nsplit); 
+    FaceIterator nfi = tri::Allocator<MeshType>::AddFaces(m,nsplit); 
+    VertexIterator nvi = tri::Allocator<MeshType>::AddVertices(m,nsplit); 
     for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) {
-      Face* fa = &*fi;
+      FaceType* fa = &*fi;
       int ea2 = -1; // border and faux face (if any)
       if (fa->FFp(0)==fa &&  fa->IsF(0) ) ea2=0;
       if (fa->FFp(1)==fa &&  fa->IsF(1) ) ea2=1;
@@ -458,12 +454,12 @@ void MakeBitQuadOnlyByRefine(Mesh &m){
         int ea1 = (ea2+2) %3;        
         
         // create new vert in center of faux edge
-        Vert *nv = &*nvi; nvi++;
+        VertexType *nv = &*nvi; nvi++;
         *nv = * fa->V0( ea2 );
         nv->P() = ( fa->V(ea2)->P() + fa->V(ea0)->P() ) /2.0;
           
         // split face: add 1 face
-        Face *fc = &*nfi; nfi++;
+        FaceType *fc = &*nfi; nfi++;
         *fc = *fa;
           
         fa->V(ea2) = fc->V(ea0) = nv ;
@@ -502,38 +498,33 @@ void MakeBitQuadOnlyByRefine(Mesh &m){
 
 // uses Catmull Clark to enforce quad only meshes
 // each old edge (but not faux) is split in two.
-template <class Mesh> 
-void MakeBitQuadOnlyByCatmullClark(Mesh &m){
-  MakeBitQuadOnlyByRefine(m);
-  MakeBitQuadOnlyByRefine(m);
+static void MakePureByCatmullClark(MeshType &m){
+  MakePureByRefine(m);
+  MakePureByRefine(m);
   // et-voilą!!!
 }
 
 // Helper funcion:
 // marks edge distance froma a given face.
 // Stops at maxDist or at the distance when a triangle is found
-template <class Mesh>
-typename Mesh::FaceType * MarkEdgeDistance(Mesh &m, typename Mesh::FaceType *f, int maxDist){
-  typedef typename Mesh::FaceType Face;
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::VertexIterator VertexIterator;
-  assert(Mesh::HasPerFaceQuality());
+static FaceType * MarkEdgeDistance(MeshType &m, FaceType *f, int maxDist){
+  assert(MeshType::HasPerFaceQuality());
   
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++)  if (!f->IsD()) {
     fi->Q()=maxDist;
   }
 
-  Face * firstTriangleFound = NULL;
+  FaceType * firstTriangleFound = NULL;
   
   f->Q() =  0;
-  std::vector<Face*> stack;
+  std::vector<FaceType*> stack;
   int stackPos=0;
   stack.push_back(f);
 
   while ( stackPos<stack.size() ) {
-    Face *f = stack[stackPos++];
+    FaceType *f = stack[stackPos++];
     for (int k=0; k<3; k++) {
-      Face *fk = f->FFp(k);
+      FaceType *fk = f->FFp(k);
       int fq = int(f->Q()) + ( ! f->IsF(k) );
       if (fk->Q()> fq && fq <= maxDist) {
         if (!fk->IsAnyF()) { firstTriangleFound = fk; maxDist = fq;}
@@ -556,13 +547,9 @@ typename Mesh::FaceType * MarkEdgeDistance(Mesh &m, typename Mesh::FaceType *f,
   
   maxdist is the maximal edge distance where to look for a companion triangle
 */
-template <class Mesh>
-int MakeBitQuadOnlyByFlipStepByStep(Mesh &m, int maxdist=10000, int restart=false){
-  typedef typename Mesh::FaceType Face;
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::VertexIterator VertexIterator;
+static int MakePureByFlipStepByStep(MeshType &m, int maxdist=10000, int restart=false){
   
-  static Face *ta, *tb; // faces to be matched into a quad
+  static FaceType *ta, *tb; // faces to be matched into a quad
   
   static int step = 0; // hack
   
@@ -595,18 +582,18 @@ if (step==0) {
     for (int k=0; k<3; k++) {
       if (tb->FFp(k) == tb) continue; // border
       
-      Face* tbk = tb->FFp(k);
+      FaceType* tbk = tb->FFp(k);
       
       if (!tbk->IsAnyF()) {done=true; marriageEdge=k; break; } // found my match
       
       int back = tb->FFi(k);
-      int faux = FauxIndex(tbk);
+      int faux = BQ::FauxIndex(tbk);
       int other = 3-back-faux;
       
       int scoreA = int(tbk->FFp(other)->Q());
       
-      Face* tbh = tbk->FFp(faux);
-      int fauxh = FauxIndex(tbh);
+      FaceType* tbh = tbk->FFp(faux);
+      int fauxh = BQ::FauxIndex(tbh);
       
       int scoreB = int(tbh->FFp( (fauxh+1)%3 )->Q());
       int scoreC = int(tbh->FFp( (fauxh+2)%3 )->Q());
@@ -623,10 +610,10 @@ if (step==0) {
     
     // use that edge to proceed
     if (mustDoFlip) {
-      FlipBitQuadDiag( *(tb->FFp(edge)) );
+      BQ::FlipDiag( *(tb->FFp(edge)) );
     }
     
-    Face* next = tb->FFp(edge)->FFp( FauxIndex(tb->FFp(edge))  );
+    FaceType* next = tb->FFp(edge)->FFp( BQ::FauxIndex(tb->FFp(edge))  );
     
     // create new edge
     next->ClearAllF();
@@ -660,12 +647,11 @@ break;
   - maxdist is the maximal edge distance where to look for a companion triangle
   - retunrs true if all triangles are merged (always, unless they are odd, or maxdist not enough).
 */
-template <class Mesh>
-bool MakeBitQuadOnlyByFlip(Mesh &m, int maxdist=10000)
+static bool MakePureByFlip(MeshType &m, int maxdist=10000)
 {
-  MakeBitQuadOnlyByFlipStepByStep(m, maxdist, true); // restart
+  MakePureByFlipStepByStep(m, maxdist, true); // restart
   int res=-1;
-  while (res==-1) res = MakeBitQuadOnlyByFlipStepByStep(m, maxdist);
+  while (res==-1) res = MakePureByFlipStepByStep(m, maxdist);
   return res==0;
 }
 
@@ -676,23 +662,22 @@ bool MakeBitQuadOnlyByFlip(Mesh &m, int maxdist=10000)
       level = 1: smarter: leaves more triangles, but makes better quality quads
       level = 2: even more so (marginally)
 */
-template <class Mesh>
-void MakeBitQuadDominant(Mesh &m, int level){
+static void MakeDominant(MeshType &m, int level){
   
-  assert(Mesh::HasPerFaceQuality());
-  assert(Mesh::HasPerFaceFlags());
+  assert(MeshType::HasPerFaceQuality());
+  assert(MeshType::HasPerFaceFlags());
   
-  typedef typename Mesh::FaceIterator FaceIterator;
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) {
     fi->ClearAllF();
     fi->Q() = 0;
   }
   
 
-  MakeQuadDominantPass<Mesh, false> (m);
-  if (level>0)  MakeQuadDominantPass<Mesh, true> (m);
-  if (level>1)  MakeQuadDominantPass<Mesh, true> (m);
-  if (level>0)  MakeQuadDominantPass<Mesh, false> (m);
+  MakeDominantPass<false> (m);
+  if (level>0)  MakeDominantPass<true> (m);
+  if (level>1)  MakeDominantPass<true> (m);
+  if (level>0)  MakeDominantPass<false> (m);
 }
 
+};
 }} // end namespace vcg::tri
diff --git a/vcg/complex/trimesh/bitquad_optimization.h b/vcg/complex/trimesh/bitquad_optimization.h
index 153cb22b..341cec21 100644
--- a/vcg/complex/trimesh/bitquad_optimization.h
+++ b/vcg/complex/trimesh/bitquad_optimization.h
@@ -1,12 +1,26 @@
 
 namespace vcg{namespace tri{
 
-// helper function: mark a quadface, setting Q at 0, and neight at .75, 0.5...
-template <class Mesh>
-void MarkFace(typename Mesh::FaceType* f, Mesh &m){
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
+template <class BQ>
+class BitQuadOptimization{
+  
+typedef typename BQ::MeshType MeshType;
 
+typedef typename MeshType::ScalarType ScalarType;
+typedef typename MeshType::CoordType CoordType;
+typedef typename MeshType::FaceType FaceType;
+typedef typename MeshType::FaceType* FaceTypeP;
+typedef typename MeshType::VertexType VertexType;
+typedef typename MeshType::FaceIterator FaceIterator;
+typedef typename MeshType::VertexIterator VertexIterator;
+
+//typedef BitQuad<MeshType> BQ; // static class to make basic quad operatins
+
+public:
+
+// helper function: mark a quadface, setting Q at 0, and neight at .75, 0.5...
+static void MarkFace(FaceType* f, MeshType &m){
+  
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) {
      fi->Q() = 1; 
   }
@@ -22,11 +36,8 @@ void MarkFace(typename Mesh::FaceType* f, Mesh &m){
 }
 
 // helper function: mark a quadface, setting Q at 0, and neight at .75, 0.5...
-template <class Mesh>
-void MarkVertex(typename Mesh::FaceType* f, int wedge, Mesh &m){
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
-  typedef typename Mesh::VertexType VertexType;
+static void MarkVertex(FaceType* f, int wedge, MeshType &m){
+
   VertexType *v = f->V(wedge);
   
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) {
@@ -36,12 +47,8 @@ void MarkVertex(typename Mesh::FaceType* f, int wedge, Mesh &m){
   
 }
 
-template <class Mesh>
-bool MarkSmallestEdge(Mesh &m, bool perform)
+static bool MarkSmallestEdge(MeshType &m, bool perform)
 {
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
-  typedef typename Mesh::ScalarType ScalarType;
   ScalarType min = std::numeric_limits<ScalarType>::max();
   
   FaceType *fa=NULL; int w=0;
@@ -64,7 +71,7 @@ bool MarkSmallestEdge(Mesh &m, bool perform)
   }
   if (fa) {
     if (perform) {
-      return CollapseQuadEdge(*fa,w,m);
+      return BQ::CollapseEdge(*fa,w,m);
     } else {
       fa->Q()=0.0;
       fa->FFp(w)->Q()=0.0;
@@ -75,12 +82,8 @@ bool MarkSmallestEdge(Mesh &m, bool perform)
 }
 
 // returns: 0 if fail. 1 if edge. 2 if diag.
-template <class Mesh>
-int MarkSmallestEdgeOrDiag(Mesh &m, typename Mesh::ScalarType edgeMult, bool perform)
+static int MarkSmallestEdgeOrDiag(MeshType &m, ScalarType edgeMult, bool perform)
 {
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
-  typedef typename Mesh::ScalarType ScalarType;
   ScalarType min = std::numeric_limits<ScalarType>::max();
   
   FaceType *fa=NULL; int w=0; bool counterDiag = false;
@@ -104,7 +107,7 @@ int MarkSmallestEdgeOrDiag(Mesh &m, typename Mesh::ScalarType edgeMult, bool per
     }
     
     if (f->IsF(k)) { // for diag faces, test counterdiag too
-      score = CounterDiag(f).Norm();
+      score = BQ::CounterDiag(f).Norm();
       if (score<min) {
         min=score; 
         fa = f;
@@ -120,12 +123,12 @@ int MarkSmallestEdgeOrDiag(Mesh &m, typename Mesh::ScalarType edgeMult, bool per
     if (perform) {
       if (fa->IsF(w)) {
         if (counterDiag) {
-          CollapseQuadCounterDiag(*fa, PosOnDiag(*fa,true), m ); return 2;
+          BQ::CollapseCounterDiag(*fa, BQ::PosOnDiag(*fa,true), m ); return 2;
         } else {
-          CollapseQuadDiag(*fa, PosOnDiag(*fa,false), m ); return 2;
+          BQ::CollapseDiag(*fa, BQ::PosOnDiag(*fa,false), m ); return 2;
         }
       } else {
-        if  (CollapseQuadEdge(*fa,w,m)) return 1;
+        if  (BQ::CollapseEdge(*fa,w,m)) return 1;
       }
     } else {
       fa->Q()=0.0;
@@ -137,12 +140,8 @@ int MarkSmallestEdgeOrDiag(Mesh &m, typename Mesh::ScalarType edgeMult, bool per
 }
 
 
-template <class Mesh>
-void MarkSmallestDiag(Mesh &m)
+static void MarkSmallestDiag(MeshType &m)
 {
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
-  typedef typename Mesh::ScalarType ScalarType;
   ScalarType min = std::numeric_limits<ScalarType>::max();
   
   FaceType *fa=NULL; 
@@ -151,13 +150,13 @@ void MarkSmallestDiag(Mesh &m)
     
     ScalarType score;
 
-    score = Diag(f).Norm();
+    score = BQ::Diag(f).Norm();
     if (score<min) {
       min=score; 
       fa = f;
     }
     
-    score = CounterDiag(f).Norm();
+    score = BQ::CounterDiag(f).Norm();
     if (score<min) {
       min=score; 
       fa = f;
@@ -166,16 +165,12 @@ void MarkSmallestDiag(Mesh &m)
   }
   if (fa) {
     fa->Q()=0.0;
-    fa->FFp(FauxIndex(fa))->Q()=0.0;
+    fa->FFp(BQ::FauxIndex(fa))->Q()=0.0;
   }
 
 }
 
-template <class Mesh>
-bool IdentifyAndCollapseSmallestDiag(Mesh &m){
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
-  typedef typename Mesh::ScalarType ScalarType;
+static bool IdentifyAndCollapseSmallestDiag(MeshType &m){
   ScalarType min = std::numeric_limits<ScalarType>::max();
   
   FaceType *fa=NULL; bool flip;
@@ -184,14 +179,14 @@ bool IdentifyAndCollapseSmallestDiag(Mesh &m){
     
     ScalarType score;
     
-    score = Diag(f).Norm();
+    score = BQ::Diag(f).Norm();
     if (score<min) {
       min=score; 
       fa = f;
       flip = false;
     }
     
-    score = CounterDiag(f).Norm();
+    score = BQ::CounterDiag(f).Norm();
     if (score<min) {
       min=score; 
       fa = f;
@@ -201,22 +196,22 @@ bool IdentifyAndCollapseSmallestDiag(Mesh &m){
   }
   if (!fa) return false;
   
-  if (TestAndRemoveDoublet(*fa,0,m)) { return true; }
-  if (TestAndRemoveDoublet(*fa,1,m)) { return true; }
-  if (TestAndRemoveDoublet(*fa,2,m)) { return true; }
-  int k = FauxIndex(fa);
-  if (TestAndRemoveDoublet( *fa->FFp(k),(fa->FFi(k)+2)%3, m )) return true;
+  if (BQ::TestAndRemoveDoublet(*fa,0,m)) { return true; }
+  if (BQ::TestAndRemoveDoublet(*fa,1,m)) { return true; }
+  if (BQ::TestAndRemoveDoublet(*fa,2,m)) { return true; }
+  int k = BQ::FauxIndex(fa);
+  if (BQ::TestAndRemoveDoublet( *fa->FFp(k),(fa->FFi(k)+2)%3, m )) return true;
 
   if (flip) {
-    if (!CheckFlipBitQuadDiag(*fa) ) {
+    if (!BQ::CheckFlipDiag(*fa) ) {
       // I can't collapse (why?)
       MarkFace(fa,m);
       return false;
     } else 
-    CollapseQuadCounterDiag(*fa, PosOnDiag(*fa,true), m );
+    BQ::CollapseCounterDiag(*fa, BQ::PosOnDiag(*fa,true), m );
   }
   else  {
-    CollapseQuadDiag(*fa, PosOnDiag(*fa,false), m );
+    BQ::CollapseDiag(*fa, BQ::PosOnDiag(*fa,false), m );
   }
   return true;
 }
@@ -228,18 +223,15 @@ seeks and removes all doublets (a pair of quads sharing two consecutive edges)
 by merging them into a single quad (thus removing one vertex and two tri faces)-
 Returns number of removed Doublets
 */
-template <class Mesh>
-int BitQuadRemoveDoublets(Mesh &m)
+static int RemoveDoublets(MeshType &m)
 {
   int res=0;
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) {
     fi->Q()=1;
     for (int k=0; k<3; k++) {
-      if ( IsDoublet(*fi,k) ){
+      if ( BQ::IsDoublet(*fi,k) ){
         res++;
-        RemoveDoublet(*fi,k,m);
+        BQ::RemoveDoublet(*fi,k,m);
         if (fi->IsD()) break; // break wedge circle, if face disappeard
       }
     }
@@ -251,14 +243,10 @@ int BitQuadRemoveDoublets(Mesh &m)
 marks (Quality=0) and approx. counts profitable vertex rotations
 (vertex rotations which make edge shorter
 */
-template <class Mesh, bool perform>
-int BitQuadMarkVertexRotations(Mesh &m)
+template <bool perform>
+static int MarkVertexRotations(MeshType &m)
 {
   int res=0;
-  typedef typename Mesh::VertexIterator VertexIterator;
-  typedef typename Mesh::VertexType VertexType;
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
   for (VertexIterator vi = m.vert.begin();  vi!=m.vert.end(); vi++) if (!vi->IsD()) vi->ClearV();
   if (!perform)
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) fi->Q()=1.0;
@@ -267,14 +255,14 @@ int BitQuadMarkVertexRotations(Mesh &m)
     
     for (int k=0; k<3; k++) {
       if (fi->V(k)->IsV()) continue;
-      if (TestBitQuadVertexRotation(*fi,k)) {
+      if (BQ::TestVertexRotation(*fi,k)) {
         res++;
         fi->V(k)->SetV();
         if (!perform) {
           res++; MarkVertex(&*fi, k, m); //fi->Q()=0;
         }
         else {
-          if (RotateBitQuadVertex(*fi, k)) res++; //fi->Q()=0;
+          if (BQ::RotateVertex(*fi, k)) res++; //fi->Q()=0;
           //if (res>1) return res; // uncomment for only one rotation
         }
       }
@@ -285,12 +273,10 @@ int BitQuadMarkVertexRotations(Mesh &m)
 
 // mark (and count) all edges that are worth rotating
 // if perform == true, actually rotate them
-template <class Mesh, bool perform>
-int BitQuadMarkEdgeRotations(Mesh &m)
+template <bool perform>
+static int MarkEdgeRotations(MeshType &m)
 {
   int count = 0;
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
   
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) fi->Q()=1;
   
@@ -299,10 +285,10 @@ int BitQuadMarkEdgeRotations(Mesh &m)
     for (int k=0; k<3; k++) {
       if (fi->IsF(k)) continue;
       if (fi->FFp(k)<= &*fi) continue; // only once per real (non faux) edge, and only for non border ones
-      int best = TestBitQuadEdgeRotation(*fi, k);
+      int best = BQ::TestEdgeRotation(*fi, k);
       if (perform) {
-        if (best==+1) if (RotateBitQuadEdge<FaceType, true>(*fi, k)) count++;
-        if (best==-1) if (RotateBitQuadEdge<FaceType,false>(*fi, k)) count++;
+        if (best==+1) if (BQ::template RotateEdge< true>(*fi, k)) count++;
+        if (best==-1) if (BQ::template RotateEdge<false>(*fi, k)) count++;
       }
       else {
         if (best!=0) { fi->Q()=0; fi->FFp(k)->Q()=0; count++; }
@@ -316,16 +302,13 @@ int BitQuadMarkEdgeRotations(Mesh &m)
 /*
 marks (Quality=0) and approx. counts doublets (a pair of quads sharing two consecutive edges)
 */
-template <class Mesh>
-int BitQuadMarkDoublets(Mesh &m)
+static int MarkDoublets(MeshType &m)
 {
   int res=0;
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) {
     fi->Q()=1;
     for (int k=0; k<3; k++) {
-      if ( IsDoublet(*fi,k) ){
+      if ( BQ::IsDoublet(*fi,k) ){
         res++;
         if (fi->IsF((k+1)%3)) res++; // counts for a quad
         fi->Q()=0;
@@ -339,16 +322,13 @@ int BitQuadMarkDoublets(Mesh &m)
 /*
 marks (Quality=0) and counts singlets (vertex B in an A-B-A-C quad)
 */
-template <class Mesh>
-int BitQuadMarkSinglets(Mesh &m)
+static int MarkSinglets(MeshType &m)
 {
   int res=0;
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) {
     fi->Q()=1;
     for (int k=0; k<3; k++) {
-      if ( IsSinglet(*fi,k) ){
+      if ( BQ::IsSinglet(*fi,k) ){
         res++;
         fi->Q()=0;
       }
@@ -361,17 +341,14 @@ int BitQuadMarkSinglets(Mesh &m)
 /*
 deletes singlets, reutrns number of
 */
-template <class Mesh>
-int BitQuadRemoveSinglets(Mesh &m)
+static int RemoveSinglets(MeshType &m)
 {
   int res=0;
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::FaceType FaceType;
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) {
     for (int k=0; k<3; k++) {
-      if ( IsSinglet(*fi,k) ){
+      if ( BQ::IsSinglet(*fi,k) ){
         res++;
-        RemoveSinglet(*fi,k,m);
+        BQ::RemoveSinglet(*fi,k,m);
         return res;
         break;
       }
@@ -383,19 +360,17 @@ int BitQuadRemoveSinglets(Mesh &m)
 
 /* returns average quad quality, and assigns it to triangle quality
 */
-template <class Mesh>
-typename Mesh::ScalarType MeasureBitQuadQuality(Mesh &m)
+static ScalarType MeasureQuality(MeshType &m)
 {
-  assert(Mesh::HasPerFaceFlags());
-  typename Mesh::ScalarType res = 0;
+  assert(MeshType::HasPerFaceFlags());
+  ScalarType res = 0;
   int div = 0;
-  typedef typename Mesh::FaceIterator FaceIterator;
   for (FaceIterator fi = m.face.begin();  fi!=m.face.end(); fi++) if (!fi->IsD()) {
     if (fi->IsAnyF()) {
       
-      typename Mesh::ScalarType q = quadQuality( &*fi, FauxIndex(&*fi) );
+      ScalarType q = BQ::quadQuality( &*fi, BQ::FauxIndex(&*fi) );
       
-      if (Mesh::HasPerFaceQuality()) fi->Q() = q;
+      if (MeshType::HasPerFaceQuality()) fi->Q() = q;
       res += q;
       div++;
     }
@@ -403,4 +378,5 @@ typename Mesh::ScalarType MeasureBitQuadQuality(Mesh &m)
   if (!div) return 0; else return res / div;
 }
 
+};
 }} // end namespace vcg::tri
diff --git a/vcg/complex/trimesh/bitquad_support.h b/vcg/complex/trimesh/bitquad_support.h
index 2f112df4..00d1af07 100644
--- a/vcg/complex/trimesh/bitquad_support.h
+++ b/vcg/complex/trimesh/bitquad_support.h
@@ -10,23 +10,23 @@
 
    [ basic operations: ]
 
-   bool IsDoublet(const Face& f, int wedge)
-   void RemoveDoublet(Face &f, int wedge, Mesh& m)
+   bool IsDoublet(const FaceType& f, int wedge)
+   void RemoveDoublet(FaceType &f, int wedge, MeshType& m)
     - identifies and removed "Doublets" (pair of quads sharing two consecutive edges)
 
-   bool IsSinglet(const Face& f, int wedge)
-   void RemoveSinglet(Face &f, int wedge, Mesh& m)
+   bool IsSinglet(const FaceType& f, int wedge)
+   void RemoveSinglet(FaceType &f, int wedge, MeshType& m)
 
-   void FlipBitQuadDiag(Face &f)
+   void FlipDiag(FaceType &f)
     - rotates the faux edge of a quad (quad only change internally)
 
-   bool RotateBitQuadEdge(Face& f, int w0a);
+   bool RotateEdge(FaceType& f, int w0a);
     - rotate a quad edge (clockwise or counterclockwise, specified via template)
    
-   bool RotateBitQuadVertex(FaceType &f, int w0)
+   bool RotateVertex(FaceType &f, int w0)
     - rotate around a quad vertex ("wind-mill" operation)
 
-   void CollapseQuadDiag(Face &f, ... p , Mesh& m)
+   void CollapseDiag(FaceType &f, ... p , MeshType& m)
     - collapses a quad on its diagonal. 
     - p identifies the pos of collapsed point 
       (as either the parametric pos on the diagonal, or a fresh coordtype)
@@ -39,10 +39,10 @@
     - (should be made into a template parameter for methods using it)
     - currently measures how squared each angle is
    
-   int FauxIndex(const Face* f);
+   int FauxIndex(const FaceType* f);
     - returns index of the only faux edge of a quad (otherwise, assert)
 
-   int CountBitPolygonInternalValency(const Face& f, int wedge)
+   int CountBitPolygonInternalValency(const FaceType& f, int wedge)
     - returns valency of vertex in terms  of polygons (quads, tris...)
     
   
@@ -55,134 +55,59 @@
 
 namespace vcg{namespace tri{
 
-// helper function:
-// cos of angle abc. This should probably go elsewhere
-template<class CoordType>
-static typename CoordType::ScalarType Cos(const CoordType &a, const CoordType &b, const CoordType &c )
-{
-  CoordType 
-    e0 = b - a,
-    e1 = b - c;
-  typename CoordType::ScalarType d =  (e0.Norm()*e1.Norm());
-  if (d==0) return 0.0;
-  return (e0*e1)/d;
-}
-
-// helper function:
-// returns quality of a quad formed by points a,b,c,d
-// quality is computed as "how squared angles are"
-template <class Coord>
-inline static typename Coord::ScalarType quadQuality(const Coord &a, const Coord &b, const Coord &c, const Coord &d){
-  typename Coord::ScalarType score = 0;
-  score += 1 - math::Abs( Cos( a,b,c) );
-  score += 1 - math::Abs( Cos( b,c,d) );
-  score += 1 - math::Abs( Cos( c,d,a) );
-  score += 1 - math::Abs( Cos( d,a,b) );
-  return score / 4;
-}
-
-// helper function:
-// returns quality of a given (potential) quad
-template <class Face>
-static typename Face::ScalarType quadQuality(Face *f, int edge){
-  
-  typedef typename Face::CoordType CoordType;
-  
-  CoordType 
-    a = f->V0(edge)->P(),
-    b = f->FFp(edge)->V2( f->FFi(edge) )->P(),
-    c = f->V1(edge)->P(),
-    d = f->V2(edge)->P();
-  
-  return quadQuality(a,b,c,d);
-
-}
-
-/**
-helper function:
-given a quad edge, retruns:
-   0 if that edge should not be rotated
-  +1 if it should be rotated clockwise (+1)
-  -1 if it should be rotated counterclockwise (-1)
-Currently an edge is rotated iff it is shortened by that rotations
-(shortcut criterion)
-*/
-template <class Face>
-int TestBitQuadEdgeRotation(const Face &f, int w0)
-{
-  const Face *fa = &f;
-  assert(! fa->IsF(w0) );
-  typename Face::ScalarType q0,q1,q2;
-  typename Face::CoordType v0,v1,v2,v3,v4,v5;
-  int w1 = (w0+1)%3;
-  int w2 = (w0+2)%3;
-  
-  v0 = fa->P(w0);
-  v3 = fa->P(w1);
-  
-  if (fa->IsF(w2) ) {
-    v1 = fa->cFFp(w2)->V2( fa->cFFi(w2) )->P();
-    v2 = fa->P(w2);
-  } else {
-    v1 = fa->P(w2);
-    v2 = fa->cFFp(w1)->V2( fa->cFFi(w1) )->P();
+/* simple geometric-interpolation mono-function class used 
+as a default template parameter to BitQuad class */
+template <class VertexType>
+class GeometricInterpolator{
+public:
+  typedef typename VertexType::ScalarType ScalarType;
+  static void Apply( const VertexType &a,  const VertexType &b, ScalarType t, VertexType &res){
+    assert (&a != &b);
+    res.P() = a.P()*(1-t) + b.P()*(t);
   }
-  
-  const Face *fb = fa->cFFp(w0);
-  w0 = fa->cFFi(w0);
-  
-  w1 = (w0+1)%3;
-  w2 = (w0+2)%3;
-  if (fb->IsF(w2) ) {
-    v4 = fb->cFFp(w2)->V2( fb->cFFi(w2) )->P();
-    v5 = fb->P(w2);
-  } else {
-    v4 = fb->P(w2);
-    v5 = fb->cFFp(w1)->V2( fb->cFFi(w1) )->P();
-  }
-  
-  /*
-  //  max overall quality criterion:
-  q0 = quadQuality(v0,v1,v2,v3) +  quadQuality(v3,v4,v5,v0); // keep as is?
-  q1 = quadQuality(v1,v2,v3,v4) +  quadQuality(v4,v5,v0,v1); // rotate CW?
-  q2 = quadQuality(v5,v0,v1,v2) +  quadQuality(v2,v3,v4,v5); // rotate CCW?
-  
-  if (q0>=q1 && q0>=q2) return 0;
-  if (q1>=q2) return 1;*/
-  
-  // min distance (shortcut criterion)
-  q0 = (v0 - v3).SquaredNorm();
-  q1 = (v1 - v4).SquaredNorm();
-  q2 = (v5 - v2).SquaredNorm();
-  if (q0<=q1 && q0<=q2) return 0;
-  if (q1<=q2) return 1;
-  return -1;
-}
+};
+
+template <
+  // first template parameter: the tri mesh (with face-edges flagged)
+  class _MeshType, 
+  // second template parameter: used to define interpolations between points
+  class Interpolator = GeometricInterpolator<typename _MeshType::VertexType> 
+>
+class BitQuad{
+public:
+
+typedef _MeshType MeshType;
+typedef typename MeshType::ScalarType ScalarType;
+typedef typename MeshType::CoordType CoordType;
+typedef typename MeshType::FaceType FaceType;
+typedef typename MeshType::FaceType* FaceTypeP;
+typedef typename MeshType::VertexType VertexType;
+typedef typename MeshType::FaceIterator FaceIterator;
+typedef typename MeshType::VertexIterator VertexIterator;
 
 
-
-template <class Face, bool verse>
-bool RotateBitQuadEdge(Face& f, int w0a){
-  Face *fa = &f;
+template <bool verse>
+static bool RotateEdge(FaceType& f, int w0a){
+  FaceType *fa = &f;
   assert(! fa->IsF(w0a) );
 
-  typename Face::VertexType *v0, *v1;
+  VertexType *v0, *v1;
   v0= fa->V0(w0a);
   v1= fa->V1(w0a);
   
   int w1a = (w0a+1)%3;
   int w2a = (w0a+2)%3;
 
-  Face *fb = fa->FFp(w0a);
+  FaceType *fb = fa->FFp(w0a);
   int w0b = fa->FFi(w0a);
   int w1b = (w0b+1)%3;
   int w2b = (w0b+2)%3;
     
   if (fa->IsF(w2a) == verse) {
-    if (!CheckFlipBitQuadDiag(*fa)) return false;
-    FlipBitQuadDiag(*fa);
+    if (!CheckFlipDiag(*fa)) return false;
+    FlipDiag(*fa);
     // recover edge index, so that (f, w0a) identifies the same edge as before
-    Face *fc = fa->FFp(FauxIndex(fa));
+    FaceType *fc = fa->FFp(FauxIndex(fa));
     for (int i=0; i<3; i++){
       if ( v0==fa->V0(i) && v1==fa->V1(i) ) w0a = i;
       if ( v0==fc->V0(i) && v1==fc->V1(i) ) { fa = fc; w0a = i; }
@@ -190,20 +115,19 @@ bool RotateBitQuadEdge(Face& f, int w0a){
   }
   
   if (fb->IsF(w2b) == verse) {
-    if (!CheckFlipBitQuadDiag(*fb)) return false;
-    FlipBitQuadDiag(*fb);
+    if (!CheckFlipDiag(*fb)) return false;
+    FlipDiag(*fb);
   }
   
   if (!CheckFlipEdge(*fa,w0a)) return false;
-  FlipBitQuadEdge(*fa,w0a);
+  FlipEdge(*fa,w0a);
   return true;
 }
 
 /* small helper function which returns the index of the only
    faux index, assuming there is exactly one (asserts out otherwise)
 */
-template <class Face>
-int FauxIndex(const Face* f){
+static int FauxIndex(const FaceType* f){
   if (f->IsF(0)) return 0;
   if (f->IsF(1)) return 1;
   assert(f->IsF(2));
@@ -211,11 +135,10 @@ int FauxIndex(const Face* f){
 }
 
 // rotates the diagonal of a quad
-template <class Face>
-void FlipBitQuadDiag(Face &f){
+static void FlipDiag(FaceType &f){
   int faux = FauxIndex(&f);
-  Face* fa = &f;
-  Face* fb = f.FFp(faux);
+  FaceType* fa = &f;
+  FaceType* fb = f.FFp(faux);
   vcg::face::FlipEdge(f, faux);
   // ripristinate faux flags
   fb->ClearAllF();
@@ -230,11 +153,9 @@ void FlipBitQuadDiag(Face &f){
 // given a vertex (i.e. a face and a wedge), 
 // this function tells us how the average edge lenght around a vertex would change
 // if that vertex is rotated
-template <class Face>
-typename Face::ScalarType AvgBitQuadEdgeLenghtVariationIfVertexRotated(const Face &f, int w0)
+static ScalarType AvgEdgeLenghtVariationIfVertexRotated(const FaceType &f, int w0)
 {
   assert(!f.IsD());
-  typedef typename Face::ScalarType ScalarType;
   
   ScalarType 
     before=0, // sum of quad edges (originating from v)
@@ -242,7 +163,7 @@ typename Face::ScalarType AvgBitQuadEdgeLenghtVariationIfVertexRotated(const Fac
   int guard = 0;
 
   // rotate arond vertex
-  const Face* pf = &f;
+  const FaceType* pf = &f;
   int pi = w0;
   int n = 0; // vertex valency
   int na = 0; 
@@ -252,11 +173,11 @@ typename Face::ScalarType AvgBitQuadEdgeLenghtVariationIfVertexRotated(const Fac
     else { before+= triEdge; n++; }
     if ( pf->IsF((pi+1)%3)) { after += CounterDiag( pf ).Norm(); na++; }
     
-    const Face *t = pf;
+    const FaceType *t = pf;
     t = pf->FFp( pi );
     if (pf == t ) return std::numeric_limits<ScalarType>::max(); // it's a mesh border! flee!
     pi = pf->cFFi( pi );
-    pi = (pi+1)%3; // Face::Next( pf->FFi( pi ) );
+    pi = (pi+1)%3; // FaceType::Next( pf->FFi( pi ) );
     pf = t;
     assert(guard++<100);
   } while (pf != &f);
@@ -265,17 +186,17 @@ typename Face::ScalarType AvgBitQuadEdgeLenghtVariationIfVertexRotated(const Fac
 }
 
 /*
-  const Face* pf = &f;
+  const FaceType* pf = &f;
   int pi = wedge;
   int res = 0, guard=0;
   do {
     if (!pf->IsAnyF()) return false; // there's a triangle!
     if (!pf->IsF(pi)) res++;
-    const Face *t = pf;
+    const FaceType *t = pf;
     t = pf->FFp( pi );
     if (pf == t ) return false;
     pi = pf->cFFi( pi );
-    pi = (pi+1)%3; // Face::Next( pf->FFi( pi ) );
+    pi = (pi+1)%3; // FaceType::Next( pf->FFi( pi ) );
     pf = t;
     assert(guard++<100);
   } while (pf != &f);
@@ -284,19 +205,16 @@ typename Face::ScalarType AvgBitQuadEdgeLenghtVariationIfVertexRotated(const Fac
 // given a vertex (i.e. a face and a wedge), 
 // this function tells us if it should be rotated or not
 // (currently, we should iff it is shortened)
-template <class Face>
-bool TestBitQuadVertexRotation(const Face &f, int w0)
+static bool TestVertexRotation(const FaceType &f, int w0)
 {
   assert(!f.IsD());
   // rotate quad IFF this way edges become shorter:
-  return AvgBitQuadEdgeLenghtVariationIfVertexRotated(f,w0)<0;
+  return AvgEdgeLenghtVariationIfVertexRotated(f,w0)<0;
 }
 
 
-template <class FaceType>
-bool RotateBitQuadVertex(FaceType &f, int w0)
+static bool RotateVertex(FaceType &f, int w0)
 {
-  typedef typename FaceType::ScalarType ScalarType;
   
   int guard = 0;
 
@@ -330,8 +248,8 @@ bool RotateBitQuadVertex(FaceType &f, int w0)
     int tmp = (pf->FFi(pi)+1)%3; pf = pf->FFp(pi); pi = tmp;  // flipF
         
     if (mustFlip) {
-      if (!CheckFlipBitQuadDiag(*lastF)) return false; // cannot flip??
-      FlipBitQuadDiag(*lastF);
+      if (!CheckFlipDiag(*lastF)) return false; // cannot flip??
+      FlipDiag(*lastF);
     }
     
   } while (pf != stopA && pf!= stopB);
@@ -354,16 +272,15 @@ bool RotateBitQuadVertex(FaceType &f, int w0)
 
 
 // flips the faux edge of a quad
-template <class Face>
-void FlipBitQuadEdge(Face &f, int k){
+static void FlipEdge(FaceType &f, int k){
   assert(!f.IsF(k));
-  Face* fa = &f;
-  Face* fb = f.FFp(k);
+  FaceType* fa = &f;
+  FaceType* fb = f.FFp(k);
   assert(fa!=fb); // else, rotating a border edge
 
   // backup prev other-quads-halves
-  Face* fa2 = fa->FFp( FauxIndex(fa) );
-  Face* fb2 = fb->FFp( FauxIndex(fb) );
+  FaceType* fa2 = fa->FFp( FauxIndex(fa) );
+  FaceType* fb2 = fb->FFp( FauxIndex(fb) );
 
   vcg::face::FlipEdge(*fa, k);
   
@@ -379,22 +296,19 @@ void FlipBitQuadEdge(Face &f, int k){
 }
 
 // check if a quad diagonal can be topologically flipped
-template <class Face>
-bool CheckFlipBitQuadDiag(Face &f){
+static bool CheckFlipDiag(FaceType &f){
   return (vcg::face::CheckFlipEdge(f, FauxIndex(&f) ) );
 }
 
 // given a face (part of a quad), returns its diagonal
-template <class Face>
-typename Face::CoordType Diag(const Face* f){
+static CoordType Diag(const FaceType* f){
   int i = FauxIndex(f);
   return f->P1( i ) - f->P0( i );
 }
 
 
 // given a face (part of a quad), returns other diagonal
-template <class Face>
-typename Face::CoordType CounterDiag(const Face* f){
+static CoordType CounterDiag(const FaceType* f){
   int i = FauxIndex(f);
   return f->cP2( i ) - f->cFFp( i )->cP2(f->cFFi(i) ) ;
 }
@@ -402,22 +316,20 @@ typename Face::CoordType CounterDiag(const Face* f){
 /* helper function:
    collapses a single face along its faux edge. 
    Updates FF adj of other edges. */
-template <class Mesh>
-void _CollapseQuadDiagHalf(typename Mesh::FaceType &f, int faux, Mesh& m)
+static void _CollapseDiagHalf(FaceType &f, int faux, MeshType& m)
 {
-  typedef typename Mesh::FaceType Face;
   int faux1 = (faux+1)%3;
   int faux2 = (faux+2)%3;
   
-  Face* fA = f.FFp( faux1 );
-  Face* fB = f.FFp( faux2 );
+  FaceType* fA = f.FFp( faux1 );
+  FaceType* fB = f.FFp( faux2 );
   int iA = f.FFi( faux1 );
   int iB = f.FFi( faux2 );
   
   if (fA==&f && fB==&f) { 
     // both non-faux edges are borders: tri-face disappears, just remove the vertex
     if (DELETE_VERTICES)
-      Allocator<Mesh>::DeleteVertex(m,*(f.V(faux2))); 
+      Allocator<MeshType>::DeleteVertex(m,*(f.V(faux2))); 
   } else {
     if (fA==&f) {
       fB->FFp(iB) = fB;  fB->FFi(iB) = iB;
@@ -432,15 +344,14 @@ void _CollapseQuadDiagHalf(typename Mesh::FaceType &f, int faux, Mesh& m)
     }
   }
   
-  Allocator<Mesh>::DeleteFace(m,f);
+  Allocator<MeshType>::DeleteFace(m,f);
 
 }
 
-template <class Mesh>
-void RemoveDoublet(typename Mesh::FaceType &f, int wedge, Mesh& m){
+static void RemoveDoublet(FaceType &f, int wedge, MeshType& m){
   if (f.IsF((wedge+1)%3) ) {
-    typename Mesh::VertexType *v = f.V(wedge);
-    FlipBitQuadDiag(f);
+    VertexType *v = f.V(wedge);
+    FlipDiag(f);
     // quick hack: recover wedge index after flip
     if (f.V(0)==v) wedge = 0;
     else if (f.V(1)==v) wedge = 1;
@@ -449,15 +360,14 @@ void RemoveDoublet(typename Mesh::FaceType &f, int wedge, Mesh& m){
       wedge = 2;
     }
   }
-  typename Mesh::ScalarType k=(f.IsF(wedge))?1:0;
-  CollapseQuadDiag(f, k, m);
-  typename Mesh::VertexType *v = f.V(wedge);
+  ScalarType k=(f.IsF(wedge))?1:0;
+  CollapseDiag(f, k, m);
+  VertexType *v = f.V(wedge);
 }
 
-template <class Mesh>
-void RemoveSinglet(typename Mesh::FaceType &f, int wedge, Mesh& m){
-  typename Mesh::FaceType *fa, *fb; // these will die
-  typename Mesh::FaceType *fc, *fd; // their former neight
+static void RemoveSinglet(FaceType &f, int wedge, MeshType& m){
+  FaceType *fa, *fb; // these will die
+  FaceType *fc, *fd; // their former neight
   fa = & f;
   fb = fa->FFp(wedge);
   int wa0 = wedge;
@@ -478,15 +388,14 @@ void RemoveSinglet(typename Mesh::FaceType &f, int wedge, Mesh& m){
   // faux status of survivors: unchanged
   assert( ! ( fc->IsF( wc) ) );
   assert( ! ( fd->IsF( wd) ) );
-  Allocator<Mesh>::DeleteFace( m,*fa );
-  Allocator<Mesh>::DeleteFace( m,*fb );
+  Allocator<MeshType>::DeleteFace( m,*fa );
+  Allocator<MeshType>::DeleteFace( m,*fb );
   if (DELETE_VERTICES)
-  Allocator<Mesh>::DeleteVertex( m,*fa->V(wedge) );
+  Allocator<MeshType>::DeleteVertex( m,*fa->V(wedge) );
 }
 
 
-template <class Mesh>
-bool TestAndRemoveDoublet(typename Mesh::FaceType &f, int wedge, Mesh& m){
+static bool TestAndRemoveDoublet(FaceType &f, int wedge, MeshType& m){
   if (IsDoublet(f,wedge)) {
      RemoveDoublet(f,wedge,m);
      return true;
@@ -494,32 +403,31 @@ bool TestAndRemoveDoublet(typename Mesh::FaceType &f, int wedge, Mesh& m){
   return false;
 }
 
-template <class Mesh>
-bool TestAndRemoveSinglet(typename Mesh::FaceType &f, int wedge, Mesh& m){
+static bool TestAndRemoveSinglet(FaceType &f, int wedge, MeshType& m){
   if (IsSinglet(f,wedge)) {
      RemoveSinglet(f,wedge,m);
      return true;
   }
   return false;
 }
-template <class Face, int verse>
-void RotateBitQuadEdge(const Face& f, int wedge){
+
+template <int verse>
+static void RotateEdge(const FaceType& f, int wedge){
 }
 
 // given a face and a wedge, counts its valency in terms of quads (and triangles)
 // uses only FF, assumes twomanyfold
 // returns -1 if border
-template <class Face>
-int CountBitPolygonInternalValency(const Face& f, int wedge){
-  const Face* pf = &f;
+static int CountBitPolygonInternalValency(const FaceType& f, int wedge){
+  const FaceType* pf = &f;
   int pi = wedge;
   int res = 0;
   do {
     if (!pf->IsF(pi)) res++;
-    const Face *t = pf;
+    const FaceType *t = pf;
     t = pf->FFp( pi );
     if (pf == t ) return -1;
-    pi = (pi+1)%3; // Face::Next( pf->FFi( pi ) );
+    pi = (pi+1)%3; // FaceType::Next( pf->FFi( pi ) );
     pf = t;
   } while (pf != &f);
   return res;
@@ -527,38 +435,36 @@ int CountBitPolygonInternalValency(const Face& f, int wedge){
 
 // given a face and a wedge, returns if it host a doubet
 // assumes tri and quad only. uses FF topology only.
-template <class Face>
-bool IsDoublet(const Face& f, int wedge){
-  const Face* pf = &f;
+static bool IsDoublet(const FaceType& f, int wedge){
+  const FaceType* pf = &f;
   int pi = wedge;
   int res = 0, guard=0;
   do {
     if (!pf->IsAnyF()) return false; // there's a triangle!
     if (!pf->IsF(pi)) res++;
-    const Face *t = pf;
+    const FaceType *t = pf;
     t = pf->FFp( pi );
     if (pf == t ) return false;
     pi = pf->cFFi( pi );
-    pi = (pi+1)%3; // Face::Next( pf->FFi( pi ) );
+    pi = (pi+1)%3; // FaceType::Next( pf->FFi( pi ) );
     pf = t;
     assert(guard++<100);
   } while (pf != &f);
   return (res == 2);
 }
 
-template <class Face>
-bool IsSinglet(const Face& f, int wedge){
-  const Face* pf = &f;
+static bool IsSinglet(const FaceType& f, int wedge){
+  const FaceType* pf = &f;
   int pi = wedge;
   int res = 0, guard=0;
   do {
     if (!pf->IsAnyF()) return false; // there's a triangle!
     if (!pf->IsF(pi)) res++;
-    const Face *t = pf;
+    const FaceType *t = pf;
     t = pf->FFp( pi );
     if (pf == t ) return false;
     pi = pf->cFFi( pi );
-    pi = (pi+1)%3; // Face::Next( pf->FFi( pi ) );
+    pi = (pi+1)%3; // FaceType::Next( pf->FFi( pi ) );
     pf = t;
     assert(guard++<100);
   } while (pf != &f);
@@ -566,32 +472,16 @@ bool IsSinglet(const Face& f, int wedge){
 }
 
 
-
-/** collapses a quad diagonal a-b 
-  forming the new vertex in between the two old vertices.
-   if k == 0, new vertex is in a
-   if k == 1, new vertex is in b
-   if k == 0.5, new vertex in the middle, etc
-*/
-template <class Mesh>
-void CollapseQuadDiag(typename Mesh::FaceType &f, typename Mesh::ScalarType k, Mesh& m){
-  typename Mesh::CoordType p;
-  int fauxa = FauxIndex(&f);
-  p = f.V(fauxa)->P()*(1-k) + f.V( (fauxa+1)%3 )->P()*(k);
-  CollapseQuadDiag(f,p,m);
-}
-
-template <class Mesh>
-bool CollapseQuadEdgeDirect(typename Mesh::FaceType &f, int w0, Mesh& m){
-  typename Mesh::FaceType * f0 = &f;
+static bool CollapseEdgeDirect(FaceType &f, int w0, MeshType& m){
+  FaceType * f0 = &f;
 
   assert( !f0->IsF(w0) );
   
-  typename Mesh::VertexType *v0, *v1;
+  VertexType *v0, *v1;
   v0 = f0->V0(w0);
   v1 = f0->V1(w0);
   
-  if (!RotateBitQuadVertex(*f0,w0)) return false;
+  if (!RotateVertex(*f0,w0)) return false;
   
   
   // quick hack: recover original wedge
@@ -603,13 +493,11 @@ bool CollapseQuadEdgeDirect(typename Mesh::FaceType &f, int w0, Mesh& m){
   assert( f0->V1(w0) == v1 );
   assert( f0->IsF(w0) );
   
-  CollapseQuadDiag(*f0,PosOnDiag(*f0,false), m);
+  CollapseDiag(*f0,PosOnDiag(*f0,false), m);
   return true;
 }
 
-template <class Mesh>
-bool CollapseQuadEdge(typename Mesh::FaceType &f, int w0, Mesh& m){
-  typedef typename Mesh::FaceType * FaceTypeP;
+static bool CollapseEdge(FaceType &f, int w0, MeshType& m){
   FaceTypeP f0 = &f;
   
   assert(!f0->IsF(w0)); // don't use this to collapse diag.
@@ -621,39 +509,69 @@ bool CollapseQuadEdge(typename Mesh::FaceType &f, int w0, Mesh& m){
   
   // choose: rotate around V0 or around V1?
   if (
-    AvgBitQuadEdgeLenghtVariationIfVertexRotated(*f0,w0)
+    AvgEdgeLenghtVariationIfVertexRotated(*f0,w0)
     <
-    AvgBitQuadEdgeLenghtVariationIfVertexRotated(*f1,w1)
-  )    return CollapseQuadEdgeDirect(*f0,w0,m); 
-  else return CollapseQuadEdgeDirect(*f1,w1,m);
+    AvgEdgeLenghtVariationIfVertexRotated(*f1,w1)
+  )    return CollapseEdgeDirect(*f0,w0,m); 
+  else return CollapseEdgeDirect(*f1,w1,m);
 }
 
 
-template <class Mesh>
-void CollapseQuadDiag(typename Mesh::FaceType &f, const typename Mesh::CoordType &p, Mesh& m){
+
+/** collapses a quad diagonal a-b 
+  forming the new vertex in between the two old vertices.
+   if k == 0, new vertex is in a
+   if k == 1, new vertex is in b
+   if k == 0.5, new vertex in the middle, etc
+*/
+
+static void CollapseCounterDiag(FaceType &f, ScalarType interpol, MeshType& m){
+  //CoordType p;
+  //int fauxa = FauxIndex(&f);
+  //p = f.V(fauxa)->P()*(1-k) + f.V( (fauxa+1)%3 )->P()*(k);
   
-  typedef typename Mesh::FaceType Face;
-  typedef typename Mesh::VertexType Vert;
-  
-  Face* fa = &f;
+  FlipDiag(f);
+  CollapseDiag(f,interpol,m);
+}
+
+/*
+static void CollapseCounterDiag(FaceType &f, ScalarType k, MeshType& m){
+  CoordType p;
+  int fauxa = FauxIndex(&f);
+  p = f.P2(fauxa)*(1-k) + f.FFp( fauxa )->P2( f.FFi( fauxa ) )*(k);
+  CollapseCounterDiag(f,p,m);
+}
+*/
+
+//static void CollapseCounterDiag(FaceType &f, const CoordType &p, MeshType& m){
+//  FlipDiag(f);
+//  CollapseDiag(f,p,m);
+//}
+
+
+static void CollapseDiag(FaceType &f, ScalarType interpol, MeshType& m){
+    
+  FaceType* fa = &f;
   int fauxa = FauxIndex(fa);
-  Face* fb = fa->FFp(fauxa);
+  FaceType* fb = fa->FFp(fauxa);
   assert (fb!=fa);
   int fauxb = FauxIndex(fb);
   
-  Vert* va = fa->V(fauxa); // va lives
-  Vert* vb = fb->V(fauxb); // vb dies
+  VertexType* va = fa->V(fauxa); // va lives
+  VertexType* vb = fb->V(fauxb); // vb dies
   
+  Interpolator::Apply( *(f.V0(fauxa)), *(f.V1(fauxa)), interpol, *va);
+
   // update FV...
   bool border = false;
   int pi = fauxb;
-  Face* pf = fb; /* pf, pi could be a Pos<Face> p(pf, pi) */
+  FaceType* pf = fb; /* pf, pi could be a Pos<FaceType> p(pf, pi) */
   // rotate around vb, (same-sense-as-face)-wise
   do {
     pf->V(pi) = va;
     
     pi=(pi+2)%3;
-    Face *t = pf->FFp(pi);
+    FaceType *t = pf->FFp(pi);
     if (t==pf) { border= true; break; }
     pi = pf->FFi(pi);
     pf = t;
@@ -662,11 +580,11 @@ void CollapseQuadDiag(typename Mesh::FaceType &f, const typename Mesh::CoordType
   // rotate around va, (counter-sense-as-face)-wise
   if (border) {
     int pi = fauxa;
-    Face* pf = fa; /* pf, pi could be a Pos<Face> p(pf, pi) */
+    FaceType* pf = fa; /* pf, pi could be a Pos<FaceType> p(pf, pi) */
     do {
       pi=(pi+1)%3;
       pf->V(pi) = va;
-      Face *t = pf->FFp(pi);
+      FaceType *t = pf->FFp(pi);
       if (t==pf) break;
       pi = pf->FFi(pi);
       pf = t;
@@ -674,39 +592,32 @@ void CollapseQuadDiag(typename Mesh::FaceType &f, const typename Mesh::CoordType
   }
     
   // update FF, delete faces
-  _CollapseQuadDiagHalf(*fb, fauxb, m);
-  _CollapseQuadDiagHalf(*fa, fauxa, m);
+  _CollapseDiagHalf(*fb, fauxb, m);
+  _CollapseDiagHalf(*fa, fauxa, m);
   
-  if (DELETE_VERTICES) Allocator<Mesh>::DeleteVertex(m,*vb);
-  va->P() = p;
+  if (DELETE_VERTICES) Allocator<MeshType>::DeleteVertex(m,*vb);
+  
+  
+
+  // for diagonals
+    
+  // for counterdiagonals
+  //Inpterpolator::Apply( *(f.V2(fauxa)), *(f.FFp( fauxa )->V2(fauxa)), interpol, va);
+  //va->P() = p;
 }      
 
 
-template <class Mesh>
-void CollapseQuadCounterDiag(typename Mesh::FaceType &f, typename Mesh::ScalarType k, Mesh& m){
-  typename Mesh::CoordType p;
-  int fauxa = FauxIndex(&f);
-  p = f.P2(fauxa)*(1-k) + f.FFp( fauxa )->P2( f.FFi( fauxa ) )*(k);
-  CollapseQuadCounterDiag(f,p,m);
-}
-
-template <class Mesh>
-void CollapseQuadCounterDiag(typename Mesh::FaceType &f, const typename Mesh::CoordType &p, Mesh& m){
-  FlipBitQuadDiag(f);
-  CollapseQuadDiag(f,p,m);
-}
 
 
 // helper function: find a good position on a diag to collapse a point
 // currently, it is point in the middle, 
 //    unless a mixed border-non border edge is collapsed, then it is an exreme
-template <class Face>
-typename Face::ScalarType PosOnDiag(const Face& f, bool counterDiag){
+static ScalarType PosOnDiag(const FaceType& f, bool counterDiag){
   bool b0, b1, b2, b3; // which side of the quads are border
   
-  const Face* fa=&f;
+  const FaceType* fa=&f;
   int ia = FauxIndex(fa);
-  const Face* fb=fa->cFFp(ia); 
+  const FaceType* fb=fa->cFFp(ia); 
   int ib = fa->cFFi(ia);
   
   b0 = fa->FFp((ia+1)%3) == fa;
@@ -715,8 +626,8 @@ typename Face::ScalarType PosOnDiag(const Face& f, bool counterDiag){
   b3 = fb->FFp((ib+2)%3) == fb;
   
   if (counterDiag) {
-    if (  (b0||b1) && !(b2||b3) ) return 0;
-    if ( !(b0||b1) &&  (b2||b3) ) return 1;
+    if (  (b0||b1) && !(b2||b3) ) return 1;
+    if ( !(b0||b1) &&  (b2||b3) ) return 0;
   } else {
     if (  (b1||b2) && !(b3||b0) ) return 0;
     if ( !(b1||b2) &&  (b3||b0) ) return 1;
@@ -725,10 +636,7 @@ typename Face::ScalarType PosOnDiag(const Face& f, bool counterDiag){
   return 0.5f;
 }
 
-template <class Mesh>
-void UpdateQualityAsBitQuadValency(Mesh& m){
-  typedef typename Mesh::FaceIterator FaceIterator;
-  typedef typename Mesh::VertexIterator VertexIterator;
+static void UpdateQualityAsValency(MeshType& m){
   for (VertexIterator vi = m.vert.begin();  vi!=m.vert.end(); vi++) if (!vi->IsD()) {
      vi->Q() = 0; 
   }
@@ -739,4 +647,106 @@ void UpdateQualityAsBitQuadValency(Mesh& m){
   }
 }
 
+private:
+  
+// helper function:
+// cos of angle abc. This should probably go elsewhere
+static ScalarType Cos(const CoordType &a, const CoordType &b, const CoordType &c )
+{
+  CoordType 
+    e0 = b - a,
+    e1 = b - c;
+  ScalarType d =  (e0.Norm()*e1.Norm());
+  if (d==0) return 0.0;
+  return (e0*e1)/d;
+}
+
+// helper function:
+// returns quality of a quad formed by points a,b,c,d
+// quality is computed as "how squared angles are"
+static ScalarType quadQuality(const CoordType &a, const CoordType &b, const CoordType &c, const CoordType &d){
+  ScalarType score = 0;
+  score += 1 - math::Abs( Cos( a,b,c) );
+  score += 1 - math::Abs( Cos( b,c,d) );
+  score += 1 - math::Abs( Cos( c,d,a) );
+  score += 1 - math::Abs( Cos( d,a,b) );
+  return score / 4;
+}
+
+// helper function:
+// returns quality of a given (potential) quad
+static ScalarType quadQuality(FaceType *f, int edge){
+  
+  CoordType 
+    a = f->V0(edge)->P(),
+    b = f->FFp(edge)->V2( f->FFi(edge) )->P(),
+    c = f->V1(edge)->P(),
+    d = f->V2(edge)->P();
+  
+  return quadQuality(a,b,c,d);
+
+}
+
+/**
+helper function:
+given a quad edge, retruns:
+   0 if that edge should not be rotated
+  +1 if it should be rotated clockwise (+1)
+  -1 if it should be rotated counterclockwise (-1)
+Currently an edge is rotated iff it is shortened by that rotations
+(shortcut criterion)
+*/
+static int TestEdgeRotation(const FaceType &f, int w0)
+{
+  const FaceType *fa = &f;
+  assert(! fa->IsF(w0) );
+  ScalarType q0,q1,q2;
+  CoordType v0,v1,v2,v3,v4,v5;
+  int w1 = (w0+1)%3;
+  int w2 = (w0+2)%3;
+  
+  v0 = fa->P(w0);
+  v3 = fa->P(w1);
+  
+  if (fa->IsF(w2) ) {
+    v1 = fa->cFFp(w2)->V2( fa->cFFi(w2) )->P();
+    v2 = fa->P(w2);
+  } else {
+    v1 = fa->P(w2);
+    v2 = fa->cFFp(w1)->V2( fa->cFFi(w1) )->P();
+  }
+  
+  const FaceType *fb = fa->cFFp(w0);
+  w0 = fa->cFFi(w0);
+  
+  w1 = (w0+1)%3;
+  w2 = (w0+2)%3;
+  if (fb->IsF(w2) ) {
+    v4 = fb->cFFp(w2)->V2( fb->cFFi(w2) )->P();
+    v5 = fb->P(w2);
+  } else {
+    v4 = fb->P(w2);
+    v5 = fb->cFFp(w1)->V2( fb->cFFi(w1) )->P();
+  }
+  
+  /*
+  //  max overall quality criterion:
+  q0 = quadQuality(v0,v1,v2,v3) +  quadQuality(v3,v4,v5,v0); // keep as is?
+  q1 = quadQuality(v1,v2,v3,v4) +  quadQuality(v4,v5,v0,v1); // rotate CW?
+  q2 = quadQuality(v5,v0,v1,v2) +  quadQuality(v2,v3,v4,v5); // rotate CCW?
+  
+  if (q0>=q1 && q0>=q2) return 0;
+  if (q1>=q2) return 1;*/
+  
+  // min distance (shortcut criterion)
+  q0 = (v0 - v3).SquaredNorm();
+  q1 = (v1 - v4).SquaredNorm();
+  q2 = (v5 - v2).SquaredNorm();
+  if (q0<=q1 && q0<=q2) return 0;
+  if (q1<=q2) return 1;
+  return -1;
+}
+
+
+};
 }} // end namespace vcg::tri