one-ring matchings are now computed in proper order on a boundary vertex

Former-commit-id: e211dc69005a7612cf17522a119a1051b968b41b

include/igl/polyvector_field_singularities_from_matchings.cpp

@@ -1,3 +1,4 @@
+#include <iostream>
 #include <igl/polyvector_field_singularities_from_matchings.h>
 #include <igl/is_border_vertex.h>
 #include <igl/vertex_triangle_adjacency.h>
@@ -22,7 +23,30 @@ void igl::polyvector_field_one_ring_matchings(const Eigen::PlainObjectBase<Deriv
   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];
+  //first, check if the vertex is on a boundary
+  //then there must be two faces that are on the boundary
+  //(other cases not supported)
+  
+  int fstart = -1;
+  int ind = 0;
+  for (int i =0; i<VF[vi].size(); ++i)
+  {
+    int fi = VF[vi][i];
+    for (int  j=0; j<3; ++j)
+      if (F(fi,(j+1)%3)==vi && TT(fi,j) == -1)
+      {
+        ind ++;
+        fstart = fi;
+        //        break;
+      }
+  }
+  if (ind >1 )
+  {
+    std::cerr<<"igl::polyvector_field_one_ring_matchings -- vertex "<<vi<< " is on an unusual boundary"<<std::endl;
+    exit(1);
+  }
+  if (fstart == -1)
+    fstart = VF[vi][0];
   int current_face = fstart;
   int i =0;
   fi[i] = current_face;
@@ -30,39 +54,47 @@ void igl::polyvector_field_one_ring_matchings(const Eigen::PlainObjectBase<Deriv
     mvi(i,j) = j;
   
   int next_face = -1;
-  while (next_face != fstart)
+  while (next_face != fstart && current_face!=-1)
   {
     // look for the vertex
     int j=-1;
     for (unsigned z=0; z<3; ++z)
       if (F(current_face,z) == vi)
+      {
         j=z;
+        break;
+      }
     assert(j!=-1);
-
+    
     next_face = TT(current_face, j);
     ++i;
-
-    // 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,k) = match_ab(current_edge,(mvi(i-1,k))%half_degree);
-    }
+    
+    if (next_face == -1)
+      mvi.row(i).setConstant(-1);
     else
     {
-      assert((E2F(current_edge,1) == current_face) &&
-             (E2F(current_edge,0) == next_face));
-      //look at match_ba
-        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);
+      // 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,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,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);
+      }
     }
     current_face = next_face;
     fi[i] = current_face;
@@ -77,17 +109,17 @@ IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
                                                                    const Eigen::PlainObjectBase<DerivedM> &match_ba,
                                                                    Eigen::PlainObjectBase<DerivedS> &singularities)
 {
-
+  
   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(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);
 }
 
@@ -104,20 +136,20 @@ IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
                                                                    const Eigen::PlainObjectBase<DerivedM> &match_ba,
                                                                    Eigen::PlainObjectBase<DerivedS> &singularities)
 {
-
+  
   int numV = V.rows();
-
+  
   std::vector<int> singularities_v;
   int half_degree = match_ab.cols();
-    for (int vi =0; vi<numV; ++vi)
-    {
-      ///check that is on border..
-      if (V_border[vi])
-        continue;
+  for (int vi =0; vi<numV; ++vi)
+  {
+    ///check that is on border..
+    if (V_border[vi])
+      continue;
     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);
-
+    
     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)
@@ -126,13 +158,13 @@ IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
       {
         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());
-
+  
   igl::list_to_matrix(singularities_v, singularities);
 }
 
@@ -175,7 +207,7 @@ IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
                                                                    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
@@ -183,13 +215,13 @@ IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
   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