improve performance for manifold TT

Former-commit-id: 73655a1758d381c206891934410eb43a01ae5771

include/igl/triangle_triangle_adjacency.cpp

@@ -1,11 +1,13 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 //
+// Copyright (C) 2018 Alec Jacobson, Marc Alexa
 // Copyright (C) 2014 Daniele Panozzo <daniele.panozzo@gmail.com>
 //
 // This Source Code Form is subject to the terms of the Mozilla Public License
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "triangle_triangle_adjacency.h"
+#include "vertex_triangle_adjacency.h"
 #include "parallel_for.h"
 #include "unique_edge_map.h"
 #include <algorithm>
@@ -37,8 +39,35 @@ IGL_INLINE void igl::triangle_triangle_adjacency(
   const Eigen::MatrixBase<DerivedF>& F,
   Eigen::PlainObjectBase<DerivedTT>& TT)
 {
-  DerivedTT TTi;
-  return triangle_triangle_adjacency(F,TT,TTi);
+  const int n = F.maxCoeff()+1;
+  typedef Eigen::Matrix<typename DerivedTT::Scalar,Eigen::Dynamic,1> VectorXI;
+  VectorXI VF,NI;
+  vertex_triangle_adjacency(F,n,VF,NI);
+  TT = DerivedTT::Constant(F.rows(),3,-1);
+  // Loop over faces
+  igl::parallel_for(F.rows(),[&](int f)
+  {
+    // Loop over corners
+    for (int k = 0; k < 3; k++)
+    {
+      int vi = F(f,k), vin = F(f,(k+1)%3);
+      // Loop over face neighbors incident on this corner
+      for (int j = NI[vi]; j < NI[vi+1]; j++)
+      {
+        int fn = VF[j];
+        // Not this face
+        if (fn != f)
+        {
+          // Face neighbor also has [vi,vin] edge
+          if (F(fn,0) == vin || F(fn,1) == vin || F(fn,2) == vin)
+          {
+            TT(f,k) = fn;
+            break;
+          }
+        }
+      }
+    }
+  });
 }
 
 template <typename DerivedF, typename TTT_type>
@@ -89,10 +118,29 @@ IGL_INLINE void igl::triangle_triangle_adjacency(
   Eigen::PlainObjectBase<DerivedTT>& TT,
   Eigen::PlainObjectBase<DerivedTTi>& TTi)
 {
-  std::vector<std::vector<int> > TTT;
-  triangle_triangle_adjacency_preprocess(F,TTT);
-  triangle_triangle_adjacency_extractTT(F,TTT,TT);
-  triangle_triangle_adjacency_extractTTi(F,TTT,TTi);
+  triangle_triangle_adjacency(F,TT);
+  TTi = DerivedTTi::Constant(TT.rows(),TT.cols(),-1);
+  //for(int f = 0; f<F.rows(); f++)
+  igl::parallel_for(F.rows(),[&](int f)
+  {
+    for(int k = 0;k<3;k++)
+    {
+      int vi = F(f,k), vj = F(f,(k+1)%3);
+      int fn = TT(f,k);
+      if(fn >= 0)
+      {
+        for(int kn = 0;kn<3;kn++)
+        {
+          int vin = F(fn,kn), vjn = F(fn,(kn+1)%3);
+          if(vi == vjn && vin == vj)
+          {
+            TTi(f,k) = kn;
+            break;
+          }
+        }
+      }
+    }
+  });
 }
 
 template <

include/igl/vertex_triangle_adjacency.cpp

@@ -6,11 +6,12 @@
 // 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 "vertex_triangle_adjacency.h"
+#include "cumsum.h"
 
 template <typename DerivedF, typename VFType, typename VFiType>
 IGL_INLINE void igl::vertex_triangle_adjacency(
   const typename DerivedF::Scalar n,
-  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedF>& F,
   std::vector<std::vector<VFType> >& VF,
   std::vector<std::vector<VFiType> >& VFi)
 {
@@ -34,28 +35,66 @@ IGL_INLINE void igl::vertex_triangle_adjacency(
 
 template <typename DerivedV, typename DerivedF, typename IndexType>
 IGL_INLINE void igl::vertex_triangle_adjacency(
-  const Eigen::PlainObjectBase<DerivedV>& V,
-  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
   std::vector<std::vector<IndexType> >& VF,
   std::vector<std::vector<IndexType> >& VFi)
 {
   return vertex_triangle_adjacency(V.rows(),F,VF,VFi);
 }
 
+template <
+  typename DerivedF,
+  typename DerivedVF,
+  typename DerivedNI>
+IGL_INLINE void igl::vertex_triangle_adjacency(
+  const Eigen::MatrixBase<DerivedF> & F,
+  const int n,
+  Eigen::PlainObjectBase<DerivedVF> & VF,
+  Eigen::PlainObjectBase<DerivedNI> & NI)
+{
+  typedef Eigen::Matrix<typename DerivedVF::Scalar,Eigen::Dynamic,1> VectorXI;
+  // vfd  #V list so that vfd(i) contains the vertex-face degree (number of
+  // faces incident on vertex i)
+  VectorXI vfd = VectorXI::Zero(n);
+  for (int i = 0; i < F.rows(); i++)
+  {
+    for (int j = 0; j < 3; j++)
+    {
+      vfd[F(i,j)]++;
+    }
+  }
+  igl::cumsum(vfd,1,NI);
+  // Prepend a zero
+  NI = (DerivedNI(n+1)<<0,NI).finished();
+  // vfd now acts as a counter
+  vfd = NI;
+
+  VF = Eigen::VectorXi(3*F.rows());
+  for (int i = 0; i < F.rows(); i++)
+  {
+    for (int j = 0; j < 3; j++)
+    {
+      VF[vfd[F(i,j)]] = i;
+      vfd[F(i,j)]++;
+    }
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // 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::PlainObjectBase<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> > > >&);
+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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<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> > > >&);
-template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, unsigned int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > >&, std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > >&);
-template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > 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> > > >&);
-template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, long, long>(Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
-template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, long, long>(Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
-template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned long, unsigned long>(Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > 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> > > >&);
-template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<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> > > >&);
+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> > > >&);
+template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, unsigned int>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > >&, std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > >&);
+template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, int>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > 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> > > >&);
+template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, long, long>(Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
+template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, long, long>(Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
+template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned long, unsigned long>(Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > 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> > > >&);
+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> > > >&);
 #ifdef WIN32
-template void igl::vertex_triangle_adjacency<class Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned __int64, unsigned __int64>(int, class Eigen::PlainObjectBase<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::PlainObjectBase<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>, 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>>>> &);
 #endif
 #endif

include/igl/vertex_triangle_adjacency.h

@@ -7,7 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_VERTEX_TRIANGLE_ADJACENCY_H
 #define IGL_VERTEX_TRIANGLE_ADJACENCY_H
-#include <igl/igl_inline.h>
+#include "igl_inline.h"
 
 #include <Eigen/Dense>
 #include <vector>
@@ -34,15 +34,34 @@ namespace igl
   template <typename DerivedF, typename VFType, typename VFiType>
   IGL_INLINE void vertex_triangle_adjacency(
     const typename DerivedF::Scalar n,
-    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedF>& F,
     std::vector<std::vector<VFType> >& VF,
     std::vector<std::vector<VFiType> >& VFi);
   template <typename DerivedV, typename DerivedF, typename IndexType>
   IGL_INLINE void vertex_triangle_adjacency(
-    const Eigen::PlainObjectBase<DerivedV>& V,
-    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::MatrixBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedF>& F,
     std::vector<std::vector<IndexType> >& VF,
     std::vector<std::vector<IndexType> >& VFi);
+  // Inputs:
+  //   F  #F by 3 list of triangle indices into some vertex list V
+  //   n  number of vertices, #V (e.g., F.maxCoeff()+1)
+  // Outputs:
+  //   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.
+  template <
+    typename DerivedF,
+    typename DerivedVF,
+    typename DerivedNI>
+  IGL_INLINE void vertex_triangle_adjacency(
+    const Eigen::MatrixBase<DerivedF> & F,
+    const int n,
+    Eigen::PlainObjectBase<DerivedVF> & VF,
+    Eigen::PlainObjectBase<DerivedNI> & NI);
 }
 
 #ifndef IGL_STATIC_LIBRARY