Merge branch 'dev' of https://github.com/libigl/libigl into dev

.appveyor.yml

@@ -1,16 +1,20 @@
 version: 1.0.{build}
 os: Visual Studio 2017
-test: off
+platform: x64
 clone_folder: C:\projects\libigl
+shallow_clone: true
 branches:
   only:
     - master
     - dev
 environment:
-  BOOST_ROOT: C:/Libraries/boost_1_65_1
+  matrix:
+  - config: Debug
+    BOOST_ROOT: C:/Libraries/boost_1_65_1
 install:
-  - git submodule update --init --recursive
   - cinstall: python
+build:
+  parallel: true
 build_script:
   - cd c:\projects\libigl
   # External libraries (boost)
@@ -32,7 +36,10 @@ build_script:
   # Tutorials and tests
   - mkdir build
   - cd build
-  - cmake -DLIBIGL_BUILD_TESTS=OFF -DLIBIGL_BUILD_TUTORIALS=ON -DLIBIGL_BUILD_PYTHON=OFF -G "Visual Studio 15 2017 Win64" ../
+  - cmake -DCMAKE_BUILD_TYPE=%config% -DLIBIGL_WITH_COMISO=OFF -DLIBIGL_BUILD_PYTHON=OFF -G "Visual Studio 15 2017 Win64" ../
   - set MSBuildLogger="C:\Program Files\AppVeyor\BuildAgent\Appveyor.MSBuildLogger.dll"
-  - set MSBuildOptions=/v:m /p:Configuration=Debug /logger:%MSBuildLogger%
+  - set MSBuildOptions=/v:m /p:Configuration=%config% /logger:%MSBuildLogger%
   - msbuild %MSBuildOptions% libigl.sln
+
+test_script:
+  - ctest -C %config% --verbose --output-on-failure -j 2

.travis.yml

@@ -60,7 +60,7 @@ script:
 # Tutorials and tests
 - mkdir build
 - cd build
-- cmake -DCMAKE_BUILD_TYPE=$CONFIG -DLIBIGL_BUILD_TESTS=ON -DLIBIGL_BUILD_TUTORIALS=ON ../
+- cmake -DCMAKE_BUILD_TYPE=$CONFIG -DLIBIGL_BUILD_TESTS=ON -DLIBIGL_BUILD_TUTORIALS=ON -DEMBREE_ISA_AVX=OFF -DEMBREE_ISA_AVX2=OFF -DEMBREE_ISA_AVX512SKX=OFF ../
 - make -j 2
 - ctest --verbose
 - ccache --show-stats

CMakeLists.txt

@@ -40,6 +40,9 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 ### conditionally compile certain modules depending on libraries found on the system
 list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_LIST_DIR}/cmake)
 
+### Set compiler flags for building the tests on Windows with Visual Studio
+include(LibiglWindows)
+
 ### Adding libIGL: choose the path to your local copy libIGL
 include(libigl)
 
@@ -48,6 +51,7 @@ if(LIBIGL_BUILD_TUTORIALS)
 endif()
 
 if(LIBIGL_BUILD_TESTS)
+	include(CTest)
 	enable_testing()
 	add_subdirectory(tests)
 endif()

cmake/LibiglDownloadExternal.cmake

@@ -58,8 +58,8 @@ endfunction()
 ## Embree
 function(igl_download_embree)
 	igl_download_project(embree
-		URL            https://github.com/embree/embree/archive/v2.17.4.tar.gz
-		URL_MD5        2038f3216b1d626e87453aee72c470e5
+		URL            https://github.com/embree/embree/archive/v3.2.3.tar.gz
+		URL_MD5        1868cda1c97d83d7a0b67b0b64b18cef
 		# GIT_REPOSITORY https://github.com/embree/embree.git
 		# GIT_TAG        cb61322db3bb7082caed21913ad14869b436fe78
 	)
@@ -77,7 +77,7 @@ endfunction()
 function(igl_download_glfw)
 	igl_download_project(glfw
 		GIT_REPOSITORY https://github.com/glfw/glfw.git
-		GIT_TAG        58cc4b2c5c2c9a245e09451437dd6f5af4d60c84
+		GIT_TAG        53c8c72c676ca97c10aedfe3d0eb4271c5b23dba
 	)
 endfunction()
 
@@ -133,11 +133,11 @@ function(igl_download_triangle)
 	)
 endfunction()
 
-## Google test
-function(igl_download_googletest)
-	igl_download_project(googletest
-		GIT_REPOSITORY https://github.com/google/googletest
-		GIT_TAG        release-1.8.1
+## Catch2
+function(igl_download_catch2)
+	igl_download_project(catch2
+		GIT_REPOSITORY https://github.com/catchorg/Catch2.git
+		GIT_TAG        03d122a35c3f5c398c43095a87bc82ed44642516
 	)
 endfunction()
 

cmake/LibiglWindows.cmake

@@ -0,0 +1,22 @@
+if(MSVC)
+	if("${MSVC_RUNTIME}" STREQUAL "")
+		set(MSVC_RUNTIME "static")
+	endif()
+	if(${MSVC_RUNTIME} STREQUAL "static")
+		message(STATUS "MSVC -> forcing use of statically-linked runtime.")
+		foreach(config ${CMAKE_CONFIGURATION_TYPES})
+			string(TOUPPER ${config} config)
+			string(REPLACE /MD /MT CMAKE_C_FLAGS_${config} "${CMAKE_C_FLAGS_${config}}")
+			string(REPLACE /MD /MT CMAKE_CXX_FLAGS_${config} "${CMAKE_CXX_FLAGS_${config}}")
+		endforeach()
+		string(REPLACE "/MDd" "/MTd" CMAKE_CXX_FLAGS_DEBUG ${CMAKE_CXX_FLAGS_DEBUG})
+	else()
+		message(STATUS "MSVC -> forcing use of dynamically-linked runtime.")
+		foreach(config ${CMAKE_CONFIGURATION_TYPES})
+			string(TOUPPER ${config} config)
+			string(REPLACE /MT /MD CMAKE_C_FLAGS_${config} "${CMAKE_C_FLAGS_${config}}")
+			string(REPLACE /MT /MD CMAKE_CXX_FLAGS_${config} "${CMAKE_CXX_FLAGS_${config}}")
+		endforeach()
+		string(REPLACE "/MTd" "/MDd" CMAKE_CXX_FLAGS_DEBUG ${CMAKE_CXX_FLAGS_DEBUG})
+	endif()
+endif()

cmake/libigl.cmake

@@ -265,39 +265,28 @@ endif()
 if(LIBIGL_WITH_EMBREE)
   set(EMBREE_DIR "${LIBIGL_EXTERNAL}/embree")
 
+  set(EMBREE_TESTING_INTENSITY 0 CACHE INT "" FORCE)
   set(EMBREE_ISPC_SUPPORT OFF CACHE BOOL " " FORCE)
   set(EMBREE_TASKING_SYSTEM "INTERNAL" CACHE BOOL " " FORCE)
   set(EMBREE_TUTORIALS OFF CACHE BOOL " " FORCE)
   set(EMBREE_MAX_ISA NONE CACHE STRINGS " " FORCE)
-  set(BUILD_TESTING OFF CACHE BOOL " " FORCE)
-
-  # set(ENABLE_INSTALLER OFF CACHE BOOL " " FORCE)
+  set(EMBREE_STATIC_LIB ON CACHE BOOL " " FORCE)
   if(MSVC)
-    # set(EMBREE_STATIC_RUNTIME OFF CACHE BOOL " " FORCE)
-    set(EMBREE_STATIC_LIB OFF CACHE BOOL " " FORCE)
-  else()
-    set(EMBREE_STATIC_LIB ON CACHE BOOL " " FORCE)
+    set(EMBREE_STATIC_RUNTIME ON CACHE BOOL " " FORCE)
   endif()
 
   if(NOT TARGET embree)
+    # TODO: Should probably save/restore the CMAKE_CXX_FLAGS_*, since embree seems to be
+    # overriding them on Windows. But well... it works for now.
     igl_download_embree()
     add_subdirectory("${EMBREE_DIR}" "embree")
   endif()
 
-  if(MSVC)
-    add_custom_target(Copy-Embree-DLL ALL # Adds a post-build event
-        COMMAND ${CMAKE_COMMAND} -E copy_if_different # which executes "cmake - E
-        $<TARGET_FILE:embree> # <--this is in-file
-        "${CMAKE_BINARY_DIR}" # <--this is out-file path
-        DEPENDS embree) # Execute after embree target has been built
-  endif()
-
   compile_igl_module("embree")
   target_link_libraries(igl_embree ${IGL_SCOPE} embree)
   target_include_directories(igl_embree ${IGL_SCOPE} ${EMBREE_DIR}/include)
-  if(NOT MSVC)
-    target_compile_definitions(igl_embree ${IGL_SCOPE} -DENABLE_STATIC_LIB)
-  endif()
+  target_compile_definitions(igl_embree ${IGL_SCOPE} -DEMBREE_STATIC_LIB)
+
 endif()
 
 ################################################################################

include/igl/bbw.cpp

@@ -118,7 +118,12 @@ IGL_INLINE bool igl::bbw(
     }
     W.col(i) = Wi;
   };
+#ifdef WIN32
+  for (int i = 0; i < m; ++i)
+    optimize_weight(i);
+#else
   parallel_for(m,optimize_weight,2);
+#endif
   if(error)
   {
     return false;

include/igl/boundary_facets.cpp

@@ -103,7 +103,7 @@ IGL_INLINE void igl::boundary_facets(
 
 template <typename DerivedT, typename DerivedF>
 IGL_INLINE void igl::boundary_facets(
-  const Eigen::PlainObjectBase<DerivedT>& T,
+  const Eigen::MatrixBase<DerivedT>& T,
   Eigen::PlainObjectBase<DerivedF>& F)
 {
   assert(T.cols() == 0 || T.cols() == 4 || T.cols() == 3);
@@ -119,7 +119,7 @@ IGL_INLINE void igl::boundary_facets(
 
 template <typename DerivedT, typename Ret>
 Ret igl::boundary_facets(
-  const Eigen::PlainObjectBase<DerivedT>& T)
+  const Eigen::MatrixBase<DerivedT>& T)
 {
   Ret F;
   igl::boundary_facets(T,F);
@@ -130,12 +130,12 @@ Ret igl::boundary_facets(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
-template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
+template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
 // generated by autoexplicit.sh
-template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
-template void igl::boundary_facets<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 void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::boundary_facets<int, int>(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
-//template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > igl::boundary_facets(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
-template Eigen::Matrix<int, -1, -1, 0, -1, -1> igl::boundary_facets<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&);
-template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
+//template Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > igl::boundary_facets(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
+template Eigen::Matrix<int, -1, -1, 0, -1, -1> igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
+template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
 #endif

include/igl/boundary_facets.h

@@ -37,12 +37,12 @@ namespace igl
   //   DerivedF  integer-value: i.e. from MatrixXi
   template <typename DerivedT, typename DerivedF>
   IGL_INLINE void boundary_facets(
-    const Eigen::PlainObjectBase<DerivedT>& T,
+    const Eigen::MatrixBase<DerivedT>& T,
     Eigen::PlainObjectBase<DerivedF>& F);
   // Same as above but returns F
   template <typename DerivedT, typename Ret>
   Ret boundary_facets(
-    const Eigen::PlainObjectBase<DerivedT>& T);
+    const Eigen::MatrixBase<DerivedT>& T);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/circulation.cpp

@@ -7,12 +7,11 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "circulation.h"
 #include "list_to_matrix.h"
+#include <cassert>
 
 IGL_INLINE std::vector<int> igl::circulation(
   const int e,
   const bool ccw,
-  const Eigen::MatrixXi & F,
-  const Eigen::MatrixXi & E,
   const Eigen::VectorXi & EMAP,
   const Eigen::MatrixXi & EF,
   const Eigen::MatrixXi & EI)
@@ -20,7 +19,8 @@ IGL_INLINE std::vector<int> igl::circulation(
   // prepare output
   std::vector<int> N;
   N.reserve(6);
-  const int m = F.rows();
+  const int m = EMAP.size()/3;
+  assert(m*3 == EMAP.size());
   const auto & step = [&](
     const int e, 
     const int ff,
@@ -59,13 +59,11 @@ IGL_INLINE std::vector<int> igl::circulation(
 IGL_INLINE void igl::circulation(
   const int e,
   const bool ccw,
-  const Eigen::MatrixXi & F,
-  const Eigen::MatrixXi & E,
   const Eigen::VectorXi & EMAP,
   const Eigen::MatrixXi & EF,
   const Eigen::MatrixXi & EI,
   Eigen::VectorXi & vN)
 {
-  std::vector<int> N = circulation(e,ccw,F,E,EMAP,EF,EI);
+  std::vector<int> N = circulation(e,ccw,EMAP,EF,EI);
   igl::list_to_matrix(N,vN);
 }

include/igl/circulation.h

@@ -20,8 +20,6 @@ namespace igl
   //   e  index into E of edge to circulate
   //   ccw  whether to _continue_ in ccw direction of edge (circulate around
   //     E(e,1))
-  //   F  #F by 3 list of face indices
-  //   E  #E by 2 list of edge indices
   //   EMAP #F*3 list of indices into E, mapping each directed edge to unique
   //     unique edge in E
   //   EF  #E by 2 list of edge flaps, EF(e,0)=f means e=(i-->j) is the edge of
@@ -30,11 +28,10 @@ namespace igl
   //   EI  #E by 2 list of edge flap corners (see above).
   // Returns list of faces touched by circulation (in cyclically order).
   //   
+  // See also: edge_flaps
   IGL_INLINE std::vector<int> circulation(
     const int e,
     const bool ccw,
-    const Eigen::MatrixXi & F,
-    const Eigen::MatrixXi & E,
     const Eigen::VectorXi & EMAP,
     const Eigen::MatrixXi & EF,
     const Eigen::MatrixXi & EI);
@@ -42,8 +39,6 @@ namespace igl
   IGL_INLINE void circulation(
     const int e,
     const bool ccw,
-    const Eigen::MatrixXi & F,
-    const Eigen::MatrixXi & E,
     const Eigen::VectorXi & EMAP,
     const Eigen::MatrixXi & EF,
     const Eigen::MatrixXi & EI,

include/igl/collapse_edge.cpp

@@ -40,7 +40,7 @@ IGL_INLINE bool igl::collapse_edge(
   }
 
   // Important to grab neighbors of d before monkeying with edges
-  const std::vector<int> nV2Fd = circulation(e,!eflip,F,E,EMAP,EF,EI);
+  const std::vector<int> nV2Fd = circulation(e,!eflip,EMAP,EF,EI);
 
   // The following implementation strongly relies on s<d
   assert(s<d && "s should be less than d");
@@ -339,8 +339,8 @@ IGL_INLINE bool igl::collapse_edge(
   Q.erase(Q.begin());
   e = p.second;
   Qit[e] = Q.end();
-  std::vector<int> N  = circulation(e, true,F,E,EMAP,EF,EI);
-  std::vector<int> Nd = circulation(e,false,F,E,EMAP,EF,EI);
+  std::vector<int> N  = circulation(e, true,EMAP,EF,EI);
+  std::vector<int> Nd = circulation(e,false,EMAP,EF,EI);
   N.insert(N.begin(),Nd.begin(),Nd.end());
   bool collapsed = true;
   if(pre_collapse(V,F,E,EMAP,EF,EI,Q,Qit,C,e))

include/igl/copyleft/cgal/delaunay_triangulation.cpp

@@ -16,7 +16,7 @@ template<
   typename DerivedV,
   typename DerivedF>
 IGL_INLINE void igl::copyleft::cgal::delaunay_triangulation(
-    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedV>& V,
     Eigen::PlainObjectBase<DerivedF>& F)
 {
   typedef typename DerivedV::Scalar Scalar;
@@ -60,3 +60,8 @@ IGL_INLINE void igl::copyleft::cgal::delaunay_triangulation(
 //  };
 }
 
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::delaunay_triangulation<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::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/copyleft/cgal/delaunay_triangulation.h

@@ -32,7 +32,7 @@ namespace igl
         typename DerivedF
         >
       IGL_INLINE void delaunay_triangulation(
-          const Eigen::PlainObjectBase<DerivedV>& V,
+          const Eigen::MatrixBase<DerivedV>& V,
           Eigen::PlainObjectBase<DerivedF>& F);
     }
   }

include/igl/copyleft/cgal/hausdorff.cpp

@@ -41,4 +41,4 @@ IGL_INLINE void igl::copyleft::cgal::hausdorff(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 template void igl::copyleft::cgal::hausdorff<Eigen::Matrix<double, -1, -1, 0, -1, -1>, CGAL::Simple_cartesian<double>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Simple_cartesian<double>, CGAL::AABB_triangle_primitive<CGAL::Simple_cartesian<double>, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > >::iterator, CGAL::Boolean_tag<false> >, CGAL::Default> > const&, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > > const&, double&, double&);
-#endif
+#endif

include/igl/copyleft/cgal/incircle.cpp

@@ -37,3 +37,9 @@ IGL_INLINE short igl::copyleft::cgal::incircle(
       throw "Invalid incircle result";
   }
 }
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template short igl::copyleft::cgal::incircle<double>(double const*, double const*, double const*, double const*);
+#endif

include/igl/copyleft/cgal/orient2D.cpp

@@ -35,3 +35,9 @@ IGL_INLINE short igl::copyleft::cgal::orient2D(
       throw "Invalid orientation";
   }
 }
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template short igl::copyleft::cgal::orient2D<double>(double const*, double const*, double const*);
+#endif

include/igl/copyleft/progressive_hulls_cost_and_placement.cpp

@@ -32,8 +32,8 @@ IGL_INLINE void igl::copyleft::progressive_hulls_cost_and_placement(
 
   assert(V.cols() == 3 && "V.cols() should be 3");
   // Gather list of unique face neighbors
-  vector<int> Nall =  circulation(e, true,F,E,EMAP,EF,EI);
-  vector<int> Nother= circulation(e,false,F,E,EMAP,EF,EI);
+  vector<int> Nall =  circulation(e, true,EMAP,EF,EI);
+  vector<int> Nother= circulation(e,false,EMAP,EF,EI);
   Nall.insert(Nall.end(),Nother.begin(),Nother.end());
   vector<int> N;
   igl::unique(Nall,N);

include/igl/cotmatrix_entries.cpp

@@ -1,6 +1,7 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 
@@ -50,6 +51,7 @@ IGL_INLINE void igl::cotmatrix_entries(
       C.resize(m,3);
       for(int i = 0;i<m;i++)
       {
+        // Alec: I'm doubtful that using l2 here is actually improving numerics.
         C(i,0) = (l2(i,1) + l2(i,2) - l2(i,0))/dblA(i)/4.0;
         C(i,1) = (l2(i,2) + l2(i,0) - l2(i,1))/dblA(i)/4.0;
         C(i,2) = (l2(i,0) + l2(i,1) - l2(i,2))/dblA(i)/4.0;
@@ -97,14 +99,49 @@ IGL_INLINE void igl::cotmatrix_entries(
   }
 }
 
+template <typename Derivedl, typename DerivedC>
+IGL_INLINE void igl::cotmatrix_entries(
+  const Eigen::MatrixBase<Derivedl>& l,
+  Eigen::PlainObjectBase<DerivedC>& C)
+{
+  using namespace Eigen;
+  const int m = l.rows();
+  assert(l.cols() == 3 && "Only triangles accepted");
+  //Compute squared Edge lengths 
+  Matrix<typename DerivedC::Scalar,Dynamic,3> l2;
+  l2 = l.array().square();
+  // Alec: It's a little annoying that there's duplicate code here. The
+  // "extrinic" version above is first computing squared edge lengths, taking
+  // the square root and calling this. We can't have a cotmatrix_entries(l,l2,C)
+  // overload because it will confuse Eigen with the cotmatrix_entries(V,F,C)
+  // overload. In the end, I'd like to be convinced that using l2 directly above
+  // is actually better numerically (or significantly faster) than just calling
+  // edge_lengths and this cotmatrix_entries(l,C);
+  //
+  // double area
+  Matrix<typename DerivedC::Scalar,Dynamic,1> dblA;
+  doublearea(l,0.,dblA);
+  // cotangents and diagonal entries for element matrices
+  // correctly divided by 4 (alec 2010)
+  C.resize(m,3);
+  for(int i = 0;i<m;i++)
+  {
+    // Alec: I'm doubtful that using l2 here is actually improving numerics.
+    C(i,0) = (l2(i,1) + l2(i,2) - l2(i,0))/dblA(i)/4.0;
+    C(i,1) = (l2(i,2) + l2(i,0) - l2(i,1))/dblA(i)/4.0;
+    C(i,2) = (l2(i,0) + l2(i,1) - l2(i,2))/dblA(i)/4.0;
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
-template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, -1, 1, -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, 1, -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::cotmatrix_entries<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
+template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::cotmatrix_entries<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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, -1, 1, -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, 1, -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::cotmatrix_entries<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> >&);
-// generated by autoexplicit.sh
 template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/cotmatrix_entries.h

@@ -1,6 +1,7 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 
@@ -28,6 +29,17 @@ namespace igl
     const Eigen::MatrixBase<DerivedV>& V,
     const Eigen::MatrixBase<DerivedF>& F,
     Eigen::PlainObjectBase<DerivedC>& C);
+  // Intrinsic version.
+  //
+  // Inputs:
+  //   l  #F by 3 list of triangle edge lengths (see edge_lengths)
+  // Outputs:
+  //   C  #F by 3 list of 1/2*cotangents corresponding angles
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  template <typename Derivedl, typename DerivedC>
+  IGL_INLINE void cotmatrix_entries(
+    const Eigen::MatrixBase<Derivedl>& l,
+    Eigen::PlainObjectBase<DerivedC>& C);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/cotmatrix_intrinsic.cpp

@@ -0,0 +1,67 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 "cotmatrix_intrinsic.h"
+#include "cotmatrix_entries.h"
+#include <iostream>
+
+template <typename Derivedl, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::cotmatrix_intrinsic(
+  const Eigen::MatrixBase<Derivedl> & l, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<Scalar>& L)
+{
+  using namespace Eigen;
+  using namespace std;
+  // Cribbed from cotmatrix
+
+  const int nverts = F.maxCoeff()+1;
+  L.resize(nverts,nverts);
+  Matrix<int,Dynamic,2> edges;
+  int simplex_size = F.cols();
+  // 3 for triangles, 4 for tets
+  assert(simplex_size == 3);
+  // This is important! it could decrease the comptuation time by a factor of 2
+  // Laplacian for a closed 2d manifold mesh will have on average 7 entries per
+  // row
+  L.reserve(10*nverts);
+  edges.resize(3,2);
+  edges << 
+    1,2,
+    2,0,
+    0,1;
+  // Gather cotangents
+  Matrix<Scalar,Dynamic,Dynamic> C;
+  cotmatrix_entries(l,C);
+  
+  vector<Triplet<Scalar> > IJV;
+  IJV.reserve(F.rows()*edges.rows()*4);
+  // Loop over triangles
+  for(int i = 0; i < F.rows(); i++)
+  {
+    // loop over edges of element
+    for(int e = 0;e<edges.rows();e++)
+    {
+      int source = F(i,edges(e,0));
+      int dest = F(i,edges(e,1));
+      IJV.push_back(Triplet<Scalar>(source,dest,C(i,e)));
+      IJV.push_back(Triplet<Scalar>(dest,source,C(i,e)));
+      IJV.push_back(Triplet<Scalar>(source,source,-C(i,e)));
+      IJV.push_back(Triplet<Scalar>(dest,dest,-C(i,e)));
+    }
+  }
+  L.setFromTriplets(IJV.begin(),IJV.end());
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::cotmatrix_intrinsic<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cotmatrix_intrinsic<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cotmatrix_intrinsic<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cotmatrix_intrinsic<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cotmatrix_intrinsic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/cotmatrix_intrinsic.h

@@ -0,0 +1,40 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_COTMATRIX_INTRINSIC_H
+#define IGL_COTMATRIX_INTRINSIC_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl 
+{
+  // Constructs the cotangent stiffness matrix (discrete laplacian) for a given
+  // mesh with faces F and edge lengths l.
+  //
+  // Inputs:
+  //   l  #F by 3 list of (half-)edge lengths
+  //   F  #F by 3 list of face indices into some (not necessarily
+  //     determined/embedable) list of vertex positions V. It is assumed #V ==
+  //     F.maxCoeff()+1
+  // Outputs:
+  //   L  #V by #V sparse Laplacian matrix
+  //
+  // See also: cotmatrix, intrinsic_delaunay_cotmatrix
+  template <typename Derivedl, typename DerivedF, typename Scalar>
+  IGL_INLINE void cotmatrix_intrinsic(
+    const Eigen::MatrixBase<Derivedl> & l, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    Eigen::SparseMatrix<Scalar>& L);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cotmatrix_intrinsic.cpp"
+#endif
+
+#endif

include/igl/cumprod.cpp

@@ -0,0 +1,76 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 "cumprod.h"
+#include <numeric>
+#include <iostream>
+
+template <typename DerivedX, typename DerivedY>
+IGL_INLINE void igl::cumprod(
+  const Eigen::MatrixBase<DerivedX > & X,
+  const int dim,
+  Eigen::PlainObjectBase<DerivedY > & Y)
+{
+  using namespace Eigen;
+  using namespace std;
+  Y.resizeLike(X);
+  // get number of columns (or rows)
+  int num_outer = (dim == 1 ? X.cols() : X.rows() );
+  // get number of rows (or columns)
+  int num_inner = (dim == 1 ? X.rows() : X.cols() );
+  // This has been optimized so that dim = 1 or 2 is roughly the same cost.
+  // (Optimizations assume ColMajor order)
+  if(dim == 1)
+  {
+#pragma omp parallel for
+    for(int o = 0;o<num_outer;o++)
+    {
+      typename DerivedX::Scalar prod = 1;
+      for(int i = 0;i<num_inner;i++)
+      {
+        prod *= X(i,o);
+        Y(i,o) = prod;
+      }
+    }
+  }else
+  {
+    for(int i = 0;i<num_inner;i++)
+    {
+      // Notice that it is *not* OK to put this above the inner loop
+      // Though here it doesn't seem to pay off...
+//#pragma omp parallel for
+      for(int o = 0;o<num_outer;o++)
+      {
+        if(i == 0)
+        {
+          Y(o,i) = X(o,i);
+        }else
+        {
+          Y(o,i) = Y(o,i-1) * X(o,i);
+        }
+      }
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::cumprod<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> >&);
+// generated by autoexplicit.sh
+template void igl::cumprod<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> >&);
+template void igl::cumprod<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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::cumprod<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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::cumprod<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&);
+template void igl::cumprod<Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1>, Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1> >(Eigen::MatrixBase<Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<unsigned long, 2, 1, 0, 2, 1> >&);
+template void igl::cumprod<Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1>, Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<unsigned long, -1, 1, 0, -1, 1> >&);
+template void igl::cumprod<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#ifdef WIN32
+template void igl::cumprod<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::cumprod<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>> &);
+#endif
+#endif

include/igl/cumprod.h

@@ -0,0 +1,44 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_CUMPROD_H
+#define IGL_CUMPROD_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Computes a cumulative product of the columns of X, like matlab's `cumprod`.
+  //
+  // Templates:
+  //   DerivedX  Type of matrix X
+  //   DerivedY  Type of matrix Y
+  // Inputs:
+  //   X  m by n Matrix to be cumulatively multiplied.
+  //   dim  dimension to take cumulative product (1 or 2)
+  // Output:
+  //   Y  m by n Matrix containing cumulative product.
+  //
+  template <typename DerivedX, typename DerivedY>
+  IGL_INLINE void cumprod(
+    const Eigen::MatrixBase<DerivedX > & X,
+    const int dim,
+    Eigen::PlainObjectBase<DerivedY > & Y);
+  //template <typename DerivedX, typename DerivedY>
+  //IGL_INLINE void cumprod(
+  //  const Eigen::MatrixBase<DerivedX > & X,
+  //  Eigen::PlainObjectBase<DerivedY > & Y);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cumprod.cpp"
+#endif
+
+#endif
+
+

include/igl/cumsum.cpp

@@ -59,6 +59,10 @@ IGL_INLINE void igl::cumsum(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+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> >&);
+// generated by autoexplicit.sh
+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> >&);
 template void igl::cumsum<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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::cumsum<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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::cumsum<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&);

include/igl/delaunay_triangulation.cpp

@@ -10,6 +10,7 @@
 #include "flip_edge.h"
 #include "lexicographic_triangulation.h"
 #include "unique_edge_map.h"
+#include "is_delaunay.h"
 
 #include <vector>
 #include <sstream>
@@ -20,9 +21,9 @@ template<
   typename InCircle,
   typename DerivedF>
 IGL_INLINE void igl::delaunay_triangulation(
-    const Eigen::PlainObjectBase<DerivedV>& V,
-    Orient2D orient2D,
-    InCircle incircle,
+    const Eigen::MatrixBase<DerivedV>& V,
+    const Orient2D orient2D,
+    const InCircle incircle,
     Eigen::PlainObjectBase<DerivedF>& F)
 {
   assert(V.cols() == 2);
@@ -36,43 +37,18 @@ IGL_INLINE void igl::delaunay_triangulation(
   }
   assert(F.cols() == 3);
 
-  Eigen::MatrixXi E;
-  Eigen::MatrixXi uE;
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> MatrixX2I;
+  MatrixX2I E,uE;
   Eigen::VectorXi EMAP;
   std::vector<std::vector<Index> > uE2E;
   igl::unique_edge_map(F, E, uE, EMAP, uE2E);
 
-  auto is_delaunay = [&V,&F,&uE2E,num_faces,&incircle](size_t uei) {
-    auto& half_edges = uE2E[uei];
-    if (half_edges.size() != 2) {
-      throw "Cannot flip non-manifold or boundary edge";
-    }
-
-    const size_t f1 = half_edges[0] % num_faces;
-    const size_t f2 = half_edges[1] % num_faces;
-    const size_t c1 = half_edges[0] / num_faces;
-    const size_t c2 = half_edges[1] / num_faces;
-    assert(c1 < 3);
-    assert(c2 < 3);
-    assert(f1 != f2);
-    const size_t v1 = F(f1, (c1+1)%3);
-    const size_t v2 = F(f1, (c1+2)%3);
-    const size_t v4 = F(f1, c1);
-    const size_t v3 = F(f2, c2);
-    const Scalar p1[] = {V(v1, 0), V(v1, 1)};
-    const Scalar p2[] = {V(v2, 0), V(v2, 1)};
-    const Scalar p3[] = {V(v3, 0), V(v3, 1)};
-    const Scalar p4[] = {V(v4, 0), V(v4, 1)};
-    auto orientation = incircle(p1, p2, p4, p3);
-    return orientation <= 0;
-  };
-
   bool all_delaunay = false;
   while(!all_delaunay) {
     all_delaunay = true;
     for (size_t i=0; i<uE2E.size(); i++) {
       if (uE2E[i].size() == 2) {
-        if (!is_delaunay(i)) {
+        if (!is_delaunay(V,F,uE2E,incircle,i)) {
           all_delaunay = false;
           flip_edge(F, E, uE, EMAP, uE2E, i);
         }
@@ -82,8 +58,9 @@ IGL_INLINE void igl::delaunay_triangulation(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-template void igl::delaunay_triangulation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, short (*)(double const*, double const*, double const*), short (*)(double const*, double const*, double const*, double const*), Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, short (*)(double const*, double const*, double const*), short (*)(double const*, double const*, double const*, double const*), Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// Explicit template instantiation
+template void igl::delaunay_triangulation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, short (*)(double const*, double const*, double const*), short (*)(double const*, double const*, double const*, double const*), Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, short (*)(double const*, double const*, double const*), short (*)(double const*, double const*, double const*, double const*), Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #ifdef WIN32
-template void igl::delaunay_triangulation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, short (*)(double const * const, double const * const, double const * const), short(*)(double const * const, double const * const, double const * const, double const * const), Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const &, short(*)(double const * const, double const * const, double const * const), short(*)(double const * const, double const * const, double const * const, double const * const), Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > &);
+template void igl::delaunay_triangulation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, short (*)(double const * const, double const * const, double const * const), short(*)(double const * const, double const * const, double const * const, double const * const), Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const &, short(*)(double const * const, double const * const, double const * const), short(*)(double const * const, double const * const, double const * const, double const * const), Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > &);
+#endif
 #endif
-#endif

include/igl/delaunay_triangulation.h

@@ -37,9 +37,9 @@ namespace igl
     typename DerivedF
     >
   IGL_INLINE void delaunay_triangulation(
-      const Eigen::PlainObjectBase<DerivedV>& V,
-      Orient2D orient2D,
-      InCircle incircle,
+      const Eigen::MatrixBase<DerivedV>& V,
+      const Orient2D orient2D,
+      const InCircle incircle,
       Eigen::PlainObjectBase<DerivedF>& F);
 }
 

include/igl/doublearea.cpp

@@ -231,13 +231,12 @@ IGL_INLINE void igl::doublearea_quad(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
+template void igl::doublearea<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::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::doublearea<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> >&);
+template void igl::doublearea<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 1, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<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::doublearea<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -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<double, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
 template void igl::doublearea<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> >&);
 template void igl::doublearea<Eigen::Matrix<float, -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<float, -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> >&);
 template void igl::doublearea<Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 1, 0, 1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 1, 0, 1, 1> >&);

include/igl/edge_collapse_is_valid.cpp

@@ -49,7 +49,7 @@ IGL_INLINE bool igl::edge_collapse_is_valid(
       const Eigen::MatrixXi & EI) 
     {
       vector<int> N,uN;
-      vector<int> V2Fe = circulation(e, ccw,F,E,EMAP,EF,EI);
+      vector<int> V2Fe = circulation(e, ccw,EMAP,EF,EI);
       for(auto f : V2Fe)
       {
         N.push_back(F(f,0));

include/igl/edge_exists_near.cpp

@@ -0,0 +1,96 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 "edge_exists_near.h"
+
+template <
+  typename DeriveduE,
+  typename DerivedEMAP,
+  typename uE2EType,
+  typename Index>
+IGL_INLINE bool igl::edge_exists_near(
+  const Eigen::MatrixBase<DeriveduE> & uE,
+  const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+  const std::vector<std::vector< uE2EType> > & uE2E,
+  const Index & a,
+  const Index & b,
+  const Index & uei)
+{
+  std::vector<Index> face_queue;
+  face_queue.reserve(32);
+  std::vector<Index> pushed;
+  // 32 is faster than 8
+  pushed.reserve(32);
+  assert(a!=b);
+  // starting with the (2) faces incident on e, consider all faces
+  // incident on edges containing either a or b.
+  //
+  // face_queue  Queue containing faces incident on exactly one of a/b
+  // Using a vector seems mildly faster
+  assert(EMAP.size()%3 == 0);
+  const Index num_faces = EMAP.size()/3;
+  const Index f1 = uE2E[uei][0]%num_faces;
+  const Index f2 = uE2E[uei][1]%num_faces;
+  // map is faster than unordered_map here, and vector + brute force
+  // is_member check is even faster
+  face_queue.push_back(f1);
+  pushed.push_back(f1);
+  face_queue.push_back(f2);
+  pushed.push_back(f2);
+  while(!face_queue.empty())
+  {
+    const Index f = face_queue.back();
+    face_queue.pop_back();
+    // consider each edge of this face
+    for(int c = 0;c<3;c++)
+    {
+      // Unique edge id
+      const Index uec = EMAP(c*num_faces+f);
+      const Index s = uE(uec,0);
+      const Index d = uE(uec,1);
+      const bool ona = s == a || d == a;
+      const bool onb = s == b || d == b;
+      // Is this the edge we're looking for?
+      if(ona && onb)
+      {
+        return true;
+      }
+      // not incident on either?
+      if(!ona && !onb)
+      {
+        continue;
+      }
+      // loop over all incident half-edges
+      for(const auto & he : uE2E[uec])
+      {
+        // face of this he
+        const Index fhe = he%num_faces;
+        bool already_pushed = false;
+        for(const auto & fp : pushed)
+        {
+          if(fp == fhe)
+          {
+            already_pushed = true;
+            break;
+          }
+        }
+        if(!already_pushed)
+        {
+          pushed.push_back(fhe);
+          face_queue.push_back(fhe);
+        }
+      }
+    }
+  }
+  return false;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::edge_exists_near<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int, int>(Eigen::MatrixBase<Eigen::Matrix<int, -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> > > > const&, int const&, int const&, int const&);
+#endif 

include/igl/edge_exists_near.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_EDGE_EXISTS_NEAR_H
+#define IGL_EDGE_EXISTS_NEAR_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // Does edge (a,b) exist in the edges of all faces incident on
+  // existing unique edge uei.
+  //
+  // Inputs:
+  //   uE  #uE by 2 list of unique undirected edges
+  //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  //     undirected edge
+  //   uE2E  #uE list of lists of indices into E of coexisting edges
+  //   E  #F*3 by 2 list of half-edges
+  //   a  1st end-point of query edge
+  //   b  2nd end-point of query edge
+  //   uei  index into uE/uE2E of unique edge
+  // Returns true if edge exists near uei.
+  //
+  // See also: unique_edge_map
+  template <
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType,
+    typename Index>
+  IGL_INLINE bool edge_exists_near(
+    const Eigen::MatrixBase<DeriveduE> & uE,
+    const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+    const std::vector<std::vector< uE2EType> > & uE2E,
+    const Index & a,
+    const Index & b,
+    const Index & uei);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "edge_exists_near.cpp"
+#endif
+#endif
+

include/igl/edge_flaps.cpp

@@ -12,14 +12,14 @@
 
 IGL_INLINE void igl::edge_flaps(
   const Eigen::MatrixXi & F,
-  const Eigen::MatrixXi & E,
+  const Eigen::MatrixXi & uE,
   const Eigen::VectorXi & EMAP,
   Eigen::MatrixXi & EF,
   Eigen::MatrixXi & EI)
 {
   // Initialize to boundary value
-  EF.setConstant(E.rows(),2,-1);
-  EI.setConstant(E.rows(),2,-1);
+  EF.setConstant(uE.rows(),2,-1);
+  EI.setConstant(uE.rows(),2,-1);
   // loop over all faces
   for(int f = 0;f<F.rows();f++)
   {
@@ -29,13 +29,13 @@ IGL_INLINE void igl::edge_flaps(
       // get edge id
       const int e = EMAP(v*F.rows()+f);
       // See if this is left or right flap w.r.t. edge orientation
-      if( F(f,(v+1)%3) == E(e,0) && F(f,(v+2)%3) == E(e,1))
+      if( F(f,(v+1)%3) == uE(e,0) && F(f,(v+2)%3) == uE(e,1))
       {
         EF(e,0) = f;
         EI(e,0) = v;
       }else
       {
-        assert(F(f,(v+1)%3) == E(e,1) && F(f,(v+2)%3) == E(e,0));
+        assert(F(f,(v+1)%3) == uE(e,1) && F(f,(v+2)%3) == uE(e,0));
         EF(e,1) = f;
         EI(e,1) = v;
       }
@@ -45,16 +45,16 @@ IGL_INLINE void igl::edge_flaps(
 
 IGL_INLINE void igl::edge_flaps(
   const Eigen::MatrixXi & F,
-  Eigen::MatrixXi & E,
+  Eigen::MatrixXi & uE,
   Eigen::VectorXi & EMAP,
   Eigen::MatrixXi & EF,
   Eigen::MatrixXi & EI)
 {
   Eigen::MatrixXi allE;
   std::vector<std::vector<int> > uE2E;
-  igl::unique_edge_map(F,allE,E,EMAP,uE2E);
+  igl::unique_edge_map(F,allE,uE,EMAP,uE2E);
   // Const-ify to call overload
-  const auto & cE = E;
+  const auto & cuE = uE;
   const auto & cEMAP = EMAP;
-  return edge_flaps(F,cE,cEMAP,EF,EI);
+  return edge_flaps(F,cuE,cEMAP,EF,EI);
 }

include/igl/edge_flaps.h

@@ -16,15 +16,17 @@ namespace igl
   //
   // Inputs:
   //   F  #F by 3 list of face indices
-  //   E  #E by 2 list of edge indices into V.
-  //   EMAP #F*3 list of indices into E, mapping each directed edge to unique
-  //     unique edge in E
+  //   uE  #uE by 2 list of edge indices into V.
+  //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  //     unique edge in uE
   // Outputs:
   //   EF  #E by 2 list of edge flaps, EF(e,0)=f means e=(i-->j) is the edge of
   //     F(f,:) opposite the vth corner, where EI(e,0)=v. Similarly EF(e,1) "
   //     e=(j->i)
   //   EI  #E by 2 list of edge flap corners (see above).
   //
+  // See also: unique_edge_map
+  //
   // TODO: This seems to be a duplicate of edge_topology.h
   // igl::edge_topology(V,F,etEV,etFE,etEF);
   // igl::edge_flaps(F,efE,efEMAP,efEF,efEI);
@@ -34,14 +36,14 @@ namespace igl
   // all(efEMAP(sub2ind(size(F),repmat(1:size(F,1),3,1)',repmat([1 2 3],size(F,1),1))) == etFE(:,[2 3 1]))
   IGL_INLINE void edge_flaps(
     const Eigen::MatrixXi & F,
-    const Eigen::MatrixXi & E,
+    const Eigen::MatrixXi & uE,
     const Eigen::VectorXi & EMAP,
     Eigen::MatrixXi & EF,
     Eigen::MatrixXi & EI);
   // Only faces as input
   IGL_INLINE void edge_flaps(
     const Eigen::MatrixXi & F,
-    Eigen::MatrixXi & E,
+    Eigen::MatrixXi & uE,
     Eigen::VectorXi & EMAP,
     Eigen::MatrixXi & EF,
     Eigen::MatrixXi & EI);

include/igl/edge_lengths.cpp

@@ -22,6 +22,8 @@ IGL_INLINE void igl::edge_lengths(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(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, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
 template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
 // generated by autoexplicit.sh
 template void igl::edge_lengths<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> >&);

include/igl/embree/EmbreeIntersector.h

@@ -21,8 +21,8 @@
 #include <Eigen/Core>
 #include <Eigen/Geometry>
 
-#include <embree2/rtcore.h>
-#include <embree2/rtcore_ray.h>
+#include <embree3/rtcore.h>
+#include <embree3/rtcore_ray.h>
 #include <iostream>
 #include <vector>
 
@@ -177,7 +177,7 @@ namespace igl
       bool initialized;
 
       inline void createRay(
-        RTCRay& ray,
+        RTCRayHit& ray,
         const Eigen::RowVector3f& origin,
         const Eigen::RowVector3f& direction,
         float tnear,
@@ -199,29 +199,28 @@ namespace igl
   {
     // Keeps track of whether the **Global** Embree intersector has been
     // initialized. This should never been done at the global scope.
-    static bool EmbreeIntersector_inited = false;
+    static RTCDevice g_device = nullptr;
   }
 }
 
 inline void igl::embree::EmbreeIntersector::global_init()
 {
-  if(!EmbreeIntersector_inited)
+  if(!g_device)
   {
-    rtcInit();
-    if(rtcGetError() != RTC_NO_ERROR)
+    g_device = rtcNewDevice (NULL);
+    if(rtcGetDeviceError (g_device) != RTC_ERROR_NONE)
       std::cerr << "Embree: An error occurred while initializing embree core!" << std::endl;
 #ifdef IGL_VERBOSE
     else
       std::cerr << "Embree: core initialized." << std::endl;
 #endif
-    EmbreeIntersector_inited = true;
   }
 }
 
 inline void igl::embree::EmbreeIntersector::global_deinit()
 {
-  EmbreeIntersector_inited = false;
-  rtcExit();
+  rtcReleaseDevice (g_device);
+  g_device = nullptr;
 }
 
 inline igl::embree::EmbreeIntersector::EmbreeIntersector()
@@ -287,43 +286,49 @@ inline void igl::embree::EmbreeIntersector::init(
     return;
   }
 
+  RTCBuildQuality buildQuality = isStatic ? RTC_BUILD_QUALITY_HIGH : RTC_BUILD_QUALITY_MEDIUM;
+
   // create a scene
-  RTCSceneFlags flags = RTC_SCENE_ROBUST | RTC_SCENE_HIGH_QUALITY;
-  if(isStatic)
-    flags = flags | RTC_SCENE_STATIC;
-  scene = rtcNewScene(flags,RTC_INTERSECT1);
+  scene = rtcNewScene(g_device);
+  rtcSetSceneFlags(scene, RTC_SCENE_FLAG_ROBUST);
+  rtcSetSceneBuildQuality(scene, buildQuality);
 
   for(int g=0;g<(int)V.size();g++)
   {
     // create triangle mesh geometry in that scene
-    geomID = rtcNewTriangleMesh(scene,RTC_GEOMETRY_STATIC,F[g]->rows(),V[g]->rows(),1);
+    RTCGeometry geom_0 = rtcNewGeometry (g_device, RTC_GEOMETRY_TYPE_TRIANGLE);
+    rtcSetGeometryBuildQuality(geom_0,buildQuality);
+    rtcSetGeometryTimeStepCount(geom_0,1);
+    geomID = rtcAttachGeometry(scene,geom_0);
+    rtcReleaseGeometry(geom_0);
 
     // fill vertex buffer
-    vertices = (Vertex*)rtcMapBuffer(scene,geomID,RTC_VERTEX_BUFFER);
+    vertices = (Vertex*)rtcSetNewGeometryBuffer(geom_0,RTC_BUFFER_TYPE_VERTEX,0,RTC_FORMAT_FLOAT3,4*sizeof(float),V[g]->rows());
     for(int i=0;i<(int)V[g]->rows();i++)
     {
       vertices[i].x = (float)V[g]->coeff(i,0);
       vertices[i].y = (float)V[g]->coeff(i,1);
       vertices[i].z = (float)V[g]->coeff(i,2);
     }
-    rtcUnmapBuffer(scene,geomID,RTC_VERTEX_BUFFER);
+    
 
     // fill triangle buffer
-    triangles = (Triangle*) rtcMapBuffer(scene,geomID,RTC_INDEX_BUFFER);
+    triangles = (Triangle*) rtcSetNewGeometryBuffer(geom_0,RTC_BUFFER_TYPE_INDEX,0,RTC_FORMAT_UINT3,3*sizeof(int),F[g]->rows());
     for(int i=0;i<(int)F[g]->rows();i++)
     {
       triangles[i].v0 = (int)F[g]->coeff(i,0);
       triangles[i].v1 = (int)F[g]->coeff(i,1);
       triangles[i].v2 = (int)F[g]->coeff(i,2);
     }
-    rtcUnmapBuffer(scene,geomID,RTC_INDEX_BUFFER);
+    
 
-    rtcSetMask(scene,geomID,masks[g]);
+    rtcSetGeometryMask(geom_0,masks[g]);
+    rtcCommitGeometry(geom_0);
   }
 
-  rtcCommit(scene);
+  rtcCommitScene(scene);
 
-  if(rtcGetError() != RTC_NO_ERROR)
+  if(rtcGetDeviceError (g_device) != RTC_ERROR_NONE)
       std::cerr << "Embree: An error occurred while initializing the provided geometry!" << endl;
 #ifdef IGL_VERBOSE
   else
@@ -342,11 +347,11 @@ igl::embree::EmbreeIntersector
 
 void igl::embree::EmbreeIntersector::deinit()
 {
-  if(EmbreeIntersector_inited && scene)
+  if(g_device && scene)
   {
-    rtcDeleteScene(scene);
+    rtcReleaseScene(scene);
 
-    if(rtcGetError() != RTC_NO_ERROR)
+    if(rtcGetDeviceError (g_device) != RTC_ERROR_NONE)
     {
         std::cerr << "Embree: An error occurred while resetting!" << std::endl;
     }
@@ -367,23 +372,31 @@ inline bool igl::embree::EmbreeIntersector::intersectRay(
   float tfar,
   int mask) const
 {
-  RTCRay ray;
+  RTCRayHit ray; // EMBREE_FIXME: use RTCRay for occlusion rays
+  ray.ray.flags = 0;
   createRay(ray, origin,direction,tnear,tfar,mask);
 
   // shot ray
-  rtcIntersect(scene,ray);
+  {
+    RTCIntersectContext context;
+    rtcInitIntersectContext(&context);
+    rtcIntersect1(scene,&context,&ray);
+    ray.hit.Ng_x = -ray.hit.Ng_x; // EMBREE_FIXME: only correct for triangles,quads, and subdivision surfaces
+    ray.hit.Ng_y = -ray.hit.Ng_y;
+    ray.hit.Ng_z = -ray.hit.Ng_z;
+  }
 #ifdef IGL_VERBOSE
-  if(rtcGetError() != RTC_NO_ERROR)
+  if(rtcGetDeviceError (g_device) != RTC_ERROR_NONE)
       std::cerr << "Embree: An error occurred while resetting!" << std::endl;
 #endif
 
-  if((unsigned)ray.geomID != RTC_INVALID_GEOMETRY_ID)
+  if((unsigned)ray.hit.geomID != RTC_INVALID_GEOMETRY_ID)
   {
-    hit.id = ray.primID;
-    hit.gid = ray.geomID;
-    hit.u = ray.u;
-    hit.v = ray.v;
-    hit.t = ray.tfar;
+    hit.id = ray.hit.primID;
+    hit.gid = ray.hit.geomID;
+    hit.u = ray.hit.u;
+    hit.v = ray.hit.v;
+    hit.t = ray.ray.tfar;
     return true;
   }
 
@@ -419,6 +432,7 @@ inline bool igl::embree::EmbreeIntersector::intersectBeam(
   const float eps= 1e-5;
 
   Eigen::RowVector3f up(0,1,0);
+  if (direction.cross(up).norm() < eps) up = Eigen::RowVector3f(1,0,0);
   Eigen::RowVector3f offset = direction.cross(up).normalized();
 
   Eigen::Matrix3f rot = Eigen::AngleAxis<float>(2*3.14159265358979/samples,direction).toRotationMatrix();
@@ -462,22 +476,30 @@ igl::embree::EmbreeIntersector
   const double eps = FLOAT_EPS;
   double min_t = tnear;
   bool large_hits_warned = false;
-  RTCRay ray;
+  RTCRayHit ray; // EMBREE_FIXME: use RTCRay for occlusion rays
+  ray.ray.flags = 0;
   createRay(ray,origin,direction,tnear,tfar,mask);
 
   while(true)
   {
-    ray.tnear = min_t;
-    ray.tfar = tfar;
-    ray.geomID = RTC_INVALID_GEOMETRY_ID;
-    ray.primID = RTC_INVALID_GEOMETRY_ID;
-    ray.instID = RTC_INVALID_GEOMETRY_ID;
+    ray.ray.tnear = min_t;
+    ray.ray.tfar = tfar;
+    ray.hit.geomID = RTC_INVALID_GEOMETRY_ID;
+    ray.hit.primID = RTC_INVALID_GEOMETRY_ID;
+    ray.hit.instID[0] = RTC_INVALID_GEOMETRY_ID;
     num_rays++;
-    rtcIntersect(scene,ray);
-    if((unsigned)ray.geomID != RTC_INVALID_GEOMETRY_ID)
+    {
+      RTCIntersectContext context;
+      rtcInitIntersectContext(&context);
+      rtcIntersect1(scene,&context,&ray);
+      ray.hit.Ng_x = -ray.hit.Ng_x; // EMBREE_FIXME: only correct for triangles,quads, and subdivision surfaces
+      ray.hit.Ng_y = -ray.hit.Ng_y;
+      ray.hit.Ng_z = -ray.hit.Ng_z;
+    }
+    if((unsigned)ray.hit.geomID != RTC_INVALID_GEOMETRY_ID)
     {
       // Hit self again, progressively advance
-      if(ray.primID == last_id0 || ray.tfar <= min_t)
+      if(ray.hit.primID == last_id0 || ray.ray.tfar <= min_t)
       {
         // push min_t a bit more
         //double t_push = pow(2.0,self_hits-4)*(hit.t<eps?eps:hit.t);
@@ -492,23 +514,23 @@ igl::embree::EmbreeIntersector
       else
       {
         Hit hit;
-        hit.id = ray.primID;
-        hit.gid = ray.geomID;
-        hit.u = ray.u;
-        hit.v = ray.v;
-        hit.t = ray.tfar;
+        hit.id = ray.hit.primID;
+        hit.gid = ray.hit.geomID;
+        hit.u = ray.hit.u;
+        hit.v = ray.hit.v;
+        hit.t = ray.ray.tfar;
         hits.push_back(hit);
 #ifdef IGL_VERBOSE
         std::cerr<<"  t: "<<hit.t<<endl;
 #endif
         // Instead of moving origin, just change min_t. That way calculations
         // all use exactly same origin values
-        min_t = ray.tfar;
+        min_t = ray.ray.tfar;
 
         // reset t_scale
         self_hits = 0;
       }
-      last_id0 = ray.primID;
+      last_id0 = ray.hit.primID;
     }
     else
       break; // no more hits
@@ -544,18 +566,26 @@ inline bool
 igl::embree::EmbreeIntersector
 ::intersectSegment(const Eigen::RowVector3f& a, const Eigen::RowVector3f& ab, Hit &hit, int mask) const
 {
-  RTCRay ray;
+  RTCRayHit ray; // EMBREE_FIXME: use RTCRay for occlusion rays
+  ray.ray.flags = 0;
   createRay(ray,a,ab,0,1.0,mask);
 
-  rtcIntersect(scene,ray);
+  {
+    RTCIntersectContext context;
+    rtcInitIntersectContext(&context);
+    rtcIntersect1(scene,&context,&ray);
+    ray.hit.Ng_x = -ray.hit.Ng_x; // EMBREE_FIXME: only correct for triangles,quads, and subdivision surfaces
+    ray.hit.Ng_y = -ray.hit.Ng_y;
+    ray.hit.Ng_z = -ray.hit.Ng_z;
+  }
 
-  if((unsigned)ray.geomID != RTC_INVALID_GEOMETRY_ID)
+  if((unsigned)ray.hit.geomID != RTC_INVALID_GEOMETRY_ID)
   {
-    hit.id = ray.primID;
-    hit.gid = ray.geomID;
-    hit.u = ray.u;
-    hit.v = ray.v;
-    hit.t = ray.tfar;
+    hit.id = ray.hit.primID;
+    hit.gid = ray.hit.geomID;
+    hit.u = ray.hit.u;
+    hit.v = ray.hit.v;
+    hit.t = ray.ray.tfar;
     return true;
   }
 
@@ -564,21 +594,23 @@ igl::embree::EmbreeIntersector
 
 inline void
 igl::embree::EmbreeIntersector
-::createRay(RTCRay& ray, const Eigen::RowVector3f& origin, const Eigen::RowVector3f& direction, float tnear, float tfar, int mask) const
+::createRay(RTCRayHit& ray, const Eigen::RowVector3f& origin, const Eigen::RowVector3f& direction, float tnear, float tfar, int mask) const
 {
-  ray.org[0] = origin[0];
-  ray.org[1] = origin[1];
-  ray.org[2] = origin[2];
-  ray.dir[0] = direction[0];
-  ray.dir[1] = direction[1];
-  ray.dir[2] = direction[2];
-  ray.tnear = tnear;
-  ray.tfar = tfar;
-  ray.geomID = RTC_INVALID_GEOMETRY_ID;
-  ray.primID = RTC_INVALID_GEOMETRY_ID;
-  ray.instID = RTC_INVALID_GEOMETRY_ID;
-  ray.mask = mask;
-  ray.time = 0.0f;
+  ray.ray.org_x = origin[0];
+  ray.ray.org_y = origin[1];
+  ray.ray.org_z = origin[2];
+  ray.ray.dir_x = direction[0];
+  ray.ray.dir_y = direction[1];
+  ray.ray.dir_z = direction[2];
+  ray.ray.tnear = tnear;
+  ray.ray.tfar = tfar;
+  ray.ray.id = RTC_INVALID_GEOMETRY_ID;
+  ray.ray.mask = mask;
+  ray.ray.time = 0.0f;
+
+  ray.hit.geomID = RTC_INVALID_GEOMETRY_ID;
+  ray.hit.instID[0] = RTC_INVALID_GEOMETRY_ID;
+  ray.hit.primID = RTC_INVALID_GEOMETRY_ID;
 }
 
 #endif //EMBREE_INTERSECTOR_H

include/igl/flip_edge.cpp

@@ -25,6 +25,14 @@ IGL_INLINE void igl::flip_edge(
   typedef typename DerivedF::Scalar Index;
   const size_t num_faces = F.rows();
   assert(F.cols() == 3);
+  // Edge to flip [v1,v2] --> [v3,v4]
+  // Before:
+  // F(f1,:) = [v1,v2,v4] // in some cyclic order
+  // F(f2,:) = [v1,v3,v2] // in some cyclic order
+  // After: 
+  // F(f1,:) = [v1,v3,v4] // in *this* order 
+  // F(f2,:) = [v2,v4,v3] // in *this* order
+  //
   //          v1                 v1
   //          /|\                / \
   //         / | \              /f1 \
@@ -149,8 +157,12 @@ IGL_INLINE void igl::flip_edge(
 
 
 #ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::flip_edge<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, unsigned long);
 template void igl::flip_edge<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<int, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -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<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, unsigned long);
 #ifdef WIN32
 template void igl::flip_edge<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<int, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -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<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, unsigned __int64);
+template void igl::flip_edge<class Eigen::Matrix<int,-1,-1,0,-1,-1>,class Eigen::Matrix<int,-1,2,0,-1,2>,class Eigen::Matrix<int,-1,2,0,-1,2>,class Eigen::Matrix<int,-1,1,0,-1,1>,int>(class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,-1,0,-1,-1> > &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,2,0,-1,2> > &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,2,0,-1,2> > &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,1,0,-1,1> > &,class std::vector<class std::vector<int,class std::allocator<int> >,class std::allocator<class std::vector<int,class std::allocator<int> > > > &,unsigned __int64);
+#endif
 #endif
-#endif

include/igl/grad.cpp

@@ -15,10 +15,12 @@
 #include "doublearea.h"
 
 template <typename DerivedV, typename DerivedF>
-IGL_INLINE void grad_tet(const Eigen::PlainObjectBase<DerivedV>&V,
-                     const Eigen::PlainObjectBase<DerivedF>&T,
-                            Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
-                            bool uniform) {
+IGL_INLINE void grad_tet(
+  const Eigen::MatrixBase<DerivedV>&V,
+  const Eigen::MatrixBase<DerivedF>&T,
+  Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+  bool uniform) 
+{
   using namespace Eigen;
   assert(T.cols() == 4);
   const int n = V.rows(); int m = T.rows();
@@ -116,15 +118,22 @@ IGL_INLINE void grad_tet(const Eigen::PlainObjectBase<DerivedV>&V,
 }
 
 template <typename DerivedV, typename DerivedF>
-IGL_INLINE void grad_tri(const Eigen::PlainObjectBase<DerivedV>&V,
-                     const Eigen::PlainObjectBase<DerivedF>&F,
-                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
-                    bool uniform)
+IGL_INLINE void grad_tri(
+  const Eigen::MatrixBase<DerivedV>&V,
+  const Eigen::MatrixBase<DerivedF>&F,
+  Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+  bool uniform)
 {
+  // Number of faces
+  const int m = F.rows();
+  // Number of vertices
+  const int nv = V.rows();
+  // Number of dimensions
+  const int dims = V.cols();
   Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3>
-    eperp21(F.rows(),3), eperp13(F.rows(),3);
+    eperp21(m,3), eperp13(m,3);
 
-  for (int i=0;i<F.rows();++i)
+  for (int i=0;i<m;++i)
   {
     // renaming indices of vertices of triangles for convenience
     int i1 = F(i,0);
@@ -132,10 +141,14 @@ IGL_INLINE void grad_tri(const Eigen::PlainObjectBase<DerivedV>&V,
     int i3 = F(i,2);
 
     // #F x 3 matrices of triangle edge vectors, named after opposite vertices
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v32 = V.row(i3) - V.row(i2);
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v13 = V.row(i1) - V.row(i3);
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v21 = V.row(i2) - V.row(i1);
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> n = v32.cross(v13);
+    typedef Eigen::Matrix<typename DerivedV::Scalar, 1, 3> RowVector3S;
+    RowVector3S v32 = RowVector3S::Zero(1,3);
+    RowVector3S v13 = RowVector3S::Zero(1,3);
+    RowVector3S v21 = RowVector3S::Zero(1,3);
+    v32.head(V.cols()) = V.row(i3) - V.row(i2);
+    v13.head(V.cols()) = V.row(i1) - V.row(i3);
+    v21.head(V.cols()) = V.row(i2) - V.row(i1);
+    RowVector3S n = v32.cross(v13);
     // area of parallelogram is twice area of triangle
     // area of parallelogram is || v1 x v2 ||
     // This does correct l2 norm of rows, so that it contains #F list of twice
@@ -174,70 +187,46 @@ IGL_INLINE void grad_tri(const Eigen::PlainObjectBase<DerivedV>&V,
     eperp13.row(i) *= norm13 / dblA;
   }
 
-  std::vector<int> rs;
-  rs.reserve(F.rows()*4*3);
-  std::vector<int> cs;
-  cs.reserve(F.rows()*4*3);
-  std::vector<double> vs;
-  vs.reserve(F.rows()*4*3);
-
-  // row indices
-  for(int r=0;r<3;r++)
+  // create sparse gradient operator matrix
+  G.resize(dims*m,nv);
+  std::vector<Eigen::Triplet<typename DerivedV::Scalar> > Gijv;
+  Gijv.reserve(4*dims*m);
+  for(int f = 0;f<F.rows();f++)
   {
-    for(int j=0;j<4;j++)
+    for(int d = 0;d<dims;d++)
     {
-      for(int i=r*F.rows();i<(r+1)*F.rows();i++) rs.push_back(i);
+      Gijv.emplace_back(f+d*m,F(f,1), eperp13(f,d));
+      Gijv.emplace_back(f+d*m,F(f,0),-eperp13(f,d));
+      Gijv.emplace_back(f+d*m,F(f,2), eperp21(f,d));
+      Gijv.emplace_back(f+d*m,F(f,0),-eperp21(f,d));
     }
   }
-
-  // column indices
-  for(int r=0;r<3;r++)
-  {
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,1));
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,2));
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));
-  }
-
-  // values
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,2));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,2));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,2));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,2));
-
-  // create sparse gradient operator matrix
-  G.resize(3*F.rows(),V.rows());
-  std::vector<Eigen::Triplet<typename DerivedV::Scalar> > triplets;
-  for (int i=0;i<(int)vs.size();++i)
-  {
-    triplets.push_back(Eigen::Triplet<typename DerivedV::Scalar>(rs[i],cs[i],vs[i]));
-  }
-  G.setFromTriplets(triplets.begin(), triplets.end());
+  G.setFromTriplets(Gijv.begin(), Gijv.end());
 }
 
 template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
-                     const Eigen::PlainObjectBase<DerivedF>&F,
-                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
-                    bool uniform)
+IGL_INLINE void igl::grad(
+  const Eigen::MatrixBase<DerivedV>&V,
+  const Eigen::MatrixBase<DerivedF>&F,
+  Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+  bool uniform)
 {
   assert(F.cols() == 3 || F.cols() == 4);
-  if (F.cols() == 3)
-    return grad_tri(V,F,G,uniform);
-  if (F.cols() == 4)
-    return grad_tet(V,F,G,uniform);
+  switch(F.cols())
+  {
+    case 3:
+      return grad_tri(V,F,G,uniform);
+    case 4:
+      return grad_tet(V,F,G,uniform);
+    default:
+      assert(false);
+  }
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template void igl::grad<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> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 0, int>&, bool);
-template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&, bool);
+// generated by autoexplicit.sh
+template void igl::grad<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -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::SparseMatrix<Eigen::Matrix<double, -1, 2, 0, -1, 2>::Scalar, 0, int>&, bool);
+template void igl::grad<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::SparseMatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 0, int>&, bool);
+template void igl::grad<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&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&, bool);
 #endif

include/igl/grad.h

@@ -5,8 +5,8 @@
 // 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/.
-#ifndef IGL_GRAD_MAT_H
-#define IGL_GRAD_MAT_H
+#ifndef IGL_GRAD_H
+#define IGL_GRAD_H
 #include "igl_inline.h"
 
 #include <Eigen/Core>
@@ -33,11 +33,12 @@ namespace igl {
   // i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
   // 90 degrees
   //
-template <typename DerivedV, typename DerivedF>
-IGL_INLINE void grad(const Eigen::PlainObjectBase<DerivedV>&V,
-                     const Eigen::PlainObjectBase<DerivedF>&F,
-                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
-                    bool uniform = false);
+  template <typename DerivedV, typename DerivedF>
+  IGL_INLINE void grad(
+    const Eigen::MatrixBase<DerivedV>&V,
+    const Eigen::MatrixBase<DerivedF>&F,
+    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+    bool uniform = false);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "grad.cpp"

include/igl/grad_intrinsic.cpp

@@ -0,0 +1,69 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 "grad_intrinsic.h"
+#include "grad.h"
+
+template <typename Derivedl, typename DerivedF, typename Gtype>
+IGL_INLINE void igl::grad_intrinsic(
+  const Eigen::MatrixBase<Derivedl>&l,
+  const Eigen::MatrixBase<DerivedF>&F,
+  Eigen::SparseMatrix<Gtype> &G)
+{
+  assert(F.cols() ==3 && "Only triangles supported");
+  // number of vertices
+  const int n = F.maxCoeff()+1;
+  // number of faces
+  const int m = F.rows();
+  typedef Eigen::Matrix<Gtype,Eigen::Dynamic,2> MatrixX2S;
+  MatrixX2S V2 = MatrixX2S::Zero(3*m,2);
+  //     1=[x,y]
+  //     /\
+  // l3 /   \ l2
+  //   /      \
+  //  /         \
+  // 2-----------3
+  //       l1
+  //
+  // x = (l2²-l1²-l3²)/(-2*l1)
+  // y = sqrt(l3² - x²)
+  //
+  //
+  // Place 3rd vertex at [l(:,1) 0]
+  V2.block(2*m,0,m,1) = l.col(0);
+  // Place second vertex at [0 0]
+  // Place third vertex at [x y]
+  V2.block(0,0,m,1) = 
+    (l.col(1).cwiseAbs2()-l.col(0).cwiseAbs2()-l.col(2).cwiseAbs2()).array()/
+    (-2.*l.col(0)).array();
+  V2.block(0,1,m,1) = 
+    (l.col(2).cwiseAbs2() - V2.block(0,0,m,1).cwiseAbs2()).array().sqrt();
+  DerivedF F2(F.rows(),F.cols());
+  std::vector<Eigen::Triplet<Gtype> > Pijv;
+  Pijv.reserve(F.size());
+  for(int f = 0;f<m;f++)
+  {
+    for(int c = 0;c<F.cols();c++)
+    {
+      F2(f,c) = f+c*m;
+      Pijv.emplace_back(F2(f,c),F(f,c),1);
+    }
+  }
+  Eigen::SparseMatrix<Gtype> P(m*3,n);
+  P.setFromTriplets(Pijv.begin(),Pijv.end());
+  Eigen::SparseMatrix<Gtype> G2;
+  grad(V2,F2,G2);
+  G = G2*P;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::grad_intrinsic<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+// generated by autoexplicit.sh
+template void igl::grad_intrinsic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/grad_intrinsic.h

@@ -0,0 +1,35 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_GRAD_INTRINSIC_H
+#define IGL_GRAD_INTRINSIC_H
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+namespace igl {
+  // GRAD_INTRINSIC Construct an intrinsic gradient operator.
+  //
+  // Inputs:
+  //  l  #F by 3 list of edge lengths
+  //  F  #F by 3 list of triangle indices into some vertex list V
+  // Outputs:
+  //  G  #F*2 by #V gradient matrix: G=[Gx;Gy] where x runs along the 23 edge and
+  //    y runs in the counter-clockwise 90° rotation.
+  template <typename Derivedl, typename DerivedF, typename Gtype>
+  IGL_INLINE void grad_intrinsic(
+    const Eigen::MatrixBase<Derivedl>&l,
+    const Eigen::MatrixBase<DerivedF>&F,
+    Eigen::SparseMatrix<Gtype> &G);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "grad_intrinsic.cpp"
+#endif
+
+#endif
+

include/igl/grid.cpp

@@ -6,6 +6,7 @@
 // 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 "grid.h"
+#include <cassert>
 
 template <
   typename Derivedres,
@@ -16,23 +17,30 @@ IGL_INLINE void igl::grid(
 {
   using namespace Eigen;
   typedef typename DerivedGV::Scalar Scalar;
-  GV.resize(res(0)*res(1)*res(2),3);
-  for(int zi = 0;zi<res(2);zi++)
+  GV.resize(res.array().prod(),res.size());
+  const auto lerp = 
+    [&res](const Scalar di, const int d)->Scalar{return di/(Scalar)(res(d)-1);};
+  int gi = 0;
+  Derivedres sub;
+  sub.resizeLike(res);
+  sub.setConstant(0);
+  for(int gi = 0;gi<GV.rows();gi++)
   {
-    const auto lerp = 
-      [&](const Scalar di, const int d)->Scalar{return di/(Scalar)(res(d)-1);};
-    const Scalar z = lerp(zi,2);
-    for(int yi = 0;yi<res(1);yi++)
+    // omg, I'm implementing addition...
+    for(int c = 0;c<res.size()-1;c++)
     {
-      const Scalar y = lerp(yi,1);
-      for(int xi = 0;xi<res(0);xi++)
+      if(sub(c)>=res(c))
       {
-        const Scalar x = lerp(xi,0);
-        const int gi = xi+res(0)*(yi + res(1)*zi);
-        GV.row(gi) = 
-          Eigen::Matrix<Scalar,1,3>(x,y,z);
+        sub(c) = 0;
+        // roll over
+        sub(c+1)++;
       }
     }
+    for(int c = 0;c<res.size();c++)
+    {
+      GV(gi,c) = lerp(sub(c),c);
+    }
+    sub(0)++;
   }
 }
 
@@ -40,6 +48,8 @@ IGL_INLINE void igl::grid(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::grid<Eigen::Matrix<int, 2, 1, 0, 2, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, 2, 1, 0, 2, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::grid<Eigen::Matrix<int, 1, 3, 1, 1, 3>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(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::grid<Eigen::Matrix<int, 3, 1, 0, 3, 1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<int, 3, 1, 0, 3, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);

include/igl/grid.h

@@ -1,6 +1,7 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 
@@ -15,9 +16,9 @@ namespace igl
   // `igl::marching_cubes`
   //
   // Inputs:
-  //   res  3-long list of number of vertices along the x y and z dimensions
+  //   res  #res list of number of vertices along each dimension
   // Outputs:
-  //   GV  res(0)*res(1)*res(2) by 3 list of mesh vertex positions.
+  //   GV  res.array().prod() by #res list of mesh vertex positions.
   //   
   template <
     typename Derivedres,

include/igl/heat_geodesics.cpp

@@ -0,0 +1,159 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 "heat_geodesics.h"
+#include "grad.h"
+#include "doublearea.h"
+#include "cotmatrix.h"
+#include "intrinsic_delaunay_cotmatrix.h"
+#include "massmatrix.h"
+#include "massmatrix_intrinsic.h"
+#include "grad_intrinsic.h"
+#include "boundary_facets.h"
+#include "unique.h"
+#include "slice.h"
+#include "avg_edge_length.h"
+
+
+template < typename DerivedV, typename DerivedF, typename Scalar >
+IGL_INLINE bool igl::heat_geodesics_precompute(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  HeatGeodesicsData<Scalar> & data)
+{
+  // default t value
+  const Scalar h = avg_edge_length(V,F);
+  const Scalar t = h*h;
+  return heat_geodesics_precompute(V,F,t,data);
+}
+
+template < typename DerivedV, typename DerivedF, typename Scalar >
+IGL_INLINE bool igl::heat_geodesics_precompute(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Scalar t,
+  HeatGeodesicsData<Scalar> & data)
+{
+  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1> VectorXS;
+  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> MatrixXS;
+  Eigen::SparseMatrix<Scalar> L,M;
+  Eigen::Matrix<Scalar,Eigen::Dynamic,3> l_intrinsic;
+  DerivedF F_intrinsic;
+  VectorXS dblA;
+  if(data.use_intrinsic_delaunay)
+  {
+    igl::intrinsic_delaunay_cotmatrix(V,F,L,l_intrinsic,F_intrinsic);
+    igl::massmatrix_intrinsic(l_intrinsic,F_intrinsic,MASSMATRIX_TYPE_DEFAULT,M);
+    igl::doublearea(l_intrinsic,0,dblA);
+    igl::grad_intrinsic(l_intrinsic,F_intrinsic,data.Grad);
+  }else
+  {
+    igl::cotmatrix(V,F,L);
+    igl::massmatrix(V,F,MASSMATRIX_TYPE_DEFAULT,M);
+    igl::doublearea(V,F,dblA);
+    igl::grad(V,F,data.Grad);
+  }
+  // div
+  assert(F.cols() == 3 && "Only triangles are supported");
+  // number of gradient components
+  data.ng = data.Grad.rows() / F.rows();
+  assert(data.ng == 3 || data.ng == 2);
+  data.Div = -0.25*data.Grad.transpose()*dblA.replicate(data.ng,1).asDiagonal();
+
+  Eigen::SparseMatrix<Scalar> Q = M - t*L;
+  Eigen::MatrixXi O;
+  igl::boundary_facets(F,O);
+  igl::unique(O,data.b);
+  {
+    Eigen::SparseMatrix<Scalar> _;
+    if(!igl::min_quad_with_fixed_precompute(
+      Q,Eigen::VectorXi(),_,true,data.Neumann))
+    {
+      return false;
+    }
+    // Only need if there's a boundary
+    if(data.b.size()>0)
+    {
+      if(!igl::min_quad_with_fixed_precompute(Q,data.b,_,true,data.Dirichlet))
+      {
+        return false;
+      }
+    }
+    const Eigen::SparseMatrix<double> Aeq = M.diagonal().transpose().sparseView();
+    if(!igl::min_quad_with_fixed_precompute(
+      (-L*0.5).eval(),Eigen::VectorXi(),Aeq,true,data.Poisson))
+    {
+      return false;
+    }
+  }
+  return true;
+}
+
+template < typename Scalar, typename Derivedgamma, typename DerivedD>
+IGL_INLINE void igl::heat_geodesics_solve(
+  const HeatGeodesicsData<Scalar> & data,
+  const Eigen::MatrixBase<Derivedgamma> & gamma,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  // number of mesh vertices
+  const int n = data.Grad.cols();
+  // Set up delta at gamma
+  DerivedD u0 = DerivedD::Zero(n,1);
+  for(int g = 0;g<gamma.size();g++)
+  {
+    u0(gamma(g)) = 1;
+  }
+  // Neumann solution
+  DerivedD u;
+  igl::min_quad_with_fixed_solve(
+    data.Neumann,u0,DerivedD(),DerivedD(),u);
+  if(data.b.size()>0)
+  {
+    // Average Dirichelt and Neumann solutions
+    DerivedD uD;
+    igl::min_quad_with_fixed_solve(
+      data.Dirichlet,u0,DerivedD::Zero(data.b.size()),DerivedD(),uD);
+    u += uD;
+    u *= 0.5;
+  }
+  DerivedD grad_u = data.Grad*u;
+  const int m = data.Grad.rows()/data.ng;
+  for(int i = 0;i<m;i++)
+  {
+    Scalar norm = 0;
+    for(int d = 0;d<data.ng;d++)
+    {
+      norm += grad_u(d*m+i)*grad_u(d*m+i);
+    }
+    norm = sqrt(norm);
+    if(norm == 0)
+    {
+      for(int d = 0;d<data.ng;d++) { grad_u(d*m+i) = 0; }
+    }else
+    {
+      for(int d = 0;d<data.ng;d++) { grad_u(d*m+i) /= norm; }
+    }
+  }
+  const DerivedD div_X = -data.Div*grad_u;
+  const DerivedD Beq = (DerivedD(1,1)<<0).finished();
+  igl::min_quad_with_fixed_solve(data.Poisson,-2.0*div_X,DerivedD(),Beq,D);
+  DerivedD Dgamma;
+  igl::slice(D,gamma,Dgamma);
+  D.array() -= Dgamma.mean();
+  if(D.mean() < 0)
+  {
+    D = -D;
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::heat_geodesics_solve<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::HeatGeodesicsData<double> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template bool igl::heat_geodesics_precompute<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, double, igl::HeatGeodesicsData<double>&);
+template bool igl::heat_geodesics_precompute<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::HeatGeodesicsData<double>&);
+#endif

include/igl/heat_geodesics.h

@@ -0,0 +1,69 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_HEAT_GEODESICS_H
+#define IGL_HEAT_GEODESICS_H
+#include "igl_inline.h"
+#include "min_quad_with_fixed.h"
+#include <Eigen/Sparse>
+#include <Eigen/Sparse>
+namespace igl
+{
+  template <typename Scalar>
+  struct HeatGeodesicsData
+  {
+    // Gradient and Divergence operators
+    Eigen::SparseMatrix<Scalar> Grad,Div;
+    // Number of gradient components
+    int ng;
+    // List of boundary vertex indices
+    Eigen::VectorXi b;
+    // Solvers for Dirichet, Neumann problems
+    min_quad_with_fixed_data<Scalar> Dirichlet,Neumann,Poisson;
+    bool use_intrinsic_delaunay = false;
+  };
+  // Precompute factorized solvers for computing a fast approximation of
+  // geodesic distances on a mesh (V,F). [Crane et al. 2013]
+  //
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by 3 list of mesh face indices into V
+  // Outputs:
+  //   data  precomputation data (see heat_geodesics_solve)
+  template < typename DerivedV, typename DerivedF, typename Scalar >
+  IGL_INLINE bool heat_geodesics_precompute(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    HeatGeodesicsData<Scalar> & data);
+  // Inputs:
+  //   t  "heat" parameter (smaller --> more accurate, less stable)
+  template < typename DerivedV, typename DerivedF, typename Scalar >
+  IGL_INLINE bool heat_geodesics_precompute(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Scalar t,
+    HeatGeodesicsData<Scalar> & data);
+  // Compute fast approximate geodesic distances using precomputed data from a
+  // set of selected source vertices (gamma)
+  //
+  // Inputs: 
+  //   data  precomputation data (see heat_geodesics_precompute)
+  //   gamma  #gamma list of indices into V of source vertices
+  // Outputs:
+  //   D  #V list of distances to gamma 
+  template < typename Scalar, typename Derivedgamma, typename DerivedD>
+  IGL_INLINE void heat_geodesics_solve(
+    const HeatGeodesicsData<Scalar> & data,
+    const Eigen::MatrixBase<Derivedgamma> & gamma,
+    Eigen::PlainObjectBase<DerivedD> & D);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "heat_geodesics.cpp"
+#endif
+
+#endif

include/igl/intrinsic_delaunay_cotmatrix.cpp

@@ -0,0 +1,54 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 "intrinsic_delaunay_cotmatrix.h"
+#include "edge_lengths.h"
+#include "intrinsic_delaunay_triangulation.h"
+#include "cotmatrix_intrinsic.h"
+#include <cassert>
+
+template <typename DerivedV, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::intrinsic_delaunay_cotmatrix(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<Scalar>& L)
+{
+  Eigen::Matrix<Scalar, Eigen::Dynamic, 3> l_intrinsic;
+  DerivedF F_intrinsic;
+  return igl::intrinsic_delaunay_cotmatrix(V,F,L,l_intrinsic,F_intrinsic);
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename Scalar,
+  typename Derivedl_intrinsic,
+  typename DerivedF_intrinsic>
+IGL_INLINE void igl::intrinsic_delaunay_cotmatrix(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<Scalar>& L,
+  Eigen::PlainObjectBase<Derivedl_intrinsic> & l_intrinsic,
+  Eigen::PlainObjectBase<DerivedF_intrinsic> & F_intrinsic)
+{
+  assert(F.cols() == 3 && "Only triangles are supported");
+  Eigen::Matrix<Scalar, Eigen::Dynamic, 3> l;
+  igl::edge_lengths(V,F,l);
+  igl::intrinsic_delaunay_triangulation(l,F,l_intrinsic,F_intrinsic);
+  igl::cotmatrix_intrinsic(l_intrinsic,F_intrinsic,L);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, -1, 3, 0, -1, 3>, 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::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, 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::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+#endif

include/igl/intrinsic_delaunay_cotmatrix.h

@@ -0,0 +1,51 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_INTRINSIC_DELAUNAY_COTMATRIX_H
+#define IGL_INTRINSIC_DELAUNAY_COTMATRIX_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // INTRINSIC_DELAUNAY_COTMATRIX Computes the discrete cotangent Laplacian of a
+  // mesh after converting it into its intrinsic Delaunay triangulation (see,
+  // e.g., [Fisher et al. 2007].
+  //
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by 3 list of mesh elements (triangles or tetrahedra)
+  // Outputs: 
+  //   L  #V by #V cotangent matrix, each row i corresponding to V(i,:)
+  //   l_intrinsic  #F by 3 list of intrinsic edge-lengths used to compute L
+  //   F_intrinsic  #F by 3 list of intrinsic face indices used to compute L
+  //
+  // See also: intrinsic_delaunay_triangulation, cotmatrix, cotmatrix_intrinsic
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename Scalar,
+    typename Derivedl_intrinsic,
+    typename DerivedF_intrinsic>
+  IGL_INLINE void intrinsic_delaunay_cotmatrix(
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    Eigen::SparseMatrix<Scalar>& L,
+    Eigen::PlainObjectBase<Derivedl_intrinsic> & l_intrinsic,
+    Eigen::PlainObjectBase<DerivedF_intrinsic> & F_intrinsic);
+  template <typename DerivedV, typename DerivedF, typename Scalar>
+  IGL_INLINE void intrinsic_delaunay_cotmatrix(
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    Eigen::SparseMatrix<Scalar>& L);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "intrinsic_delaunay_cotmatrix.cpp"
+#endif
+
+#endif

include/igl/intrinsic_delaunay_triangulation.cpp

@@ -0,0 +1,199 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 "intrinsic_delaunay_triangulation.h"
+#include "is_intrinsic_delaunay.h"
+#include "tan_half_angle.h"
+#include "unique_edge_map.h"
+#include "flip_edge.h"
+#include "EPS.h"
+#include "matlab_format.h"
+#include <iostream>
+#include <queue>
+#include <map>
+
+template <
+  typename Derivedl_in,
+  typename DerivedF_in,
+  typename Derivedl,
+  typename DerivedF>
+IGL_INLINE void igl::intrinsic_delaunay_triangulation(
+  const Eigen::MatrixBase<Derivedl_in> & l_in,
+  const Eigen::MatrixBase<DerivedF_in> & F_in,
+  Eigen::PlainObjectBase<Derivedl> & l,
+  Eigen::PlainObjectBase<DerivedF> & F)
+{
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> MatrixX2I;
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> VectorXI;
+  MatrixX2I E,uE;
+  VectorXI EMAP;
+  std::vector<std::vector<typename DerivedF::Scalar> > uE2E;
+  return intrinsic_delaunay_triangulation(l_in,F_in,l,F,E,uE,EMAP,uE2E);
+}
+
+template <
+  typename Derivedl_in,
+  typename DerivedF_in,
+  typename Derivedl,
+  typename DerivedF,
+  typename DerivedE,
+  typename DeriveduE,
+  typename DerivedEMAP,
+  typename uE2EType>
+IGL_INLINE void igl::intrinsic_delaunay_triangulation(
+  const Eigen::MatrixBase<Derivedl_in> & l_in,
+  const Eigen::MatrixBase<DerivedF_in> & F_in,
+  Eigen::PlainObjectBase<Derivedl> & l,
+  Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DeriveduE> & uE,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+  std::vector<std::vector<uE2EType> > & uE2E)
+{
+  igl::unique_edge_map(F_in, E, uE, EMAP, uE2E);
+  // We're going to work in place
+  l = l_in;
+  F = F_in;
+  typedef typename DerivedF::Scalar Index;
+  typedef typename Derivedl::Scalar Scalar;
+  const Index num_faces = F.rows();
+
+  // Vector is faster than queue...
+  std::vector<Index> Q;
+  Q.reserve(uE2E.size());
+  for (size_t uei=0; uei<uE2E.size(); uei++) 
+  {
+    Q.push_back(uei);
+  }
+  // I tried using a "delaunay_since = iter" flag to avoid duplicates, but there
+  // was no speed up.
+
+  while(!Q.empty())
+  {
+#ifdef IGL_INTRINSIC_DELAUNAY_TRIANGULATION_DEBUG
+    // Expensive sanity check
+    {
+      Eigen::Matrix<bool,Eigen::Dynamic,1> inQ(uE2E.size(),1);
+      inQ.setConstant(false);
+      for(const auto uei : Q)
+      {
+        inQ(uei) = true;
+      }
+      for (Index uei=0; uei<uE2E.size(); uei++) 
+      {
+        if(!inQ(uei) && !is_intrinsic_delaunay(l,uE2E,num_faces,uei))
+        {
+          std::cout<<"  "<<uei<<" might never be fixed!"<<std::endl;
+        }
+      }
+    }
+#endif
+
+    const Index uei = Q.back();
+    Q.pop_back();
+    if (uE2E[uei].size() == 2) 
+    {
+      if(!is_intrinsic_delaunay(l,uE2E,num_faces,uei)) 
+      {
+        // update l just before flipping edge
+        //      .        //
+        //     /|\       //
+        //   a/ | \d     //
+        //   /  e  \     //
+        //  /   |   \    //
+        // .----|-f--.   //
+        //  \   |   /    //
+        //   \  |  /     //
+        //   b\α|δ/c     //
+        //     \|/       //
+        //      .        //
+        // Annotated from flip_edge:
+        // Edge to flip [v1,v2] --> [v3,v4]
+        // Before:
+        // F(f1,:) = [v1,v2,v4] // in some cyclic order
+        // F(f2,:) = [v1,v3,v2] // in some cyclic order
+        // After: 
+        // F(f1,:) = [v1,v3,v4] // in *this* order 
+        // F(f2,:) = [v2,v4,v3] // in *this* order
+        //
+        //          v1                 v1
+        //          /|\                / \
+        //        c/ | \b            c/f1 \b
+        //     v3 /f2|f1\ v4  =>  v3 /__f__\ v4
+        //        \  e  /            \ f2  /
+        //        d\ | /a            d\   /a
+        //          \|/                \ /
+        //          v2                 v2
+        //
+        // Compute intrinsic length of oppposite edge
+        assert(uE2E[uei].size() == 2 && "edge should have 2 incident faces");
+        const Index f1 = uE2E[uei][0]%num_faces;
+        const Index f2 = uE2E[uei][1]%num_faces;
+        const Index c1 = uE2E[uei][0]/num_faces;
+        const Index c2 = uE2E[uei][1]/num_faces;
+        assert(c1 < 3);
+        assert(c2 < 3);
+        assert(f1 != f2);
+        const Index v1 = F(f1, (c1+1)%3);
+        const Index v2 = F(f1, (c1+2)%3);
+        const Index v4 = F(f1, c1);
+        const Index v3 = F(f2, c2);
+        assert(F(f2, (c2+2)%3) == v1);
+        assert(F(f2, (c2+1)%3) == v2);
+        assert( std::abs(l(f1,c1)-l(f2,c2)) < igl::EPS<Scalar>() );
+        const Scalar e = l(f1,c1);
+        const Scalar a = l(f1,(c1+1)%3);
+        const Scalar b = l(f1,(c1+2)%3);
+        const Scalar c = l(f2,(c2+1)%3);
+        const Scalar d = l(f2,(c2+2)%3);
+        // tan(α/2)
+        const Scalar tan_a_2= tan_half_angle(a,b,e);
+        // tan(δ/2)
+        const Scalar tan_d_2 = tan_half_angle(d,e,c);
+        // tan((α+δ)/2)
+        const Scalar tan_a_d_2 = (tan_a_2 + tan_d_2)/(1.0-tan_a_2*tan_d_2);
+        // cos(α+δ)
+        const Scalar cos_a_d = 
+          (1.0 - tan_a_d_2*tan_a_d_2)/(1.0+tan_a_d_2*tan_a_d_2);
+        const Scalar f = sqrt(b*b + c*c - 2.0*b*c*cos_a_d);
+        l(f1,0) = f;
+        l(f1,1) = b;
+        l(f1,2) = c;
+        l(f2,0) = f;
+        l(f2,1) = d;
+        l(f2,2) = a;
+        // Important to grab these indices _before_ calling flip_edges (they
+        // will be correct after)
+        const size_t e_24 = f1 + ((c1 + 1) % 3) * num_faces;
+        const size_t e_41 = f1 + ((c1 + 2) % 3) * num_faces;
+        const size_t e_13 = f2 + ((c2 + 1) % 3) * num_faces;
+        const size_t e_32 = f2 + ((c2 + 2) % 3) * num_faces;
+        const size_t ue_24 = EMAP(e_24);
+        const size_t ue_41 = EMAP(e_41);
+        const size_t ue_13 = EMAP(e_13);
+        const size_t ue_32 = EMAP(e_32);
+        flip_edge(F, E, uE, EMAP, uE2E, uei);
+        Q.push_back(ue_24);
+        Q.push_back(ue_41);
+        Q.push_back(ue_13);
+        Q.push_back(ue_32);
+      }
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_triangulation<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_triangulation<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_triangulation<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<int, -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>, int>(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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -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<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_triangulation<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<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/intrinsic_delaunay_triangulation.h

@@ -0,0 +1,79 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_INTRINSIC_DELAUNAY_TRIANGULATION_H
+#define IGL_INTRINSIC_DELAUNAY_TRIANGULATION_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // INTRINSIC_DELAUNAY_TRIANGULATION Flip edges _intrinsically_ until all are
+  // "intrinsic Delaunay". See "An algorithm for the construction of intrinsic
+  // delaunay triangulations with applications to digital geometry processing"
+  // [Fisher et al. 2007].
+  //
+  // Inputs:
+  //   l_in  #F_in by 3 list of edge lengths (see edge_lengths)
+  //   F_in  #F_in by 3 list of face indices into some unspecified vertex list V
+  // Outputs:
+  //   l  #F by 3 list of edge lengths
+  //   F  #F by 3 list of new face indices. Note: Combinatorially F may contain
+  //     non-manifold edges, duplicate faces and self-loops (e.g., an edge [1,1]
+  //     or a face [1,1,1]). However, the *intrinsic geometry* is still
+  //     well-defined and correct. See [Fisher et al. 2007] Figure 3 and 2nd to
+  //     last paragraph of 1st page. Since F may be "non-eddge-manifold" in the
+  //     usual combinatorial sense, it may be useful to call the more verbose
+  //     overload below if disentangling edges will be necessary later on.
+  //     Calling unique_edge_map on this F will give a _different_ result than
+  //     those outputs.
+  //
+  // See also: is_intrinsic_delaunay
+  template <
+    typename Derivedl_in,
+    typename DerivedF_in,
+    typename Derivedl,
+    typename DerivedF>
+  IGL_INLINE void intrinsic_delaunay_triangulation(
+    const Eigen::MatrixBase<Derivedl_in> & l_in,
+    const Eigen::MatrixBase<DerivedF_in> & F_in,
+    Eigen::PlainObjectBase<Derivedl> & l,
+    Eigen::PlainObjectBase<DerivedF> & F);
+  // Outputs:
+  //   E  #F*3 by 2 list of all directed edges, such that E.row(f+#F*c) is the
+  //     edge opposite F(f,c)
+  //   uE  #uE by 2 list of unique undirected edges
+  //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  //     undirected edge
+  //   uE2E  #uE list of lists of indices into E of coexisting edges
+  //
+  // See also: unique_edge_map
+  template <
+    typename Derivedl_in,
+    typename DerivedF_in,
+    typename Derivedl,
+    typename DerivedF,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType>
+  IGL_INLINE void intrinsic_delaunay_triangulation(
+    const Eigen::MatrixBase<Derivedl_in> & l_in,
+    const Eigen::MatrixBase<DerivedF_in> & F_in,
+    Eigen::PlainObjectBase<Derivedl> & l,
+    Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DeriveduE> & uE,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    std::vector<std::vector<uE2EType> > & uE2E);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "intrinsic_delaunay_triangulation.cpp"
+#endif
+
+#endif

include/igl/is_border_vertex.cpp

@@ -10,14 +10,13 @@
 
 #include "triangle_triangle_adjacency.h"
 
-template <typename DerivedV, typename DerivedF>
+template <typename DerivedF>
 IGL_INLINE std::vector<bool> igl::is_border_vertex(
-    const Eigen::PlainObjectBase<DerivedV> &V,
-    const Eigen::PlainObjectBase<DerivedF> &F)
+  const Eigen::MatrixBase<DerivedF> &F)
 {
   DerivedF FF;
   igl::triangle_triangle_adjacency(F,FF);
-  std::vector<bool> ret(V.rows());
+  std::vector<bool> ret(F.maxCoeff()+1);
   for(unsigned i=0; i<ret.size();++i)
     ret[i] = false;
 
@@ -31,9 +30,17 @@ IGL_INLINE std::vector<bool> igl::is_border_vertex(
   return ret;
 }
 
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE std::vector<bool> igl::is_border_vertex(
+  const Eigen::MatrixBase<DerivedV> &/*V*/,
+  const Eigen::MatrixBase<DerivedF> &F)
+{
+  return igl::is_border_vertex(F);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
-template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<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> > const&);
-template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<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> > const&);
+template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<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&);
+template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<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&);
+template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<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&);
 #endif

include/igl/is_border_vertex.h

@@ -22,13 +22,16 @@ namespace igl
   //   F  #F by 3 list of triangle indices
   // Returns #V vector of bools revealing whether vertices are on boundary
   //
-  // Known Bugs: - does not depend on V
-  //             - assumes mesh is edge manifold
+  // Known Bugs: - assumes mesh is edge manifold
   // 
+  template <typename DerivedF>
+  IGL_INLINE std::vector<bool> is_border_vertex(
+   const Eigen::MatrixBase<DerivedF> &F);
+  // Deprecated:
   template <typename DerivedV, typename DerivedF>
   IGL_INLINE std::vector<bool> is_border_vertex(
-   const Eigen::PlainObjectBase<DerivedV> &V,
-   const Eigen::PlainObjectBase<DerivedF> &F);
+   const Eigen::MatrixBase<DerivedV> &V,
+   const Eigen::MatrixBase<DerivedF> &F);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/is_delaunay.cpp

@@ -0,0 +1,105 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 "is_delaunay.h"
+#include "unique_edge_map.h"
+#include <cassert>
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedD>
+IGL_INLINE void igl::is_delaunay(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  typedef typename DerivedV::Scalar Scalar;
+  // Should use Shewchuk's predicates instead.
+  const auto float_incircle = [](
+    const Scalar pa[2],
+    const Scalar pb[2],
+    const Scalar pc[2],
+    const Scalar pd[2])->short
+  {
+    // I acknowledge that I am cating to double
+    const Eigen::Matrix3d A = (Eigen::Matrix3d(3,3)<<
+      pa[0]-pd[0], pa[1]-pd[1],(pa[0]-pd[0])*(pa[0]-pd[0])+(pa[1]-pd[1])*(pa[1]-pd[1]),
+      pb[0]-pd[0], pb[1]-pd[1],(pb[0]-pd[0])*(pb[0]-pd[0])+(pb[1]-pd[1])*(pb[1]-pd[1]),
+      pc[0]-pd[0], pc[1]-pd[1],(pc[0]-pd[0])*(pc[0]-pd[0])+(pc[1]-pd[1])*(pc[1]-pd[1])
+      ).finished();
+    const Scalar detA = A.determinant();
+    return (Scalar(0) < detA) - (detA < Scalar(0));
+  };
+
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> MatrixX2I;
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> VectorXI;
+  MatrixX2I E,uE;
+  VectorXI EMAP;
+  std::vector<std::vector<typename DerivedF::Scalar> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+  const int num_faces = F.rows();
+  D.setConstant(F.rows(),F.cols(),false);
+  // loop over all unique edges
+  for(int ue = 0;ue < uE2E.size(); ue++)
+  {
+    const bool ue_is_d = is_delaunay(V,F,uE2E,float_incircle,ue);
+    // Set for all instances
+    for(int e = 0;e<uE2E[ue].size();e++)
+    {
+      D( uE2E[ue][e]%num_faces, uE2E[ue][e]/num_faces) = ue_is_d;
+    }
+  }
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename uE2EType,
+  typename InCircle,
+  typename ueiType>
+IGL_INLINE bool igl::is_delaunay(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const std::vector<std::vector<uE2EType> > & uE2E,
+  const InCircle incircle,
+  const ueiType uei)
+{
+  if(uE2E[uei].size() == 1) return true;
+  if(uE2E[uei].size() > 2) return false;
+  const int num_faces = F.rows();
+  typedef typename DerivedV::Scalar Scalar;
+  const auto& half_edges = uE2E[uei];
+  assert((half_edges.size() == 2)  && "uE2E[uei].size() should be 2");
+  const size_t f1 = half_edges[0] % num_faces;
+  const size_t f2 = half_edges[1] % num_faces;
+  const size_t c1 = half_edges[0] / num_faces;
+  const size_t c2 = half_edges[1] / num_faces;
+  assert(c1 < 3);
+  assert(c2 < 3);
+  assert(f1 != f2);
+  const size_t v1 = F(f1, (c1+1)%3);
+  const size_t v2 = F(f1, (c1+2)%3);
+  const size_t v4 = F(f1, c1);
+  const size_t v3 = F(f2, c2);
+  const Scalar p1[] = {V(v1, 0), V(v1, 1)};
+  const Scalar p2[] = {V(v2, 0), V(v2, 1)};
+  const Scalar p3[] = {V(v3, 0), V(v3, 1)};
+  const Scalar p4[] = {V(v4, 0), V(v4, 1)};
+  auto orientation = incircle(p1, p2, p4, p3);
+  return orientation <= 0;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::is_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -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::PlainObjectBase<Eigen::Matrix<bool, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::is_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -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<bool, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template bool igl::is_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, short (*)(double const*, double const*, double const*, double const*), unsigned long>(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> > > > const&, short (*)(double const*, double const*, double const*, double const*), unsigned long);
+#endif

include/igl/is_delaunay.h

@@ -0,0 +1,64 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_IS_DELAUNAY_H
+#define IGL_IS_DELAUNAY_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // IS_DELAUNAY Determine if each edge in the mesh (V,F) is Delaunay.
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   F  #F by 3 list of triangles indices
+  // Outputs:
+  //   D  #F by 3 list of bools revealing whether edges corresponding 23 31 12
+  //     are locally Delaunay. Boundary edges are by definition Delaunay.
+  //     Non-Manifold edges are by definition not Delaunay.
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedD>
+  IGL_INLINE void is_delaunay(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedD> & D);
+  // Determine whether a single edge is Delaunay using a provided (extrinsic) incirle
+  // test.
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   F  #F by 3 list of triangles indices
+  //   uE2E  #uE list of lists of indices into E of coexisting edges (see
+  //     unique_edge_map)
+  //   incircle  A functor such that incircle(pa, pb, pc, pd) returns
+  //               1 if pd is on the positive size of circumcirle of (pa,pb,pc)
+  //              -1 if pd is on the positive size of circumcirle of (pa,pb,pc)
+  //               0 if pd is cocircular with pa, pb, pc.
+  //               (see delaunay_triangulation)
+  //   uei  index into uE2E of edge to check
+  // Returns true iff edge is Delaunay
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename uE2EType,
+    typename InCircle,
+    typename ueiType>
+  IGL_INLINE bool is_delaunay(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const std::vector<std::vector<uE2EType> > & uE2E,
+    const InCircle incircle,
+    const ueiType uei);
+
+}
+#ifndef IGL_STATIC_LIBRARY
+#include "is_delaunay.cpp"
+#endif
+#endif

include/igl/is_intrinsic_delaunay.cpp

@@ -0,0 +1,147 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 "is_intrinsic_delaunay.h"
+#include "unique_edge_map.h"
+#include "tan_half_angle.h"
+#include "EPS.h"
+#include "matlab_format.h"
+#include "cotmatrix_entries.h"
+#include <iostream>
+#include <cassert>
+template <
+  typename Derivedl,
+  typename DerivedF,
+  typename DerivedD>
+IGL_INLINE void igl::is_intrinsic_delaunay(
+  const Eigen::MatrixBase<Derivedl> & l,
+  const Eigen::MatrixBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> MatrixX2I;
+  typedef Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> VectorXI;
+  MatrixX2I E,uE;
+  VectorXI EMAP;
+  std::vector<std::vector<typename DerivedF::Scalar> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+  return is_intrinsic_delaunay(l,F,uE2E,D);
+}
+
+template <
+  typename Derivedl,
+  typename DerivedF,
+  typename uE2EType,
+  typename DerivedD>
+IGL_INLINE void igl::is_intrinsic_delaunay(
+  const Eigen::MatrixBase<Derivedl> & l,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const std::vector<std::vector<uE2EType> > & uE2E,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  const int num_faces = F.rows();
+  D.setConstant(F.rows(),F.cols(),false);
+  // loop over all unique edges
+  for(int ue = 0;ue < uE2E.size(); ue++)
+  {
+    const bool ue_is_d = is_intrinsic_delaunay(l,uE2E,num_faces,ue);
+    // Set for all instances
+    for(int e = 0;e<uE2E[ue].size();e++)
+    {
+      D( uE2E[ue][e]%num_faces, uE2E[ue][e]/num_faces) = ue_is_d;
+    }
+  }
+};
+
+template <
+  typename Derivedl,
+  typename uE2EType,
+  typename Index>
+IGL_INLINE bool igl::is_intrinsic_delaunay(
+  const Eigen::MatrixBase<Derivedl> & l,
+  const std::vector<std::vector<uE2EType> > & uE2E,
+  const Index num_faces,
+  const Index uei)
+{
+  if(uE2E[uei].size() == 1) return true;
+  if(uE2E[uei].size() > 2) return false;
+  typedef typename Derivedl::Scalar Scalar;
+
+  const auto cot_alpha = [](
+    const Scalar & a,
+    const Scalar & b,
+    const Scalar & c)->Scalar
+  {
+    // Fisher 2007
+    const Scalar t = tan_half_angle(a,b,c);
+    return (1.0-t*t)/(2*t);
+  };
+
+  //      2        //
+  //     /|\       //
+  //   a/ | \d     //
+  //   /  e  \     //
+  //  /   |   \    //
+  //0.α---|-f-β.3  //
+  //  \   |   /    //
+  //   \  |  /     //
+  //   b\ | /c     //
+  //     \|/       //
+  //      .        //
+  //      1
+
+  // Fisher 2007
+  assert(uE2E[uei].size() == 2 && "edge should have 2 incident faces");
+  const Index he1 = uE2E[uei][0];
+  const Index he2 = uE2E[uei][1];
+  const Index f1 = he1%num_faces;
+  const Index c1 = he1/num_faces;
+  const Index f2 = he2%num_faces;
+  const Index c2 = he2/num_faces;
+  assert( std::abs(l(f1,c1)-l(f2,c2)) < igl::EPS<Scalar>());
+  const Scalar e = l(f1,c1);
+  const Scalar a = l(f1,(c1+1)%3);
+  const Scalar b = l(f1,(c1+2)%3);
+  const Scalar c = l(f2,(c2+1)%3);
+  const Scalar d = l(f2,(c2+2)%3);
+  const Scalar w = cot_alpha(e,a,b) + cot_alpha(e,c,d);
+
+  //// Test
+  //{
+  //  Eigen::MatrixXd l(2,3);
+  //  l<<e,a,b,
+  //     d,e,c;
+  //  //Eigen::MatrixXi F(2,3);
+  //  //F<<0,1,2,
+  //  //   1,3,2;
+  //  Eigen::MatrixXd C;
+  //  cotmatrix_entries(l,C);
+  //  if(std::abs(w-(2.0*(C(0,0)+C(1,1))))>1e-10)
+  //  {
+  //    std::cout<<matlab_format(l,"l")<<std::endl;
+  //    std::cout<<w<<std::endl;
+  //    std::cout<<(2.0*(C(0,0)+C(1,1)))<<std::endl;
+  //    std::cout<<w-(2.0*(C(0,0)+C(1,1)))<<std::endl;
+  //    std::cout<<std::endl;
+  //  }
+  //}
+
+  return w >= 0;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, 3, 0, -1, 3>, int, int>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, int, int);
+// generated by autoexplicit.sh
+template bool igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, int, int>(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&, int, int);
+// generated by autoexplicit.sh
+template void igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, Eigen::Matrix<bool, -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&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -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<bool, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::is_intrinsic_delaunay<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -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::PlainObjectBase<Eigen::Matrix<bool, -1, 3, 0, -1, 3> >&);
+#endif

include/igl/is_intrinsic_delaunay.h

@@ -0,0 +1,69 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_IS_INTRINSIC_DELAUNAY_H
+#define IGL_IS_INTRINSIC_DELAUNAY_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // IS_INTRINSIC_DELAUNAY Determine if each edge in the mesh (V,F) is Delaunay.
+  //
+  // Inputs:
+  //   l  #l by dim list of edge lengths
+  //   F  #F by 3 list of triangles indices
+  // Outputs:
+  //   D  #F by 3 list of bools revealing whether edges corresponding 23 31 12
+  //     are locally Delaunay. Boundary edges are by definition Delaunay.
+  //     Non-Manifold edges are by definition not Delaunay.
+  template <
+    typename Derivedl,
+    typename DerivedF,
+    typename DerivedD>
+  IGL_INLINE void is_intrinsic_delaunay(
+    const Eigen::MatrixBase<Derivedl> & l,
+    const Eigen::MatrixBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedD> & D);
+  // Inputs:
+  //   uE2E  #uE list of lists mapping unique edges to (half-)edges
+  template <
+    typename Derivedl,
+    typename DerivedF,
+    typename uE2EType,
+    typename DerivedD>
+  IGL_INLINE void is_intrinsic_delaunay(
+    const Eigen::MatrixBase<Derivedl> & l,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const std::vector<std::vector<uE2EType> > & uE2E,
+    Eigen::PlainObjectBase<DerivedD> & D);
+  // Determine whether a single edge is Delaunay using a provided (extrinsic) incirle
+  // test.
+  //
+  // Inputs:
+  //   l  #l by dim list of edge lengths
+  //   uE2E  #uE list of lists of indices into E of coexisting edges (see
+  //     unique_edge_map)
+  //   num_faces  number of faces (==#F)
+  //   uei  index into uE2E of edge to check
+  // Returns true iff edge is Delaunay
+  template <
+    typename Derivedl,
+    typename uE2EType,
+    typename Index>
+  IGL_INLINE bool is_intrinsic_delaunay(
+    const Eigen::MatrixBase<Derivedl> & l,
+    const std::vector<std::vector<uE2EType> > & uE2E,
+    const Index num_faces,
+    const Index uei);
+
+}
+#ifndef IGL_STATIC_LIBRARY
+#include "is_intrinsic_delaunay.cpp"
+#endif
+#endif
+

include/igl/lexicographic_triangulation.cpp

@@ -19,7 +19,7 @@ template<
   typename DerivedF
   >
 IGL_INLINE void igl::lexicographic_triangulation(
-    const Eigen::PlainObjectBase<DerivedP>& P,
+    const Eigen::MatrixBase<DerivedP>& P,
     Orient2D orient2D,
     Eigen::PlainObjectBase<DerivedF>& F)
 {
@@ -128,5 +128,5 @@ IGL_INLINE void igl::lexicographic_triangulation(
 
 
 #ifdef IGL_STATIC_LIBRARY
-template void igl::lexicographic_triangulation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, short (*)(double const*, double const*, double const*), Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, short (*)(double const*, double const*, double const*), Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-#endif
+template void igl::lexicographic_triangulation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, short (*)(double const*, double const*, double const*), Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, short (*)(double const*, double const*, double const*), Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/lexicographic_triangulation.h

@@ -33,7 +33,7 @@ namespace igl
     typename DerivedF
     >
   IGL_INLINE void lexicographic_triangulation(
-      const Eigen::PlainObjectBase<DerivedP>& P,
+      const Eigen::MatrixBase<DerivedP>& P,
       Orient2D orient2D,
       Eigen::PlainObjectBase<DerivedF>& F);
 }

include/igl/massmatrix.cpp

@@ -6,6 +6,8 @@
 // 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 "massmatrix.h"
+#include "massmatrix_intrinsic.h"
+#include "edge_lengths.h"
 #include "normalize_row_sums.h"
 #include "sparse.h"
 #include "doublearea.h"
@@ -37,92 +39,18 @@ IGL_INLINE void igl::massmatrix(
   // Not yet supported
   assert(type!=MASSMATRIX_TYPE_FULL);
 
-  Matrix<int,Dynamic,1> MI;
-  Matrix<int,Dynamic,1> MJ;
-  Matrix<Scalar,Dynamic,1> MV;
   if(simplex_size == 3)
   {
     // Triangles
     // edge lengths numbered same as opposite vertices
-    Matrix<Scalar,Dynamic,3> l(m,3);
-    // loop over faces
-    for(int i = 0;i<m;i++)
-    {
-      l(i,0) = (V.row(F(i,1))-V.row(F(i,2))).norm();
-      l(i,1) = (V.row(F(i,2))-V.row(F(i,0))).norm();
-      l(i,2) = (V.row(F(i,0))-V.row(F(i,1))).norm();
-    }
-    Matrix<Scalar,Dynamic,1> dblA;
-    doublearea(l,0.,dblA);
-
-    switch(eff_type)
-    {
-      case MASSMATRIX_TYPE_BARYCENTRIC:
-        // diagonal entries for each face corner
-        MI.resize(m*3,1); MJ.resize(m*3,1); MV.resize(m*3,1);
-        MI.block(0*m,0,m,1) = F.col(0);
-        MI.block(1*m,0,m,1) = F.col(1);
-        MI.block(2*m,0,m,1) = F.col(2);
-        MJ = MI;
-        repmat(dblA,3,1,MV);
-        MV.array() /= 6.0;
-        break;
-      case MASSMATRIX_TYPE_VORONOI:
-        {
-          // diagonal entries for each face corner
-          // http://www.alecjacobson.com/weblog/?p=874
-          MI.resize(m*3,1); MJ.resize(m*3,1); MV.resize(m*3,1);
-          MI.block(0*m,0,m,1) = F.col(0);
-          MI.block(1*m,0,m,1) = F.col(1);
-          MI.block(2*m,0,m,1) = F.col(2);
-          MJ = MI;
-
-          // Holy shit this needs to be cleaned up and optimized
-          Matrix<Scalar,Dynamic,3> cosines(m,3);
-          cosines.col(0) = 
-            (l.col(2).array().pow(2)+l.col(1).array().pow(2)-l.col(0).array().pow(2))/(l.col(1).array()*l.col(2).array()*2.0);
-          cosines.col(1) = 
-            (l.col(0).array().pow(2)+l.col(2).array().pow(2)-l.col(1).array().pow(2))/(l.col(2).array()*l.col(0).array()*2.0);
-          cosines.col(2) = 
-            (l.col(1).array().pow(2)+l.col(0).array().pow(2)-l.col(2).array().pow(2))/(l.col(0).array()*l.col(1).array()*2.0);
-          Matrix<Scalar,Dynamic,3> barycentric = cosines.array() * l.array();
-          normalize_row_sums(barycentric,barycentric);
-          Matrix<Scalar,Dynamic,3> partial = barycentric;
-          partial.col(0).array() *= dblA.array() * 0.5;
-          partial.col(1).array() *= dblA.array() * 0.5;
-          partial.col(2).array() *= dblA.array() * 0.5;
-          Matrix<Scalar,Dynamic,3> quads(partial.rows(),partial.cols());
-          quads.col(0) = (partial.col(1)+partial.col(2))*0.5;
-          quads.col(1) = (partial.col(2)+partial.col(0))*0.5;
-          quads.col(2) = (partial.col(0)+partial.col(1))*0.5;
-
-          quads.col(0) = (cosines.col(0).array()<0).select( 0.25*dblA,quads.col(0));
-          quads.col(1) = (cosines.col(0).array()<0).select(0.125*dblA,quads.col(1));
-          quads.col(2) = (cosines.col(0).array()<0).select(0.125*dblA,quads.col(2));
-
-          quads.col(0) = (cosines.col(1).array()<0).select(0.125*dblA,quads.col(0));
-          quads.col(1) = (cosines.col(1).array()<0).select(0.25*dblA,quads.col(1));
-          quads.col(2) = (cosines.col(1).array()<0).select(0.125*dblA,quads.col(2));
-
-          quads.col(0) = (cosines.col(2).array()<0).select(0.125*dblA,quads.col(0));
-          quads.col(1) = (cosines.col(2).array()<0).select(0.125*dblA,quads.col(1));
-          quads.col(2) = (cosines.col(2).array()<0).select( 0.25*dblA,quads.col(2));
-
-          MV.block(0*m,0,m,1) = quads.col(0);
-          MV.block(1*m,0,m,1) = quads.col(1);
-          MV.block(2*m,0,m,1) = quads.col(2);
-          
-          break;
-        }
-      case MASSMATRIX_TYPE_FULL:
-        assert(false && "Implementation incomplete");
-        break;
-      default:
-        assert(false && "Unknown Mass matrix eff_type");
-    }
-
+    Matrix<Scalar,Dynamic,3> l;
+    igl::edge_lengths(V,F,l);
+    return massmatrix_intrinsic(l,F,type,M);
   }else if(simplex_size == 4)
   {
+    Matrix<int,Dynamic,1> MI;
+    Matrix<int,Dynamic,1> MJ;
+    Matrix<Scalar,Dynamic,1> MV;
     assert(V.cols() == 3);
     assert(eff_type == MASSMATRIX_TYPE_BARYCENTRIC);
     MI.resize(m*4,1); MJ.resize(m*4,1); MV.resize(m*4,1);
@@ -145,12 +73,12 @@ IGL_INLINE void igl::massmatrix(
       MV(i+2*m) = v/4.0;
       MV(i+3*m) = v/4.0;
     }
+    sparse(MI,MJ,MV,n,n,M);
   }else
   {
     // Unsupported simplex size
     assert(false && "Unsupported simplex size");
   }
-  sparse(MI,MJ,MV,n,n,M);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/massmatrix.h

@@ -5,11 +5,10 @@
 // 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/.
-#ifndef IGL_MASSMATRIX_TYPE_H
-#define IGL_MASSMATRIX_TYPE_H
+#ifndef IGL_MASSMATRIX_H
+#define IGL_MASSMATRIX_H
 #include "igl_inline.h"
 
-#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 

include/igl/massmatrix_intrinsic.cpp

@@ -0,0 +1,128 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@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 "massmatrix_intrinsic.h"
+#include "edge_lengths.h"
+#include "normalize_row_sums.h"
+#include "sparse.h"
+#include "doublearea.h"
+#include "repmat.h"
+#include <Eigen/Geometry>
+#include <iostream>
+#include <cassert>
+
+template <typename Derivedl, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::massmatrix_intrinsic(
+  const Eigen::MatrixBase<Derivedl> & l, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  const MassMatrixType type,
+  Eigen::SparseMatrix<Scalar>& M)
+{
+  const int n = F.maxCoeff()+1;
+  return massmatrix_intrinsic(l,F,type,n,M);
+}
+
+template <typename Derivedl, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::massmatrix_intrinsic(
+  const Eigen::MatrixBase<Derivedl> & l, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  const MassMatrixType type,
+  const int n,
+  Eigen::SparseMatrix<Scalar>& M)
+{
+  using namespace Eigen;
+  using namespace std;
+  MassMatrixType eff_type = type;
+  const int m = F.rows();
+  const int simplex_size = F.cols();
+  // Use voronoi of for triangles by default, otherwise barycentric
+  if(type == MASSMATRIX_TYPE_DEFAULT)
+  {
+    eff_type = (simplex_size == 3?MASSMATRIX_TYPE_VORONOI:MASSMATRIX_TYPE_BARYCENTRIC);
+  }
+  assert(F.cols() == 3 && "only triangles supported");
+  Matrix<Scalar,Dynamic,1> dblA;
+  doublearea(l,0.,dblA);
+  Matrix<int,Dynamic,1> MI;
+  Matrix<int,Dynamic,1> MJ;
+  Matrix<Scalar,Dynamic,1> MV;
+
+  switch(eff_type)
+  {
+    case MASSMATRIX_TYPE_BARYCENTRIC:
+      // diagonal entries for each face corner
+      MI.resize(m*3,1); MJ.resize(m*3,1); MV.resize(m*3,1);
+      MI.block(0*m,0,m,1) = F.col(0);
+      MI.block(1*m,0,m,1) = F.col(1);
+      MI.block(2*m,0,m,1) = F.col(2);
+      MJ = MI;
+      repmat(dblA,3,1,MV);
+      MV.array() /= 6.0;
+      break;
+    case MASSMATRIX_TYPE_VORONOI:
+      {
+        // diagonal entries for each face corner
+        // http://www.alecjacobson.com/weblog/?p=874
+        MI.resize(m*3,1); MJ.resize(m*3,1); MV.resize(m*3,1);
+        MI.block(0*m,0,m,1) = F.col(0);
+        MI.block(1*m,0,m,1) = F.col(1);
+        MI.block(2*m,0,m,1) = F.col(2);
+        MJ = MI;
+
+        // Holy shit this needs to be cleaned up and optimized
+        Matrix<Scalar,Dynamic,3> cosines(m,3);
+        cosines.col(0) = 
+          (l.col(2).array().pow(2)+l.col(1).array().pow(2)-l.col(0).array().pow(2))/(l.col(1).array()*l.col(2).array()*2.0);
+        cosines.col(1) = 
+          (l.col(0).array().pow(2)+l.col(2).array().pow(2)-l.col(1).array().pow(2))/(l.col(2).array()*l.col(0).array()*2.0);
+        cosines.col(2) = 
+          (l.col(1).array().pow(2)+l.col(0).array().pow(2)-l.col(2).array().pow(2))/(l.col(0).array()*l.col(1).array()*2.0);
+        Matrix<Scalar,Dynamic,3> barycentric = cosines.array() * l.array();
+        normalize_row_sums(barycentric,barycentric);
+        Matrix<Scalar,Dynamic,3> partial = barycentric;
+        partial.col(0).array() *= dblA.array() * 0.5;
+        partial.col(1).array() *= dblA.array() * 0.5;
+        partial.col(2).array() *= dblA.array() * 0.5;
+        Matrix<Scalar,Dynamic,3> quads(partial.rows(),partial.cols());
+        quads.col(0) = (partial.col(1)+partial.col(2))*0.5;
+        quads.col(1) = (partial.col(2)+partial.col(0))*0.5;
+        quads.col(2) = (partial.col(0)+partial.col(1))*0.5;
+
+        quads.col(0) = (cosines.col(0).array()<0).select( 0.25*dblA,quads.col(0));
+        quads.col(1) = (cosines.col(0).array()<0).select(0.125*dblA,quads.col(1));
+        quads.col(2) = (cosines.col(0).array()<0).select(0.125*dblA,quads.col(2));
+
+        quads.col(0) = (cosines.col(1).array()<0).select(0.125*dblA,quads.col(0));
+        quads.col(1) = (cosines.col(1).array()<0).select(0.25*dblA,quads.col(1));
+        quads.col(2) = (cosines.col(1).array()<0).select(0.125*dblA,quads.col(2));
+
+        quads.col(0) = (cosines.col(2).array()<0).select(0.125*dblA,quads.col(0));
+        quads.col(1) = (cosines.col(2).array()<0).select(0.125*dblA,quads.col(1));
+        quads.col(2) = (cosines.col(2).array()<0).select( 0.25*dblA,quads.col(2));
+
+        MV.block(0*m,0,m,1) = quads.col(0);
+        MV.block(1*m,0,m,1) = quads.col(1);
+        MV.block(2*m,0,m,1) = quads.col(2);
+        
+        break;
+      }
+    case MASSMATRIX_TYPE_FULL:
+      assert(false && "Implementation incomplete");
+      break;
+    default:
+      assert(false && "Unknown Mass matrix eff_type");
+  }
+  sparse(MI,MJ,MV,n,n,M);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::massmatrix_intrinsic<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::massmatrix_intrinsic<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::massmatrix_intrinsic<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::massmatrix_intrinsic<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::MassMatrixType, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/massmatrix_intrinsic.h

@@ -0,0 +1,56 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_MASSMATRIX_INTRINSIC_H
+#define IGL_MASSMATRIX_INTRINSIC_H
+#include "igl_inline.h"
+#include "massmatrix.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl 
+{
+
+  // Constructs the mass (area) matrix for a given mesh (V,F).
+  //
+  // Inputs:
+  //   l  #l by simplex_size list of mesh edge lengths
+  //   F  #F by simplex_size list of mesh elements (triangles or tetrahedra)
+  //   type  one of the following ints:
+  //     MASSMATRIX_TYPE_BARYCENTRIC  barycentric
+  //     MASSMATRIX_TYPE_VORONOI voronoi-hybrid {default}
+  //     MASSMATRIX_TYPE_FULL full {not implemented}
+  // Outputs: 
+  //   M  #V by #V mass matrix
+  //
+  // See also: adjacency_matrix
+  //
+  template <typename Derivedl, typename DerivedF, typename Scalar>
+  IGL_INLINE void massmatrix_intrinsic(
+    const Eigen::MatrixBase<Derivedl> & l, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    const MassMatrixType type,
+    Eigen::SparseMatrix<Scalar>& M);
+  // Inputs:
+  //   n  number of vertices (>= F.maxCoeff()+1)
+  template <typename Derivedl, typename DerivedF, typename Scalar>
+  IGL_INLINE void massmatrix_intrinsic(
+    const Eigen::MatrixBase<Derivedl> & l, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    const MassMatrixType type,
+    const int n,
+    Eigen::SparseMatrix<Scalar>& M);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "massmatrix_intrinsic.cpp"
+#endif
+
+#endif
+
+

include/igl/matlab_format.cpp

@@ -114,6 +114,12 @@ IGL_INLINE const std::string igl::matlab_format(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<bool, int>, Eigen::Matrix<bool, -1, 3, 0, -1, 3> const> > const igl::matlab_format<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<bool, int>, Eigen::Matrix<bool, -1, 3, 0, -1, 3> const> >(Eigen::DenseBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<bool, int>, Eigen::Matrix<bool, -1, 3, 0, -1, 3> const> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> const> const> > const igl::matlab_format<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> const> const> >(Eigen::DenseBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> const> const> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
+template Eigen::WithFormat<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const> > const igl::matlab_format<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const> >(Eigen::DenseBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<int>, Eigen::ArrayWrapper<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+// generated by autoexplicit.sh
 template Eigen::WithFormat<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const igl::matlab_format<Eigen::Matrix<float, 1, 3, 1, 1, 3> >(Eigen::DenseBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 // generated by autoexplicit.sh
 template Eigen::WithFormat<Eigen::Matrix<int, 4, 1, 0, 4, 1> > const igl::matlab_format<Eigen::Matrix<int, 4, 1, 0, 4, 1> >(Eigen::DenseBase<Eigen::Matrix<int, 4, 1, 0, 4, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
@@ -152,4 +158,5 @@ template Eigen::WithFormat<Eigen::Matrix<double, 3, 2, 0, 3, 2> > const igl::mat
 template Eigen::WithFormat<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const igl::matlab_format<Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 template Eigen::WithFormat<Eigen::Matrix<int, 2, 2, 0, 2, 2> > const igl::matlab_format<Eigen::Matrix<int, 2, 2, 0, 2, 2> >(Eigen::DenseBase<Eigen::Matrix<int, 2, 2, 0, 2, 2> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 template Eigen::WithFormat<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const igl::matlab_format<Eigen::Matrix<float, 4, 4, 0, 4, 4> >(Eigen::DenseBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
+template Eigen::WithFormat<Eigen::Matrix<bool, -1, 1, 0, -1, 1> > const igl::matlab_format<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >);
 #endif

include/igl/min_quad_with_fixed.cpp

@@ -582,7 +582,8 @@ IGL_INLINE bool igl::min_quad_with_fixed(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
+template bool igl::min_quad_with_fixed_solve<double, Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -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> >&);
+template bool igl::min_quad_with_fixed_solve<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, Eigen::Matrix<double, -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> >&);
 template bool igl::min_quad_with_fixed<double, 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::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, 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::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template bool igl::min_quad_with_fixed_solve<double, 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::Matrix<double, -1, 1, 0, -1, 1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -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> >&);
 template bool igl::min_quad_with_fixed_precompute<double, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int> const&, bool, igl::min_quad_with_fixed_data<double>&);

include/igl/next_filename.h

@@ -17,7 +17,7 @@ namespace igl
   // Inputs:
   //   prefix  path to containing dir and filename prefix
   //   zeros number of leading zeros as if digit printed with printf
-  //   suffix  suffix of filename and extension (should included dot)
+  //   suffix  suffix of filename and extension (should include dot)
   // Outputs:
   //   next  path to next file
   // Returns true if found, false if exceeding range in zeros

include/igl/oriented_facets.h

@@ -12,14 +12,15 @@
 namespace igl
 {
   // ORIENTED_FACETS Determines all "directed
-  // [facets](https://en.wikipedia.org/wiki/Simplex#Elements)" of a given set
-  // of simplicial elements. For a manifold triangle mesh, this computes all
+  // [facets](https://en.wikipedia.org/wiki/Simplex#Elements)" of a given set of
+  // simplicial elements. For a manifold triangle mesh, this computes all
   // half-edges. For a manifold tetrahedral mesh, this computes all half-faces.
   //
   // Inputs:
-  //   F  #F by simplex_size list of simplices
+  //   F  #F by simplex_size  list of simplices
   // Outputs:
-  //   E  #E by simplex_size-1  list of facets
+  //   E  #E by simplex_size-1  list of facets, such that E.row(f+#F*c) is the
+  //     facet opposite F(f,c)
   //
   // Note: this is not the same as igl::edges because this includes every
   // directed edge including repeats (meaning interior edges on a surface will

include/igl/per_face_normals.cpp

@@ -105,6 +105,8 @@ IGL_INLINE void igl::per_face_normals_stable(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// Nonsense template
+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

include/igl/principal_curvature.cpp

@@ -938,4 +938,4 @@ template void igl::principal_curvature<Eigen::Matrix<double, -1, -1, 0, -1, -1>,
 template void igl::principal_curvature<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::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, unsigned int, bool);
 template void igl::principal_curvature<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::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> >&, 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> >&, unsigned int, bool);
 template void igl::principal_curvature<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::Matrix<double, -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&, 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> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, std::vector<int, std::allocator<int> >&, unsigned int, bool);
-#endif
+#endif

include/igl/scaf.cpp

@@ -704,4 +704,4 @@ IGL_INLINE Eigen::MatrixXd igl::scaf_solve(SCAFData &s, int iter_num)
 }
 
 #ifdef IGL_STATIC_LIBRARY
-#endif
+#endif

include/igl/simplify_polyhedron.cpp

@@ -55,7 +55,7 @@ IGL_INLINE void igl::simplify_polyhedron(
       const auto vi = E(e,positive);
       const auto vj = E(e,!positive);
       p = V.row(vj);
-      std::vector<int> faces = igl::circulation(e,positive,F,E,EMAP,EF,EI);
+      std::vector<int> faces = igl::circulation(e,positive,EMAP,EF,EI);
       cost = 0;
       for(auto f : faces)
       {

include/igl/slice.cpp

@@ -248,6 +248,8 @@ IGL_INLINE DerivedX igl::slice(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::slice<Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_multiple_op<double>, Eigen::Matrix<double, -1, 1, 0, -1, 1> const> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+// generated by autoexplicit.sh
 template void igl::slice<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
 // generated by autoexplicit.sh
 template void igl::slice<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);

include/igl/sort.cpp

@@ -317,6 +317,8 @@ if(!ascending)
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/squared_edge_lengths.cpp

@@ -73,6 +73,8 @@ IGL_INLINE void igl::squared_edge_lengths(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(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, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
 template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
 // generated by autoexplicit.sh
 template void igl::squared_edge_lengths<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> >&);

include/igl/tan_half_angle.cpp

@@ -0,0 +1,37 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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 "tan_half_angle.h"
+#include <cmath>
+template < typename Scalar>
+IGL_INLINE Scalar igl::tan_half_angle(
+  const Scalar & a,
+  const Scalar & b,
+  const Scalar & c)
+{
+  //      .
+  //     /|
+  //   c/ |
+  //   /  |
+  //  /   |
+  // .α   | a
+  //  \   |
+  //   \  |
+  //   b\ |  
+  //     \| 
+  //
+  // tan(α/2)
+  // Fisher 2007
+  return sqrt(((a-b+c)*(a+b-c))/((a+b+c)*(-a+b+c)));
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template double igl::tan_half_angle<double>(double const&, double const&, double const&);
+#endif

include/igl/tan_half_angle.h

@@ -0,0 +1,35 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_TAN_HALF_ANGLE_H
+#define IGL_TAN_HALF_ANGLE_H
+#include "igl_inline.h"
+namespace igl
+{
+  // TAN_HALF_ANGLE Compute the tangent of half of the angle opposite the side
+  // with length a, in a triangle with side lengths (a,b,c).
+  //
+  // Inputs:
+  //   a  scalar edge length of first side of triangle
+  //   b  scalar edge length of second side of triangle
+  //   c  scalar edge length of third side of triangle
+  // Returns tangent of half of the angle opposite side with length a
+  //
+  // See also: is_intrinsic_delaunay
+  template < typename Scalar>
+  IGL_INLINE Scalar tan_half_angle(
+    const Scalar & a,
+    const Scalar & b,
+    const Scalar & c);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "tan_half_angle.cpp"
+#endif
+
+#endif
+

include/igl/topological_hole_fill.cpp

@@ -52,4 +52,4 @@ IGL_INLINE void igl::topological_hole_fill(
 
 #ifdef IGL_STATIC_LIBRARY
 template void igl::topological_hole_fill<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::vector<int, std::allocator<int> > >(Eigen::MatrixBase<Eigen::Matrix<int, -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> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-#endif
+#endif

include/igl/triangulated_grid.cpp

@@ -0,0 +1,53 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@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 "triangulated_grid.h"
+#include "grid.h"
+#include <cassert>
+
+template <
+  typename XType,
+  typename YType,
+  typename DerivedGV,
+  typename DerivedGF>
+IGL_INLINE void igl::triangulated_grid(
+  const XType & nx,
+  const YType & ny,
+  Eigen::PlainObjectBase<DerivedGV> & GV,
+  Eigen::PlainObjectBase<DerivedGF> & GF)
+{
+  using namespace Eigen;
+  Eigen::Matrix<XType,2,1> res(nx,ny);
+  igl::grid(res,GV);
+  GF.resize((nx-1)*(ny-1)*2,3);
+  for(int y = 0;y<ny-1;y++)
+  {
+    for(int x = 0;x<nx-1;x++)
+    {
+      // index of southwest corner
+      const XType sw = (x  +nx*(y+0));
+      const XType se = (x+1+nx*(y+0));
+      const XType ne = (x+1+nx*(y+1));
+      const XType nw = (x  +nx*(y+1));
+      // Index of first triangle in this square
+      const XType gf = 2*(x+(nx-1)*y);
+      GF(gf+0,0) = sw;
+      GF(gf+0,1) = se;
+      GF(gf+0,2) = nw;
+      GF(gf+1,0) = se;
+      GF(gf+1,1) = ne;
+      GF(gf+1,2) = nw;
+    }
+  }
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::triangulated_grid<int, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(int const&, int const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/triangulated_grid.h

@@ -0,0 +1,38 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Alec Jacobson <alecjacobson@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/.
+#ifndef IGL_TRIANGULATEGRID_H
+#define IGL_TRIANGULATEGRID_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Create a regular grid of elements (only 2D supported, currently)
+  // Vertex position order is compatible with `igl::grid`
+  //
+  // Inputs:
+  //   nx  number of vertices in the x direction
+  //   ny  number of vertices in the y direction
+  // Outputs:
+  //   GV  nx*ny by 2 list of mesh vertex positions.
+  //   GF  2*(nx-1)*(ny-1) by 3  list of triangle indices
+  template <
+    typename XType,
+    typename YType,
+    typename DerivedGV,
+    typename DerivedGF>
+  IGL_INLINE void triangulated_grid(
+    const XType & nx,
+    const YType & ny,
+    Eigen::PlainObjectBase<DerivedGV> & GV,
+    Eigen::PlainObjectBase<DerivedGF> & GF);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "triangulated_grid.cpp"
+#endif
+#endif 
+

include/igl/unique_edge_map.cpp

@@ -58,6 +58,7 @@ template void igl::unique_edge_map<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::M
 template void igl::unique_edge_map<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
 template void igl::unique_edge_map<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<int, -1, 1, 0, -1, 1>, unsigned long>(Eigen::MatrixBase<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > >&);
+template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
 
 #ifdef WIN32
 template void igl::unique_edge_map<class Eigen::Matrix<int, -1, 3, 0, -1, 3>, class Eigen::Matrix<int, -1, 2, 0, -1, 2>, class Eigen::Matrix<int, -1, 2, 0, -1, 2>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, __int64>(class Eigen::MatrixBase<class Eigen::Matrix<int, -1, 3, 0, -1, 3> > const &, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, 2, 0, -1, 2> > &, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, 2, 0, -1, 2> > &, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> > &, class std::vector<class std::vector<__int64, class std::allocator<__int64> >, class std::allocator<class std::vector<__int64, class std::allocator<__int64> > > > &);

include/igl/unique_edge_map.h

@@ -18,7 +18,8 @@ namespace igl
   // Inputs:
   //   F  #F by 3  list of simplices
   // Outputs:
-  //   E  #F*3 by 2 list of all of directed edges
+  //   E  #F*3 by 2 list of all directed edges, such that E.row(f+#F*c) is the
+  //     edge opposite F(f,c)
   //   uE  #uE by 2 list of unique undirected edges
   //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
   //     undirected edge

include/igl/volume.cpp

@@ -105,8 +105,12 @@ IGL_INLINE void igl::volume(
        (192.*u*v*w);
   }
 }
+
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// Nonsense template
+namespace igl{ template<> void volume<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::volume<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh

tests/CMakeLists.txt

@@ -1,55 +1,86 @@
 cmake_minimum_required(VERSION 3.1)
 project(libigl_tests)
 
-list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/../cmake)		
-		
-### Adding libIGL: choose the path to your local copy libIGL		
-if(NOT TARGET igl_common)		
-  include(libigl)		
-else()		
-  include(LibiglDownloadExternal)		
+list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/../cmake)
+
+### Adding libIGL: choose the path to your local copy libIGL
+if(NOT TARGET igl::core)
+  include(libigl)
+else()
+  include(LibiglDownloadExternal)
 endif()
 
 ### Download data
 igl_download_test_data()
 set(IGL_TEST_DATA ${LIBIGL_EXTERNAL}/../tests/data)
 
-### Download Google unit test framework.
-igl_download_googletest()
+### Download Catch2 unit test framework
+igl_download_catch2()
+list(APPEND CMAKE_MODULE_PATH ${LIBIGL_EXTERNAL}/catch2/contrib)
+
+
+# Add catch2
+add_subdirectory(${LIBIGL_EXTERNAL}/catch2 catch2)
 
-SET(TEST_ROOT_DIR ${CMAKE_CURRENT_LIST_DIR})
-INCLUDE_DIRECTORIES(${TEST_ROOT_DIR})
 
-# Set TEST_DIR definition
-ADD_DEFINITIONS(-DLIBIGL_DATA_DIR="${IGL_TEST_DATA}")
+# Create test executable
+add_executable(libigl_tests main.cpp test_common.h)
+target_link_libraries(libigl_tests PUBLIC igl::core Catch2::Catch2)
+target_include_directories(libigl_tests PUBLIC ${CMAKE_CURRENT_LIST_DIR})
+
+# Set DATA_DIR definition
+set(DATA_DIR "${CMAKE_CURRENT_SOURCE_DIR}/data/")
+target_compile_definitions(libigl_tests PUBLIC -DLIBIGL_DATA_DIR="${IGL_TEST_DATA}")
 
-# Add googletest googlemock support
-ADD_SUBDIRECTORY(
-  ${LIBIGL_EXTERNAL}/googletest/googlemock
-  ${CMAKE_CURRENT_BINARY_DIR}/gtest)
-SET(GTEST_BOTH_LIBRARIES gtest gmock)
-INCLUDE_DIRECTORIES(${gmock_SOURCE_DIR})
-INCLUDE_DIRECTORIES(${gmock_SOURCE_DIR}/include)
-INCLUDE_DIRECTORIES(${gtest_SOURCE_DIR})
-INCLUDE_DIRECTORIES(${gtest_SOURCE_DIR}/include)
 
 # Process code in each subdirectories: add in decreasing order of complexity
 # (last added will run first and those should be the fastest tests)
-IF(LIBIGL_WITH_MOSEK)
-  ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl/mosek)
-ENDIF()
+if(LIBIGL_WITH_MOSEK)
+  file(GLOB TEST_SRC_FILES ./include/igl/mosek/*.cpp)
+  file(GLOB TEST_INC_FILES ./include/igl/mosek/*.h ./include/igl/mosek/*.inl)
+  target_sources(libigl_tests PRIVATE ${TEST_SRC_FILES} ${TEST_INC_FILES})
+
+  target_link_libraries(libigl_tests PUBLIC igl::mosek)
+endif()
+
+if(LIBIGL_WITH_CGAL)
+  file(GLOB TEST_SRC_FILES ./include/igl/copyleft/boolean/*.cpp ./include/igl/copyleft/cgal/*.cpp)
+  file(GLOB TEST_INC_FILES ./include/igl/copyleft/boolean/*.h ./include/igl/copyleft/cgal/*.h ./include/igl/copyleft/boolean/*.inl ./include/igl/copyleft/cgal/*.inl)
+  target_sources(libigl_tests PRIVATE ${TEST_SRC_FILES} ${TEST_INC_FILES})
+
+  target_link_libraries(libigl_tests PUBLIC igl::cgal)
+endif()
 
-IF(LIBIGL_WITH_CGAL)
-  ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl/copyleft/boolean)
-  ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl/copyleft/cgal)
-ENDIF()
+if(LIBIGL_WITH_TETGEN)
+  file(GLOB TEST_SRC_FILES ./include/igl/copyleft/tetgen/*.cpp)
+  file(GLOB TEST_INC_FILES ./include/igl/copyleft/tetgen/*.h ./include/igl/copyleft/tetgen/*.inl)
+  target_sources(libigl_tests PRIVATE ${TEST_SRC_FILES} ${TEST_INC_FILES})
+
+  target_link_libraries(libigl_tests PUBLIC igl::tetgen)
+endif()
+
+if(LIBIGL_WITH_COMISO)
+  file(GLOB TEST_SRC_FILES ./include/igl/copyleft/comiso/*.cpp)
+  file(GLOB TEST_INC_FILES ./include/igl/copyleft/comiso/*.h ./include/igl/copyleft/comiso/*.inl)
+  target_sources(libigl_tests PRIVATE ${TEST_SRC_FILES} ${TEST_INC_FILES})
+
+  target_link_libraries(libigl_tests PUBLIC igl::comiso)
+endif()
+
+if(LIBIGL_WITH_EMBREE)
+  file(GLOB TEST_SRC_FILES ./include/igl/embree/*.cpp)
+  file(GLOB TEST_INC_FILES ./include/igl/embree/*.h ./include/igl/embree/*.inl)
+  target_sources(libigl_tests PRIVATE ${TEST_SRC_FILES} ${TEST_INC_FILES})
+
+  target_link_libraries(libigl_tests PUBLIC igl::embree)
+endif()
 
-IF(LIBIGL_WITH_TETGEN)
-  ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl/copyleft/tetgen)
-ENDIF()
+file(GLOB TEST_SRC_FILES ./include/igl/*.cpp)
+file(GLOB TEST_INC_FILES ./include/igl/*.h ./include/igl/*.inl)
+target_sources(libigl_tests PRIVATE ${TEST_SRC_FILES} ${TEST_INC_FILES})
 
-IF(LIBIGL_WITH_COMISO)
-  ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl/copyleft/comiso)
-ENDIF()
 
-ADD_SUBDIRECTORY(${TEST_ROOT_DIR}/include/igl)
+# Register tests
+set(PARSE_CATCH_TESTS_ADD_TO_CONFIGURE_DEPENDS ON)
+include(Catch)
+catch_discover_tests(libigl_tests)

tests/include/igl/CMakeLists.txt

@@ -1,6 +0,0 @@
-file(GLOB TEST_SRC_FILES *.cpp main.cpp)
-file(GLOB TEST_INC_FILES *.h *.inl)
-
-add_executable(igl_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
-target_link_libraries(igl_tests igl::core gtest_main)
-add_test(NAME run_igl_tests COMMAND igl_tests)

tests/include/igl/avg_edge_length.cpp

@@ -3,7 +3,7 @@
 #include <iostream>
 
 
-TEST(avg_edge_length, cube)
+TEST_CASE("avg_edge_length: cube", "[igl]")
 {
   //The allowed error for this test
   const double epsilon = 1e-15;
@@ -31,32 +31,32 @@ TEST(avg_edge_length, cube)
   double avg;
 
   avg = igl::avg_edge_length(V,F);
-  ASSERT_NEAR((12.*sqrt(side_sq) + 6.*sqrt(diag_sq))/(12.+6.), avg, epsilon);
+  REQUIRE (avg == Approx ((12.*sqrt(side_sq) + 6.*sqrt(diag_sq))/(12.+6.)).margin( epsilon));
 
   //Check the regular tetrahedron
   avg = igl::avg_edge_length(V,F_tet);
-  ASSERT_NEAR(sqrt(diag_sq), avg, epsilon);
+  REQUIRE (avg == Approx (sqrt(diag_sq)).margin( epsilon));
 
 
   //Scale the cube to have huge sides
   side_sq = huge_scale * huge_scale;  //squared lenght of a side
   diag_sq = 2.0 * side_sq;  //squared lenght of a diagonal
   avg = igl::avg_edge_length(V_huge,F);
-  ASSERT_NEAR((12.*sqrt(side_sq) + 6.*sqrt(diag_sq))/(12.+6.), avg, epsilon*huge_scale);
+  REQUIRE (avg == Approx ((12.*sqrt(side_sq) + 6.*sqrt(diag_sq))/(12.+6.)).margin( epsilon*huge_scale));
 
   //Check the equilateral triangles
   avg = igl::avg_edge_length(V_huge,F_tet);
-  ASSERT_NEAR(sqrt(diag_sq), avg, epsilon*huge_scale);
+  REQUIRE (avg == Approx (sqrt(diag_sq)).margin( epsilon*huge_scale));
 
 
   //Scale the cube to have tiny sides
   side_sq = tiny_scale * tiny_scale;  //squared lenght of a side
   diag_sq = 2.0 * side_sq;  //squared lenght of a diagonal
   avg = igl::avg_edge_length(V_tiny,F);
-  ASSERT_NEAR((12.*sqrt(side_sq) + 6.*sqrt(diag_sq))/(12.+6.), avg, epsilon*tiny_scale);
+  REQUIRE (avg == Approx ((12.*sqrt(side_sq) + 6.*sqrt(diag_sq))/(12.+6.)).margin( epsilon*tiny_scale));
 
   //Check the regular tetrahedron
   avg = igl::avg_edge_length(V_tiny,F_tet);
-  ASSERT_NEAR(sqrt(diag_sq), avg, epsilon*tiny_scale);
+  REQUIRE (avg == Approx (sqrt(diag_sq)).margin( epsilon*tiny_scale));
 
 }

tests/include/igl/bbw.cpp

@@ -1,4 +1,3 @@
-#include <gtest/gtest.h>
 #include <test_common.h>
 #include <igl/boundary_conditions.h>
 #include <igl/readMESH.h>
@@ -6,7 +5,7 @@
 #include <igl/readTGF.h>
 #include <igl/bbw.h>
 
-TEST(bbw, decimated_knight)
+TEST_CASE("bbw: decimated_knight", "[igl]")
 {
   Eigen::MatrixXd V,C;
   Eigen::MatrixXi T,F,E;
@@ -22,6 +21,6 @@ TEST(bbw, decimated_knight)
   params.active_set_params.max_iter = 100;
   igl::bbw(V,T,b,bc,params,Was);
   // igl::writeDMAT("decimated-knight-as.dmat",Was);
-  ASSERT_LT( (Was-W_groundtruth).array().abs().maxCoeff() ,1e-4);
+  REQUIRE (1e-4 > (Was-W_groundtruth).array().abs().maxCoeff());
 }
 

tests/include/igl/boundary_loop.cpp

@@ -4,7 +4,7 @@
 #include <algorithm>
 #include <iostream>
 
-TEST(boundary_loop, cube)
+TEST_CASE("boundary_loop: cube", "[igl]")
 {
   Eigen::MatrixXd V;
   Eigen::MatrixXi F;
@@ -16,10 +16,10 @@ TEST(boundary_loop, cube)
   igl::boundary_loop(F, boundary);
 
   //The cube has no boundary
-  ASSERT_EQ(0, boundary.size());
+  REQUIRE (boundary.size() == 0);
 }
 
-TEST(boundary_loop, bunny)
+TEST_CASE("boundary_loop: bunny", "[igl]")
 {
   Eigen::MatrixXd V;
   Eigen::MatrixXi F;
@@ -31,7 +31,7 @@ TEST(boundary_loop, bunny)
   igl::boundary_loop(F, boundaries);
 
   //Compare our result with known results taken from meshlab
-  ASSERT_EQ(5, boundaries.size());
+  REQUIRE (boundaries.size() == 5);
 
   //Compute min, max and sum of boundaries
   size_t boundaryMin=9999999;
@@ -45,9 +45,9 @@ TEST(boundary_loop, bunny)
   }
 
   //Total boundary has 223 vertex
-  ASSERT_EQ(223, boundarySum);
+  REQUIRE (boundarySum == 223);
   //Largest loop has 80 vertex
-  ASSERT_EQ(80, boundaryMax);
+  REQUIRE (boundaryMax == 80);
   //Smallest loop has 22 vertex
-  ASSERT_EQ(22, boundaryMin);
+  REQUIRE (boundaryMin == 22);
 }

tests/include/igl/circulation.cpp

@@ -0,0 +1,69 @@
+#include <test_common.h>
+#include <igl/circulation.h>
+#include <igl/edge_flaps.h>
+#include <igl/unique_edge_map.h>
+#include <igl/matlab_format.h>
+
+TEST_CASE("circulation: single_edge", "[igl]")
+{
+  //       7    
+  //     /₆|₇\
+  //   4 - 5 - 6
+  //   |₂/₃|₄\₅|
+  //   1 - 2 - 3
+  //     \₀|₁/  
+  //       0     
+  const Eigen::MatrixXi F = (Eigen::MatrixXi(8,3)<<
+    0,2,1,
+    0,3,2,
+    1,5,4,
+    1,2,5,
+    2,3,5,
+    3,6,5,
+    4,5,7,
+    5,6,7).finished();
+  Eigen::MatrixXi E,uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<int> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+  Eigen::MatrixXi EI,EF;
+  {
+    const auto & cuE = uE;
+    const auto & cEMAP = EMAP;
+    igl::edge_flaps(F,cuE,cEMAP,EF,EI);
+  }
+  // Find (2,5) in uE
+  int ei = 0;
+  bool flip = false;
+  for(;ei<E.rows();ei++)
+  {
+    if(uE(ei,0) == 2 && uE(ei,1) == 5){flip=false;break;}
+    if(uE(ei,1) == 2 && uE(ei,0) == 5){flip=true;break;}
+  }
+  Eigen::VectorXi Nccw;
+  igl::circulation(ei,!flip,EMAP,EF,EI,Nccw);
+  Eigen::VectorXi Nccwgt = 
+    (Eigen::VectorXi(6)<<
+     4,
+     5,
+     7,
+     6,
+     2,
+     3).finished();
+  Eigen::VectorXi Ncwgt = 
+    (Eigen::VectorXi(4)<<
+     4,
+     1,
+     0,
+     3).finished();
+  if(flip)
+  {
+    Nccwgt = Nccwgt.reverse().eval();
+    Ncwgt = Ncwgt.reverse().eval();
+  }
+  test_common::assert_eq(Nccw,Nccwgt);
+  Eigen::VectorXi Ncw;
+  igl::circulation(ei, flip,EMAP,EF,EI,Ncw);
+  test_common::assert_eq(Ncw,Ncwgt);
+}
+

tests/include/igl/copyleft/boolean/CMakeLists.txt

@@ -1,7 +0,0 @@
-# Check if CGAL is available
-file(GLOB TEST_SRC_FILES *.cpp main.cpp)
-file(GLOB TEST_INC_FILES *.h *.inl)
-
-add_executable(igl_boolean_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
-target_link_libraries(igl_boolean_tests igl::core igl::cgal gtest_main)
-add_test(NAME run_igl_boolean_tests COMMAND igl_boolean_tests)

tests/include/igl/copyleft/boolean/main.cpp

@@ -1,6 +0,0 @@
-#include <gtest/gtest.h>
-
-int main(int argc, char **argv) {
-      ::testing::InitGoogleTest(&argc, argv);
-        return RUN_ALL_TESTS();
-}

tests/include/igl/copyleft/boolean/mesh_boolean.cpp

@@ -17,7 +17,7 @@ namespace mesh_boolean_test {
         Eigen::MatrixXi Eb;
         igl::exterior_edges(F, Eb);
 
-        ASSERT_EQ(0, Eb.rows());
+        REQUIRE (Eb.rows() == 0);
     }
 
     template<typename DerivedV, typename DerivedF>
@@ -25,8 +25,8 @@ namespace mesh_boolean_test {
             const Eigen::PlainObjectBase<DerivedV>& V,
             const Eigen::PlainObjectBase<DerivedF>& F) {
         Eigen::MatrixXi B;
-        ASSERT_TRUE(igl::is_vertex_manifold(F, B));
-        ASSERT_TRUE(igl::is_edge_manifold(F));
+        REQUIRE (igl::is_vertex_manifold(F, B));
+        REQUIRE (igl::is_edge_manifold(F));
     }
 
     template<typename DerivedV, typename DerivedF>
@@ -47,10 +47,10 @@ namespace mesh_boolean_test {
 
         const int num_edges = uE.rows();
         const int euler = num_vertices - num_edges + num_faces;
-        ASSERT_EQ(euler, 2 - 2 * genus);
+        REQUIRE (2 - 2 * genus == euler);
     }
 
-TEST(MeshBoolean, TwoCubes) {
+TEST_CASE("MeshBoolean: TwoCubes", "[igl/copyleft/boolean]")
     Eigen::MatrixXd V1;
     Eigen::MatrixXi F1;
     test_common::load_mesh("two-boxes-bad-self-union.ply", V1, F1);
@@ -70,7 +70,7 @@ TEST(MeshBoolean, TwoCubes) {
     assert_genus_eq(Vo, Fo, 0);
 }
 
-TEST(MeshBoolean, MinusTest) {
+TEST_CASE("MeshBoolean: MinusTest", "[igl/copyleft/boolean]")
     // Many thanks to Eric Yao for submitting this test case.
     Eigen::MatrixXd V1, V2, Vo;
     Eigen::MatrixXi F1, F2, Fo;
@@ -86,7 +86,7 @@ TEST(MeshBoolean, MinusTest) {
     assert_genus_eq(Vo, Fo, 1);
 }
 
-TEST(MeshBoolean, IntersectWithSelf) {
+TEST_CASE("MeshBoolean: IntersectWithSelf", "[igl/copyleft/boolean]")
     Eigen::MatrixXd V1, Vo;
     Eigen::MatrixXi F1, Fo;
     test_common::load_mesh("cube.obj", V1, F1);

tests/include/igl/copyleft/cgal/CMakeLists.txt

@@ -1,6 +0,0 @@
-file(GLOB TEST_SRC_FILES *.cpp main.cpp)
-file(GLOB TEST_INC_FILES *.h *.inl)
-
-add_executable(igl_cgal_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
-target_link_libraries(igl_cgal_tests igl::core igl::cgal gtest_main)
-add_test(NAME run_igl_cgal_tests COMMAND igl_cgal_tests)

tests/include/igl/copyleft/cgal/CSGTree.cpp

@@ -2,7 +2,7 @@
 
 #include <igl/copyleft/cgal/CSGTree.h>
 
-TEST(CSGTree, extrusion) {
+TEST_CASE("CSGTree: extrusion", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V;
     Eigen::MatrixXi F;
     test_common::load_mesh("extrusion.obj", V, F);
@@ -12,6 +12,6 @@ TEST(CSGTree, extrusion) {
     Eigen::MatrixXd V2 = inter.cast_V<Eigen::MatrixXd>();
     Eigen::MatrixXi F2 = inter.F();
 
-    ASSERT_EQ(V.rows(), V2.rows());
-    ASSERT_EQ(F.rows(), F2.rows());
+    REQUIRE (V2.rows() == V.rows());
+    REQUIRE (F2.rows() == F.rows());
 }

tests/include/igl/copyleft/cgal/delaunay_triangulation.cpp

@@ -0,0 +1,26 @@
+#include <test_common.h>
+#include <igl/copyleft/cgal/delaunay_triangulation.h>
+#include <igl/unique_simplices.h>
+#include <igl/matlab_format.h>
+
+TEST_CASE("igl_copyleft_cgal_delaunay_triangulation: two_triangles", "[igl/copyleft/cgal]")
+{
+  const Eigen::MatrixXd V = 
+    (Eigen::MatrixXd(4,2)<<
+     0,10,
+     1,0,
+     1,20,
+     2,10).finished();
+  Eigen::MatrixXi F;
+  igl::copyleft::cgal::delaunay_triangulation(V,F);
+  // Ground truth
+  Eigen::MatrixXi Fgt = (Eigen::MatrixXi(2,3)<<0,1,3,0,3,2).finished();
+  REQUIRE (2 == F.rows());
+  Eigen::MatrixXi Fu;
+  Eigen::VectorXi IA,IC;
+  igl::unique_simplices(
+    (Eigen::MatrixXi(4,3)<<F,Fgt).finished(),
+    Fu,IA,IC);
+  // Now new faces w.r.t. ground truth
+  REQUIRE (2 == Fu.rows());
+}

tests/include/igl/copyleft/cgal/hausdorff.cpp

@@ -5,7 +5,7 @@
 #include <igl/copyleft/cgal/point_mesh_squared_distance.h>
 #include <igl/upsample.h>
 
-TEST(hausdorff, knightVScheburashka) 
+TEST_CASE("hausdorff: knightVScheburashka", "[igl/copyleft/cgal]")
 {
   Eigen::MatrixXd VA,VB;
   Eigen::MatrixXi FA,FB;
@@ -46,7 +46,7 @@ TEST(hausdorff, knightVScheburashka)
         {
           u4 = std::min(u4,U[j](i-u));
         }
-        ASSERT_LE(u4,U[j-1](i/4));
+        REQUIRE (U[j-1](i/4) >= u4);
       }
     }
     break;

tests/include/igl/copyleft/cgal/main.cpp

@@ -1,6 +0,0 @@
-#include <gtest/gtest.h>
-
-int main(int argc, char **argv) {
-      ::testing::InitGoogleTest(&argc, argv);
-        return RUN_ALL_TESTS();
-}

tests/include/igl/copyleft/cgal/order_facets_around_edges.cpp

@@ -25,7 +25,7 @@ void assert_consistently_oriented(size_t num_faces,
         const std::vector<int>& expected_face_order,
         const std::vector<int>& e_order) {
     const size_t num_items = expected_face_order.size();
-    ASSERT_EQ(num_items, e_order.size());
+    REQUIRE (e_order.size() == num_items);
 
     std::vector<int> order(num_items);
     std::transform(e_order.begin(), e_order.end(), order.begin(),
@@ -33,7 +33,7 @@ void assert_consistently_oriented(size_t num_faces,
 
     size_t ref_start = index_of(order, expected_face_order[0]);
     for (size_t i=0; i<num_items; i++) {
-        ASSERT_EQ(expected_face_order[i], order[(ref_start + i) % num_items]);
+        REQUIRE (order[(ref_start + i) % num_items] == expected_face_order[i]);
     }
 }
 
@@ -83,7 +83,7 @@ void assert_order(
     }
 }
 
-TEST(copyleft_cgal_order_facets_around_edges, Simple) {
+TEST_CASE("copyleft_cgal_order_facets_around_edges: Simple", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V(4, 3);
     V << 0.0, 0.0, 0.0,
          1.0, 0.0, 0.0,
@@ -96,7 +96,7 @@ TEST(copyleft_cgal_order_facets_around_edges, Simple) {
     assert_order(V, F, 1, 2, {0, 1});
 }
 
-TEST(copyleft_cgal_order_facets_around_edges, TripletFaces) {
+TEST_CASE("copyleft_cgal_order_facets_around_edges: TripletFaces", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V(5, 3);
     V << 0.0, 0.0, 0.0,
          1.0, 0.0, 0.0,
@@ -111,7 +111,7 @@ TEST(copyleft_cgal_order_facets_around_edges, TripletFaces) {
     assert_order(V, F, 1, 2, {0, 1, 2});
 }
 
-TEST(copyleft_cgal_order_facets_around_edges, DuplicatedFaces) {
+TEST_CASE("copyleft_cgal_order_facets_around_edges: DuplicatedFaces", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V(5, 3);
     V << 0.0, 0.0, 0.0,
          1.0, 0.0, 0.0,
@@ -127,7 +127,7 @@ TEST(copyleft_cgal_order_facets_around_edges, DuplicatedFaces) {
     assert_order(V, F, 1, 2, {0, 1, 3, 2});
 }
 
-TEST(copyleft_cgal_order_facets_around_edges, MultipleDuplicatedFaces) {
+TEST_CASE("copyleft_cgal_order_facets_around_edges: MultipleDuplicatedFaces", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V(5, 3);
     V << 0.0, 0.0, 0.0,
          1.0, 0.0, 0.0,
@@ -145,7 +145,7 @@ TEST(copyleft_cgal_order_facets_around_edges, MultipleDuplicatedFaces) {
     assert_order(V, F, 1, 2, {1, 0, 2, 3, 5, 4});
 }
 
-TEST(copyleft_cgal_order_facets_around_edges, Debug) {
+TEST_CASE("copyleft_cgal_order_facets_around_edges: Debug", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V(5, 3);
     V <<
         -44.3205080756887781, 4.22994972382184579e-15, 75,
@@ -162,7 +162,7 @@ TEST(copyleft_cgal_order_facets_around_edges, Debug) {
     assert_order(V, F, 1, 2, {0, 2, 1});
 }
 
-TEST(copyleft_cgal_order_facets_around_edges, Debug2) {
+TEST_CASE("copyleft_cgal_order_facets_around_edges: Debug2", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V(5, 3);
     V <<
         -22.160254037844382, 38.3826859021798441, 75,
@@ -178,7 +178,7 @@ TEST(copyleft_cgal_order_facets_around_edges, Debug2) {
     assert_order(V, F, 1, 2, {1, 0, 2});
 }
 
-TEST(copyleft_cgal_order_facets_around_edges, NormalSensitivity) {
+TEST_CASE("copyleft_cgal_order_facets_around_edges: NormalSensitivity", "[igl/copyleft/cgal]")
     // This example shows that epsilon difference in normal vectors could
     // results in very different ordering of facets.
 

tests/include/igl/copyleft/cgal/outer_facet.cpp

@@ -16,7 +16,7 @@ void assert_outer_facet_is_correct(
     // Todo.
 }
 
-TEST(OuterFacet, Simple) {
+TEST_CASE("OuterFacet: Simple", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V;
     Eigen::MatrixXi F;
     test_common::load_mesh("cube.obj", V, F);
@@ -30,11 +30,11 @@ TEST(OuterFacet, Simple) {
     bool flipped;
     igl::copyleft::cgal::outer_facet(V, F, I, fid, flipped);
 
-    ASSERT_LT(fid, num_faces);
-    ASSERT_FALSE(flipped);
+    REQUIRE (num_faces > fid);
+    REQUIRE (!flipped);
 }
 
-TEST(OuterFacet, DuplicatedOppositeFaces) {
+TEST_CASE("OuterFacet: DuplicatedOppositeFaces", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V;
     Eigen::MatrixXi F1;
     test_common::load_mesh("cube.obj", V, F1);
@@ -52,11 +52,11 @@ TEST(OuterFacet, DuplicatedOppositeFaces) {
     bool flipped;
     igl::copyleft::cgal::outer_facet(V, F, I, fid, flipped);
 
-    ASSERT_LT(fid, F.rows());
-    ASSERT_FALSE(flipped);
+    REQUIRE (F.rows() > fid);
+    REQUIRE (!flipped);
 }
 
-TEST(OuterFacet, FullyDegnerated) {
+TEST_CASE("OuterFacet: FullyDegnerated", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V;
     Eigen::MatrixXi F;
     test_common::load_mesh("degenerated.obj", V, F);
@@ -68,11 +68,11 @@ TEST(OuterFacet, FullyDegnerated) {
     bool flipped;
     igl::copyleft::cgal::outer_facet(V, F, I, fid, flipped);
 
-    ASSERT_LT(fid, F.rows());
-    ASSERT_FALSE(flipped);
+    REQUIRE (F.rows() > fid);
+    REQUIRE (!flipped);
 }
 
-TEST(OuterFacet, InvertedNormal) {
+TEST_CASE("OuterFacet: InvertedNormal", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V;
     Eigen::MatrixXi F;
     test_common::load_mesh("cube.obj", V, F);
@@ -85,11 +85,11 @@ TEST(OuterFacet, InvertedNormal) {
     bool flipped;
     igl::copyleft::cgal::outer_facet(V, F, I, fid, flipped);
 
-    ASSERT_LT(fid, F.rows());
-    ASSERT_TRUE(flipped);
+    REQUIRE (F.rows() > fid);
+    REQUIRE (flipped);
 }
 
-TEST(OuterFacet, SliverTet) {
+TEST_CASE("OuterFacet: SliverTet", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V;
     Eigen::MatrixXi F;
     test_common::load_mesh("sliver_tet.ply", V, F);
@@ -101,8 +101,8 @@ TEST(OuterFacet, SliverTet) {
     bool flipped;
     igl::copyleft::cgal::outer_facet(V, F, I, fid, flipped);
 
-    ASSERT_LT(fid, F.rows());
-    ASSERT_FALSE(flipped);
+    REQUIRE (F.rows() > fid);
+    REQUIRE (!flipped);
 }
 
 }

tests/include/igl/copyleft/cgal/outer_hull.cpp

@@ -3,7 +3,7 @@
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <igl/copyleft/cgal/outer_hull.h>
 
-TEST(OuterHull, CubeWithFold) {
+TEST_CASE("OuterHull: CubeWithFold", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V;
     Eigen::MatrixXi F;
     test_common::load_mesh("cube_with_fold.ply", V, F);

tests/include/igl/copyleft/cgal/peel_outer_hull_layers.cpp

@@ -11,12 +11,12 @@
 
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
-TEST(copyleft_cgal_peel_outer_hull_layers, TwoCubes) {
+TEST_CASE("copyleft_cgal_peel_outer_hull_layers: TwoCubes", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V;
     Eigen::MatrixXi F;
     test_common::load_mesh("two-boxes-bad-self-union.ply", V, F);
-    ASSERT_EQ(486, V.rows());
-    ASSERT_EQ(708, F.rows());
+    REQUIRE (V.rows() == 486);
+    REQUIRE (F.rows() == 708);
 
     typedef CGAL::Exact_predicates_exact_constructions_kernel K;
     typedef K::FT Scalar;
@@ -45,12 +45,12 @@ TEST(copyleft_cgal_peel_outer_hull_layers, TwoCubes) {
             vertices.data(), [](Scalar v) { return CGAL::to_double(v); });
     igl::writeOBJ("debug.obj", vertices, Ft);
 
-    ASSERT_EQ(num_faces, I.rows());
-    ASSERT_EQ(0, I.minCoeff());
-    ASSERT_EQ(1, I.maxCoeff());
+    REQUIRE (I.rows() == num_faces);
+    REQUIRE (I.minCoeff() == 0);
+    REQUIRE (I.maxCoeff() == 1);
 }
 
-TEST(PeelOuterHullLayers, CubeWithFold) {
+TEST_CASE("PeelOuterHullLayers: CubeWithFold", "[igl/copyleft/cgal]")
     Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> V;
     Eigen::MatrixXi F;
     test_common::load_mesh("cube_with_fold.ply", V, F);

tests/include/igl/copyleft/cgal/points_inside_component.cpp

@@ -5,7 +5,7 @@
 
 namespace PointInsideComponentHelper {
 
-TEST(PointInsideComponent, simple) {
+TEST_CASE("PointInsideComponent: simple", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V1;
     Eigen::MatrixXi F1;
     test_common::load_mesh("cube.obj", V1, F1);
@@ -17,14 +17,14 @@ TEST(PointInsideComponent, simple) {
          0.0, 0.0, 1.0;
     Eigen::VectorXi inside;
 
-    EXPECT_NO_THROW(igl::copyleft::cgal::points_inside_component(V1, F1, P, inside));
-    ASSERT_EQ(1, inside[0]);
-    ASSERT_EQ(0, inside[1]);
-    ASSERT_EQ(0, inside[2]);
-    ASSERT_EQ(0, inside[3]);
+    CHECK_NOTHROW (igl::copyleft::cgal::points_inside_component(V1, F1, P, inside));
+    REQUIRE (inside[0] == 1);
+    REQUIRE (inside[1] == 0);
+    REQUIRE (inside[2] == 0);
+    REQUIRE (inside[3] == 0);
 }
 
-TEST(PointInsideComponent, near_boundary) {
+TEST_CASE("PointInsideComponent: near_boundary", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V1;
     Eigen::MatrixXi F1;
     test_common::load_mesh("cube.obj", V1, F1);
@@ -39,16 +39,16 @@ TEST(PointInsideComponent, near_boundary) {
          0.0, 0.0, 0.5 - EPS;
 
     Eigen::VectorXi inside;
-    EXPECT_NO_THROW(igl::copyleft::cgal::points_inside_component(V1, F1, P, inside));
-    ASSERT_EQ(0, inside[0]);
-    ASSERT_EQ(0, inside[1]);
-    ASSERT_EQ(0, inside[2]);
-    ASSERT_EQ(1, inside[3]);
-    ASSERT_EQ(1, inside[4]);
-    ASSERT_EQ(1, inside[5]);
+    CHECK_NOTHROW (igl::copyleft::cgal::points_inside_component(V1, F1, P, inside));
+    REQUIRE (inside[0] == 0);
+    REQUIRE (inside[1] == 0);
+    REQUIRE (inside[2] == 0);
+    REQUIRE (inside[3] == 1);
+    REQUIRE (inside[4] == 1);
+    REQUIRE (inside[5] == 1);
 }
 
-TEST(PointInsideComponent, near_corner) {
+TEST_CASE("PointInsideComponent: near_corner", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V1;
     Eigen::MatrixXi F1;
     test_common::load_mesh("cube.obj", V1, F1);
@@ -65,8 +65,8 @@ TEST(PointInsideComponent, near_corner) {
             -0.5 - EPS,-0.5 - EPS,-0.5 - EPS;
 
     Eigen::VectorXi inside;
-    EXPECT_NO_THROW(igl::copyleft::cgal::points_inside_component(V1, F1, P_out, inside));
-    ASSERT_TRUE((inside.array()==0).all());
+    CHECK_NOTHROW (igl::copyleft::cgal::points_inside_component(V1, F1, P_out, inside));
+    REQUIRE ((inside.array()==0).all());
 
     Eigen::MatrixXd P_in(8, 3);
     P_in << 0.5 - EPS, 0.5 - EPS, 0.5 - EPS,
@@ -77,8 +77,8 @@ TEST(PointInsideComponent, near_corner) {
            -0.5 + EPS, 0.5 - EPS,-0.5 + EPS,
             0.5 - EPS,-0.5 + EPS,-0.5 + EPS,
            -0.5 + EPS,-0.5 + EPS,-0.5 + EPS;
-    EXPECT_NO_THROW(igl::copyleft::cgal::points_inside_component(V1, F1, P_in, inside));
-    ASSERT_TRUE((inside.array()==1).all());
+    CHECK_NOTHROW (igl::copyleft::cgal::points_inside_component(V1, F1, P_in, inside));
+    REQUIRE ((inside.array()==1).all());
 }
 
 }

tests/include/igl/copyleft/cgal/remesh_self_intersections.cpp

@@ -6,7 +6,7 @@
 
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
-TEST(RemeshSelfIntersections, CubeWithFold) {
+TEST_CASE("RemeshSelfIntersections: CubeWithFold", "[igl/copyleft/cgal]")
     Eigen::MatrixXd V;
     Eigen::MatrixXi F;
     test_common::load_mesh("cube_with_fold.ply", V, F);

tests/include/igl/copyleft/comiso/CMakeLists.txt

@@ -1,6 +0,0 @@
-file(GLOB TEST_SRC_FILES *.cpp main.cpp)
-file(GLOB TEST_INC_FILES *.h *.inl)
-
-add_executable(igl_comiso_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
-target_link_libraries(igl_comiso_tests igl::core igl::comiso gtest_main)
-add_test(NAME run_igl_comiso_tests COMMAND igl_comiso_tests)