polygon mesh support

include/igl/per_corner_normals.cpp

@@ -101,6 +101,97 @@ IGL_INLINE void igl::per_corner_normals(
     }
   }
 }
+
+#include "doublearea.h"
+
+template <
+  typename DerivedV, 
+  typename DerivedI, 
+  typename DerivedC, 
+  typename DerivedN,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedJ,
+  typename DerivedNN>
+IGL_INLINE void igl::per_corner_normals(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedC> & C,
+  const typename DerivedV::Scalar corner_threshold,
+  Eigen::PlainObjectBase<DerivedN>  & N,
+  Eigen::PlainObjectBase<DerivedVV> & VV,
+  Eigen::PlainObjectBase<DerivedFF> & FF,
+  Eigen::PlainObjectBase<DerivedJ>  & J,
+  Eigen::PlainObjectBase<DerivedNN> & NN)
+{
+  const Eigen::Index m = C.size()-1;
+  typedef Eigen::Index Index;
+  Eigen::MatrixXd FN;
+  per_face_normals(V,I,C,FN,VV,FF,J);
+  typedef typename DerivedN::Scalar Scalar;
+  Eigen::Matrix<Scalar,Eigen::Dynamic,1> AA;
+  doublearea(VV,FF,AA);
+  // VF[i](j) = p means p is the jth face incident on vertex i
+  // to-do micro-optimization to avoid vector<vector>
+  std::vector<std::vector<Eigen::Index>> VF(V.rows());
+  for(Eigen::Index p = 0;p<m;p++)
+  {
+    // number of faces/vertices in this simple polygon
+    const Index np = C(p+1)-C(p);
+    for(Eigen::Index i = 0;i<np;i++)
+    {
+      VF[I(C(p)+i)].push_back(p);
+    }
+  }
+
+  N.resize(I.rows(),3);
+  for(Eigen::Index p = 0;p<m;p++)
+  {
+    // number of faces/vertices in this simple polygon
+    const Index np = C(p+1)-C(p);
+    Eigen::Matrix<Scalar,3,1> fn = FN.row(p);
+    for(Eigen::Index i = 0;i<np;i++)
+    {
+      N.row(C(p)+i).setZero();
+      // Loop over faces sharing this vertex
+      for(const auto & n : VF[I(C(p)+i)])
+      {
+        Eigen::Matrix<Scalar,3,1> ifn = AA(n) * FN.row(n);
+        // dot product between face's normal and other face's normal
+        Scalar dp = fn.dot(ifn);
+        if(dp > cos(corner_threshold*igl::PI/180))
+        {
+          // add to running sum
+          N.row(C(p)+i) += ifn;
+        }
+      }
+      N.row(C(p)+i).normalize();
+    }
+  }
+
+  // Relies on order of FF 
+  NN.resize(FF.rows()*3,3);
+  {
+    Eigen::Index k = 0;
+    for(Eigen::Index p = 0;p<m;p++)
+    {
+      // number of faces/vertices in this simple polygon
+      const Index np = C(p+1)-C(p);
+      for(Eigen::Index i = 0;i<np;i++)
+      {
+        assert(FF(k,0) == I(C(p)+((i+0)%np)));
+        assert(FF(k,1) == I(C(p)+((i+1)%np)));
+        assert(FF(k,2) == V.rows()+p);
+        NN.row(k*3+0) = N.row(C(p)+((i+0)%np));
+        NN.row(k*3+1) = N.row(C(p)+((i+1)%np));
+        NN.row(k*3+2) = FN.row(p);
+        k++;
+      }
+    }
+    assert(k == FF.rows());
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh

include/igl/per_corner_normals.h

@@ -55,6 +55,38 @@ namespace igl
     const std::vector<std::vector<IndexType> >& VF,
     const double corner_threshold,
     Eigen::PlainObjectBase<DerivedCN> & CN);
+  // Inputs:
+  //   V  #V by 3 list of mesh vertex positions
+  //   I  #I vectorized list of polygon corner indices into rows of some matrix V
+  //   C  #polygons+1 list of cumulative polygon sizes so that C(i+1)-C(i) = size of
+  //     the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the indices of
+  //     the ith polygon
+  //   corner_threshold  threshold in degrees on sharp angles
+  // Outputs:
+  //   N  #I by 3 list of per corner normals
+  //   VV  #I+#polygons by 3 list of auxiliary triangle mesh vertex positions
+  //   FF  #I by 3 list of triangle indices into rows of VV
+  //   J  #I list of indices into original polygons
+  //   NN  #FF by 3 list of normals for each auxiliary triangle
+  template <
+    typename DerivedV, 
+    typename DerivedI, 
+    typename DerivedC, 
+    typename DerivedN,
+    typename DerivedVV,
+    typename DerivedFF,
+    typename DerivedJ,
+    typename DerivedNN>
+  IGL_INLINE void per_corner_normals(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedI> & I,
+    const Eigen::MatrixBase<DerivedC> & C,
+    const typename DerivedV::Scalar corner_threshold,
+    Eigen::PlainObjectBase<DerivedN>  & N,
+    Eigen::PlainObjectBase<DerivedVV> & VV,
+    Eigen::PlainObjectBase<DerivedFF> & FF,
+    Eigen::PlainObjectBase<DerivedJ>  & J,
+    Eigen::PlainObjectBase<DerivedNN> & NN);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/per_face_normals.cpp

@@ -103,21 +103,76 @@ IGL_INLINE void igl::per_face_normals_stable(
 
 }
 
+#include "cotmatrix.h"
+
+template <
+  typename DerivedV,
+  typename DerivedI,
+  typename DerivedC,
+  typename DerivedN,
+  typename DerivedVV,
+  typename DerivedFF,
+  typename DerivedJ>
+IGL_INLINE void igl::per_face_normals(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedN> & N,
+  Eigen::PlainObjectBase<DerivedVV> & VV,
+  Eigen::PlainObjectBase<DerivedFF> & FF,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  assert(V.cols() == 3);
+  typedef Eigen::Index Index;
+  typedef typename DerivedN::Scalar Scalar;
+  // Use Bunge et al. algorithm in igl::cotmatrix to insert a point for each
+  // polygon which minimizes squared area.
+  {
+    Eigen::SparseMatrix<Scalar> _1,_2,P;
+    igl::cotmatrix(V,I,C,_1,_2,P);
+    VV = P*V;
+  }
+  // number of polygons
+  const Eigen::Index m = C.size()-1;
+  N.resize(m,3);
+  FF.resize(C(m),3);
+  J.resize(C(m));
+  {
+    Eigen::Index k = 0;
+    for(Eigen::Index p = 0;p<m;p++)
+    {
+      N.row(p).setZero();
+      // number of faces/vertices in this simple polygon
+      const Index np = C(p+1)-C(p);
+      for(Eigen::Index i = 0;i<np;i++)
+      {
+        FF.row(k) << 
+          I(C(p)+((i+0)%np)),
+          I(C(p)+((i+1)%np)),
+          V.rows()+p;
+        J(k) = p;
+        k++;
+        typedef Eigen::Matrix<Scalar,1,3> V3;
+        N.row(p) +=
+          V3(VV.row(I(C(p)+((i+0)%np)))-VV.row(V.rows()+p)).cross(
+          V3(VV.row(I(C(p)+((i+1)%np)))-VV.row(V.rows()+p)));
+      }
+      // normalize to take average
+      N.row(p) /= N.row(p).stableNorm();
+    }
+    assert(k == FF.rows());
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// Nonsense template
+// Nonsense template. Where'd this come from? AABB nonsense?
 namespace igl{template<> void per_face_normals<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&){} }
-// generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
-// generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
 template void igl::per_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::per_face_normals<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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);

include/igl/per_face_normals.h

@@ -42,6 +42,34 @@ namespace igl
     const Eigen::MatrixBase<DerivedV>& V,
     const Eigen::MatrixBase<DerivedF>& F,
     Eigen::PlainObjectBase<DerivedN> & N);
+  // Inputs:
+  //   V  #V by 3 list of mesh vertex positions
+  //   I  #I vectorized list of polygon corner indices into rows of some matrix V
+  //   C  #polygons+1 list of cumulative polygon sizes so that C(i+1)-C(i) = size of
+  //     the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the indices of
+  //     the ith polygon
+  //   corner_threshold  threshold in degrees on sharp angles
+  // Outputs:
+  //   N  #polygons by 3 list of per face normals
+  //   VV  #I+#polygons by 3 list of auxiliary triangle mesh vertex positions
+  //   FF  #I by 3 list of triangle indices into rows of VV
+  //   J  #I list of indices into original polygons
+  template <
+    typename DerivedV,
+    typename DerivedI,
+    typename DerivedC,
+    typename DerivedN,
+    typename DerivedVV,
+    typename DerivedFF,
+    typename DerivedJ>
+  IGL_INLINE void per_face_normals(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedI> & I,
+    const Eigen::MatrixBase<DerivedC> & C,
+    Eigen::PlainObjectBase<DerivedN> & N,
+    Eigen::PlainObjectBase<DerivedVV> & VV,
+    Eigen::PlainObjectBase<DerivedFF> & FF,
+    Eigen::PlainObjectBase<DerivedJ> & J);
 }
 
 #ifndef IGL_STATIC_LIBRARY