olga's recent changes

Former-commit-id: 2aafc90507d5829a018bc0dded5ac158690a0c49

include/igl/hsv_to_rgb.cpp

@@ -68,4 +68,5 @@ template void igl::hsv_to_rgb<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::M
 template void igl::hsv_to_rgb<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
 template void igl::hsv_to_rgb<Eigen::Matrix<unsigned char, 64, 3, 1, 64, 3>, Eigen::Matrix<unsigned char, 64, 3, 1, 64, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<unsigned char, 64, 3, 1, 64, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned char, 64, 3, 1, 64, 3> >&);
 template void igl::hsv_to_rgb<Eigen::Matrix<float, 64, 3, 1, 64, 3>, Eigen::Matrix<float, 64, 3, 1, 64, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, 64, 3, 1, 64, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 64, 3, 1, 64, 3> >&);
+template void igl::hsv_to_rgb<double>(double const*, double*);
 #endif

include/igl/integrable_polyvector_fields.cpp

@@ -190,8 +190,7 @@ makeFieldCCW(Eigen::MatrixXd &sol3D)
   {
     //take all 4 vectors (including opposites) and pick two that are in ccw order
     all << sol3D.row(fi), -sol3D.row(fi);
-	Eigen::VectorXi inv_order_unused;
-    igl::sort_vectors_ccw(all, FN.row(fi).eval(), order, true, t, false, inv_order_unused);
+	  igl::sort_vectors_ccw(all, FN.row(fi).eval(), order, t);
     //if we are in a constrained face, we need to make sure that the first vector is always the same vector as in the constraints
     if(is_constrained_face[fi])
     {

include/igl/is_vertex_manifold.cpp

@@ -97,4 +97,5 @@ IGL_INLINE bool igl::is_vertex_manifold(
 
 #ifdef IGL_STATIC_LIBRARY
 template bool igl::is_vertex_manifold<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template bool igl::is_vertex_manifold<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/matlab_format.cpp

@@ -54,7 +54,7 @@ igl::matlab_format(
   if(!name.empty())
   {
     prefix = name + "IJV = ";
-    suffix = "\n"+name + " = sparse("+name+"IJV(:,1),"+name+"IJV(:,2),"+name+"IJV(:,3));";
+    suffix = "\n"+name + " = sparse("+name+"IJV(:,1),"+name+"IJV(:,2),"+name+"IJV(:,3),"+std::to_string(S.rows())+","+std::to_string(S.cols())+" );";
   }
   return STR(""<<
     SIJV.format(Eigen::IOFormat(

include/igl/polyvector_field_matchings.cpp

@@ -13,21 +13,36 @@ IGL_INLINE void igl::polyvector_field_matching(
                                                const Eigen::PlainObjectBase<DerivedV>& e,
                                                bool match_with_curl,
                                                Eigen::PlainObjectBase<DerivedM>& mab,
-                                               Eigen::PlainObjectBase<DerivedM>& mba)
+                                               Eigen::PlainObjectBase<DerivedM>& mba,
+                                               bool is_symmetric)
 {
   // make sure the matching preserve ccw order of the vectors across the edge
   // 1) order vectors in a, ccw  e.g. (0,1,2,3)_a not ccw --> (0,3,2,1)_a ccw
   // 2) order vectors in b, ccw  e.g. (0,1,2,3)_b not ccw --> (0,2,1,3)_b ccw
   // 3) the vectors in b that match the ordered vectors in a (in this case  (0,3,2,1)_a ) must be a circular shift of the ccw ordered vectors in b  - so we need to explicitely check only these matchings to find the best ccw one, there are N of them
   int hN = _ua.cols()/3;
-  int N = 2*hN;
-  Eigen::Matrix<typename DerivedS::Scalar,1,Eigen::Dynamic> ua (1,N*3); ua <<_ua, -_ua;
-  Eigen::Matrix<typename DerivedS::Scalar,1,Eigen::Dynamic> ub (1,N*3); ub <<_ub, -_ub;
+  int N;
+  if (is_symmetric)
+    N = 2*hN;
+  else
+    N = hN;
+
+  Eigen::Matrix<typename DerivedS::Scalar,1,Eigen::Dynamic> ua (1,N*3);
+  Eigen::Matrix<typename DerivedS::Scalar,1,Eigen::Dynamic> ub (1,N*3);
+  if (is_symmetric)
+  {
+    ua <<_ua, -_ua;
+    ub <<_ub, -_ub;
+  }
+  else
+  {
+    ua =_ua;
+    ub =_ub;
+  }
+
   Eigen::Matrix<typename DerivedM::Scalar,Eigen::Dynamic,1> order_a, order_b;
-  Eigen::Matrix<typename DerivedS::Scalar, 1, Eigen::Dynamic> sorted_unused;
-  Eigen::VectorXi inv_order_unused;
-  igl::sort_vectors_ccw(ua, na, order_a,false,sorted_unused,false,inv_order_unused);
-  igl::sort_vectors_ccw(ub, nb, order_b,false,sorted_unused,false,inv_order_unused);
+  igl::sort_vectors_ccw(ua, na, order_a);
+  igl::sort_vectors_ccw(ub, nb, order_b);
 
   //checking all possible circshifts of order_b as matches for order_a
   Eigen::Matrix<typename DerivedM::Scalar,Eigen::Dynamic,Eigen::Dynamic> all_matches(N,N);
@@ -71,14 +86,16 @@ IGL_INLINE void igl::polyvector_field_matching(
   for (int i=0; i< N; ++i)
     mba[mab[i]] = i;
 
+  if (is_symmetric)
+  {
   mab = mab.head(hN);
   mba = mba.head(hN);
-
+}
 }
 
 
-template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedM, typename DerivedC>
-IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
+template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedM>
+IGL_INLINE void igl::polyvector_field_matchings(
                                                                      const Eigen::PlainObjectBase<DerivedS>& sol3D,
                                                                      const Eigen::PlainObjectBase<DerivedV>&V,
                                                                      const Eigen::PlainObjectBase<DerivedF>&F,
@@ -86,9 +103,9 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
                                                                      const Eigen::PlainObjectBase<DerivedV>& FN,
                                                                      const Eigen::MatrixXi &E2F,
                                                                      bool match_with_curl,
+                                                bool is_symmetric,
                                                                      Eigen::PlainObjectBase<DerivedM>& match_ab,
-                                                                     Eigen::PlainObjectBase<DerivedM>& match_ba,
-                                                                     Eigen::PlainObjectBase<DerivedC>& curl)
+                                                Eigen::PlainObjectBase<DerivedM>& match_ba)
 {
   int numEdges = E.rows();
   int half_degree = sol3D.cols()/3;
@@ -101,11 +118,9 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
       isBorderEdge[i] = 1;
   }
 
-  curl.setZero(numEdges,1);
   match_ab.setZero(numEdges, half_degree);
   match_ba.setZero(numEdges, half_degree);
 
-  typename DerivedC::Scalar meanCurl = 0;
   for (int ei=0; ei<numEdges; ++ei)
   {
     if (isBorderEdge[ei])
@@ -124,18 +139,69 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
                                    ce,
                                    match_with_curl,
                                    mab,
-                                   mba);
+                                   mba,
+                                   is_symmetric);
 
     match_ab.row(ei) = mab;
     match_ba.row(ei) = mba;
+  }
+}
+
+
+template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedM, typename DerivedC>
+IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
+                                                                     const Eigen::PlainObjectBase<DerivedS>& sol3D,
+                                                                     const Eigen::PlainObjectBase<DerivedV>&V,
+                                                                     const Eigen::PlainObjectBase<DerivedF>&F,
+                                                                     const Eigen::PlainObjectBase<DerivedE>&E,
+                                                                     const Eigen::PlainObjectBase<DerivedV>& FN,
+                                                                     const Eigen::MatrixXi &E2F,
+                                                                     bool match_with_curl,
+                                                                     bool is_symmetric,
+                                                                     Eigen::PlainObjectBase<DerivedM>& match_ab,
+                                                                     Eigen::PlainObjectBase<DerivedM>& match_ba,
+                                                                     Eigen::PlainObjectBase<DerivedC>& curl)
+{
+  int numEdges = E.rows();
+  int half_degree = sol3D.cols()/3;
+
+  Eigen::VectorXi isBorderEdge;
+  isBorderEdge.setZero(numEdges,1);
+  for(unsigned i=0; i<numEdges; ++i)
+  {
+    if ((E2F(i,0) == -1) || ((E2F(i,1) == -1)))
+      isBorderEdge[i] = 1;
+  }
+
+  igl::polyvector_field_matchings(sol3D, V, F, E, FN, E2F, match_with_curl, is_symmetric, match_ab, match_ba);
+  curl.setZero(numEdges,1);
+  typename DerivedC::Scalar meanCurl = 0;
+  for (int ei=0; ei<numEdges; ++ei)
+  {
+    if (isBorderEdge[ei])
+      continue;
+    // the two faces of the flap
+    int a = E2F(ei,0);
+    int b = E2F(ei,1);
+
+    const Eigen::Matrix<typename DerivedM::Scalar, 1, Eigen::Dynamic> &mab = match_ab.row(ei);
+    const Eigen::Matrix<typename DerivedM::Scalar, 1, Eigen::Dynamic> &mba = match_ba.row(ei);
+
     Eigen::Matrix<typename DerivedS::Scalar, 1, Eigen::Dynamic> matched;
     matched.resize(1, 3*half_degree);
     for (int i = 0; i<half_degree; ++i)
     {
-      int sign = (mab[i]<half_degree)?1:-1;
-      matched.segment(i*3, 3) = sign*sol3D.row(b).segment((mab[i]%half_degree)*3, 3);
+      int sign = 1;
+      int m = mab[i];
+      if (is_symmetric)
+      {
+        sign = (mab[i]<half_degree)?1:-1;
+        m = m%half_degree;
+    }
+      matched.segment(i*3, 3) = sign*sol3D.row(b).segment(m*3, 3);
     }
 
+    Eigen::Matrix<typename DerivedV::Scalar, 1, Eigen::Dynamic> ce = (V.row(E(ei,1))-V.row(E(ei,0))).normalized().template cast<typename DerivedV::Scalar>();
     typename DerivedC::Scalar avgCurl = 0;
     for (int i = 0; i<half_degree; ++i)
       avgCurl += fabs(sol3D.row(a).segment(i*3, 3).dot(ce) - matched.segment(i*3, 3).dot(ce));
@@ -153,14 +219,13 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
   return meanCurl;
 }
 
-
-
 template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedM, typename DerivedC>
 IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
                                                                      const Eigen::PlainObjectBase<DerivedS>& sol3D,
                                                                      const Eigen::PlainObjectBase<DerivedV>&V,
                                                                      const Eigen::PlainObjectBase<DerivedF>&F,
                                                                      bool match_with_curl,
+                                                                     bool is_symmetric,
                                                                      Eigen::PlainObjectBase<DerivedM>& match_ab,
                                                                      Eigen::PlainObjectBase<DerivedM>& match_ba,
                                                                      Eigen::PlainObjectBase<DerivedC>& curl)
@@ -172,12 +237,31 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
   DerivedV FN;
   igl::per_face_normals(V,F,FN);
 
-  return igl::polyvector_field_matchings(sol3D, V, F, E, FN, E2F, match_with_curl, match_ab, match_ba, curl);
+  return igl::polyvector_field_matchings(sol3D, V, F, E, FN, E2F, match_with_curl, is_symmetric, match_ab, match_ba, curl);
+}
+
+template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedM>
+IGL_INLINE void igl::polyvector_field_matchings(
+                                                const Eigen::PlainObjectBase<DerivedS>& sol3D,
+                                                const Eigen::PlainObjectBase<DerivedV>&V,
+                                                const Eigen::PlainObjectBase<DerivedF>&F,
+                                                bool match_with_curl,
+                                                bool is_symmetric,
+                                                Eigen::PlainObjectBase<DerivedM>& match_ab,
+                                                Eigen::PlainObjectBase<DerivedM>& match_ba)
+{
+  Eigen::MatrixXi E, E2F, F2E;
+  igl::edge_topology(V,F,E,F2E,E2F);
+
+  Eigen::PlainObjectBase<DerivedV> FN;
+  igl::per_face_normals(V,F,FN);
+
+  igl::polyvector_field_matchings(sol3D, V, F, E, FN, E2F, match_with_curl, is_symmetric, match_ab, match_ba);
 }
 
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template Eigen::Matrix<float, -1, 1, 0, -1, 1>::Scalar igl::polyvector_field_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
-template Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar igl::polyvector_field_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::polyvector_field_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar igl::polyvector_field_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/polyvector_field_matchings.h

@@ -25,6 +25,10 @@ namespace igl {
   //   e                1 by 3, the vector corresponding to the shared edge between a and b
   //   match_with_curl  boolean flag, determines whether a curl or a smoothness matching will
   //                    be computed
+  //   is_symmetric     boolean flag, determines whether the input vector set field is symmetric(
+  //                    =consists of pairwise collinear vectors in each set, in which case only one
+  //                    of the vectors in the pair is stored) or not, i.e. the set contains all the vectors
+  // )
   // Outputs:
   //   mab              1 by N row vector, describing the matching a->b (i.e. vector #i of the
   //                    vector set in a is matched to vector #mab[i] in b)
@@ -40,7 +44,8 @@ namespace igl {
                                             const Eigen::PlainObjectBase<DerivedV>& e,
                                             bool match_with_curl,
                                             Eigen::PlainObjectBase<DerivedM>& mab,
-                                            Eigen::PlainObjectBase<DerivedM>& mba);
+                                            Eigen::PlainObjectBase<DerivedM>& mba,
+                                            bool is_symmetric);
 
 
   // Given a mesh and a vector set field consisting of unordered N-vector sets defined
@@ -58,6 +63,9 @@ namespace igl {
   //                    via igl::edge_topology)
   //   match_with_curl  boolean flag, determines whether curl or smoothness matchings will
   //                    be computed
+  //   is_symmetric     boolean flag, determines whether the input vector set field is symmetric(
+  //                    =consists of pairwise collinear vectors in each set, in which case only one
+  //                    of the vectors in the pair is stored) or not, i.e. the set contains all the vectors
   // Outputs:
   //   match_ab         #E by N matrix, describing for each edge the matching a->b, where a
   //                    and b are the faces adjacent to the edge (i.e. vector #i of
@@ -78,6 +86,7 @@ namespace igl {
                                                                   const Eigen::PlainObjectBase<DerivedV>& FN,
                                                                   const Eigen::MatrixXi &E2F,
                                                                   bool match_with_curl,
+                                                                  bool is_symmetric,
                                                                   Eigen::PlainObjectBase<DerivedM>& match_ab,
                                                                   Eigen::PlainObjectBase<DerivedM>& match_ba,
                                                                   Eigen::PlainObjectBase<DerivedC>& curl);
@@ -90,10 +99,34 @@ namespace igl {
                                                                   const Eigen::PlainObjectBase<DerivedV>&V,
                                                                   const Eigen::PlainObjectBase<DerivedF>&F,
                                                                   bool match_with_curl,
+                                                                  bool is_symmetric,
                                                                   Eigen::PlainObjectBase<DerivedM>& match_ab,
                                                                   Eigen::PlainObjectBase<DerivedM>& match_ba,
                                                                   Eigen::PlainObjectBase<DerivedC>& curl);
 
+  //Wrappers with no curl output
+  template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedM>
+  IGL_INLINE void polyvector_field_matchings(
+                                             const Eigen::PlainObjectBase<DerivedS>& sol3D,
+                                             const Eigen::PlainObjectBase<DerivedV>&V,
+                                             const Eigen::PlainObjectBase<DerivedF>&F,
+                                             bool match_with_curl,
+                                             bool is_symmetric,
+                                             Eigen::PlainObjectBase<DerivedM>& match_ab,
+                                             Eigen::PlainObjectBase<DerivedM>& match_ba);
+  template <typename DerivedS, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedM>
+  IGL_INLINE void polyvector_field_matchings(
+                                             const Eigen::PlainObjectBase<DerivedS>& sol3D,
+                                             const Eigen::PlainObjectBase<DerivedV>&V,
+                                             const Eigen::PlainObjectBase<DerivedF>&F,
+                                             const Eigen::PlainObjectBase<DerivedE>&E,
+                                             const Eigen::PlainObjectBase<DerivedV>& FN,
+                                             const Eigen::MatrixXi &E2F,
+                                             bool match_with_curl,
+                                             bool is_symmetric,
+                                             Eigen::PlainObjectBase<DerivedM>& match_ab,
+                                             Eigen::PlainObjectBase<DerivedM>& match_ba);
+  
 };
 
 

include/igl/polyvector_field_poisson_reconstruction.cpp

@@ -9,14 +9,12 @@
 
 #include <Eigen/Sparse>
 
-template <typename DerivedV, typename DerivedF, typename DerivedSF, typename DerivedS, typename DerivedE>
-IGL_INLINE double igl::polyvector_field_poisson_reconstruction(
+template <typename DerivedV, typename DerivedF, typename DerivedSF, typename DerivedS>
+IGL_INLINE void igl::polyvector_field_poisson_reconstruction(
   const Eigen::PlainObjectBase<DerivedV> &Vcut,
   const Eigen::PlainObjectBase<DerivedF> &Fcut,
   const Eigen::PlainObjectBase<DerivedS> &sol3D_combed,
-  Eigen::PlainObjectBase<DerivedSF> &scalars,
-  Eigen::PlainObjectBase<DerivedS> &sol3D_recon,
-  Eigen::PlainObjectBase<DerivedE> &max_error )
+                                                               Eigen::PlainObjectBase<DerivedSF> &scalars)
   {
     Eigen::SparseMatrix<typename DerivedV::Scalar> gradMatrix;
     igl::grad(Vcut, Fcut, gradMatrix);
@@ -34,8 +32,6 @@ IGL_INLINE double igl::polyvector_field_poisson_reconstruction(
 
     int half_degree = sol3D_combed.cols()/3;
 
-    sol3D_recon.setZero(sol3D_combed.rows(),sol3D_combed.cols());
-
     int numF = Fcut.rows();
     scalars.setZero(Vcut.rows(),half_degree);
 
@@ -70,8 +66,28 @@ IGL_INLINE double igl::polyvector_field_poisson_reconstruction(
       Eigen::VectorXi index = i*Eigen::VectorXi::Ones(Iu.rows(),1);
       igl::slice_into(yu, Iu, index, scalars);scalars(Ik[0],i)=xk[0];
     }
+}
+
+template <typename DerivedV, typename DerivedF, typename DerivedSF, typename DerivedS, typename DerivedE>
+IGL_INLINE double igl::polyvector_field_poisson_reconstruction(
+                                                               const Eigen::PlainObjectBase<DerivedV> &Vcut,
+                                                               const Eigen::PlainObjectBase<DerivedF> &Fcut,
+                                                               const Eigen::PlainObjectBase<DerivedS> &sol3D_combed,
+                                                               Eigen::PlainObjectBase<DerivedSF> &scalars,
+                                                               Eigen::PlainObjectBase<DerivedS> &sol3D_recon,
+                                                               Eigen::PlainObjectBase<DerivedE> &max_error )
+{
+  
+  igl::polyvector_field_poisson_reconstruction(Vcut, Fcut, sol3D_combed, scalars);
+
+  Eigen::SparseMatrix<typename DerivedV::Scalar> gradMatrix;
+  igl::grad(Vcut, Fcut, gradMatrix);
+  int numF = Fcut.rows();
+  int half_degree = sol3D_combed.cols()/3;
 
     //    evaluate gradient of found scalar function
+  sol3D_recon.setZero(sol3D_combed.rows(),sol3D_combed.cols());
+  
     for (int i =0; i<half_degree; ++i)
     {
       Eigen::VectorXd vec_poisson = gradMatrix*scalars.col(i);
@@ -91,7 +107,9 @@ IGL_INLINE double igl::polyvector_field_poisson_reconstruction(
     return max_error.mean();
   }
 
+
   #ifdef IGL_STATIC_LIBRARY
   // Explicit template specialization
   template double igl::polyvector_field_poisson_reconstruction<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
+template double igl::polyvector_field_poisson_reconstruction<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
   #endif

include/igl/polyvector_field_poisson_reconstruction.h

@@ -41,6 +41,13 @@ namespace igl {
                                                             Eigen::PlainObjectBase<DerivedSF> &scalars,
                                                             Eigen::PlainObjectBase<DerivedS> &sol3D_recon,
                                                             Eigen::PlainObjectBase<DerivedE> &max_error );
+  //Wrappers with less output
+  template <typename DerivedV, typename DerivedF, typename DerivedSF, typename DerivedS>
+  IGL_INLINE void polyvector_field_poisson_reconstruction(
+                                                            const Eigen::PlainObjectBase<DerivedV> &Vcut,
+                                                            const Eigen::PlainObjectBase<DerivedF> &Fcut,
+                                                            const Eigen::PlainObjectBase<DerivedS> &sol3D_combed,
+                                                            Eigen::PlainObjectBase<DerivedSF> &scalars);
 
 
 };

include/igl/polyvector_field_singularities_from_matchings.cpp

@@ -15,19 +15,20 @@ void igl::polyvector_field_one_ring_matchings(const Eigen::PlainObjectBase<Deriv
                                               const Eigen::PlainObjectBase<DerivedM> &match_ab,
                                               const Eigen::PlainObjectBase<DerivedM> &match_ba,
                                               const int vi,
-                                              const int vector_to_match,
-                                              Eigen::VectorXi &mvi,
+                                              Eigen::MatrixXi &mvi,
                                               Eigen::VectorXi &fi)
 {
   int half_degree = match_ab.cols();
-  mvi.resize(VF[vi].size()+1,1);
+  mvi.resize(VF[vi].size()+1,half_degree);
   fi.resize(VF[vi].size()+1,1);
   //start from one face
-  const int &fstart = VF[vi][0];
+  const int fstart = VF[vi][0];
   int current_face = fstart;
   int i =0;
-  mvi[i] = vector_to_match;
   fi[i] = current_face;
+  for (int j=0; j<half_degree; ++j)
+    mvi(i,j) = j;
+  
   int next_face = -1;
   while (next_face != fstart)
   {
@@ -44,23 +45,25 @@ void igl::polyvector_field_one_ring_matchings(const Eigen::PlainObjectBase<Deriv
     // look at the edge between the two faces
     const int &current_edge = F2E(current_face,j);
 
+    for (int k=0; k<half_degree; ++k)
+    {
     // check its orientation to determine whether match_ab or match_ba should be used
     if ((E2F(current_edge,0) == current_face) &&
         (E2F(current_edge,1) == next_face) )
     {
       //look at match_ab
-      mvi[i] = match_ab(current_edge,(mvi[i-1])%half_degree);
+        mvi(i,k) = match_ab(current_edge,(mvi(i-1,k))%half_degree);
     }
     else
     {
       assert((E2F(current_edge,1) == current_face) &&
              (E2F(current_edge,0) == next_face));
       //look at match_ba
-      mvi[i] = match_ba(current_edge,(mvi[i-1])%half_degree);
+        mvi(i,k) = match_ba(current_edge,(mvi(i-1,k))%half_degree);
+    }
+      if (mvi(i-1,k)>=half_degree)
+        mvi(i,k) = (mvi(i,k)+half_degree)%(2*half_degree);
     }
-    if (mvi[i-1]>=half_degree)
-      mvi[i] = (mvi[i]+half_degree)%(2*half_degree);
-
     current_face = next_face;
     fi[i] = current_face;
   }
@@ -106,22 +109,26 @@ IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
 
   std::vector<int> singularities_v;
   int half_degree = match_ab.cols();
-  //pick one of the vectors to check for singularities
-  for (int vector_to_match = 0; vector_to_match < half_degree; ++vector_to_match)
-  {
     for (int vi =0; vi<numV; ++vi)
     {
       ///check that is on border..
       if (V_border[vi])
         continue;
-      Eigen::VectorXi mvi,fi;
-
-      igl::polyvector_field_one_ring_matchings(V, F, VF, E2F, F2E, TT, match_ab, match_ba, vi, vector_to_match, mvi, fi);
+    Eigen::VectorXi fi;
+    Eigen::MatrixXi mvi;
+    igl::polyvector_field_one_ring_matchings(V, F, VF, E2F, F2E, TT, match_ab, match_ba, vi, mvi, fi);
 
-      if(mvi.tail(1) != mvi.head(1))
+    int num = fi.size();
+    //pick one of the vectors to check for singularities
+    for (int vector_to_match = 0; vector_to_match < half_degree; ++vector_to_match)
+    {
+      if(mvi(num-1,vector_to_match) != mvi(0,vector_to_match))
+      {
         singularities_v.push_back(vi);
+        break;
     }
   }
+  }
   std::sort(singularities_v.begin(), singularities_v.end());
   auto last = std::unique(singularities_v.begin(), singularities_v.end());
   singularities_v.erase(last, singularities_v.end());
@@ -129,8 +136,67 @@ IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
   igl::list_to_matrix(singularities_v, singularities);
 }
 
+
+template <typename DerivedV, typename DerivedF, typename DerivedM, typename DerivedS>
+IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
+                                                                   const Eigen::PlainObjectBase<DerivedV> &V,
+                                                                   const Eigen::PlainObjectBase<DerivedF> &F,
+                                                                   const Eigen::PlainObjectBase<DerivedM> &match_ab,
+                                                                   const Eigen::PlainObjectBase<DerivedM> &match_ba,
+                                                                   Eigen::PlainObjectBase<DerivedS> &singularities,
+                                                                   Eigen::PlainObjectBase<DerivedS> &singularity_indices)
+{
+  
+  std::vector<bool> V_border = igl::is_border_vertex(V,F);
+  std::vector<std::vector<int> > VF, VFi;
+  igl::vertex_triangle_adjacency(V,F,VF,VFi);
+  
+  Eigen::MatrixXi TT, TTi;
+  igl::triangle_triangle_adjacency(V,F,TT,TTi);
+  
+  Eigen::MatrixXi E, E2F, F2E;
+  igl::edge_topology(V,F,E,F2E,E2F);
+  
+  igl::polyvector_field_singularities_from_matchings(V, F, V_border, VF, TT, E2F, F2E, match_ab, match_ba, singularities, singularity_indices);
+}
+
+template <typename DerivedV, typename DerivedF, typename DerivedM, typename VFType, typename DerivedS>
+IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
+                                                                   const Eigen::PlainObjectBase<DerivedV> &V,
+                                                                   const Eigen::PlainObjectBase<DerivedF> &F,
+                                                                   const std::vector<bool> &V_border,
+                                                                   const std::vector<std::vector<VFType> > &VF,
+                                                                   const Eigen::MatrixXi &TT,
+                                                                   const Eigen::MatrixXi &E2F,
+                                                                   const Eigen::MatrixXi &F2E,
+                                                                   const Eigen::PlainObjectBase<DerivedM> &match_ab,
+                                                                   const Eigen::PlainObjectBase<DerivedM> &match_ba,
+                                                                   Eigen::PlainObjectBase<DerivedS> &singularities,
+                                                                   Eigen::PlainObjectBase<DerivedS> &singularity_indices)
+{
+  igl::polyvector_field_singularities_from_matchings(V, F, V_border, VF, TT, E2F, F2E, match_ab, match_ba, singularities);
+
+  singularity_indices.setZero(singularities.size(), 1);
+  
+  //get index from first vector only
+  int vector_to_match = 0;
+  for (int i =0; i<singularities.size(); ++i)
+  {
+    int vi = singularities[i];
+
+    Eigen::VectorXi mvi,fi;
+    igl::polyvector_field_one_ring_matchings(V, F, VF, E2F, F2E, TT, match_ab, match_ba, vi, vector_to_match, mvi, fi);
+    
+    singularity_indices[i] = (mvi.tail(1)[0] - vector_to_match);
+  }
+
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::polyvector_field_singularities_from_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<bool, std::allocator<bool> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::polyvector_field_singularities_from_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::polyvector_field_one_ring_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::Matrix<int, -1, -1, 0, -1, -1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
+template void igl::polyvector_field_singularities_from_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::__1::vector<bool, std::__1::allocator<bool> > const&, std::__1::vector<std::__1::vector<int, std::__1::allocator<int> >, std::__1::allocator<std::__1::vector<int, std::__1::allocator<int> > > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+//template void igl::polyvector_field_singularities_from_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<bool, std::allocator<bool> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+//template void igl::polyvector_field_singularities_from_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+//template void igl::polyvector_field_singularities_from_matchings<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/polyvector_field_singularities_from_matchings.h

@@ -17,7 +17,7 @@ namespace igl {
   
   // We are given a polyvector field on a mesh (with N vectors per face), and matchings between the vector sets
   // across all non-boundary edges.  The function computes, for the one-ring
-  // neighborhood of a given vertex, and for the selected vector of the vector set,
+  // neighborhood of a given vertex, and for each of the vectors of the vector set,
   // the sequence of the vectors that match to it around the one-ring. If the vector that
   // we land on by following the matchings is not the original vector that we started from,
   // the vertex is a singularity.
@@ -40,7 +40,6 @@ namespace igl {
   //                    and b are the faces adjacent to the edge (i.e. vector #mba[i] of
   //                    the vector set in a is matched to vector #i in b)
   //   vi               the selected one ring
-  //   vector_id        the selected vector of the polyvector
   //
   // Output:
   //   mvi              #numOneRingFaces by 1 list of the indices of the sequentially matching vectors
@@ -58,8 +57,7 @@ namespace igl {
                                            const Eigen::PlainObjectBase<DerivedM> &match_ab,
                                            const Eigen::PlainObjectBase<DerivedM> &match_ba,
                                            const int vi,
-                                           const int vector_id,
-                                           Eigen::VectorXi &mvi,
+                                           Eigen::MatrixXi &mvi,
                                            Eigen::VectorXi &fi);
   
   
@@ -114,6 +112,31 @@ namespace igl {
                                                                 Eigen::PlainObjectBase<DerivedS> &singularities);
   
   
+  // Same pair as above but also returns singularity indices
+  template <typename DerivedV, typename DerivedF, typename DerivedM, typename DerivedS>
+  IGL_INLINE void polyvector_field_singularities_from_matchings(
+                                                                const Eigen::PlainObjectBase<DerivedV> &V,
+                                                                const Eigen::PlainObjectBase<DerivedF> &F,
+                                                                const Eigen::PlainObjectBase<DerivedM> &match_ab,
+                                                                const Eigen::PlainObjectBase<DerivedM> &match_ba,
+                                                                Eigen::PlainObjectBase<DerivedS> &singularities,
+                                                                Eigen::PlainObjectBase<DerivedS> &singularity_indices);
+  
+  template <typename DerivedV, typename DerivedF, typename DerivedM, typename VFType, typename DerivedS>
+  IGL_INLINE void polyvector_field_singularities_from_matchings(
+                                                                const Eigen::PlainObjectBase<DerivedV> &V,
+                                                                const Eigen::PlainObjectBase<DerivedF> &F,
+                                                                const std::vector<bool> &V_border,
+                                                                const std::vector<std::vector<VFType> > &VF,
+                                                                const Eigen::MatrixXi &TT,
+                                                                const Eigen::MatrixXi &E2F,
+                                                                const Eigen::MatrixXi &F2E,
+                                                                const Eigen::PlainObjectBase<DerivedM> &match_ab,
+                                                                const Eigen::PlainObjectBase<DerivedM> &match_ba,
+                                                                Eigen::PlainObjectBase<DerivedS> &singularities,
+                                                                Eigen::PlainObjectBase<DerivedS> &singularity_indices);
+
+  
   
 };
 

include/igl/sort_vectors_ccw.cpp

@@ -6,14 +6,9 @@ template <typename DerivedS, typename DerivedI>
 IGL_INLINE void igl::sort_vectors_ccw(
   const Eigen::PlainObjectBase<DerivedS>& P,
   const Eigen::PlainObjectBase<DerivedS>& N,
-  Eigen::PlainObjectBase<DerivedI> &order,
-  const bool do_sorted,
-  Eigen::PlainObjectBase<DerivedS> &sorted,
-  const bool do_inv_order,
-  Eigen::PlainObjectBase<DerivedI> &inv_order)
+  Eigen::PlainObjectBase<DerivedI> &order)
 {
   int half_degree = P.cols()/3;
-
   //local frame
   Eigen::Matrix<typename DerivedS::Scalar,1,3> e1 = P.head(3).normalized();
   Eigen::Matrix<typename DerivedS::Scalar,1,3> e3 = N.normalized();
@@ -39,14 +34,31 @@ IGL_INLINE void igl::sort_vectors_ccw(
       order[i] = order[i+1];
     order(half_degree-1) = temp;
   }
-  if (do_sorted)
+}
+
+template <typename DerivedS, typename DerivedI>
+IGL_INLINE void igl::sort_vectors_ccw(
+  const Eigen::PlainObjectBase<DerivedS>& P,
+  const Eigen::PlainObjectBase<DerivedS>& N,
+  Eigen::PlainObjectBase<DerivedI> &order,
+  Eigen::PlainObjectBase<DerivedS> &sorted)
   {
+  int half_degree = P.cols()/3;
+  igl::sort_vectors_ccw(P,N,order);
     sorted.resize(1,half_degree*3);
     for (int i=0; i<half_degree; ++i)
       sorted.segment(i*3,3) = P.segment(order[i]*3,3);
   }
-  if (do_inv_order)
+
+template <typename DerivedS, typename DerivedI>
+IGL_INLINE void igl::sort_vectors_ccw(
+  const Eigen::PlainObjectBase<DerivedS>& P,
+  const Eigen::PlainObjectBase<DerivedS>& N,
+  Eigen::PlainObjectBase<DerivedI> &order,
+  Eigen::PlainObjectBase<DerivedI> &inv_order)
   {
+  int half_degree = P.cols()/3;
+  igl::sort_vectors_ccw(P,N,order);
     inv_order.resize(half_degree,1);
     for (int i=0; i<half_degree; ++i)
     {
@@ -60,9 +72,25 @@ IGL_INLINE void igl::sort_vectors_ccw(
     assert(inv_order[0] == 0);
   }
 
+template <typename DerivedS, typename DerivedI>
+IGL_INLINE void igl::sort_vectors_ccw(
+  const Eigen::PlainObjectBase<DerivedS>& P,
+  const Eigen::PlainObjectBase<DerivedS>& N,
+  Eigen::PlainObjectBase<DerivedI> &order,
+  Eigen::PlainObjectBase<DerivedS> &sorted,
+  Eigen::PlainObjectBase<DerivedI> &inv_order)
+{
+  int half_degree = P.cols()/3;
+
+  igl::sort_vectors_ccw(P,N,order,inv_order);
+
+  sorted.resize(1,half_degree*3);
+  for (int i=0; i<half_degree; ++i)
+    sorted.segment(i*3,3) = P.segment(order[i]*3,3);
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::sort_vectors_ccw<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> >&, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::sort_vectors_ccw<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::sort_vectors_ccw<Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> >&);
 #endif

include/igl/sort_vectors_ccw.h

@@ -21,8 +21,6 @@ namespace igl {
   // Inputs:
   //   P               1 by 3N row vector of the vectors to be sorted, stacked horizontally
   //   N               #1 by 3 normal of the plane where the vectors lie
-  //   do_sorted       boolean flag, determines whether to return the sorted vector set
-  //   do_inv_order    boolean flag, determines whether to return the "inverse" set of indices
   // Output:
   //   order           N by 1 order of the vectors (indices of the unordered vectors into
   //                   the ordered vector set)
@@ -35,10 +33,30 @@ namespace igl {
                                    const Eigen::PlainObjectBase<DerivedS>& P,
                                    const Eigen::PlainObjectBase<DerivedS>& N,
                                    Eigen::PlainObjectBase<DerivedI> &order,
-                                   const bool do_sorted,
                                    Eigen::PlainObjectBase<DerivedS> &sorted,
-                                   const bool do_inv_order,
                                    Eigen::PlainObjectBase<DerivedI> &inv_order);
+
+   template <typename DerivedS, typename DerivedI>
+   IGL_INLINE void sort_vectors_ccw(
+                                    const Eigen::PlainObjectBase<DerivedS>& P,
+                                    const Eigen::PlainObjectBase<DerivedS>& N,
+                                    Eigen::PlainObjectBase<DerivedI> &order,
+                                    Eigen::PlainObjectBase<DerivedS> &sorted);
+
+    template <typename DerivedS, typename DerivedI>
+    IGL_INLINE void sort_vectors_ccw(
+                                     const Eigen::PlainObjectBase<DerivedS>& P,
+                                     const Eigen::PlainObjectBase<DerivedS>& N,
+                                     Eigen::PlainObjectBase<DerivedI> &order,
+                                     Eigen::PlainObjectBase<DerivedI> &inv_order);
+
+
+     template <typename DerivedS, typename DerivedI>
+     IGL_INLINE void sort_vectors_ccw(
+                                      const Eigen::PlainObjectBase<DerivedS>& P,
+                                      const Eigen::PlainObjectBase<DerivedS>& N,
+                                      Eigen::PlainObjectBase<DerivedI> &order);
+
 };