* clean up on miq code

* added generic rotate_vectors


Former-commit-id: a58d8cec628ae0f189ba7b612ee8725b6c57cbc1

include/igl/comiso/miq.h

@@ -6,19 +6,32 @@
 
 namespace igl
 {
-  // Creates a quad mesh from a triangular mesh and a set of two directions
-  // per face, using the algorithm described in the paper
+  // Global seamless parametrization aligned with a given per-face jacobian (PD1,PD2).
+  // The algorithm is based on
   // "Mixed-Integer Quadrangulation" by D. Bommes, H. Zimmer, L. Kobbelt
   // ACM SIGGRAPH 2009, Article No. 77 (http://dl.acm.org/citation.cfm?id=1531383)
 
   // Inputs:
-  //   Vin        #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F          #F by 4 eigen Matrix of face (quad) indices
-  //   maxIter    maximum numbers of iterations
-  //   threshold  minimum allowed threshold for non-planarity
+  //   V              #V by 3 list of mesh vertex 3D positions
+  //   F              #F by 3 list of faces indices in V
+  //   PD1            #V by 3 first line of the Jacobian per triangle
+  //   PD2            #V by 3 second line of the Jacobian per triangle
+  //                  (optional, if empty it will be a vector in the tangent plane orthogonal to PD1)
+  //   scale          global scaling for the gradient (controls the quads resolution)
+  //   stiffness      weight for the stiffness iterations
+  //   direct_round   greedily round all integer variables at once (greatly improves optimization speed but lowers quality)
+  //   iter           stiffness iterations (0 = no stiffness)
+  //   local_iter     number of local iterations for the integer rounding
+  //   do_round       enables the integer rounding (disabling it could be useful for debugging)
+  //   round_vertices id of additional vertices that should be snapped to integer coordinates
+  //   hard_features  #H by 2 list of pairs of vertices that belongs to edges that should be snapped to integer coordinates
+  //
   // Output:
-  //   Vout       #V by 3 eigen Matrix of planar mesh vertex 3D positions
+  //   UV             #UV by 2 list of vertices in 2D
+  //   FUV            #FUV by 3 list of face indices in UV
   //
+  // TODO: rename the parameters name in the cpp consistenly
+  //       improve the handling of hard_features, right now it might fail in difficult cases
 
   template <typename DerivedV, typename DerivedF, typename DerivedU>
   IGL_INLINE void miq(const Eigen::PlainObjectBase<DerivedV> &V,
@@ -27,13 +40,23 @@ namespace igl
                                               const Eigen::PlainObjectBase<DerivedV> &PD2,
                                               Eigen::PlainObjectBase<DerivedU> &UV,
                                               Eigen::PlainObjectBase<DerivedF> &FUV,
-                                              double GradientSize = 30.0,
-                                              double Stiffness = 5.0,
-                                              bool DirectRound = false,
+                                              double scale = 30.0,
+                                              double stiffness = 5.0,
+                                              bool direct_round = false,
                                               int iter = 5,
-                                              int localIter = 5, bool DoRound = true,
-                                              std::vector<int> roundVertices = std::vector<int>(),
-                                              std::vector<std::vector<int> > hardFeatures = std::vector<std::vector<int> >());
+                                              int local_iter = 5,
+                                              bool DoRound = true,
+                                              std::vector<int> round_vertices = std::vector<int>(),
+                                              std::vector<std::vector<int> > hard_features = std::vector<std::vector<int> >());
+
+  // Helper function that allows to directly provided pre-combed bisectors for an already cut mesh
+  // Additional input:
+  // PD1_combed, PD2_combed  :   #F by 3 combed jacobian
+  // BIS1_combed, BIS2_combed:   #F by 3 pre combed bi-sectors
+  // MMatch:                     #F by 3 list of per-corner integer PI/2 rotations
+  // Singular:                   #V list of flag that denotes if a vertex is singular or not
+  // SingularDegree:             #V list of flag that denotes the degree of the singularity
+  // Seams:                      #F by 3 list of per-corner flag that denotes seams
 
   template <typename DerivedV, typename DerivedF, typename DerivedU>
   IGL_INLINE void miq(const Eigen::PlainObjectBase<DerivedV> &V,
@@ -42,10 +65,10 @@ namespace igl
                                               const Eigen::PlainObjectBase<DerivedV> &PD2_combed,
                                               const Eigen::PlainObjectBase<DerivedV> &BIS1_combed,
                                               const Eigen::PlainObjectBase<DerivedV> &BIS2_combed,
-                                              const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_MMatch,
-                                              const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_Singular,
-                                              const Eigen::Matrix<int, Eigen::Dynamic, 1> &Handle_SingularDegree,
-                                              const Eigen::Matrix<int, Eigen::Dynamic, 3> &Handle_Seams,
+                                              const Eigen::Matrix<int, Eigen::Dynamic, 3> &MMatch,
+                                              const Eigen::Matrix<int, Eigen::Dynamic, 1> &Singular,
+                                              const Eigen::Matrix<int, Eigen::Dynamic, 1> &SingularDegree,
+                                              const Eigen::Matrix<int, Eigen::Dynamic, 3> &Seams,
                                               Eigen::PlainObjectBase<DerivedU> &UV,
                                               Eigen::PlainObjectBase<DerivedF> &FUV,
                                               double GradientSize = 30.0,

include/igl/rotate_vectors.cpp

@@ -0,0 +1,35 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "rotate_vectors.h"
+IGL_INLINE Eigen::MatrixXd igl::rotate_vectors(
+                    const Eigen::MatrixXd& V,
+                    const Eigen::VectorXd& A,
+                    const Eigen::MatrixXd& B1,
+                    const Eigen::MatrixXd& B2)
+{
+  Eigen::MatrixXd RV(V.rows(),V.cols());
+
+  for (unsigned i=0; i<V.rows();++i)
+  {
+    // project onto the tangent plane and convert to angle
+    double a = atan2(B2.row(i).dot(V.row(i)),B1.row(i).dot(V.row(i)));
+
+    // rotate
+    a += A.size() == 1 ? A(0) : A(1);
+
+    // move it back to global coordinates
+    RV.row(i) = cos(a) * B1.row(i) + sin(a) * B2.row(i);
+  }
+
+  return RV;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif

include/igl/rotate_vectors.h

@@ -0,0 +1,37 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef IGL_ROTATE_VECTORS_H
+#define IGL_ROTATE_VECTORS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Rotate the vectors V by A radiants on the tangent plane spanned by B1 and B2
+  // Inputs:
+  //   V     #V by 3 eigen Matrix of vectors
+  //   A     #V eigen vector of rotation angles or a single angle to be applied to all vectors
+  //   B1    #V by 3 eigen Matrix of base vector 1
+  //   B2    #V by 3 eigen Matrix of base vector 2
+  //
+  // Output:
+  //   Returns the rotated vectors
+  //
+  IGL_INLINE Eigen::MatrixXd rotate_vectors(
+                                            const Eigen::MatrixXd& V,
+                                            const Eigen::VectorXd& A,
+                                            const Eigen::MatrixXd& B1,
+                                            const Eigen::MatrixXd& B2);
+
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "rotate_vectors.cpp"
+#endif
+
+#endif

tutorial/505_MIQ/main.cpp

@@ -1,12 +1,26 @@
-#define IGL_HEADER_ONLY
-#include <igl/readOBJ.h>
+#include <igl/readOFF.h>
 #include <igl/viewer/Viewer.h>
 #include <igl/comiso/miq.h>
 #include <igl/barycenter.h>
 #include <igl/avg_edge_length.h>
+#include <igl/comiso/nrosy.h>
 #include <sstream>
+#include <igl/rotate_vectors.h>
 
 
+Eigen::VectorXi Seams;
+
+// Cuts
+Eigen::VectorXi C;
+
+// Singularities
+Eigen::VectorXd S;
+
+// Cross field
+Eigen::MatrixXd X;
+Eigen::MatrixXd X2;
+
+// Create a texture that hides the integer translation in the parametrization
 void line_texture(Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic> &texture_R,
                   Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic> &texture_G,
                   Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic> &texture_B)
@@ -26,88 +40,47 @@ void line_texture(Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic> &texture_R,
     texture_B = texture_R;
   }
 
-
-bool readPolyVf(const char *fname,
-                Eigen::VectorXi &isConstrained,
-                std::vector<Eigen::MatrixXd> &polyVF)
-{
-  FILE *fp = fopen(fname,"r");
-  if (!fp)
-    return false;
-  int degree, numF;
-  if (fscanf(fp,"%d %d", &degree, &numF) !=2)
-    return false;
-  polyVF.resize(degree, Eigen::MatrixXd::Zero(numF, 3));
-  isConstrained.setZero(numF,1);
-  int vali; float u0,u1,u2;
-  for (int i = 0; i<numF; ++i)
-  {
-    if (fscanf(fp,"%d", &vali)!=1)
-      return false;
-    isConstrained[i] = vali;
-    for (int j = 0; j<degree; ++j)
-    {
-      if (fscanf(fp,"%g %g %g", &u0, &u1, &u2) !=3)
-        return false;
-      polyVF[j](i,0) = u0;
-      polyVF[j](i,1) = u1;
-      polyVF[j](i,2) = u2;
-    }
-  }
-  fclose(fp);
-  return true;
-}
-
-void writePolyVf(const char *fname,
-                 const Eigen::VectorXi &isConstrained,
-                 const std::vector<Eigen::MatrixXd> &polyVF)
-{
-  int numF = polyVF[0].rows();
-  int degree = polyVF.size();
-  FILE *fp = fopen(fname,"w");
-  fprintf(fp,"%d %d\n", degree,numF);
-  for (int i = 0; i<numF; ++i)
-  {
-    fprintf(fp,"%d ", isConstrained[i]);
-    for (int j = 0; j<degree; ++j)
-      fprintf(fp,"%.15g %.15g %.15g ", polyVF[j](i,0), polyVF[j](i,1), polyVF[j](i,2));
-    fprintf(fp, "\n");
-  }
-  fclose(fp);
-
-}
-
-
 int main(int argc, char *argv[])
 {
+  using namespace Eigen;
+
   Eigen::MatrixXd V;
   Eigen::MatrixXi F;
-  // Load a mesh in OFF format
-  igl::readOBJ("../shared/lilium.obj", V, F);
-
-
-  // Load a frame field
-  Eigen::VectorXi isConstrained;
-  std::vector<Eigen::MatrixXd> polyVF;
-  readPolyVf("../shared/lilium.crossfield", isConstrained, polyVF);
 
+  // Load a mesh in OFF format
+  igl::readOFF("../shared/3holes.off", V, F);
+
+  // Contrain one face
+  VectorXi b(1);
+  b << 0;
+  MatrixXd bc(1,3);
+  bc << 1, 0, 0;
+
+  // Create a smooth 4-RoSy field
+  igl::nrosy(V,F,b,bc,VectorXi(),VectorXd(),MatrixXd(),4,0.5,X,S);
+
+  // Find the the orthogonal vector
+  MatrixXd B1,B2,B3;
+  igl::local_basis(V,F,B1,B2,B3);
+  X2 = igl::rotate_vectors(X, VectorXd::Constant(1,M_PI/2), B1, B2);
+  
   Eigen::MatrixXd UV;
   Eigen::MatrixXi FUV;
 
-  double gradientSize = 50;
-  double quadIter = 0;
+  double gradient_size = 50;
+  double iter = 0;
   double stiffness = 5.0;
-  bool directRound = 1;
+  bool direct_round = 0;
   igl::miq(V,
            F,
-           polyVF[0],
-           polyVF[1],
+           X,
+           X2,
            UV,
            FUV,
-           gradientSize,
+           gradient_size,
            stiffness,
-           directRound,
-           quadIter);
+           direct_round,
+           iter);
 
 
   // Face barycenters
@@ -121,8 +94,8 @@ int main(int argc, char *argv[])
   viewer.set_mesh(V, F);
 
   // Plot the field
-  viewer.add_edges (MF, MF+scale*polyVF[0],Eigen::RowVector3d(1,0,1));
-  viewer.add_edges (MF, MF+scale*polyVF[1],Eigen::RowVector3d(1,0,1));
+  viewer.add_edges (MF, MF+scale*X ,Eigen::RowVector3d(1,0,1));
+  viewer.add_edges (MF, MF+scale*X2,Eigen::RowVector3d(1,0,1));
   viewer.set_uv(UV,FUV);
   viewer.options.show_texture = true;