From 4fe40ecef515fab19378562738b2e354fd451133 Mon Sep 17 00:00:00 2001
From: nicopietroni <nico.pietroni@isti.cnr.it>
Date: Thu, 11 Sep 2014 20:42:57 +0000
Subject: [PATCH] completely restructured the whole class, much more robust

---
 wrap/miq/quadrangulator.h | 794 ++++++++++++++++++++++++--------------
 1 file changed, 503 insertions(+), 291 deletions(-)

diff --git a/wrap/miq/quadrangulator.h b/wrap/miq/quadrangulator.h
index 180c3b19..ff958ce0 100644
--- a/wrap/miq/quadrangulator.h
+++ b/wrap/miq/quadrangulator.h
@@ -8,32 +8,11 @@
 #include <vcg/complex/algorithms/refine.h>
 #include <vcg/complex/algorithms/smooth.h>
 #include <vcg/complex/algorithms/clean.h>
+#include <vcg/complex/algorithms/update/bounding.h>
+#include <wrap/io_trimesh/export.h>
+#include <vcg/complex/algorithms/update/texture.h>
 
-
-template <class MeshType>
-inline void ExtractVertex(const MeshType & srcMesh,
-                          const typename MeshType::FaceType & f,
-                          int whichWedge,
-                          const MeshType &dstMesh,
-                          typename MeshType::VertexType & v)
-{
-    (void)srcMesh;
-    (void)dstMesh;
-
-    //v.P() = f.cP(whichWedge);
-    v.ImportData(*f.cV(whichWedge));
-    v.T() = f.cWT(whichWedge);
-}
-
-template <class MeshType>
-inline bool CompareVertex(const MeshType & m,
-                          const typename MeshType::VertexType & vA,
-                          const typename MeshType::VertexType & vB)
-{
-    (void)m;
-    return ((vA.cT() == vB.cT())&&(vA.cP()==vB.cP()));
-}
-
+#define precisionQ 0.0000000001
 
 template <class TriMesh,class PolyMesh>
 class Quadrangulator
@@ -51,68 +30,265 @@ public:
     typedef typename PolyMesh::ScalarType PolyScalarType;
 
 
-    ///the set of all edges that belongs to integer lines
-    std::set<std::pair<TriFaceType*,int> > IntegerEdges;
+    struct InterpolationInfo
+    {
+        CoordType Pos3D;
+        vcg::Point2<ScalarType> PosUV;
+        ScalarType alpha;
+        bool to_split;
 
-    ///the set of all integer vertices and the other vertices on integer lines which is connectes to
-    std::map<TriVertexType*,std::vector<TriVertexType*> > IntegerLineAdj;
-    ///the set of integer vertices
-    std::set<TriVertexType*> IntegerVertices;
-    ///temporary polygons
-    std::vector<std::vector<TriVertexType *> > polygons;
+        InterpolationInfo()
+        {
+            Pos3D=CoordType(0,0,0);
+            PosUV=vcg::Point2<ScalarType>(0,0);
+            to_split=false;
+            alpha=-1;
+        }
+    };
 
-    ///drawing debug structures
-    std::vector<std::pair<TriFaceType*,int> > IntegerLines;
-    std::vector<TriVertexType*> IntegerVertex;
+    //the interpolation map that is saved once to be univoque per edge
+    typedef std::pair<CoordType,CoordType > KeyEdgeType;
+
+    std::map<KeyEdgeType,InterpolationInfo> InterpMap;
+
+    //ScalarType UVtolerance;
 
 private:
-    static bool ToSplit(const vcg::Point2<ScalarType> &uv0,
-                        const vcg::Point2<ScalarType> &uv1,
-                        int Dir,
-                        int IntegerLine,
-                        ScalarType &alpha,
-                        ScalarType tolerance=0.0001)
+
+    bool ToSplit(const vcg::Point2<ScalarType> &uv0,
+                 const vcg::Point2<ScalarType> &uv1,
+                 int Dir,
+                 ScalarType &alpha)
     {
-        ScalarType lineF=(ScalarType)IntegerLine;
-        ScalarType val0=std::min(uv0.V(Dir),uv1.V(Dir));
-        ScalarType val1=std::max(uv0.V(Dir),uv1.V(Dir));
-        if (lineF<(val0+tolerance))return false;
-        if (lineF>(val1-tolerance))return false;
-        ScalarType dist=fabs(uv0.V(Dir)-uv1.V(Dir));
-        if (dist<tolerance) return false;
-        alpha=(1.0-fabs(uv0.V(Dir)-lineF)/dist);
+        ScalarType val0=uv0.V(Dir);
+        ScalarType val1=uv1.V(Dir);
+        int IntegerLine0=floor(val0);
+        int IntegerLine1=floor(val1);
+        if (IntegerLine0==IntegerLine1)
+            return false;//no integer line pass throught the edge
+
+        bool swapped=false;
+        if (IntegerLine0>IntegerLine1)
+        {
+            std::swap(IntegerLine0,IntegerLine1);
+            std::swap(val0,val1);
+            assert(val1>=val0);
+            swapped=true;
+        }
+
+        //then get the first if extist that overcome the threshold
+        int IntegerSplit=IntegerLine0+1;
+        bool found=false;
+        ScalarType dist1,dist0;
+        for (int i=IntegerSplit;i<=IntegerLine1;i++)
+        {
+            dist1=fabs(val1-IntegerSplit);
+            dist0=fabs(val0-IntegerSplit);
+
+//            if ((dist0>=UVtolerance)&&
+//                (dist1>=UVtolerance))
+            if ((val0!=IntegerSplit)&&
+                (val1!=IntegerSplit))
+            {
+                found=true;
+                break;
+            }
+            IntegerSplit++;
+        }
+        if (!found)return false;
+
+        //have to check distance also in opposite direction
+        ScalarType lenght=val1-val0;
+        assert(lenght>=0);
+        //alpha=1.0-(dist/lenght);
+        alpha=(dist1/lenght);
+        if (swapped)alpha=1-alpha;
+        assert((alpha>0)&&(alpha<1));
         return true;
     }
 
-    ///return true if the edge has to be splitted,
-    ///by considering the tolerance to the closest integer
-    static bool ToSplit(const vcg::face::Pos<TriFaceType> &ep,
-                        ScalarType &alpha,
-                        ScalarType factor=1.0,
-                        ScalarType tolerance=0.0001)
+    void RoundInitial(TriMesh &to_split)
     {
-        //TriFaceType *f=ep.f;
-        //int z=ep.z;
-        TriVertexType* v0=ep.f->V(ep.z);
-        TriVertexType* v1=ep.f->V1(ep.z);
-        vcg::Point2<ScalarType> uv0=v0->T().P()*factor;
-        vcg::Point2<ScalarType> uv1=v1->T().P()*factor;
-
-        ///then test integer for each direction
-        for (int dir=0;dir<2;dir++)
+        ScalarType minTolerance=precisionQ;
+        //first add all eddge
+        for (int i=0;i<to_split.face.size();i++)
         {
-            int Int0=std::min((int)uv0.V(dir),(int)uv1.V(dir));
-            int Int1=std::max((int)uv0.V(dir),(int)uv1.V(dir));
-
-            for (int i=Int0;i<=Int1;i++)
+            TriFaceType *f=&to_split.face[i];
+            for (int j =0;j<3;j++)
             {
-                bool to_split=ToSplit(uv0,uv1,dir,i,alpha,tolerance);
-                if (to_split)return true;
+                vcg::Point2<ScalarType> UV=f->WT(j).P();
+
+                int int0=floor(UV.X()+0.5);
+                int int1=floor(UV.Y()+0.5);
+
+                ScalarType diff0=(fabs(UV.X()-(ScalarType)int0));
+                ScalarType diff1=(fabs(UV.Y()-(ScalarType)int1));
+
+                if (diff0<minTolerance)
+                    UV.X()=(ScalarType)int0;
+                if (diff1<minTolerance)
+                    UV.Y()=(ScalarType)int1;
+
+                f->WT(j).P()=UV;
             }
         }
-        return false;
     }
 
+    void RoundSplits(TriMesh &to_split,int dir)
+    {
+        ScalarType minTolerance=precisionQ;
+        //first add all eddge
+        for (int i=0;i<to_split.face.size();i++)
+        {
+            TriFaceType *f=&to_split.face[i];
+            for (int j =0;j<3;j++)
+            {
+                CoordType p0=f->P0(j);
+                CoordType p1=f->P1(j);
+                KeyEdgeType k(p0,p1);
+                assert(InterpMap.count(k)==1);
+                if (!InterpMap[k].to_split)continue;
+                //then get the intepolated value
+                vcg::Point2<ScalarType> UV=InterpMap[k].PosUV;
+
+                int int0=floor(UV.X()+0.5);
+                int int1=floor(UV.Y()+0.5);
+
+                ScalarType diff0=(fabs(UV.X()-(ScalarType)int0));
+                ScalarType diff1=(fabs(UV.Y()-(ScalarType)int1));
+
+                if (diff0<minTolerance)
+                    UV.X()=(ScalarType)int0;
+                if (diff1<minTolerance)
+                    UV.Y()=(ScalarType)int1;
+
+                InterpMap[k].PosUV=UV;
+            }
+        }
+    }
+
+    void InitSplitMap(TriMesh &to_split,
+                      int dir)
+    {
+       assert((dir==0)||(dir==1));
+       InterpMap.clear();
+       //printf("direction %d\n",dir );
+       //first add all eddge
+       for (int i=0;i<to_split.face.size();i++)
+       {
+           TriFaceType *f=&to_split.face[i];
+           for (int j =0;j<3;j++)
+           {
+               CoordType p0=f->P0(j);
+               CoordType p1=f->P1(j);
+               vcg::Point2<ScalarType> Uv0=f->V0(j)->T().P();
+               vcg::Point2<ScalarType> Uv1=f->V1(j)->T().P();
+               KeyEdgeType k(p0,p1);
+//               printf("p0 (%5.5f,%5.5f,%5.5f) p1(%5.5f,%5.5f,%5.5f) \n",p0.X(),p0.Y(),p0.Z(),p1.X(),p1.Y(),p1.Z());
+//               printf("uv0 (%5.5f,%5.5f) uv1(%5.5f,%5.5f) \n",Uv0.X(),Uv0.Y(),Uv1.X(),Uv1.Y());
+//               fflush(stdout);
+               assert(InterpMap.count(k)==0);
+               InterpMap[k]=InterpolationInfo();
+           }
+       }
+
+       //then set the ones to be splitted
+       for (int i=0;i<to_split.face.size();i++)
+       {
+           TriFaceType *f=&to_split.face[i];
+           for (int j =0;j<3;j++)
+           {
+               CoordType p0=f->P0(j);
+               CoordType p1=f->P1(j);
+               vcg::Point2<ScalarType> uv0=f->V0(j)->T().P();
+               vcg::Point2<ScalarType> uv1=f->V1(j)->T().P();
+
+               ScalarType alpha;
+               if (!ToSplit(uv0,uv1,dir,alpha))continue;
+
+               KeyEdgeType k(p0,p1);
+               assert(InterpMap.count(k)==1);
+               InterpMap[k].Pos3D=p0*alpha+p1*(1-alpha);
+               InterpMap[k].PosUV=uv0*alpha+uv1*(1-alpha);
+               InterpMap[k].to_split=true;
+               InterpMap[k].alpha=alpha;
+           }
+       }
+
+       //then make them coherent
+       for (int i=0;i<to_split.face.size();i++)
+       {
+           TriFaceType *f=&to_split.face[i];
+           for (int j =0;j<3;j++)
+           {
+               CoordType p0=f->P0(j);
+               CoordType p1=f->P1(j);
+               vcg::Point2<ScalarType> uv0=f->V0(j)->T().P();
+               vcg::Point2<ScalarType> uv1=f->V1(j)->T().P();
+//               if (p0>p1)continue; //only one verse of coherence
+
+               KeyEdgeType k0(p0,p1);
+               assert(InterpMap.count(k0)==1);//there should be already in the
+                                              //table and it should be coherent
+
+               KeyEdgeType k1(p1,p0);
+               if(InterpMap.count(k1)==0)continue;//REAL border, no need for update
+
+               bool to_split0=InterpMap[k0].to_split;
+               bool to_split1=InterpMap[k1].to_split;
+
+               //the find all possible cases
+               if ((!to_split0)&&(!to_split1))continue;
+
+               if ((to_split0)&&(to_split1))
+               {
+                   CoordType Pos3D=InterpMap[k1].Pos3D;
+                   InterpMap[k0].Pos3D=Pos3D;
+
+                   //check if need to make coherent also the UV Position
+                   //skip the fake border and do the rest
+                   bool IsBorderFF=(f->FFp(j)==f);
+
+                   if (!IsBorderFF) //in this case they should have same UVs
+                       InterpMap[k0].PosUV=InterpMap[k1].PosUV;
+                   else
+                   {
+                       ScalarType alpha=InterpMap[k1].alpha;
+                       assert((alpha>=0)&&(alpha<=1));
+                       alpha=1-alpha;
+                       InterpMap[k0].PosUV=alpha*uv0+(1-alpha)*uv1;
+                       InterpMap[k0].alpha=alpha;
+                   }
+
+               }
+               else
+               if ((!to_split0)&&(to_split1))
+               {
+                   CoordType Pos3D=InterpMap[k1].Pos3D;
+                   InterpMap[k0].Pos3D=Pos3D;
+
+                   //check if need to make coherent also the UV Position
+                   //skip the fake border and do the rest
+                   bool IsBorderFF=(f->FFp(j)==f);
+
+                   InterpMap[k0].to_split=true;
+
+                   if (!IsBorderFF) //in this case they should have same UVs
+                       InterpMap[k0].PosUV=InterpMap[k1].PosUV;
+                   else //recalculate , it pass across a seam
+                   {
+                       ScalarType alpha=InterpMap[k1].alpha;
+                       assert((alpha>=0)&&(alpha<=1));
+                       alpha=1-alpha;
+                       InterpMap[k0].PosUV=alpha*uv0+(1-alpha)*uv1;
+                       InterpMap[k0].alpha=alpha;
+                   }
+              }
+           }
+       }
+
+       RoundSplits(to_split,dir);
+    }
 
     // Basic subdivision class
     // This class must provide methods for finding the position of the newly created vertices
@@ -125,41 +301,38 @@ private:
         typedef typename MESH_TYPE::FaceType FaceType;
         typedef typename MESH_TYPE::CoordType CoordType;
 
-        ScalarType factor;
-        ScalarType tolerance;
-        ScalarType alpha;
+        std::map<KeyEdgeType,InterpolationInfo> *MapEdge;
 
         void operator()(typename MESH_TYPE::VertexType &nv,
                         vcg::face::Pos<typename MESH_TYPE::FaceType>  ep)
         {
-            bool to_split=ToSplit(ep,alpha,factor,tolerance);
+            VertexType* v0=ep.f->V0(ep.z);
+            VertexType* v1=ep.f->V1(ep.z);
+            assert(v0!=v1);
+
+            CoordType p0=v0->P();
+            CoordType p1=v1->P();
+            assert(p0!=p1);
+
+            KeyEdgeType k(p0,p1);
+            bool found=(MapEdge->count(k)==1);
+            assert(found);
+            bool to_split=(*MapEdge)[k].to_split;
             assert(to_split);
 
-            ///get the value on which the edge must be splitted
-            VertexType* v0=ep.f->V(ep.z);
-            VertexType* v1=ep.f->V1(ep.z);
-
-            nv.P()= v0->P()*alpha+v1->P()*(1.0-alpha);
+            //get the value on which the edge must be splitted
+            nv.P()= (*MapEdge)[k].Pos3D;
             //nv.N()= v0->N()*alpha+v1->N()*(1.0-alpha);
-            nv.T().P()=v0->T().P()*alpha+v1->T().P()*(1.0-alpha);
+            nv.T().P()=(*MapEdge)[k].PosUV;
         }
 
         vcg::TexCoord2<ScalarType> WedgeInterp(vcg::TexCoord2<ScalarType> &t0,
                                                vcg::TexCoord2<ScalarType> &t1)
         {
-            vcg::TexCoord2<ScalarType> tmp;
-//            if (t0.n() != t1.n())
-//                cerr << "Failed assertion: Quadrangulator::WedgeInterp1" << endl;
-//            // assert(t0.n()== t1.n()); TODO put back
-            tmp.n()=t0.n();
-            // assert(alpha>=0); TODO put back
-            if (alpha<0)
-                cerr << "Failed assertion: Quadrangulator::WedgeInterp2" << endl;
-            tmp.t()=(alpha*t0.t()+(1.0-alpha)*t1.t());
-            return tmp;
+            return (vcg::TexCoord2<ScalarType>(0,0));
         }
 
-        SplitMidPoint(){alpha=-1;}
+        SplitMidPoint(std::map<KeyEdgeType,InterpolationInfo> *_MapEdge){MapEdge=_MapEdge;}
     };
 
     template <class MESH_TYPE>
@@ -169,41 +342,56 @@ private:
         typedef typename MESH_TYPE::FaceType FaceType;
         typedef typename MESH_TYPE::ScalarType ScalarType;
 
+        std::map<KeyEdgeType,InterpolationInfo> *MapEdge;
+
     public:
-        ScalarType factor;
-        ScalarType tolerance;
 
         bool operator()(vcg::face::Pos<typename MESH_TYPE::FaceType> ep) const
         {
-            ScalarType alpha;
-            return(ToSplit(ep,alpha,factor,tolerance));
+            VertexType* v0=ep.f->V0(ep.z);
+            VertexType* v1=ep.f->V1(ep.z);
+            assert(v0!=v1);
+
+            CoordType p0=v0->P();
+            CoordType p1=v1->P();
+            assert(p0!=p1);
+
+            KeyEdgeType k(p0,p1);
+            bool found=(MapEdge->count(k)==1);
+            assert(found);
+            bool to_split=(*MapEdge)[k].to_split;
+            return(to_split);
         }
+
+        EdgePredicate(std::map<KeyEdgeType,InterpolationInfo> *_MapEdge){MapEdge=_MapEdge;}
     };
 
-    void SplitTris(TriMesh &to_split,
-                   ScalarType factor=1.0,
-                   ScalarType tolerance=0.0001)
+    void SplitTrisDir(TriMesh &to_split,
+                       int dir)
     {
         bool done=true;
-        SplitMidPoint<TriMesh> splMd;
-        EdgePredicate<TriMesh> eP;
-
-        splMd.tolerance=tolerance;
-        splMd.factor=factor;
-        eP.tolerance=tolerance;
-        eP.factor=factor;
-
+        //int step=0;
         while (done)
+        {
+            printf("Number of Vertices %d \n",to_split.vn);
+            fflush(stdout);
+
+            InitSplitMap(to_split,dir);
+
+            SplitMidPoint<TriMesh> splMd(&InterpMap);
+            EdgePredicate<TriMesh> eP(&InterpMap);
+
             done=vcg::tri::RefineE<TriMesh,SplitMidPoint<TriMesh>,EdgePredicate<TriMesh> >(to_split,splMd,eP);
 
-        for (unsigned int i=0;i<to_split.face.size();i++)
-            for (int j=0;j<3;j++) to_split.face[i].WT(j).P()=to_split.face[i].V(j)->T().P();
+        }
+        printf("Number of Vertices %d \n",to_split.vn);
+        fflush(stdout);
+         fflush(stdout);
     }
 
 
     bool IsOnIntegerLine(vcg::Point2<ScalarType> uv0,
-                         vcg::Point2<ScalarType> uv1,
-                         ScalarType tolerance=0.0001)
+                         vcg::Point2<ScalarType> uv1)
     {
         for (int dir=0;dir<2;dir++)
         {
@@ -212,210 +400,208 @@ private:
             int integer0=floor(uv0.V(dir)+0.5);
             int integer1=floor(uv1.V(dir)+0.5);
             if (integer0!=integer1)continue;
-            if ((fabs(val0-(ScalarType)integer0))>tolerance)continue;
-            if ((fabs(val1-(ScalarType)integer1))>tolerance)continue;
+//            if ((fabs(val0-(ScalarType)integer0))>=UVtolerance)continue;
+//            if ((fabs(val1-(ScalarType)integer1))>=UVtolerance)continue;
+            if (val0!=(ScalarType)floor(val0))continue;
+            if (val1!=(ScalarType)floor(val1))continue;
             return true;
         }
         return false;
     }
 
     bool IsOnIntegerVertex(vcg::Point2<ScalarType> uv,
-                           ScalarType tolerance=0.0001)
+                           bool IsB)
     {
+        int onIntegerL=0;
         for (int dir=0;dir<2;dir++)
         {
             ScalarType val0=uv.V(dir);
             int integer0=floor(val0+0.5);
-            if ((fabs(val0-(ScalarType)integer0))>tolerance)return false;
+            //if ((fabs(val0-(ScalarType)integer0))<UVtolerance)onIntegerL++;
+            if (val0==(ScalarType)floor(val0))onIntegerL++;
         }
-        return true;
+        if ((IsB)&&(onIntegerL>0))return true;
+        return (onIntegerL==2);
     }
 
-    void InitIntegerVectors()
-    {
-        IntegerLines=std::vector<std::pair<TriFaceType*,int> >(IntegerEdges.begin(),IntegerEdges.end());
-        IntegerVertex=std::vector<TriVertexType* > (IntegerVertices.begin(),IntegerVertices.end());
-    }
 
-    void EraseIntegerEdge(const vcg::face::Pos<TriFaceType> &ep)
+    void InitIntegerEdgesVert(TriMesh &Tmesh)
     {
-        std::pair<TriFaceType*,int> edge(ep.F(),ep.E());
-        assert(IntegerEdges.count(edge)!=0);
-        IntegerEdges.erase(edge);
-    }
+        //IntegerEdges.clear();
+        vcg::tri::UpdateFlags<TriMesh>::FaceSetF(Tmesh);
+        vcg::tri::UpdateFlags<TriMesh>::FaceClearS(Tmesh);
+        vcg::tri::UpdateFlags<TriMesh>::VertexClearS(Tmesh);
 
-    void EraseIntegerEdge(const std::vector<std::pair<TriFaceType*,int> > &to_erase)
-    {
-        for (unsigned int i=0;i<to_erase.size();i++)
-            IntegerEdges.erase(to_erase[i]);
-    }
-
-    void TestIntegerEdges()
-    {
-        typedef typename std::pair<TriFaceType*,int> pair_type;
-        typedef typename std::vector< pair_type > vect_type;
-        typename vect_type::iterator IteIntl;
-        for (IteIntl=IntegerLines.begin();
-             IteIntl!=IntegerLines.end();
-             IteIntl++)
-        {
-            int E=(*IteIntl).second;
-            TriFaceType *F=(*IteIntl).first;
-            TriFaceType *F1=F->FFp(E);
-            if (F==F1) continue;
-            int E1=F->FFi(E);
-            std::pair<TriFaceType*,int> curr_edge(F1,E1);
-            assert(IntegerEdges.count(curr_edge)!=0);
-        }
-    }
-
-    void InitIntegerEdgesVert(TriMesh &Tmesh,
-                              ScalarType factor=1.0,
-                              ScalarType tolerance=0.0001)
-    {
-        IntegerEdges.clear();
         for (unsigned int i=0;i<Tmesh.face.size();i++)
         {
             TriFaceType *f=&Tmesh.face[i];
             if (f->IsD())continue;
             for (int j=0;j<3;j++)
             {
-                TriFaceType *f1=f->FFp(j);
-                int e1=f->FFi(j);
                 bool IsBorder=f->IsB(j);
-                TriVertexType *v0=f->V0(j);
-                TriVertexType *v1=f->V1(j);
-                vcg::Point2<ScalarType> uv0=f->WT(j).P()*factor;
-                vcg::Point2<ScalarType> uv1=f->WT((j+1)%3).P()*factor;
-                if (IsOnIntegerLine(uv0,uv1,tolerance)||IsBorder)
+                if (IsBorder)
+                    f->ClearF(j);
+                else
                 {
-                    //IntegerEdges.insert(std::pair<TriVertexType*,TriVertexType*>(v0,v1));
-                    IntegerEdges.insert(std::pair<TriFaceType*,int>(f,j));
-                    if (!IsBorder)
-                        IntegerEdges.insert(std::pair<TriFaceType*,int>(f1,e1));
-                    else
+                    vcg::Point2<ScalarType> uv0=f->WT(j).P();
+                    vcg::Point2<ScalarType> uv1=f->WT((j+1)%3).P();
+
+                    if (IsOnIntegerLine(uv0,uv1))
                     {
-                        IntegerVertices.insert(v0);
-                        IntegerVertices.insert(v1);
+                        f->ClearF(j);
+                        TriFaceType *f1=f->FFp(j);
+                        int z=f->FFi(j);
+                        assert(f1!=f);
+                        f1->ClearF(z);
                     }
                 }
-                if (IsOnIntegerVertex(uv0))
-                    IntegerVertices.insert(v0);
 
-                if (IsOnIntegerVertex(uv1))
-                    IntegerVertices.insert(v1);
+                bool BorderV=f->V(j)->IsB();
 
+                if (IsOnIntegerVertex(f->WT(j).P(),BorderV))
+                    f->V(j)->SetS();
             }
         }
-        //InitIntegerNeigh(Tmesh);
-        InitIntegerVectors();
-        TestIntegerEdges();
     }
 
-
-    ///return the first and the last edge
-    ///following an integer line
-    ///until if reach anothe integer edge
-    bool OneIntegerStep(vcg::face::Pos<TriFaceType> &ep)
+    short int AlignmentEdge(TriFaceType *f,
+                      int edge_index)
     {
-        TriFaceType *f_init=ep.f;
-        TriFaceType *currF=f_init;
-        //int edge_init=ep.z;
-        //ep.V()=f_init->V(edge_init);
-        TriVertexType* v_init=ep.V();
-        bool complete_turn=false;
-        do
-        {
-            ep.FlipE();
-            ///see if found an integer vert
-            currF=ep.F();
-            int currE=ep.E();
-            assert((currE>=0)&&(currE<=4));
-
-            std::pair<TriFaceType*,int> curr_edge(currF,currE);
-
-            if (IntegerEdges.count(curr_edge)!=0)
-            {
-                ///go to the other side
-                ep.FlipV();
-                assert(ep.V()!=v_init);
-                return true;
-            }
-            ep.FlipF();
-            ///see if there's a border
-            bool jumped=(currF==ep.F());
-            if (jumped)
-                return false;
-            ///test the complete turn
-            complete_turn=(ep.F()==f_init);
-        }while (!complete_turn);
-        return false;
+        vcg::Point2<ScalarType> uv0=f->WT(edge_index).P();
+        vcg::Point2<ScalarType> uv1=f->WT((edge_index+1)%3).P();
+        if (uv0.X()==uv1.X())return 0;
+        if (uv0.Y()==uv1.Y())return 1;
+        return -1;
     }
 
-    ///find a polygon starting from half edge ep, return true if found
-    bool FindPolygon(vcg::face::Pos<TriFaceType> &ep,
-                     std::vector<TriVertexType*> &poly)
+    void FindPolygon(vcg::face::Pos<TriFaceType> &currPos,
+                     std::vector<TriVertexType *> &poly,
+                     std::vector<short int> &UVpoly)
     {
-
+        currPos.F()->SetV();
+        currPos.F()->C()=vcg::Color4b(255,0,0,255);
         poly.clear();
-        TriVertexType* v_init=ep.V();
-        ///it must start from an integer vert
-        assert(IntegerVertices.count(v_init)!=0);
-        poly.push_back(v_init);
-        std::vector<std::pair<TriFaceType*,int> > to_erase;
-        to_erase.push_back(std::pair<TriFaceType*,int>(ep.F(),ep.E()));
+        assert(currPos.V()->IsS());
+        TriVertexType *v_init=currPos.V();
+        poly.push_back(currPos.V());
+
+        //retrieve UV
+        int indexV0=currPos.E();
+
+        short int Align=AlignmentEdge(currPos.F(),currPos.E());
+
+        std::vector<short int> TempUVpoly;
+        TempUVpoly.push_back(Align);
+
         do
         {
-            bool done=OneIntegerStep(ep);
-            if (!done)
+            currPos.NextNotFaux();
+            currPos.F()->SetV();
+            currPos.F()->C()=vcg::Color4b(255,0,0,255);
+
+            if ((currPos.V()->IsS())&&(currPos.V()!=v_init))
             {
-                EraseIntegerEdge(to_erase);
-                return false;
+                poly.push_back(currPos.V());
+
+                short int Align=AlignmentEdge(currPos.F(),currPos.E());
+
+                TempUVpoly.push_back(Align);
             }
-            to_erase.push_back(std::pair<TriFaceType*,int>(ep.F(),ep.E()));
-            TriVertexType* v_curr=ep.V();
-            if ((IntegerVertices.count(v_curr)!=0)&&
-                    (v_curr!=v_init))
-                poly.push_back(v_curr);
-        }while(ep.V()!=v_init);
-        EraseIntegerEdge(to_erase);
-        return true;
+
+        }while (currPos.V()!=v_init);
+
+        //then shift the order of UV by one
+        //to be consistent with edge ordering
+        int size=TempUVpoly.size();
+        for (int i=0;i<size;i++)
+            UVpoly.push_back(TempUVpoly[(i+1)%size]);
     }
 
     void FindPolygons(TriMesh &Tmesh,
-                      std::vector<std::vector<TriVertexType *> > &polygons)
+                      std::vector<std::vector<TriVertexType *> > &polygons,
+                      std::vector<std::vector<short int> > &UV)
     {
-        //int limit=2;
+        vcg::tri::UpdateFlags<TriMesh>::FaceClearV(Tmesh);
         for (unsigned int i=0;i<Tmesh.face.size();i++)
         {
             TriFaceType * f=&Tmesh.face[i];
+            if (f->IsV())continue;
+
             for (int j=0;j<3;j++)
             {
                 TriVertexType* v0=f->V0(j);
-                //TriVertexType* v1=f->V1(j);
+                if (!v0->IsS())continue;
+                if (f->IsF(j))continue;
+
+                vcg::face::Pos<TriFaceType> startPos(f,j);
 
-                //std::pair<TriVertexType*,TriVertexType*> edge(v0,v1);*/
-                std::pair<TriFaceType*,int> edge(f,j);
-                if (IntegerEdges.count(edge)==0)continue;///edge already used or not integer
-                if (IntegerVertices.count(v0)==0)continue;	///must start from integer vert
-                ///create the pos
-                vcg::face::Pos<TriFaceType> ep(f,j);
                 std::vector<TriVertexType *> poly;
+                std::vector< short int> UVpoly;
 
-                bool found=FindPolygon(ep,poly);
-                if (found)
+                FindPolygon(startPos,poly,UVpoly);
+
+                if (poly.size()>2)
                 {
-                    std::reverse(poly.begin(),poly.end());///REVERSE ORDER
+                    assert(poly.size()==UVpoly.size());
                     polygons.push_back(poly);
+                    UV.push_back(UVpoly);
                 }
+                //only one polygon per initial face
+                break;
             }
         }
-
     }
 
-    void InitVertexQuadMesh(TriMesh &Tmesh)
+    //FUNCTIONS NEEDED BY "UV WEDGE TO VERTEX" FILTER
+    static void ExtractVertex(const TriMesh & srcMesh,
+                              const TriFaceType & f,
+                              int whichWedge,
+                              const TriMesh & dstMesh,
+                              TriVertexType & v)
     {
-        FindPolygons(Tmesh,polygons);
+        (void)srcMesh;
+        (void)dstMesh;
+        // This is done to preserve every single perVertex property
+        // perVextex Texture Coordinate is instead obtained from perWedge one.
+        v.ImportData(*f.cV(whichWedge));
+        v.T() = f.cWT(whichWedge);
+    }
+
+    static bool CompareVertex(const TriMesh & m,
+                              TriVertexType & vA,
+                              TriVertexType & vB)
+    {
+        (void)m;
+        return (vA.cT() == vB.cT());
+    }
+
+    void ConvertWTtoVT(TriMesh &Tmesh)
+    {
+        int vn = Tmesh.vn;
+        vcg::tri::AttributeSeam::SplitVertex(Tmesh, ExtractVertex, CompareVertex);
+        vcg::tri::UpdateTopology<TriMesh>::FaceFace(Tmesh);
+       // vcg::tri::UpdateFlags<TriMesh>::FaceBorderFromFF(Tmesh);
+    }
+
+    void ConvertVTtoWT(TriMesh &Tmesh)
+    {
+        vcg::tri::UpdateTexture<TriMesh>::WedgeTexFromVertexTex(Tmesh);
+        vcg::tri::Clean<TriMesh>::RemoveDuplicateVertex(Tmesh);
+    }
+
+    void ReupdateMesh(TriMesh &Tmesh)
+    {
+        vcg::tri::UpdateNormal<TriMesh>::PerFaceNormalized(Tmesh);	 // update Normals
+        vcg::tri::UpdateNormal<TriMesh>::PerVertexNormalized(Tmesh);// update Normals
+        //compact the mesh
+        vcg::tri::Allocator<TriMesh>::CompactVertexVector(Tmesh);
+        vcg::tri::Allocator<TriMesh>::CompactFaceVector(Tmesh);
+        vcg::tri::UpdateTopology<TriMesh>::FaceFace(Tmesh);			 // update Topology
+        vcg::tri::UpdateTopology<TriMesh>::TestFaceFace(Tmesh);		 //and test it
+        //set flags
+        vcg::tri::UpdateFlags<TriMesh>::VertexClearV(Tmesh);
+        vcg::tri::UpdateFlags<TriMesh>::FaceBorderFromFF(Tmesh);
+        vcg::tri::UpdateFlags<TriMesh>::VertexBorderFromFace(Tmesh);
     }
 
 public:
@@ -423,7 +609,9 @@ public:
 
     void TestIsProper(TriMesh &Tmesh)
     {
-        ///test manifoldness
+
+
+        //test manifoldness
         int test=vcg::tri::Clean<TriMesh>::CountNonManifoldVertexFF(Tmesh);
         //assert(test==0);
         if (test != 0)
@@ -450,59 +638,80 @@ public:
         }
     }
 
+
+
     void Quadrangulate(TriMesh &Tmesh,
                        PolyMesh &Pmesh,
-                       ScalarType factor=1.0,
-                       ScalarType tolerance=0.000001)
+                       std::vector< std::vector< short int> > &UV)
     {
         TestIsProper(Tmesh);
-        vcg::tri::AttributeSeam::SplitVertex(Tmesh, ExtractVertex<TriMesh>, CompareVertex<TriMesh>);
+
+        RoundInitial(Tmesh);
+
+        //UVtolerance=tolerance;
+
+        //split to per vert
+        ConvertWTtoVT(Tmesh);
+
+
         vcg::tri::Allocator<TriMesh>::CompactVertexVector(Tmesh);
         vcg::tri::Allocator<TriMesh>::CompactFaceVector(Tmesh);
         vcg::tri::UpdateTopology<TriMesh>::FaceFace(Tmesh);
         (void)Pmesh;
+        //TestIsProper(Tmesh);
+
+        //then split the tris along X
+        SplitTrisDir(Tmesh,0);
+        SplitTrisDir(Tmesh,1);
+
+        //merge back the mesh and WT coords
+        ConvertVTtoWT(Tmesh);
+
+        //CleanMesh(Pmesh);
+
+        //update properties of the mesh
+        ReupdateMesh(Tmesh);
+
+        //test manifoldness
         TestIsProper(Tmesh);
 
-        ///then split the tris
-        SplitTris(Tmesh,factor,tolerance);
-        ///join the vertices back!
-        //ScalarType EPS=(ScalarType)0.00000001;
-        ScalarType EPS=(ScalarType)0.000001;
-        vcg::tri::Clean<TriMesh>::MergeCloseVertex(Tmesh,EPS);
+        InitIntegerEdgesVert(Tmesh);
 
-        vcg::tri::UpdateNormal<TriMesh>::PerFaceNormalized(Tmesh);	 // update Normals
-        vcg::tri::UpdateNormal<TriMesh>::PerVertexNormalized(Tmesh);// update Normals
-        ///compact the mesh
-        vcg::tri::Allocator<TriMesh>::CompactVertexVector(Tmesh);
-        vcg::tri::Allocator<TriMesh>::CompactFaceVector(Tmesh);
-        vcg::tri::UpdateTopology<TriMesh>::FaceFace(Tmesh);			 // update Topology
-        vcg::tri::UpdateTopology<TriMesh>::TestFaceFace(Tmesh);		 //and test it
-        ///set flags
-        vcg::tri::UpdateFlags<TriMesh>::VertexClearV(Tmesh);
-        vcg::tri::UpdateFlags<TriMesh>::FaceBorderFromFF(Tmesh);
-        vcg::tri::UpdateFlags<TriMesh>::VertexBorderFromFace(Tmesh);
-        ///test manifoldness
-        TestIsProper(Tmesh);
+//        int E3=vcg::tri::io::ExporterPLY<TriMesh>::Save(Tmesh,"./maremma_hane.ply",
+//                                                        vcg::tri::io::Mask::IOM_FACEFLAGS|
+//                                                        vcg::tri::io::Mask::IOM_VERTFLAGS);
 
-        vcg::tri::UpdateFlags<TriMesh>::VertexClearV(Tmesh);
+        for (int i=0;i<Tmesh.face.size();i++)
+            Tmesh.face[i].C()=vcg::Color4b(255,255,255,255);
 
-        InitIntegerEdgesVert(Tmesh,factor,tolerance);
-        InitVertexQuadMesh(Tmesh);
+        std::vector<std::vector<TriVertexType *> > polygons;
+        FindPolygons(Tmesh,polygons,UV);
 
-        ///then add to the polygonal mesh
+        //then add to the polygonal mesh
         Pmesh.Clear();
-        ///first create vertices
-        vcg::tri::Allocator<PolyMesh>::AddVertices(Pmesh,IntegerVertex.size());
-        std::map<TriVertexType*,int> VertMap;
-        for(unsigned int i=0;i<IntegerVertex.size();i++)
+
+        int numV=vcg::tri::UpdateSelection<TriMesh>::VertexCount(Tmesh);
+
+        //first create vertices
+        vcg::tri::Allocator<PolyMesh>::AddVertices(Pmesh,numV);
+
+        std::map<CoordType,int> VertMap;
+        int index=0;
+        for(unsigned int i=0;i<Tmesh.vert.size();i++)
         {
-            CoordType pos=IntegerVertex[i]->P();
-            CoordType norm=IntegerVertex[i]->N();
-            Pmesh.vert[i].P()=typename PolyMesh::CoordType(pos.X(),pos.Y(),pos.Z());
-            Pmesh.vert[i].N()=typename PolyMesh::CoordType(norm.X(),norm.Y(),norm.Z());
-            VertMap[IntegerVertex[i]]=i;
+            if (!Tmesh.vert[i].IsS())continue;
+
+            CoordType pos=Tmesh.vert[i].P();
+            CoordType norm=Tmesh.vert[i].N();
+            vcg::Point2<ScalarType> UV=Tmesh.vert[i].T().P();
+            Pmesh.vert[index].P()=typename PolyMesh::CoordType(pos.X(),pos.Y(),pos.Z());
+            Pmesh.vert[index].N()=typename PolyMesh::CoordType(norm.X(),norm.Y(),norm.Z());
+            Pmesh.vert[index].T().P()=UV;
+            VertMap[pos]=index;
+            index++;
         }
-        ///then add polygonal mesh
+
+        //then add polygonal mesh
         vcg::tri::Allocator<PolyMesh>::AddFaces(Pmesh,polygons.size());
         for (unsigned int i=0;i<polygons.size();i++)
         {
@@ -510,11 +719,14 @@ public:
             Pmesh.face[i].Alloc(size);
             for (int j=0;j<size;j++)
             {
-                TriVertexType* v=polygons[i][j];
-                int index=VertMap[v];
+                CoordType pos=(polygons[i][j])->P();
+                assert(VertMap.count(pos)==1);
+                int index=VertMap[pos];
                 Pmesh.face[i].V(j)=&(Pmesh.vert[index]);
             }
         }
+
+        //E3=vcg::tri::io::ExporterOBJ<PolyMesh>::Save(Pmesh,"./maremma_hane.obj",0);
     }
 };