diff --git a/vcg/complex/trimesh/bitquad_support.h b/vcg/complex/trimesh/bitquad_support.h
index 531a1093..e6c49d11 100644
--- a/vcg/complex/trimesh/bitquad_support.h
+++ b/vcg/complex/trimesh/bitquad_support.h
@@ -1,8 +1,10 @@
 #ifndef VCG_BITQUAD_SUPPORT
 #define VCG_BITQUAD_SUPPORT
 #include <vector>
+#include <set>
 #include <vcg/complex/trimesh/subset.h>
 #include <vcg/simplex/face/jumping_pos.h>
+#include <vcg/simplex/face/topology.h>
 #include <vcg/space/planar_polygon_tessellation.h>
 
 /** BIT-QUAD creation support:
@@ -91,6 +93,7 @@ typedef typename MeshType::FaceType* FaceTypeP;
 typedef typename MeshType::VertexType VertexType;
 typedef typename MeshType::FaceIterator FaceIterator;
 typedef typename MeshType::VertexIterator VertexIterator;
+typedef typename MeshType::VertexPointer VertexPointer;
 
 class Pos{
   FaceType *f;
@@ -1126,7 +1129,51 @@ static ScalarType Cos(const CoordType &a, const CoordType &b, const CoordType &c
   if (d==0) return 0.0;
   return (e0*e1)/d;
 }
+public:
+/**
+  Generic quad triangulation function.
+  It take in input 4 vertex pointrs and rotate them so that a simple fan triangulation is Ok.
+  It uses geometric criteria for avoiding bad shaped triangles, and folds
+  and it use an internal set of already created diagonal to avoid the creation of non manifold situations.
+  At the begin you shoud call this function with an empty vector to reset the set of existing diagonals.
+  */
+static void QuadTriangulate(std::vector<VertexPointer> &q)
+{
+  typedef typename std::set<std::pair<VertexPointer,VertexPointer> > diagSetType;
+  static diagSetType diagSet; // the set of already created diagonals
+  if(q.size()!=4)
+  {
+    diagSet.clear();
+    return;
+  }
+  const CoordType &P0=q[0]->cP();
+  const CoordType &P1=q[1]->cP();
+  const CoordType &P2=q[2]->cP();
+  const CoordType &P3=q[3]->cP();
 
+  CoordType N00 = Normal(P0,P1,P2);
+  CoordType N01 = Normal(P0,P2,P3);
+  CoordType N10 = Normal(P1,P2,P3);
+  CoordType N11 = Normal(P1,P3,P0);
+
+  ScalarType Angle0Rad=Angle(N00,N01);
+  ScalarType Angle1Rad=Angle(N10,N11);
+
+  // QualityRadii is inradius/circumradius; bad when close to zero.
+  // swap diagonal if the worst triangle improve.
+  bool qualityImprove = std::min(QualityRadii(P0,P1,P2),QualityRadii(P0,P2,P3)) < std::min(QualityRadii(P1,P2,P3),QualityRadii(P1,P3,P0));
+  bool swapCauseFlip = (Angle1Rad > M_PI/2.0) && (Angle0Rad <M_PI/2.0);
+
+  if ( qualityImprove && ! swapCauseFlip)
+         std::rotate(q.begin(), q.begin()+1, q.end());
+
+  std::pair<typename diagSetType::iterator,bool> res;
+  if(q[0]<q[2]) res= diagSet.insert(std::make_pair(q[0],q[2]));
+  else res= diagSet.insert(std::make_pair(q[2],q[0]));
+
+  if(!res.second) // res.second is false if an element with the same value existed; in that case rotate again!
+    std::rotate(q.begin(), q.begin()+1, q.end());
+}
 };
 }} // end namespace vcg::tri
 
diff --git a/wrap/io_trimesh/import_off.h b/wrap/io_trimesh/import_off.h
index 610a40ef..d340c0e1 100644
--- a/wrap/io_trimesh/import_off.h
+++ b/wrap/io_trimesh/import_off.h
@@ -29,6 +29,7 @@
 #include <assert.h>
 #include <vcg/space/color4.h>
 #include<vcg/complex/trimesh/allocate.h>
+#include<vcg/complex/trimesh/bitquad_support.h>
 #include <wrap/callback.h>
 #include <wrap/io_trimesh/io_mask.h>
 
@@ -371,6 +372,12 @@ namespace vcg
 
 					Allocator<MESH_TYPE>::AddFaces(mesh, nFaces);
 					unsigned int f0=0;
+
+
+          // Initial call to the QuadTriangulate with an empty vector to just reset the static set of existing diagonals
+          std::vector<VertexPointer> qtmp;
+          BitQuad<MESH_TYPE>::QuadTriangulate(qtmp);
+
 					for (unsigned int f=0; f < nFaces; f++)
 					{
 						f0 = f;
@@ -419,30 +426,13 @@ namespace vcg
 								k++;
 							}
                             if(vert_per_face==4)
-                            {   // To well triangulate the quad:
-                                // if the quad is flat and convex we use the shortest diag ,
-                                // else we should use the diag that makes the smallest                                
-                                
-                                const CoordType &P0=mesh.vert[vertIndices[0]].cP();
-                                const CoordType &P1=mesh.vert[vertIndices[1]].cP();
-                                const CoordType &P2=mesh.vert[vertIndices[2]].cP();
-                                const CoordType &P3=mesh.vert[vertIndices[3]].cP();
-
-                                CoordType N00 = Normal(P0,P1,P2);
-                                CoordType N01 = Normal(P0,P2,P3);
-                                CoordType N10 = Normal(P1,P2,P3);
-                                CoordType N11 = Normal(P1,P3,P0);
-
-                                ScalarType Angle0Rad=Angle(N00,N01);
-                                ScalarType Angle1Rad=Angle(N10,N11);
-
-                                // QualityRadii is inradius/circumradius; bad when close to zero. 
-                                // swap diagonal if the worst triangle improve. 
-                                bool qualityImprove = std::min(QualityRadii(P0,P1,P2),QualityRadii(P0,P2,P3)) < std::min(QualityRadii(P1,P2,P3),QualityRadii(P1,P3,P0));
-                                bool swapCauseFlip = (Angle1Rad > M_PI/2.0) && (Angle0Rad <M_PI/2.0);
-                                if(qualityImprove && ! swapCauseFlip)
-                                        std::rotate(vertIndices.begin(), vertIndices.begin()+1, vertIndices.end());
-
+                            {   // To well triangulate use the bitquad support function that reorders vertex for a simple fan
+                                std::vector<VertexPointer> q(4);
+                                for(int qqi=0;qqi<4;++qqi)
+                                  q[qqi]=& mesh.vert[vertIndices[qqi]];
+                                BitQuad<MESH_TYPE>::QuadTriangulate(q);
+                                for(int qqi=0;qqi<4;++qqi)
+                                  vertIndices[qqi] = q[qqi]- & mesh.vert[0];
                             }
                             // standard fan triangulation (we hope the polygon is convex...)
                             for (int j=0; j<=vert_per_face-3; j++)