diff --git a/vcg/complex/algorithms/create/platonic.h b/vcg/complex/algorithms/create/platonic.h
index 8309c0e1..9d654e52 100644
--- a/vcg/complex/algorithms/create/platonic.h
+++ b/vcg/complex/algorithms/create/platonic.h
@@ -522,6 +522,47 @@ void Cone( MeshType& in,
         }
 }
 
+template <class MeshType>
+void OrientedCone(MeshType & m,
+                  const typename MeshType::CoordType origin,
+                  const typename MeshType::CoordType end,
+                  const typename MeshType::ScalarType r1,
+                  const typename MeshType::ScalarType r2,
+                  const int SubDiv = 36  )
+{
+  typedef typename MeshType::ScalarType ScalarType;
+  typedef typename MeshType::CoordType CoordType;
+  typedef Matrix44<typename MeshType::ScalarType> Matrix44x;
+  Cone(m,r1,r2,Distance(origin,end),SubDiv);
+
+//  tri::UpdatePosition<MeshType>::Translate(m,CoordType(0,1,0));
+//  tri::UpdatePosition<MeshType>::Scale(m,CoordType(1,0.5f,1));
+//  tri::UpdatePosition<MeshType>::Scale(m,CoordType(xScale,1.0f,yScale));
+
+//  float height = Distance(origin,end);
+//  tri::UpdatePosition<MeshType>::Scale(m,CoordType(radius,height,radius));
+  CoordType norm = end-origin;
+  ScalarType angleRad = Angle(CoordType(0,1,0),norm);
+  const ScalarType Delta= 0.000000001;
+  Matrix44x rotM;
+  if (fabs(angleRad)<Delta)
+      rotM.SetIdentity();
+  else
+  if (fabs(angleRad-M_PI)<Delta)
+  {
+      CoordType axis = CoordType(0,0,1)^norm;
+      rotM.SetRotateRad(angleRad,axis);
+  }
+  else
+  {
+    CoordType axis = CoordType(0,1,0)^norm;
+    rotM.SetRotateRad(angleRad,axis);
+  }
+  tri::UpdatePosition<MeshType>::Matrix(m,rotM);
+  tri::UpdatePosition<MeshType>::Translate(m,origin);
+
+}
+
 
 template <class MeshType >
 void Box(MeshType &in, const typename MeshType::BoxType & bb )
@@ -1000,9 +1041,21 @@ void OrientedEllipticPrism(MeshType & m, const typename MeshType::CoordType orig
   tri::UpdatePosition<MeshType>::Scale(m,CoordType(radius,height,radius));
   CoordType norm = end-origin;
   ScalarType angleRad = Angle(CoordType(0,1,0),norm);
-  CoordType axis = CoordType(0,1,0)^norm;
+  const ScalarType Delta= 0.000000001;
   Matrix44x rotM;
-  rotM.SetRotateRad(angleRad,axis);
+  if (fabs(angleRad)<Delta)
+      rotM.SetIdentity();
+  else
+  if (fabs(angleRad-M_PI)<Delta)
+  {
+      CoordType axis = CoordType(0,0,1)^norm;
+      rotM.SetRotateRad(angleRad,axis);
+  }
+  else
+  {
+    CoordType axis = CoordType(0,1,0)^norm;
+    rotM.SetRotateRad(angleRad,axis);
+  }
   tri::UpdatePosition<MeshType>::Matrix(m,rotM);
   tri::UpdatePosition<MeshType>::Translate(m,origin);
 
diff --git a/vcg/complex/algorithms/curve_on_manifold.h b/vcg/complex/algorithms/curve_on_manifold.h
index c791667e..713e5498 100644
--- a/vcg/complex/algorithms/curve_on_manifold.h
+++ b/vcg/complex/algorithms/curve_on_manifold.h
@@ -201,7 +201,7 @@ public:
       {
         VertexPointer v0 = FindVertexSnap(f0,ip0);
         VertexPointer v1 = FindVertexSnap(f1,ip1);
-        assert(v1>0 && v0>0 && v0!=v1);
+        assert(v1 != NULL && v1 != NULL && v0 != v1);
         FacePointer ff0,ff1;
         int e0,e1;
         ret &= face::FindSharedFaces<FaceType>(v0,v1,ff0,ff1,e0,e1);
diff --git a/vcg/complex/algorithms/dual_meshing.h b/vcg/complex/algorithms/dual_meshing.h
index 8d15673a..1a0f8ede 100644
--- a/vcg/complex/algorithms/dual_meshing.h
+++ b/vcg/complex/algorithms/dual_meshing.h
@@ -27,6 +27,7 @@
 #include <vcg/complex/algorithms/update/topology.h>
 #include <vcg/simplex/face/pos.h>
 #include <vcg/complex/algorithms/clean.h>
+#include <vcg/complex/algorithms/polygonal_algorithms.h>
 
 namespace vcg {
 namespace tri {
@@ -61,65 +62,8 @@ class DualMeshing
             vertSeq.push_back(indexV);
         }
 
-        if (startV.IsB())
-        {
-            vcg::face::Pos<FaceType> firstPos=posVec[0];
-            firstPos.FlipE();
-            assert(firstPos.IsBorder());
-
-            int indexVt0=vcg::tri::Index(primal,firstPos.V());
-            int indexVt1=vcg::tri::Index(primal,firstPos.VFlip());
-            std::pair<int,int> key(std::min(indexVt0,indexVt1),
-                                   std::max(indexVt0,indexVt1));
-            assert(EdgeMap.count(key)>0);
-            int indexV0=EdgeMap[key];
-
-            vcg::face::Pos<FaceType> lastPos=posVec.back();
-            assert(lastPos.IsBorder());
-
-            indexVt0=vcg::tri::Index(primal,lastPos.V());
-            indexVt1=vcg::tri::Index(primal,lastPos.VFlip());
-            key=std::pair<int,int> (std::min(indexVt0,indexVt1),
-                                    std::max(indexVt0,indexVt1));
-            assert(EdgeMap.count(key)>0);
-            int indexV1=EdgeMap[key];
-
-            vertSeq.push_back(indexV1);
-            vertSeq.push_back(indexV0);
-        }
     }
 
-    static void CreateBorderEdgeVert(PolyMeshType &primal,
-                                     PolyMeshType &dual,
-                                     std::map<std::pair<int,int>, int> &VertMap)
-    {
-        VertMap.clear();
-        vcg::tri::UpdateFlags<PolyMeshType>::VertexClearB(primal);
-
-        for (size_t i=0;i<primal.face.size();i++)
-            for (int j=0;j<primal.face[i].VN();j++)
-            {
-                int edge_size=primal.face[i].VN();
-                FaceType *nextF=primal.face[i].cFFp(j);
-
-                if (nextF!=&primal.face[i])continue;
-
-                VertexType *v0=primal.face[i].V(j);
-                VertexType *v1=primal.face[i].V((j+1)%edge_size);
-
-                v0->SetB();
-                v1->SetB();
-
-                int V0Index=vcg::tri::Index(primal,v0);
-                int V1Index=vcg::tri::Index(primal,v1);
-                CoordType pos=(v0->P()+v1->P())/2;
-                vcg::tri::Allocator<PolyMeshType>::AddVertex(dual,pos);
-                std::pair<int,int> key(std::min(V0Index,V1Index),
-                                       std::max(V0Index,V1Index));
-
-                VertMap[key]=dual.vert.size()-1;
-            }
-    }
 
     static void CreateFaceVert(PolyMeshType &primal,
                                PolyMeshType &dual,
@@ -137,7 +81,6 @@ class DualMeshing
             int num=0;
             if (snapBorder)//search for border edge
             {
-                std::vector<CoordType> BorderPos;
                 for (int j=0;j<primal.face[i].VN();j++)
                 {
                    if (!primal.face[i].V(j)->IsB())continue;
@@ -145,7 +88,7 @@ class DualMeshing
                     num++;
                 }
                 if (num>0)
-                pos/=num;
+                    pos/=num;
             }
 
             if (num==0)
@@ -153,15 +96,22 @@ class DualMeshing
                 for (int j=0;j<primal.face[i].VN();j++)
                 {
                     pos+=primal.face[i].V(j)->P();
-                    int indexV=vcg::tri::Index(primal,primal.face[i].V(j));
-                    if (VertFace[indexV]!=-1)continue;
-                    VertFace[indexV]=i;
                 }
                 pos/=(ScalarType)primal.face[i].VN();
             }
             vcg::tri::Allocator<PolyMeshType>::AddVertex(dual,pos);
             VertMap[i]=dual.vert.size()-1;
+            dual.vert.back().Q()=i;
         }
+
+        //then initialize VF first face
+        for (size_t i=0;i<primal.face.size();i++)
+            for (int j=0;j<primal.face[i].VN();j++)
+            {
+                int indexV=vcg::tri::Index(primal,primal.face[i].V(j));
+                if (VertFace[indexV]!=-1)continue;
+                VertFace[indexV]=i;
+            }
     }
 
 public:
@@ -177,16 +127,18 @@ public:
         vcg::tri::RequireFFAdjacency(primal);
 
         vcg::tri::UpdateTopology<PolyMeshType>::FaceFace(primal);
+        vcg::tri::UpdateFlags<PolyMeshType>::VertexBorderFromFaceAdj(primal);
 
         std::map<std::pair<int,int>, int> VertEdgeMap;
-        CreateBorderEdgeVert(primal,dual,VertEdgeMap);
 
+        std::cout<<"Creating Dual Vertices"<<std::endl;
         std::vector<int> VertFaceMap,VertFace;
         CreateFaceVert(primal,dual,VertFaceMap,VertFace,snapBorder);
 
+        std::cout<<"Creating Dual Faces"<<std::endl;
         for (size_t i=0;i<primal.vert.size();i++)
         {
-            if ((snapBorder)&&(primal.vert[i].IsB()))continue;
+            if (primal.vert[i].IsB())continue;
 
             FaceType *firstF=&primal.face[VertFace[i]];
             std::vector<int> VertSeq;
@@ -195,9 +147,11 @@ public:
             dual.face.back().Alloc(VertSeq.size());
             for (size_t j=0;j<VertSeq.size();j++)
                 dual.face.back().V(j)=&dual.vert[VertSeq[j]];
+            dual.face.back().Q()=i;
         }
-
         vcg::tri::Clean<PolyMeshType>::RemoveUnreferencedVertex(dual);
+        //finally remove valence 1 vertices on the border
+        vcg::PolygonalAlgorithm<PolyMeshType>::RemoveValence2Vertices(dual);
     }
 
 };
diff --git a/wrap/igl/miq_parametrization.h b/wrap/igl/miq_parametrization.h
index 3003f523..6ed928fa 100644
--- a/wrap/igl/miq_parametrization.h
+++ b/wrap/igl/miq_parametrization.h
@@ -74,6 +74,8 @@ public:
         bool crease_as_feature;
         //true if roound selected vert
         bool round_selected;
+        //the anisotropy in MIQ sense (see paper)
+        double miqAnisotropy;
 
         MIQParameters()
         {
@@ -89,6 +91,7 @@ public:
             Ndir=4;
             crease_thr=0.2;
             hexaLine=false;
+            miqAnisotropy=1;
         }
     };
 
@@ -165,7 +168,7 @@ private:
 
         igl::copyleft::comiso::miq(V,F,X1,X2,UV,FUV,MiqP.gradient,MiqP.stiffness,MiqP.directRound,
                  MiqP.stiffness_iter,MiqP.local_iter,MiqP.doRound,MiqP.round_singularities,
-                 extra_round,hard_features);
+                 extra_round,hard_features,MiqP.miqAnisotropy);
 
         // then copy UV
         for (size_t i=0;i<trimesh.face.size();i++)
@@ -336,8 +339,8 @@ public:
     static void MIQParametrize(MeshType &trimesh,
                                MIQParameters &MiqP)
     {
-        if (MiqP.crease_as_feature)
-            SetCreases(trimesh,MiqP.crease_thr);
+//        if (MiqP.crease_as_feature)
+//            SetCreases(trimesh,MiqP.crease_thr);
 
         if (MiqP.Ndir==4)
             CrossFieldParam(trimesh,MiqP);