10-100x speedup on per corner normals + caching adjacency (#1745)

* improve documentation

* 10-100x speed up for per-corner normal computation. Cache adjacency.

* windows confused about templates

* windows error?

* rm commented code

* template for windows

* better documentation

include/igl/cumsum.cpp

@@ -81,6 +81,8 @@ IGL_INLINE void igl::cumsum(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::cumsum<Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 template void igl::cumsum<Eigen::Matrix<float, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
 template void igl::cumsum<Eigen::Matrix<double, 4, 1, 0, 4, 1>, Eigen::Matrix<double, 4, 1, 0, 4, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> >&);
 template void igl::cumsum<Eigen::Matrix<double, 1, 4, 1, 1, 4>, Eigen::Matrix<double, 1, 4, 1, 1, 4> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 4, 1, 1, 4> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 4, 1, 1, 4> >&);
@@ -94,5 +96,6 @@ template void igl::cumsum<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int
 #ifdef WIN32
 template void igl::cumsum<class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>, class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>>(class Eigen::MatrixBase<class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<unsigned __int64, -1, 1, 0, -1, 1>> &);
 template void igl::cumsum<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>, class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>>(class Eigen::MatrixBase<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<unsigned __int64, 2, 1, 0, 2, 1>> &);
+template void igl::cumsum<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> >(class Eigen::MatrixBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> > const&, int, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> >&);
 #endif
-#endif
+#endif

include/igl/per_corner_normals.cpp

@@ -1,108 +1,167 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 //
+// Copyright (C) 2021 Alec Jacobson <[email protected]>
 // Copyright (C) 2013 Alec Jacobson <[email protected]>
 //
 // 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 "per_corner_normals.h"
-
 #include "vertex_triangle_adjacency.h"
 #include "per_face_normals.h"
 #include "PI.h"
+#include "parallel_for.h"
 #include "doublearea.h"
-
-template <typename DerivedV, typename DerivedF, typename DerivedCN>
-IGL_INLINE void igl::per_corner_normals(
-  const Eigen::MatrixBase<DerivedV>& V,
-  const Eigen::MatrixBase<DerivedF>& F,
-  const double corner_threshold,
-  Eigen::PlainObjectBase<DerivedCN> & CN)
-{
-  using namespace Eigen;
-  using namespace std;
-  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> FN;
-  per_face_normals(V,F,FN);
-  vector<vector<int> > VF,VFi;
-  vertex_triangle_adjacency(V.rows(),F,VF,VFi);
-  return per_corner_normals(V,F,FN,VF,corner_threshold,CN);
-}
+#include <Eigen/Geometry>
 
 template <
   typename DerivedV,
   typename DerivedF,
-  typename DerivedFN,
   typename DerivedCN>
 IGL_INLINE void igl::per_corner_normals(
-  const Eigen::MatrixBase<DerivedV>& V,
-  const Eigen::MatrixBase<DerivedF>& F,
-  const Eigen::MatrixBase<DerivedFN>& FN,
-  const double corner_threshold,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const typename DerivedV::Scalar corner_threshold_degrees,
   Eigen::PlainObjectBase<DerivedCN> & CN)
 {
-  using namespace Eigen;
-  using namespace std;
-  vector<vector<int> > VF,VFi;
-  vertex_triangle_adjacency(V.rows(),F,VF,VFi);
-  return per_corner_normals(V,F,FN,VF,corner_threshold,CN);
+  Eigen::Matrix<Eigen::Index,Eigen::Dynamic,1> VF,NK;
+  vertex_triangle_adjacency(F,V.rows(),VF,NK);
+  return per_corner_normals(V,F,corner_threshold_degrees,VF,NK,CN);
 }
 
 template <
   typename DerivedV,
   typename DerivedF,
-  typename DerivedFN,
-  typename IndexType,
+  typename DerivedVF,
+  typename DerivedNI,
   typename DerivedCN>
 IGL_INLINE void igl::per_corner_normals(
-  const Eigen::MatrixBase<DerivedV>& /*V*/,
-  const Eigen::MatrixBase<DerivedF>& F,
-  const Eigen::MatrixBase<DerivedFN>& FN,
-  const std::vector<std::vector<IndexType> >& VF,
-  const double corner_threshold,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const typename DerivedV::Scalar corner_threshold_degrees,
+  const Eigen::MatrixBase<DerivedVF> & VF,
+  const Eigen::MatrixBase<DerivedNI> & NI,
   Eigen::PlainObjectBase<DerivedCN> & CN)
 {
-  using namespace Eigen;
-  using namespace std;
+  typedef typename DerivedV::Scalar Scalar;
+  typedef Eigen::Index Index;
+  // unit normals
+  Eigen::Matrix<Scalar,Eigen::Dynamic,3,Eigen::RowMajor> FN(F.rows(),3);
+  // face areas
+  Eigen::Matrix<Scalar,Eigen::Dynamic,1> FA(F.rows());
+  igl::parallel_for(F.rows(),[&](const Index f)
+  {
+    const Eigen::Matrix<Scalar,1,3> v10 = V.row(F(f,1))-V.row(F(f,0));
+    const Eigen::Matrix<Scalar,1,3> v20 = V.row(F(f,2))-V.row(F(f,0));
+    const Eigen::Matrix<Scalar,1,3> n = v10.cross(v20);
+    const Scalar a = n.norm();
+    FA(f) = a;
+    FN.row(f) = n/a;
+  },10000);
 
   // number of faces
-  const int m = F.rows();
+  const Index m = F.rows();
   // valence of faces
-  const int n = F.cols();
+  const Index n = F.cols();
+  assert(n == 3);
 
   // initialize output to ***zero***
   CN.setZero(m*n,3);
 
+  const Scalar cos_thresh = cos(corner_threshold_degrees*igl::PI/180);
   // loop over faces
-  for(size_t i = 0;int(i)<m;i++)
+  //for(Index i = 0;i<m;i++)
+  igl::parallel_for(F.rows(),[&](const Index i)
   {
     // Normal of this face
-    Eigen::Matrix<typename DerivedV::Scalar,3,1> fn = FN.row(i);
+    const auto & fnhat = FN.row(i);
     // loop over corners
-    for(size_t j = 0;int(j)<n;j++)
+    for(Index j = 0;j<n;j++)
     {
-      const std::vector<IndexType> &incident_faces = VF[F(i,j)];
-      // loop over faces sharing vertex of this corner
-      for(int k = 0;k<(int)incident_faces.size();k++)
+      const auto & v = F(i,j);
+      for(int k = NI[v]; k<NI[v+1]; k++)
       {
-        Eigen::Matrix<typename DerivedV::Scalar,3,1> ifn = FN.row(incident_faces[k]);
+        const auto & ifn = FN.row(VF[k]);
         // dot product between face's normal and other face's normal
-        double dp = fn.dot(ifn);
+        const Scalar dp = fnhat.dot(ifn);
         // if difference in normal is slight then add to average
-        if(dp > cos(corner_threshold*PI/180))
+        if(dp > cos_thresh)
         {
           // add to running sum
-          CN.row(i*n+j) += ifn;
-        // else ignore
-        }else
-        {
+          CN.row(i*n+j) += ifn*FA(VF[k]);
         }
       }
       // normalize to take average
       CN.row(i*n+j).normalize();
     }
-  }
+  },10000);
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedCI,
+  typename DerivedCC,
+  typename DerivedCN>
+IGL_INLINE void igl::per_corner_normals(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedCI> & CI,
+  const Eigen::MatrixBase<DerivedCC> & CC,
+  Eigen::PlainObjectBase<DerivedCN> & CN)
+{
+  typedef typename DerivedV::Scalar Scalar;
+  typedef Eigen::Index Index;
+  assert(CC.rows() == F.rows()*3+1);
+  // area weighted normals
+  Eigen::Matrix<Scalar,Eigen::Dynamic,3,Eigen::RowMajor> FN(F.rows(),3);
+  //for(Index f = 0;f<F.rows();f++)
+  igl::parallel_for(F.rows(),[&](const Index f)
+  {
+    const Eigen::Matrix<Scalar,1,3> v10 = V.row(F(f,1))-V.row(F(f,0));
+    const Eigen::Matrix<Scalar,1,3> v20 = V.row(F(f,2))-V.row(F(f,0));
+    FN.row(f) = v10.cross(v20);
+  },10000);
+
+  // number of faces
+  const Index m = F.rows();
+  // valence of faces
+  const Index n = F.cols();
+  assert(n == 3);
+
+  // initialize output to ***zero***
+  CN.setZero(m*n,3);
+  // loop over faces
+  igl::parallel_for(m*n,[&](const Index ci)
+  {
+    for(int k = CC(ci); k<CC(ci+1); k++)
+    {
+      // add to running sum
+      const auto cfk = CI(k);
+      CN.row(ci) += FN.row(cfk);
+    }
+    // normalize to take average
+    CN.row(ci).normalize();
+  },10000);
 }
 
+template <typename DerivedNV, typename DerivedNF, typename DerivedCN>
+IGL_INLINE void igl::per_corner_normals(
+  const Eigen::MatrixBase<DerivedNV> & NV,
+  const Eigen::MatrixBase<DerivedNF> & NF,
+  Eigen::PlainObjectBase<DerivedCN> & CN)
+{
+  const auto m = NF.rows();
+  const auto nc = NF.cols();
+  CN.resize(m*nc,3);
+  for(Eigen::Index i = 0;i<m;i++)
+  {
+    for(Eigen::Index c = 0;c<nc;c++)
+    {
+      CN.row(i*nc+c) = NV.row(NF(i,c));
+    }
+  }
+}
 
 template <
   typename DerivedV, 
@@ -117,7 +176,7 @@ 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,
+  const typename DerivedV::Scalar corner_threshold_degrees,
   Eigen::PlainObjectBase<DerivedN>  & N,
   Eigen::PlainObjectBase<DerivedVV> & VV,
   Eigen::PlainObjectBase<DerivedFF> & FF,
@@ -159,7 +218,7 @@ IGL_INLINE void igl::per_corner_normals(
         Eigen::Matrix<Scalar,3,1> ifn = 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))
+        if(dp > cos(corner_threshold_degrees*igl::PI/180))
         {
           // add to running sum
           N.row(C(p)+i) += AA(n) * ifn;
@@ -194,12 +253,12 @@ IGL_INLINE void igl::per_corner_normals(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
+template void igl::per_corner_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::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_corner_normals<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::Matrix<double, -1, 3, 0, -1, 3> >(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<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::per_corner_normals<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::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::per_corner_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::per_corner_normals<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::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::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<int, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, 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> >&);
-// generated by autoexplicit.sh
-template void igl::per_corner_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&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::per_corner_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&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::per_corner_normals<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
-template void igl::per_corner_normals<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
-template void igl::per_corner_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>, int, 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&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// Clang thinks this is the same as the one below, but Windows doesn't?
+// template void igl::per_corner_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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::per_corner_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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, double, Eigen::PlainObjectBase< Eigen::Matrix<double, -1, -1, 0, -1, -1> > &);
 #endif

include/igl/per_corner_normals.h

@@ -1,5 +1,6 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
+// Copyright (C) 2021 Alec Jacobson <[email protected]>
 // Copyright (C) 2013 Alec Jacobson <[email protected]>
 // 
 // This Source Code Form is subject to the terms of the Mozilla Public License 
@@ -13,47 +14,79 @@
 
 namespace igl
 {
-  // Compute vertex normals via vertex position list, face list
+  // Compute per corner normals for a triangle mesh by computing the
+  // area-weighted average of normals at incident faces whose normals deviate
+  // less than the provided threshold.
+  //
   // Inputs:
-  //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigen Matrix of face (triangle) indices
-  //   corner_threshold  threshold in degrees on sharp angles
-  // Output:
-  //   CN  #F*3 by 3 eigen Matrix of mesh vertex 3D normals, where the normal
-  //     for corner F(i,j) is at CN(i*3+j,:) 
-  template <typename DerivedV, typename DerivedF, typename DerivedCN>
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of mesh triangle indices into V
+  //   corner_threshold_degrees  threshold in degrees on sharp angles
+  // Outputs:
+  //   CN  #F*3 by 3 list of mesh vertex 3D normals, where the normal
+  //     for corner F(i,j) is at CN.row(i*3+j)
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedCN>
   IGL_INLINE void per_corner_normals(
-    const Eigen::MatrixBase<DerivedV>& V,
-    const Eigen::MatrixBase<DerivedF>& F,
-    const double corner_threshold,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const typename DerivedV::Scalar corner_threshold_degrees,
     Eigen::PlainObjectBase<DerivedCN> & CN);
-  // Other Inputs:
-  //   FN  #F by 3 eigen Matrix of face normals
+  // Inputs:
+  //   VF  3*#F list  List of faces indice on each vertex, so that VF(NI(i)+j) =
+  //     f, means that face f is the jth face (in no particular order) incident
+  //     on vertex i.
+  //   NI  #V+1 list  cumulative sum of vertex-triangle degrees with a
+  //     preceeding zero. "How many faces" have been seen before visiting this
+  //     vertex and its incident faces.
+  //
+  // See also: vertex_triangle_adjacency
   template <
-    typename DerivedV, 
-    typename DerivedF, 
-    typename DerivedFN, 
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedVF,
+    typename DerivedNI,
     typename DerivedCN>
   IGL_INLINE void per_corner_normals(
-    const Eigen::MatrixBase<DerivedV>& V,
-    const Eigen::MatrixBase<DerivedF>& F,
-    const Eigen::MatrixBase<DerivedFN>& FN,
-    const double corner_threshold,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const typename DerivedV::Scalar corner_threshold_degrees,
+    const Eigen::MatrixBase<DerivedVF> & VF,
+    const Eigen::MatrixBase<DerivedNI> & NI,
     Eigen::PlainObjectBase<DerivedCN> & CN);
-  // Other Inputs:
-  //   VF  map from vertices to list of incident faces
+  // Inputs:
+  //   CI  #CI list of face neighbors as indices into rows of F
+  //   CC  3*#F+1 list of cumulative sizes so that CC(i*3+j+1) - CC(i*3+j) is
+  //     the number of faces considered smoothly incident on corner at F(i,j)
+  //
+  // See also smooth_corner_adjacency
   template <
-    typename DerivedV, 
-    typename DerivedF, 
-    typename DerivedFN, 
-    typename IndexType,
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedCI,
+    typename DerivedCC,
     typename DerivedCN>
   IGL_INLINE void per_corner_normals(
-    const Eigen::MatrixBase<DerivedV>& V,
-    const Eigen::MatrixBase<DerivedF>& F,
-    const Eigen::MatrixBase<DerivedFN>& FN,
-    const std::vector<std::vector<IndexType> >& VF,
-    const double corner_threshold,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedCI> & CI,
+    const Eigen::MatrixBase<DerivedCC> & CC,
+    Eigen::PlainObjectBase<DerivedCN> & CN);
+  // Given indexed normals (e.g., read from a .obj file), explode into
+  // per-corner normals (e.g., as expected by igl::opengl::ViewerData)
+  //
+  // Inputs:
+  //   NV  #NV by 3 list of index normal vectors
+  //   NF  #F by nc list of indices into rows of NV
+  // Outputs
+  //   CN  #F*nc by 3 list of per-corner normals so that 
+  //     CN.row(i*nc+c) = NV.row(NF(i,c))
+  template <typename DerivedNV, typename DerivedNF, typename DerivedCN>
+  IGL_INLINE void per_corner_normals(
+    const Eigen::MatrixBase<DerivedNV> & NV,
+    const Eigen::MatrixBase<DerivedNF> & NF,
     Eigen::PlainObjectBase<DerivedCN> & CN);
   // Inputs:
   //   V  #V by 3 list of mesh vertex positions
@@ -81,7 +114,7 @@ namespace igl
     const Eigen::MatrixBase<DerivedV> & V,
     const Eigen::MatrixBase<DerivedI> & I,
     const Eigen::MatrixBase<DerivedC> & C,
-    const typename DerivedV::Scalar corner_threshold,
+    const typename DerivedV::Scalar corner_threshold_degrees,
     Eigen::PlainObjectBase<DerivedN>  & N,
     Eigen::PlainObjectBase<DerivedVV> & VV,
     Eigen::PlainObjectBase<DerivedFF> & FF,

include/igl/smooth_corner_adjacency.cpp

@@ -0,0 +1,148 @@
+#include "smooth_corner_adjacency.h"
+#include "vertex_triangle_adjacency.h"
+#include "matlab_format.h"
+#include "parallel_for.h"
+#include "unzip_corners.h"
+#include <iostream>
+
+void igl::smooth_corner_adjacency(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const double corner_threshold_radians,
+  Eigen::VectorXi & CI,
+  Eigen::VectorXi & CC)
+{
+  typedef double Scalar;
+  typedef Eigen::Index Index;
+  Eigen::Matrix<Index,Eigen::Dynamic,1> VF,NI;
+  igl::vertex_triangle_adjacency(F,V.rows(),VF,NI);
+  // unit normals
+  Eigen::Matrix<Scalar,Eigen::Dynamic,3,Eigen::RowMajor> FN(F.rows(),3);
+  igl::parallel_for(F.rows(),[&](const Index f)
+  {
+    const Eigen::Matrix<Scalar,1,3> v10 = V.row(F(f,1))-V.row(F(f,0));
+    const Eigen::Matrix<Scalar,1,3> v20 = V.row(F(f,2))-V.row(F(f,0));
+    const Eigen::Matrix<Scalar,1,3> n = v10.cross(v20);
+    const Scalar a = n.norm();
+    FN.row(f) = n/a;
+  },10000);
+
+  // number of faces
+  const Index m = F.rows();
+  // valence of faces
+  const Index n = F.cols();
+  assert(n == 3);
+
+  CI.resize(m*n*8);
+  CI.setConstant(-1);
+  Index ncc = 0;
+  Index ci = -1;
+  // assumes that ci is strictly increasing and we're appending to CI
+  const auto append_CI = [&](Index nf)
+  {
+    // make room
+    if(ncc >= CI.size()) { CI.conservativeResize(CI.size()*2+1); }
+    CI(ncc++) = nf;
+    CC(ci+1)++;
+  };
+  CC.resize(m*3+1);
+  CC.setConstant(-1);
+  CC(0) = 0;
+
+  const Scalar cos_thresh = cos(corner_threshold_radians);
+  // parallelizing this probably requires map-reduce
+  for(Index i = 0;i<m;i++)
+  {
+    // Normal of this face
+    const auto & fnhat = FN.row(i);
+    // loop over corners
+    for(Index j = 0;j<n;j++)
+    {
+      // increment ci
+      ci++;
+      assert(ci == i*n+j);
+      CC(ci+1) = CC(ci);
+      const auto & v = F(i,j);
+      for(int k = NI(v); k<NI(v+1); k++)
+      {
+        const int nf = VF(k);
+        const auto & ifn = FN.row(nf);
+        // dot product between face's normal and other face's normal
+        const Scalar dp = fnhat.dot(ifn);
+        // if difference in normal is slight then add to average
+        if(dp > cos_thresh)
+        {
+          append_CI(nf);
+        }
+      }
+    }
+  }
+  CI.conservativeResize(ncc);
+}
+
+
+void igl::smooth_corner_adjacency(
+  const Eigen::MatrixXi & FV,
+  const Eigen::MatrixXi & FN,
+  Eigen::VectorXi & CI,
+  Eigen::VectorXi & CC)
+{
+  typedef double Scalar;
+  typedef Eigen::Index Index;
+  assert(FV.rows() == FN.rows());
+  assert(FV.cols() == 3);
+  assert(FN.cols() == 3);
+  Eigen::VectorXi J;
+  Index nu = -1;
+  {
+    Eigen::MatrixXi U;
+    Eigen::MatrixXi _;
+    igl::unzip_corners<const Eigen::MatrixXi>({FV,FN},U,_,J);
+    nu = U.rows();
+    assert(J.maxCoeff() == nu-1);
+  }
+  // could use linear arrays here if every becomes bottleneck
+  std::vector<std::vector<Index>> U2F(nu);
+  const Index m = FV.rows();
+  for(Index j = 0;j<3;j++)
+  {
+    for(Index i = 0;i<m;i++)
+    {
+      // unzip_corners uses convention i+j*#F
+      const Index ci = i+j*m;
+      U2F[J(ci)].emplace_back(i);
+    }
+  }
+  CC.resize(m*3+1);
+  CC.setConstant(-1);
+  CC(0) = 0;
+  CI.resize(m*3*8);
+  int ncc = 0;
+  // assumes that ci is strictly increasing and we're appending to CI
+  const auto append_CI = [&CI,&CC,&ncc](Index ci, Index nf)
+  {
+    // make room
+    if(ncc >= CI.size()) { CI.conservativeResize(CI.size()*2+1); }
+    CI(ncc++) = nf;
+    CC(ci+1)++;
+  };
+  for(Index i = 0;i<m;i++)
+  {
+    for(Index j = 0;j<3;j++)
+    {
+      const Index J_ci = i+j*m;
+      // CI,CC uses convention i*3 + j
+      const Index C_ci = i*3+j;
+      CC(C_ci+1) = CC(C_ci);
+      for(const auto & nf : U2F[J(J_ci)])
+      {
+        append_CI(C_ci,nf);
+      }
+    }
+  }
+  CI.conservativeResize(ncc);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiations
+#endif

include/igl/smooth_corner_adjacency.h

@@ -0,0 +1,49 @@
+#ifndef IGL_SMOOTH_CORNER_ADJACENCY_H
+#define IGL_SMOOTH_CORNER_ADJACENCY_H
+
+#include <Eigen/Core>
+#include "igl_inline.h"
+
+namespace igl
+{
+  // Determine the corner-to-face adjacency relationship that can be used for
+  // computing crease-aware per-corner normals.
+  //
+  // Inputs:
+  //   V  #V by 3 list of mesh vertex positions
+  //   F  #F by 3 list of triangle mesh indices into rows of V
+  //   corner_threshold_radians  dihedral angle considered non-smooth (in
+  //     radians)
+  // Outputs:
+  //   CI  #CI list of face neighbors as indices into rows of F
+  //   CC  3*#F+1 list of cumulative sizes so that CC(i*3+j+1) - CC(i*3+j) is
+  //     the number of faces considered smoothly incident on corner at F(i,j)
+  void smooth_corner_adjacency(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const double corner_threshold_radians,
+    Eigen::VectorXi & CI,
+    Eigen::VectorXi & CC);
+  // Determine the effective corner-to-face adjacency relationship implied by a
+  // set of indexed vertex positions (FV) and normals (FV) (e.g., those read in
+  // from a .obj file).
+  //
+  // Inputs:
+  //   FV  #F by 3 list of triangle mesh indices into rows of some V
+  //   FN  #F by 3 list of triangle mesh indices into rows of some N
+  // Outputs:
+  //   CI  #CI list of face neighbors as indices into rows of F
+  //   CC  3*#F+1 list of cumulative sizes so that CC(i*3+j+1) - CC(i*3+j) is
+  //     the number of faces considered smoothly incident on corner at F(i,j)
+  void smooth_corner_adjacency(
+    const Eigen::MatrixXi & FV,
+    const Eigen::MatrixXi & FN,
+    Eigen::VectorXi & CI,
+    Eigen::VectorXi & CC);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "smooth_corner_adjacency.cpp"
+#endif
+
+#endif

include/igl/unzip_corners.h

@@ -32,6 +32,7 @@ namespace igl
   //   U  #U by #A list of indices into each attribute for each unique mesh
   //     vertex: U(v,a) is the attribute index of vertex v in attribute a.
   //   G  #F by 3 list of triangle indices into U
+  //   J  #F*3 by 1 list of indices so that A[](i,j) = U.row(i+j*#F)
   // Example:
   //   [V,F,TC,FTC] = readOBJ('~/Downloads/kiwis/kiwi.obj');
   //   [U,G] = unzip_corners(cat(3,F,FTC));

include/igl/vertex_triangle_adjacency.cpp

@@ -84,6 +84,8 @@ IGL_INLINE void igl::vertex_triangle_adjacency(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, 3, 1, -1, 3>, unsigned long, unsigned long>(Eigen::Matrix<int, -1, 3, 1, -1, 3>::Scalar, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > >&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > >&);
 // generated by autoexplicit.sh
 template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
@@ -100,5 +102,6 @@ template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, 3, 0, -1, 3>
 template void igl::vertex_triangle_adjacency<class Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned __int64, unsigned __int64>(int, class Eigen::MatrixBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>> const &, class std::vector<class std::vector<unsigned __int64, class std::allocator<unsigned __int64>>, class std::allocator<class std::vector<unsigned __int64, class std::allocator<unsigned __int64>>>> &, class std::vector<class std::vector<unsigned __int64, class std::allocator<unsigned __int64>>, class std::allocator<class std::vector<unsigned __int64, class std::allocator<unsigned __int64>>>> &);
 template void igl::vertex_triangle_adjacency<class Eigen::Matrix<int, -1, 3, 1, -1, 3>, unsigned __int64, unsigned __int64>(int, class Eigen::MatrixBase<class Eigen::Matrix<int, -1, 3, 1, -1, 3>> const &, class std::vector<class std::vector<unsigned __int64, class std::allocator<unsigned __int64>>, class std::allocator<class std::vector<unsigned __int64, class std::allocator<unsigned __int64>>>> &, class std::vector<class std::vector<unsigned __int64, class std::allocator<unsigned __int64>>, class std::allocator<class std::vector<unsigned __int64, class std::allocator<unsigned __int64>>>> &);
 template void igl::vertex_triangle_adjacency<class Eigen::Matrix<int,-1,-1,0,-1,-1>,__int64,__int64>(int,class Eigen::MatrixBase<class Eigen::Matrix<int,-1,-1,0,-1,-1> > const &,class std::vector<class std::vector<__int64,class std::allocator<__int64> >,class std::allocator<class std::vector<__int64,class std::allocator<__int64> > > > &,class std::vector<class std::vector<__int64,class std::allocator<__int64> >,class std::allocator<class std::vector<__int64,class std::allocator<__int64> > > > &);
+template void igl::vertex_triangle_adjacency<class Eigen::Matrix<int,-1,-1,0,-1,-1>,class Eigen::Matrix<__int64,-1,1,0,-1,1>,class Eigen::Matrix<__int64,-1,1,0,-1,1> >(class Eigen::MatrixBase<class Eigen::Matrix<int,-1,-1,0,-1,-1> > const &,int,class Eigen::PlainObjectBase<class Eigen::Matrix<__int64,-1,1,0,-1,1> > &,class Eigen::PlainObjectBase<class Eigen::Matrix<__int64,-1,1,0,-1,1> > &);
 #endif
 #endif