diff --git a/wrap/miq/MIQ.h b/wrap/miq/MIQ.h
index 430f4363..f5185328 100644
--- a/wrap/miq/MIQ.h
+++ b/wrap/miq/MIQ.h
@@ -1,11 +1,6 @@
 #ifndef __MIQ__
 #define __MIQ__
 
-
-#define SIZEQUADS 512
-#define V_SIZE 1
-#define SIZEPARA 1024
-
 #include <iostream>
 #include <vector>
 #include "mesh_type.h"
@@ -14,24 +9,15 @@
 #include "param_stats.h"
 #include "seams_initializer.h"
 #include "vertex_indexing.h"
-
+#include "stiffening.h"
+#include <vcg/complex/algorithms/clean.h>
 #include <wrap/io_trimesh/import.h>
 #include <wrap/io_trimesh/export.h>
 
 #define USECOMISO
 
-template <class ScalarType>
-ScalarType Gauss(ScalarType &value)
-{
-    const ScalarType E_NEPER=2.71828;
-    const ScalarType den=sqrt(2.0*M_PI);
-    ScalarType exponent=-0.5*pow(value,2);
-    ScalarType res=((1.0/den)*pow(E_NEPER,exponent));
-    return res;
-}
-
-template <class MeshType,class PolyMesh>
-class MIQ{
+template <class MeshType>
+class MIQ_parametrization{
 
 public:
   typename MeshType::template PerFaceAttributeHandle<float> Handle_Stiffness;
@@ -39,20 +25,18 @@ public:
   // Different stiffening mode
   enum StiffMode{NO_STIFF = 0,GAUSSIAN = 1,ITERATIVE = 2};
 
-  // Init
-  MIQ(MeshType &_mesh):mesh(_mesh),PSolver(mesh){};
-
-
   // Parametrize the mesh
-  void Solve(StiffMode stiffMode, double Stiffness = 5.0,
-             double GradientSize = 30.0, bool DirectRound = false,
-             int iter = 5, int localIter = 5, bool DoRound = true)
+  static void DoParameterize(MeshType &mesh,StiffMode stiffMode,
+                             double Stiffness = 5.0,double GradientSize = 30.0,
+                             bool DirectRound = false,int iter = 5,
+                             int localIter = 5, bool DoRound = true)
   {
+      PoissonSolver<MeshType> PSolver(mesh);
       if (mesh.fn==0)return;
-      InitDefaultStiffening();
+      StiffeningInitializer<MeshType>::InitDefaultStiffening(mesh);
       if (stiffMode==GAUSSIAN)
       {
-        AddGaussStiffening(Stiffness);
+        StiffeningInitializer<MeshType>::AddGaussStiffening(mesh,Stiffness);
         PSolver.SolvePoisson(GradientSize,1.f,DirectRound,localIter,DoRound);
       }
       else
@@ -62,9 +46,7 @@ public:
           {
             PSolver.SolvePoisson(GradientSize,1.f,DirectRound,localIter,DoRound);
             int nflips=NumFlips(mesh);
-            ScaleGLtoInt();
-            bool folded=updateStiffening(GradientSize);
-            ScaleInttoGL();
+            bool folded=StiffeningInitializer<MeshType>::updateStiffeningJacobianDistorsion(mesh,GradientSize);
             printf("ITERATION %d FLIPS %d \n",i,nflips);
             if (!folded)break;
           }
@@ -78,203 +60,43 @@ public:
       printf("**** END OPTIMIZING #FLIPS %d  ****\n",nflips);
       fflush(stdout);
       SelectFlippedFaces(mesh);
-      UpdateUVBox();
+
   }
 
-  // Generate a quad mesh starting from the parametrization
-  void Quadrangulate(PolyMesh &poly,double factor = 1)
-  {
-      Quad.Quadrangulate(mesh,poly,factor);
-  }
-
-  // Bounding box of the param domain
-  vcg::Box2<double> UVBox;
-
-  void ScaleGLtoInt()
-  {
-      double factor=(double)SIZEPARA/(double)SIZEQUADS;
-      for (unsigned int i=0;i<mesh.face.size();i++)
-      {
-          if (mesh.face[i].IsD())continue;
-          mesh.face[i].WT(0).P()*=factor;
-          mesh.face[i].WT(1).P()*=factor;
-          mesh.face[i].WT(2).P()*=factor;
-      }
-  }
-
-  void ScaleInttoGL()
-  {
-      double factor=(double)SIZEQUADS/(double)SIZEPARA;
-      for (unsigned int i=0;i<mesh.face.size();i++)
-      {
-          if (mesh.face[i].IsD())continue;
-          mesh.face[i].WT(0).P()*=factor;
-          mesh.face[i].WT(1).P()*=factor;
-          mesh.face[i].WT(2).P()*=factor;
-      }
-  }
-
-  void UpdateUVBox()
-  {
-      UVBox.SetNull();
-      for (unsigned int i=0;i<mesh.face.size();i++)
-      {
-          UVBox.Add((mesh.face[i].WT(0).P()));
-          UVBox.Add((mesh.face[i].WT(1).P()));
-          UVBox.Add((mesh.face[i].WT(2).P()));
-      }
-  }
-
-  // Quadrangulator
-  Quadrangulator<CMesh,PolyMesh> Quad;
-
-  void colorByStiffening(typename MeshType::ScalarType MaxVal=16)
-  {
-      bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
-      assert(hasStiffness);
-      for (unsigned int i=0;i<mesh.face.size();i++)
-      {
-          //CMesh::ScalarType val=MaxVal-mesh.face[i].stiffening+1;
-          CMesh::ScalarType val=MaxVal-Handle_Stiffness[i]+1;
-          if (val<1)val=1;
-          mesh.face[i].C()=vcg::Color4b::ColorRamp(1.0,MaxVal,val);
-      }
-  }
-
-private:
-
-  void AddStiffening(typename MeshType::ScalarType C,int radius=4)
-  {
-      bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
-      if(!hasStiffness)
-          Handle_Stiffness=vcg::tri::Allocator<CMesh>::AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));
-
-      bool hasSingular = vcg::tri::HasPerVertexAttribute(mesh,std::string("Singular"));
-      assert(hasSingular);
-      CMesh::PerVertexAttributeHandle<bool> Handle_Singular;
-      Handle_Singular=vcg::tri::Allocator<CMesh>::GetPerVertexAttribute<bool>(mesh,std::string("Singular"));
-
-      std::vector<CMesh::VertexType*> to_stiff;
-      for(unsigned int i=0;i<mesh.vert.size();i++)
-      {
-          CMesh::VertexType *v=&mesh.vert[i];
-          if (v->IsD())continue;
-          //if (!v->IsSingular())continue;
-          if (!Handle_Singular[v])continue;
-          to_stiff.push_back(v);
-      }
-      for(unsigned int i=0;i<mesh.face.size();i++)
-      {
-
-          CMesh::FaceType *f=&(mesh.face[i]);
-          if (f->IsD())continue;
-          if (!f->IsV())continue;
-          to_stiff.push_back(f->V(0));
-          to_stiff.push_back(f->V(1));
-          to_stiff.push_back(f->V(2));
-      }
-      std::sort(to_stiff.begin(),to_stiff.end());
-      std::vector<CMesh::VertexType*>::iterator new_end=std::unique(to_stiff.begin(),to_stiff.end());
-      int dist=distance(to_stiff.begin(),new_end);
-      to_stiff.resize(dist);
-      for (unsigned int i=0;i<to_stiff.size();i++)
-      {
-          CMesh::VertexType *v=to_stiff[i];
-          for (int r=0;r<radius;r++)
-          {
-              CMesh::ScalarType stiffVal=((CMesh::ScalarType)r)/(CMesh::ScalarType)radius;//((ScalarType)(radius-r))/(ScalarType)radius;
-              stiffVal*=3.0;
-              stiffVal=Gauss(stiffVal)/0.4;
-              stiffVal=1+(stiffVal*C);
-              std::vector<CMesh::FaceType*> ring;
-              //mesh.GetConnectedFaces(v,r,ring);
-              VFExtendedStarVF(v,r,ring);
-              ///then set stiffening
-              for (unsigned int k=0;k<ring.size();k++)
-              {
-                  CMesh::FaceType* f=ring[k];
-                  //if (f->stiffening<stiffVal)
-                  //    f->stiffening=stiffVal;
-                  if (Handle_Stiffness[f]<stiffVal)
-                      Handle_Stiffness[f]=stiffVal;
-              }
-          }
-      }
-  }
-
-
-  bool updateStiffening(typename MeshType::ScalarType grad_size)
-  {
-      bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
-      if(!hasStiffness)
-          Handle_Stiffness=vcg::tri::Allocator<CMesh>::AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));
-
-      bool flipped = NumFlips(mesh)>0;
-      //if (h == 0.0)
-      //    return flipped;
-      //
-      //assert(h != 0.0);
-
-      if (!flipped)
-          return false;
-      CMesh::ScalarType maxL=0;
-      CMesh::ScalarType maxD=0;
-      if (flipped)
-      {
-          const double c = 1.0;
-          const double d = 5.0;
-
-          for (unsigned int i = 0; i < mesh.face.size(); ++i)
-          {
-              CMesh::ScalarType dist=Distortion(mesh.face[i],grad_size);
-              if (dist>maxD)maxD=dist;
-              CMesh::ScalarType absLap=fabs(LaplaceDistortion(mesh.face[i], grad_size));
-              if (absLap>maxL)maxL=absLap;
-              CMesh::ScalarType stiffDelta = std::min(c * absLap, d);
-              //mesh.face[i].stiffening+=stiffDelta;
-              Handle_Stiffness[i]+=stiffDelta;
-          }
-      }
-      printf("Maximum Distorsion %4.4f \n",maxD);
-      printf("Maximum Laplacian %4.4f \n",maxL);
-      return flipped;
-  }
-
-  void InitDefaultStiffening()
-  {
-      bool hasStiffness = vcg::tri::HasPerFaceAttribute(mesh,std::string("Stiffness"));
-      if(!hasStiffness)
-          Handle_Stiffness=vcg::tri::Allocator<CMesh>::AddPerFaceAttribute<float>(mesh,std::string("Stiffness"));
-
-      for(unsigned int i=0;i<mesh.face.size();i++)
-      {
-          CMesh::FaceType *f=&(mesh.face[i]);
-          //f->stiffening=1;
-          Handle_Stiffness[f]=1;
-      }
-  }
-
-  void AddGaussStiffening(typename MeshType::ScalarType C)
-  {
-      int radius=floor(C);
-      if (C<4)radius=4;
-      AddStiffening(C,radius);
-  }
 public:
-  // Mesh class
-  MeshType &mesh;
 
-  // Quad mesh class
-  //CMesh quadmesh;
+  static bool IsValid(MeshType &mesh)
+  {
+      int n_comp=vcg::tri::Clean<MeshType>::CountConnectedComponents(mesh);
+      int non_manifE=vcg::tri::Clean<MeshType>::CountNonManifoldEdgeFF(mesh);
+      int non_manifV=vcg::tri::Clean<MeshType>::CountNonManifoldVertexFF(mesh);
+      return ((n_comp==1)&&(non_manifE==0)&&(non_manifV==0));
+  }
+  static void InitSeamsSing(MeshType &mesh,
+                            bool orient_globally,
+                            bool initMM,
+                             bool initCuts)
+  {
+      SeamsInitializer<MeshType> SInit;
+      SInit.Init(&mesh,orient_globally,initMM,initCuts);
+  }
 
-  ///seams initializer
-  SeamsInitializer<MeshType> SInit;
+  static void Parametrize(MeshType &mesh,StiffMode stiffMode,
+                          double Stiffness = 5.0,
+                          double GradientSize = 30.0,
+                          bool DirectRound = false,int iter = 5,
+                          int localIter = 5, bool DoRound = true)
+  {
+      VertexIndexing<MeshType> VInd;
 
-  // Vertex indexing class used for the solver
-  VertexIndexing<MeshType> VInd;
+      VInd.Init(&mesh);
+      VInd.InitMapping();
+      VInd.InitFaceIntegerVal();
+      VInd.InitSeamInfo();
+
+      DoParameterize(mesh,stiffMode,Stiffness,GradientSize,DirectRound,iter,localIter , DoRound);
+  }
 
-  // Poisson solver
-  PoissonSolver<MeshType> PSolver;
 
 };