Merge branch 'master' of https://github.com/libigl/libigl

.github/workflows/continuous.yml

@@ -25,13 +25,13 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-18.04, macos-latest]
+        os: [ubuntu-20.04, macos-latest]
         config: [Release]
         static: [ON, OFF]
         include:
           - os: macos-latest
             name: macOS
-          - os: ubuntu-18.04
+          - os: ubuntu-20.04
             name: Linux
     env:
       LIBIGL_NUM_THREADS: 1  # See https://github.com/libigl/libigl/pull/996

.github/workflows/nightly.yml

@@ -23,47 +23,47 @@ jobs:
       fail-fast: false
       matrix:
         name: [
-          ubuntu-18.04-gcc-7,
-          ubuntu-18.04-gcc-8,
-          ubuntu-18.04-gcc-9,
-          ubuntu-18.04-clang-7,
-          ubuntu-18.04-clang-8,
-          ubuntu-18.04-clang-9,
+          ubuntu-20.04-gcc-8,
+          ubuntu-20.04-gcc-9,
+          ubuntu-20.04-gcc-10,
+          ubuntu-20.04-clang-8,
+          ubuntu-20.04-clang-9,
+          ubuntu-20.04-clang-10,
           macOS-latest,
         ]
         config: [Debug, Release]
         static: [ON, OFF]
         include:
-          - name: ubuntu-18.04-gcc-7
-            os: ubuntu-18.04
-            compiler: gcc
-            version: "7"
-
-          - name: ubuntu-18.04-gcc-8
-            os: ubuntu-18.04
+          - name: ubuntu-20.04-gcc-8
+            os: ubuntu-20.04
             compiler: gcc
             version: "8"
 
-          - name: ubuntu-18.04-gcc-9
-            os: ubuntu-18.04
+          - name: ubuntu-20.04-gcc-9
+            os: ubuntu-20.04
             compiler: gcc
             version: "9"
 
-          - name: ubuntu-18.04-clang-7
-            os: ubuntu-18.04
-            compiler: clang
-            version: "7"
+          - name: ubuntu-20.04-gcc-10
+            os: ubuntu-20.04
+            compiler: gcc
+            version: "10"
 
-          - name: ubuntu-18.04-clang-8
-            os: ubuntu-18.04
+          - name: ubuntu-20.04-clang-8
+            os: ubuntu-20.04
             compiler: clang
             version: "8"
 
-          - name: ubuntu-18.04-clang-9
-            os: ubuntu-18.04
+          - name: ubuntu-20.04-clang-9
+            os: ubuntu-20.04
             compiler: clang
             version: "9"
 
+          - name: ubuntu-20.04-clang-10
+            os: ubuntu-20.04
+            compiler: clang
+            version: "10"
+
           - name: macOS-latest
             os: macOS-latest
 
@@ -84,10 +84,6 @@ jobs:
       - name: Dependencies (Linux)
         if: runner.os == 'Linux'
         run: |
-            # LLVM 9 is not in Bionic's repositories so we add the official LLVM repository.
-            if [ "${{ matrix.compiler }}" = "clang" ] && [ "${{ matrix.version }}" = "9" ]; then
-              sudo add-apt-repository "deb http://apt.llvm.org/bionic/ llvm-toolchain-bionic-9 main"
-            fi
             sudo apt-get update
 
             if [ "${{ matrix.compiler }}" = "gcc" ]; then

cmake/LibiglDownloadExternal.cmake

@@ -67,7 +67,7 @@ endfunction()
 set(LIBIGL_EIGEN_VERSION 3.3.7 CACHE STRING "Default version of Eigen used by libigl.")
 function(igl_download_eigen)
 	igl_download_project(eigen
-		GIT_REPOSITORY https://github.com/eigenteam/eigen-git-mirror.git
+		GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
 		GIT_TAG        ${LIBIGL_EIGEN_VERSION}
 		${LIBIGL_BRANCH_OPTIONS}
 	)
@@ -157,7 +157,7 @@ endfunction()
 function(igl_download_triangle)
 	igl_download_project(triangle
 		GIT_REPOSITORY https://github.com/libigl/triangle.git
-		GIT_TAG        d284c4a843efac043c310f5fa640b17cf7d96170
+		GIT_TAG        5a70326574b34d6a51d9eaf6a9f78813657ee108
 	)
 endfunction()
 
@@ -173,7 +173,7 @@ endfunction()
 function(igl_download_predicates)
 	igl_download_project(predicates
 		GIT_REPOSITORY https://github.com/libigl/libigl-predicates.git
-		GIT_TAG        5a1d2194ec114bff51d5a33230586cafb83adc86
+		GIT_TAG        488242fa2b1f98a9c5bd1441297fb4a99a6a9ae4
 	)
 endfunction()
 
@@ -196,4 +196,3 @@ function(igl_download_tutorial_data)
 		GIT_TAG        37d4e836054c9c2d2125a817c489ed8e07cd56fc
 	)
 endfunction()
-

include/igl/exact_geodesic.cpp

@@ -1027,7 +1027,6 @@ inline void Mesh::build_adjacencies()
 
 inline bool Mesh::verify()		//verifies connectivity of the mesh and prints some debug info
 {
-	std::cout << std::endl;
 	// make sure that all vertices are mentioned at least once.
 	// though the loose vertex is not a bug, it most likely indicates that something is wrong with the mesh
 	std::vector<bool> map(m_vertices.size(), false);

include/igl/fast_winding_number.cpp

@@ -434,21 +434,37 @@ IGL_INLINE void igl::fast_winding_number(
   assert(Q.cols() == 3 && "Q should be 3D");
   W.resize(Q.rows(),1);
   igl::parallel_for(Q.rows(),[&](int p)
-    {
+  {
     FastWindingNumber::HDK_Sample::UT_Vector3T<float>Qp;
-          Qp[0] = Q(p,0);
-          Qp[1] = Q(p,1);
-          Qp[2] = Q(p,2);
-      W(p) = fwn_bvh.ut_solid_angle.computeSolidAngle(
-        Qp,
-        accuracy_scale) 
-        / (4.0*igl::PI);
-    },1000);
+    Qp[0] = Q(p,0);
+    Qp[1] = Q(p,1);
+    Qp[2] = Q(p,2);
+    W(p) = fwn_bvh.ut_solid_angle.computeSolidAngle(Qp,accuracy_scale) / (4.0*igl::PI);
+  },1000);
+}
+
+template <typename Derivedp>
+IGL_INLINE typename Derivedp::Scalar igl::fast_winding_number(
+  const FastWindingNumberBVH & fwn_bvh,
+  const float accuracy_scale,
+  const Eigen::MatrixBase<Derivedp> & p)
+{
+  assert(p.cols() == 3 && "p should be 3D");
+
+  FastWindingNumber::HDK_Sample::UT_Vector3T<float>Qp;
+  Qp[0] = p(0,0);
+  Qp[1] = p(0,1);
+  Qp[2] = p(0,2);
+
+  typename Derivedp::Scalar w = fwn_bvh.ut_solid_angle.computeSolidAngle(Qp,accuracy_scale) / (4.0*igl::PI);
+
+  return w;
 }
 
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template Eigen::Matrix<float, -1, -1, 0, -1, -1>::Scalar igl::fast_winding_number<Eigen::Matrix<float, -1, -1, 0, -1, -1> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&);
 // generated by autoexplicit.sh
 template void igl::fast_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::fast_winding_number<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::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::Matrix<double, -1, 1, 0, -1, 1> >(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&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::MatrixBase<Eigen::Matrix<int, -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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
@@ -457,4 +473,10 @@ template void igl::fast_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>,
 template void igl::fast_winding_number<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&, int, igl::FastWindingNumberBVH&);
 template void igl::fast_winding_number<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
 template void igl::fast_winding_number<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, igl::FastWindingNumberBVH&);
+
+// tom did this manually. Unsure how to generate otherwise... sorry.
+template Eigen::Matrix<float, 1, 3, 1, 1, 3>::Scalar igl::fast_winding_number<Eigen::Matrix<float, 1, 3, 1, 1, 3> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&);
+template void igl::fast_winding_number<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -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&, int, igl::FastWindingNumberBVH&);
+template void igl::fast_winding_number<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, int, igl::FastWindingNumberBVH&);
+template Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, 1, 3, 1, 1, 3> const>::Scalar igl::fast_winding_number<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, 1, 3, 1, 1, 3> const> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, 1, 3, 1, 1, 3> const> > const&);
 #endif

include/igl/fast_winding_number.h

@@ -206,6 +206,19 @@ namespace igl
     const float accuracy_scale,
     const Eigen::MatrixBase<DerivedQ> & Q,
     Eigen::PlainObjectBase<DerivedW> & W);
+  // After precomputation, compute winding number at a a single point
+  //
+  // Inputs:
+  //   fwn_bvh  Precomputed bounding volume hierarchy
+  //   accuracy_scale  parameter controlling accuracy (e.g., 2)
+  //   p single position
+  // Outputs:
+  //   w  winding number of this point
+  template <typename Derivedp>
+  IGL_INLINE typename Derivedp::Scalar fast_winding_number(
+    const FastWindingNumberBVH & fwn_bvh,
+    const float accuracy_scale,
+    const Eigen::MatrixBase<Derivedp> & p);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "fast_winding_number.cpp"

include/igl/file_utils.cpp

@@ -0,0 +1,29 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Jérémie Dumas <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "file_utils.h"
+
+#include <cstring>
+
+namespace igl {
+
+IGL_INLINE void read_file_binary(FILE *fp,
+                                 std::vector<uint8_t> &fileBufferBytes) {
+  if (!ferror(fp)) {
+    fseek(fp, 0, SEEK_END);
+    size_t sizeBytes = ftell(fp);
+    fseek(fp, 0, SEEK_SET);
+    fileBufferBytes.resize(sizeBytes);
+
+    if (fread((char *)fileBufferBytes.data(), 1, sizeBytes, fp) == sizeBytes)
+      return;
+  }
+
+  throw std::runtime_error("error reading from file");
+}
+
+}

include/igl/file_utils.h

@@ -0,0 +1,62 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Jérémie Dumas <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_FILE_UTILS_H
+#define IGL_FILE_UTILS_H
+
+#include "igl_inline.h"
+
+#include <streambuf>
+#include <istream>
+#include <string>
+#include <vector>
+
+namespace igl {
+
+IGL_INLINE void read_file_binary(FILE *fp,
+                                 std::vector<uint8_t> &fileBufferBytes);
+
+struct file_memory_buffer : public std::streambuf {
+  char *p_start{nullptr};
+  char *p_end{nullptr};
+  size_t size;
+
+  file_memory_buffer(char const *first_elem, size_t size)
+      : p_start(const_cast<char *>(first_elem)), p_end(p_start + size),
+        size(size) {
+    setg(p_start, p_start, p_end);
+  }
+
+  pos_type seekoff(off_type off, std::ios_base::seekdir dir,
+                   std::ios_base::openmode which) override {
+    if (dir == std::ios_base::cur) {
+      gbump(static_cast<int>(off));
+    } else {
+      setg(p_start, (dir == std::ios_base::beg ? p_start : p_end) + off, p_end);
+    }
+    return gptr() - p_start;
+  }
+
+  pos_type seekpos(pos_type pos, std::ios_base::openmode which) override {
+    return seekoff(pos, std::ios_base::beg, which);
+  }
+};
+
+struct file_memory_stream : virtual file_memory_buffer, public std::istream {
+  file_memory_stream(char const *first_elem, size_t size)
+      : file_memory_buffer(first_elem, size), std::istream(
+                                                  static_cast<std::streambuf *>(
+                                                      this)) {}
+};
+
+} // namespace igl
+
+#ifndef IGL_STATIC_LIBRARY
+#include "file_utils.cpp"
+#endif
+
+#endif

include/igl/list_to_matrix.cpp

@@ -52,6 +52,39 @@ IGL_INLINE bool igl::list_to_matrix(const std::vector<std::vector<T > > & V,Eige
   return true;
 }
 
+template <typename T, size_t N, typename Derived>
+IGL_INLINE bool igl::list_to_matrix(const std::vector<std::array<T, N> > & V,Eigen::PlainObjectBase<Derived>& M)
+{
+  // number of rows
+  int m = V.size();
+  if(m == 0)
+  {
+    M.resize(
+            Derived::RowsAtCompileTime>=0?Derived::RowsAtCompileTime:0
+            ,
+            Derived::ColsAtCompileTime>=0?Derived::ColsAtCompileTime:0
+            );
+    return true;
+  }
+  // number of columns
+  int n = static_cast<int>(N);
+  assert(n != -1);
+  // Resize output
+  M.resize(m,n);
+
+  // Loop over rows
+  for(int i = 0;i<m;i++)
+  {
+    // Loop over cols
+    for(int j = 0;j<n;j++)
+    {
+      M(i,j) = V[i][j];
+    }
+  }
+
+  return true;
+}
+
 template <typename T, typename Derived>
 IGL_INLINE bool igl::list_to_matrix(
   const std::vector<std::vector<T > > & V,
@@ -172,6 +205,18 @@ template bool igl::list_to_matrix<int, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std
 template bool igl::list_to_matrix<unsigned long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(std::vector<unsigned long, std::allocator<unsigned long> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 template bool igl::list_to_matrix<bool, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(std::vector<bool, std::allocator<bool> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
 
+template bool igl::list_to_matrix<double, 3ul, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::vector<std::array<double, 3ul>, std::allocator<std::array<double, 3ul> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::list_to_matrix<double, 3ul, Eigen::Matrix<double, -1, -1, 1, -1, -1> >(std::vector<std::array<double, 3ul>, std::allocator<std::array<double, 3ul> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> >&);
+template bool igl::list_to_matrix<double, 3ul, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(std::vector<std::array<double, 3ul>, std::allocator<std::array<double, 3ul> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template bool igl::list_to_matrix<double, 3ul, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(std::vector<std::array<double, 3ul>, std::allocator<std::array<double, 3ul> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
+template bool igl::list_to_matrix<float, 3ul, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(std::vector<std::array<float, 3ul>, std::allocator<std::array<float, 3ul> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
+template bool igl::list_to_matrix<float, 3ul, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(std::vector<std::array<float, 3ul>, std::allocator<std::array<float, 3ul> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
+template bool igl::list_to_matrix<float, 3ul, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(std::vector<std::array<float, 3ul>, std::allocator<std::array<float, 3ul> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
+template bool igl::list_to_matrix<int, 3ul, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<std::array<int, 3ul>, std::allocator<std::array<int, 3ul> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template bool igl::list_to_matrix<int, 3ul, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::vector<std::array<int, 3ul>, std::allocator<std::array<int, 3ul> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template bool igl::list_to_matrix<int, 3ul, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::vector<std::array<int, 3ul>, std::allocator<std::array<int, 3ul> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
+template bool igl::list_to_matrix<unsigned int, 3ul, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(std::vector<std::array<unsigned int, 3ul>, std::allocator<std::array<unsigned int, 3ul> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
+
 #ifdef WIN32
 template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::vector<double, std::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::vector<double, std::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);

include/igl/list_to_matrix.h

@@ -9,6 +9,7 @@
 #define IGL_LIST_TO_MATRIX_H
 #include "igl_inline.h"
 #include <vector>
+#include <array>
 #include <Eigen/Core>
 
 namespace igl
@@ -28,6 +29,12 @@ namespace igl
   IGL_INLINE bool list_to_matrix(
     const std::vector<std::vector<T > > & V,
     Eigen::PlainObjectBase<Derived>& M);
+
+  template <typename T, size_t N, typename Derived>
+  IGL_INLINE bool list_to_matrix(
+    const std::vector<std::array<T, N> > & V,
+    Eigen::PlainObjectBase<Derived>& M);
+
   // Convert a list of row vectors of `n` or less to a matrix and pad on
   // the right with `padding`:
   //

include/igl/opengl/MeshGL.cpp

@@ -373,14 +373,12 @@ R"(#version 330
     in float offset;
     uniform mat4 view;
     uniform mat4 proj;
-    out int vPosition;
     out int vCharacter;
     out float vOffset;
     void main()
     {
       vCharacter = int(character);
       vOffset = offset;
-      vPosition = gl_VertexID;
       gl_Position = proj * view * vec4(position, 1.0);
     }
 )";
@@ -397,7 +395,6 @@ R"(#version 150 core
     uniform vec2 RenderSize;
     uniform vec2 RenderOrigin;
     uniform float TextShiftFactor;
-    in int vPosition[1];
     in int vCharacter[1];
     in float vOffset[1];
     void main()

include/igl/orientable_patches.cpp

@@ -102,4 +102,6 @@ IGL_INLINE void igl::orientable_patches(
 // Explicit template instantiation
 template void igl::orientable_patches<Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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, 1, 0, -1, 1> >&, Eigen::SparseMatrix<int, 0, int>&);
 template void igl::orientable_patches<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::orientable_patches<Eigen::Matrix<int, -1, 3, 0, -1, 3>, 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::SparseMatrix<int, 0, int>&);
+template void igl::orientable_patches<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/per_vertex_point_to_plane_quadrics.cpp

@@ -143,7 +143,7 @@ IGL_INLINE void igl::per_vertex_point_to_plane_quadrics(
       assert(N.rows() == ev.size()-2);
       Eigen::MatrixXd S(N.rows()+1,ev.size());
       S<<ev,N;
-      Quadric boundary_edge_quadric = subspace_quadric(p,S,length);
+      Quadric boundary_edge_quadric = subspace_quadric(p,S,length*length);
       for(int c = 0;c<3;c++)
       {
         if(c != infinite_corner)

include/igl/readPLY.cpp

@@ -6,6 +6,7 @@
 #include <Eigen/Core>
 
 #include "tinyply.h"
+#include "file_utils.h"
 
 
 namespace igl
@@ -137,54 +138,6 @@ IGL_INLINE bool tinyply_tristrips_to_faces(
   }
 }
 
-
-inline void read_file_binary(FILE *fp, std::vector<uint8_t>& fileBufferBytes)
-{
-    if (!ferror(fp))
-    {
-        fseek(fp,0,SEEK_END);
-        size_t sizeBytes = ftell(fp);
-        fseek(fp,0,SEEK_SET);
-        fileBufferBytes.resize(sizeBytes);
-
-        if(fread((char*)fileBufferBytes.data(), 1, sizeBytes,fp)==sizeBytes)
-          return;
-    }
-
-    throw std::runtime_error("error reading from file");
-}
-
-struct ply_memory_buffer : public std::streambuf
-{
-    char * p_start {nullptr};
-    char * p_end {nullptr};
-    size_t size;
-
-    ply_memory_buffer(char const * first_elem, size_t size)
-        : p_start(const_cast<char*>(first_elem)), p_end(p_start + size), size(size)
-    {
-        setg(p_start, p_start, p_end);
-    }
-
-    pos_type seekoff(off_type off, std::ios_base::seekdir dir, std::ios_base::openmode which) override
-    {
-        if (dir == std::ios_base::cur) gbump(static_cast<int>(off));
-        else setg(p_start, (dir == std::ios_base::beg ? p_start : p_end) + off, p_end);
-        return gptr() - p_start;
-    }
-
-    pos_type seekpos(pos_type pos, std::ios_base::openmode which) override
-    {
-        return seekoff(pos, std::ios_base::beg, which);
-    }
-};
-
-struct ply_memory_stream : virtual ply_memory_buffer, public std::istream
-{
-    ply_memory_stream(char const * first_elem, size_t size)
-        : ply_memory_buffer(first_elem, size), std::istream(static_cast<std::streambuf*>(this)) {}
-};
-
 template <
   typename DerivedV,
   typename DerivedF,
@@ -218,7 +171,7 @@ IGL_INLINE bool readPLY(
   {
     std::vector<uint8_t> fileBufferBytes;
     read_file_binary(fp,fileBufferBytes);
-    ply_memory_stream stream((char*)fileBufferBytes.data(), fileBufferBytes.size());
+    file_memory_stream stream((char*)fileBufferBytes.data(), fileBufferBytes.size());
     return readPLY(stream,V,F,E,N,UV,VD,Vheader,FD,Fheader,ED,Eheader,comments);
   }
   catch(const std::exception& e)

include/igl/readSTL.cpp

@@ -1,290 +1,324 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2014 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// Copyright (C) 2018 Qingnan Zhou <[email protected]>
+// Copyright (C) 2020 Jérémie Dumas <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "readSTL.h"
 #include "list_to_matrix.h"
+#include "string_utils.h"
+#include "file_utils.h"
 
 #include <iostream>
+
+namespace igl {
+
 template <typename DerivedV, typename DerivedF, typename DerivedN>
-IGL_INLINE bool igl::readSTL(
-  const std::string & filename,
-  Eigen::PlainObjectBase<DerivedV> & V,
-  Eigen::PlainObjectBase<DerivedF> & F,
-  Eigen::PlainObjectBase<DerivedN> & N)
-{
-  using namespace std;
-  vector<vector<typename DerivedV::Scalar> > vV;
-  vector<vector<typename DerivedN::Scalar> > vN;
-  vector<vector<typename DerivedF::Scalar> > vF;
-  if(!readSTL(filename,vV,vF,vN))
-  {
+IGL_INLINE bool readSTL(std::istream &input,
+                        Eigen::PlainObjectBase<DerivedV> &V,
+                        Eigen::PlainObjectBase<DerivedF> &F,
+                        Eigen::PlainObjectBase<DerivedN> &N) {
+  std::vector<std::array<typename DerivedV::Scalar, 3>> vV;
+  std::vector<std::array<typename DerivedN::Scalar, 3>> vN;
+  std::vector<std::array<typename DerivedF::Scalar, 3>> vF;
+  if (!readSTL(input, vV, vF, vN)) {
     return false;
   }
 
-  if(!list_to_matrix(vV,V))
-  {
+  if (!list_to_matrix(vV, V)) {
     return false;
   }
 
-  if(!list_to_matrix(vF,F))
-  {
+  if (!list_to_matrix(vF, F)) {
     return false;
   }
 
-  if(!list_to_matrix(vN,N))
-  {
+  if (!list_to_matrix(vN, N)) {
     return false;
   }
   return true;
 }
 
-template <typename TypeV, typename TypeF, typename TypeN>
-IGL_INLINE bool igl::readSTL(
-  const std::string & filename,
-  std::vector<std::vector<TypeV> > & V,
-  std::vector<std::vector<TypeF> > & F,
-  std::vector<std::vector<TypeN> > & N)
-{
-  using namespace std;
-  // Should test for ascii
-
-  // Open file, and check for error
-  FILE * stl_file = fopen(filename.c_str(),"rb");
-  if(NULL==stl_file)
-  {
-    fprintf(stderr,"IOError: %s could not be opened...\n",
-            filename.c_str());
-    return false;
-  }
-  return readSTL(stl_file,V,F,N);
-}
+IGL_INLINE bool is_stl_binary(std::istream &input) {
+  std::streampos start_pos = input.tellg();
 
-template <typename TypeV, typename TypeF, typename TypeN>
-IGL_INLINE bool igl::readSTL(
-  FILE * stl_file, 
-  std::vector<std::vector<TypeV> > & V,
-  std::vector<std::vector<TypeF> > & F,
-  std::vector<std::vector<TypeN> > & N)
-{
-  using namespace std;
-  //stl_file = freopen(NULL,"rb",stl_file);
-  if(NULL==stl_file)
-  {
-    fprintf(stderr,"IOError: stl file could not be reopened as binary (1) ...\n");
+  constexpr size_t HEADER_SIZE = 80;
+  char header[HEADER_SIZE];
+  input.read(header, HEADER_SIZE);
+  if (!starts_with(header, "solid")) {
+    input.seekg(start_pos);
+    return true;
+  }
+  if (!input.good()) {
+    input.seekg(start_pos);
     return false;
   }
 
-  V.clear();
-  F.clear();
-  N.clear();
+  // Check if filesize matches the number of faces claimed.
+  char buf[4];
+  input.read(buf, 4);
+  size_t num_faces = *reinterpret_cast<uint32_t *>(buf);
+  input.seekg(0, input.end);
+  size_t file_size = input.tellg();
 
+  input.seekg(start_pos);
 
-  // Specifically 80 character header
-  char header[80];
-  char solid[80] = {0};
-  bool is_ascii = true;
-  if(fread(header,1,80,stl_file) != 80)
-  {
-    cerr<<"IOError: too short (1)."<<endl;
-    goto close_false;
+  if (file_size == 80 + 4 + (4 * 12 + 2) * num_faces) {
+    return true;
+  } else {
+    return false;
   }
-  sscanf(header,"%79s",solid);
-  if(string("solid") != solid)
-  {
-    // definitely **not** ascii 
-    is_ascii = false;
-  }else
-  {
-    // might still be binary
-    char buf[4];
-    if(fread(buf,1,4,stl_file) != 4)
-    {
-      cerr<<"IOError: too short (3)."<<endl;
-      goto close_false;
-    }
-    size_t num_faces = *reinterpret_cast<unsigned int*>(buf);
-    fseek(stl_file,0,SEEK_END);
-    int file_size = ftell(stl_file);
-    if(file_size == 80 + 4 + (4*12 + 2) * num_faces)
-    {
-      is_ascii = false;
-    }else
-    {
-      is_ascii = true;
-    }
+}
+
+template <typename TypeV, typename TypeF, typename TypeN>
+IGL_INLINE bool read_stl_ascii(std::istream &input,
+                               std::vector<std::array<TypeV, 3>> &V,
+                               std::vector<std::array<TypeF, 3>> &F,
+                               std::vector<std::array<TypeN, 3>> &N) {
+  constexpr size_t LINE_SIZE = 256;
+  char line[LINE_SIZE];
+  bool success = true;
+
+  if (!input) {
+    throw std::runtime_error("Failed to open file");
   }
 
-  if(is_ascii)
-  {
-    // Rewind to end of header
-    //stl_file = fopen(filename.c_str(),"r");
-    //stl_file = freopen(NULL,"r",stl_file);
-    fseek(stl_file, 0, SEEK_SET);
-    if(NULL==stl_file)
-    {
-      fprintf(stderr,"IOError: stl file could not be reopened as ascii ...\n");
+  // skip header line.
+  input.getline(line, LINE_SIZE);
+
+  auto parse_ascii_normal = [&N](const char *line) {
+    double x, y, z;
+    size_t n = sscanf(line, " facet normal %lf %lf %lf", &x, &y, &z);
+    assert(n == 3);
+    if (n != 3) {
       return false;
     }
-    // Read 80 header
-    // Eat file name
-#ifndef IGL_LINE_MAX
-#  define IGL_LINE_MAX 2048
-#endif
-    char name[IGL_LINE_MAX];
-    if(NULL==fgets(name,IGL_LINE_MAX,stl_file))
-    {
-      cerr<<"IOError: ascii too short (2)."<<endl;
-      goto close_false;
+
+    N.push_back({{static_cast<TypeN>(x), static_cast<TypeN>(y),
+                  static_cast<TypeN>(z)}});
+    return true;
+  };
+
+  auto parse_ascii_vertex = [&V](const char *line) {
+    double x, y, z;
+    size_t n = sscanf(line, " vertex %lf %lf %lf", &x, &y, &z);
+    assert(n == 3);
+    if (n != 3) {
+      return false;
     }
-    // ascii
-    while(true)
-    {
-      int ret;
-      char facet[IGL_LINE_MAX],normal[IGL_LINE_MAX];
-      vector<TypeN > n(3);
-      double nd[3];
-      ret = fscanf(stl_file,"%s %s %lg %lg %lg",facet,normal,nd,nd+1,nd+2);
-      if(string("endsolid") == facet)
-      {
-        break;
-      }
-      if(ret != 5 || 
-          !(string("facet") == facet || 
-          string("faced") == facet) ||
-          string("normal") != normal)
-      {
-        cout<<"facet: "<<facet<<endl;
-        cout<<"normal: "<<normal<<endl;
-        cerr<<"IOError: bad format (1)."<<endl;
-        goto close_false;
-      }
-      // copy casts to Type
-      n[0] = nd[0]; n[1] = nd[1]; n[2] = nd[2];
-      N.push_back(n);
-      char outer[IGL_LINE_MAX], loop[IGL_LINE_MAX];
-      ret = fscanf(stl_file,"%s %s",outer,loop);
-      if(ret != 2 || string("outer") != outer || string("loop") != loop)
-      {
-        cerr<<"IOError: bad format (2)."<<endl;
-        goto close_false;
+
+    V.push_back({{static_cast<TypeV>(x), static_cast<TypeV>(y),
+                  static_cast<TypeV>(z)}});
+    return true;
+  };
+
+  auto parse_ascii_facet = [&parse_ascii_vertex, &parse_ascii_normal](std::istream &fin) {
+    constexpr size_t LINE_SIZE = 256;
+    constexpr size_t WORD_SIZE = 128;
+    char line[LINE_SIZE];
+    char first_word[WORD_SIZE];
+    const char *face_begin = "facet";
+    const char *face_end = "endfacet";
+    const char *loop_begin = "outer";
+    const char *loop_end = "endloop";
+    const char *vertex_flag = "vertex";
+
+    bool reading_facet = false;
+    bool reading_loop = false;
+    bool success = true;
+    size_t num_vts = 0;
+    while (!fin.eof()) {
+      fin.getline(line, LINE_SIZE);
+      size_t n = sscanf(line, " %s", first_word);
+      if (n == 0)
+        continue;
+      if (starts_with(first_word, face_begin)) {
+        success = parse_ascii_normal(line);
+        assert(success);
+        reading_facet = true;
+      } else if (starts_with(first_word, face_end)) {
+        assert(reading_facet);
+        reading_facet = false;
+      } else if (starts_with(first_word, loop_begin)) {
+        reading_loop = true;
+      } else if (starts_with(first_word, loop_end)) {
+        assert(reading_loop);
+        reading_loop = false;
+      } else if (starts_with(first_word, vertex_flag)) {
+        assert(reading_facet);
+        assert(reading_loop);
+        success = parse_ascii_vertex(line);
+        assert(success);
+        num_vts += 1;
       }
-      vector<TypeF> f;
-      while(true)
-      {
-        char word[IGL_LINE_MAX];
-        int ret = fscanf(stl_file,"%s",word);
-        if(ret == 1 && string("endloop") == word)
-        {
-          break;
-        }else if(ret == 1 && string("vertex") == word)
-        {
-          vector<TypeV> v(3);
-          double vd[3];
-          int ret = fscanf(stl_file,"%lg %lg %lg",vd,vd+1,vd+2);
-          if(ret != 3)
-          {
-            cerr<<"IOError: bad format (3)."<<endl;
-            goto close_false;
-          }
-          f.push_back(V.size());
-          // copy casts to Type
-          v[0] = vd[0]; v[1] = vd[1]; v[2] = vd[2];
-          V.push_back(v);
-        }else
-        {
-          cerr<<"IOError: bad format (4)."<<endl;
-          goto close_false;
-        }
+      if (!success) {
+        return false;
       }
-      F.push_back(f);
-      char endfacet[IGL_LINE_MAX];
-      ret = fscanf(stl_file,"%s",endfacet);
-      if(ret != 1 || string("endfacet") != endfacet)
-      {
-        cerr<<"IOError: bad format (5)."<<endl;
-        goto close_false;
+      if (!reading_facet) {
+        break;
       }
     }
-    // read endfacet
-    goto close_true;
-  }else
-  {
-    // Binary
-    //stl_file = freopen(NULL,"rb",stl_file);
-    fseek(stl_file, 0, SEEK_SET);
-    // Read 80 header
-    char header[80];
-    if(fread(header,sizeof(char),80,stl_file)!=80)
-    {
-      cerr<<"IOError: bad format (6)."<<endl;
-      goto close_false;
+    if (num_vts == 0) {
+      return true;
     }
-    // Read number of triangles
-    unsigned int num_tri;
-    if(fread(&num_tri,sizeof(unsigned int),1,stl_file)!=1)
-    {
-      cerr<<"IOError: bad format (7)."<<endl;
-      goto close_false;
+    assert(num_vts == 3);
+    if (num_vts != 3) {
+      std::cerr << "Warning: mesh contain face made of " << num_vts
+                << " vertices" << std::endl;
+      return false;
     }
-    V.resize(num_tri*3,vector<TypeV >(3,0));
-    N.resize(num_tri,vector<TypeN >(3,0));
-    F.resize(num_tri,vector<TypeF >(3,0));
-    for(int t = 0;t<(int)num_tri;t++)
-    {
-      // Read normal
-      float n[3];
-      if(fread(n,sizeof(float),3,stl_file)!=3)
-      {
-        cerr<<"IOError: bad format (8)."<<endl;
-        goto close_false;
-      }
-      // Read each vertex
-      for(int c = 0;c<3;c++)
-      {
-        F[t][c] = 3*t+c;
-        N[t][c] = n[c];
-        float v[3];
-        if(fread(v,sizeof(float),3,stl_file)!=3)
-        {
-          cerr<<"IOError: bad format (9)."<<endl;
-          goto close_false;
-        }
-        V[3*t+c][0] = v[0];
-        V[3*t+c][1] = v[1];
-        V[3*t+c][2] = v[2];
-      }
-      // Read attribute size
-      unsigned short att_count;
-      if(fread(&att_count,sizeof(unsigned short),1,stl_file)!=1)
-      {
-        cerr<<"IOError: bad format (10)."<<endl;
-        goto close_false;
+    return true;
+  };
+
+  while (!input.eof()) {
+    success = parse_ascii_facet(input);
+    if (!success) {
+      return false;
+    }
+  }
+
+  F.resize(V.size() / 3);
+    for (size_t f = 0; f < F.size(); ++f) {
+    auto v = static_cast<TypeF>(f * 3);
+    F[f] = {{v, v + 1, v + 2}};
+  }
+  return success;
+}
+
+template <typename TypeV, typename TypeF, typename TypeN>
+IGL_INLINE bool read_stl_binary(std::istream &input,
+                                std::vector<std::array<TypeV, 3>> &V,
+                                std::vector<std::array<TypeF, 3>> &F,
+                                std::vector<std::array<TypeN, 3>> &N) {
+  if (!input) {
+    throw std::runtime_error("Failed to open file");
+  }
+
+  constexpr size_t FLOAT_SIZE = sizeof(float);
+  static_assert(FLOAT_SIZE == 4, "float type is not 4 bytes");
+  constexpr size_t LINE_SIZE = 256;
+  char buf[LINE_SIZE];
+
+  // 80 bytes header, no data significance.
+  input.read(buf, 80);
+  if (!input.good()) {
+    throw std::runtime_error("Unable to parse STL header.");
+  }
+
+  input.read(buf, 4);
+  const size_t num_faces = *reinterpret_cast<uint32_t *>(buf);
+  if (!input.good()) {
+    throw std::runtime_error("Unable to parse STL number of faces.");
+  }
+
+  for (size_t i = 0; i < num_faces; i++) {
+    // Parse normal
+    input.read(buf, FLOAT_SIZE * 3);
+    auto nx = static_cast<TypeN>(*reinterpret_cast<float *>(buf));
+    auto ny = static_cast<TypeN>(*reinterpret_cast<float *>(buf + FLOAT_SIZE));
+    auto nz =
+        static_cast<TypeN>(*reinterpret_cast<float *>(buf + FLOAT_SIZE * 2));
+    assert(input.good());
+
+    // vertex 1
+    input.read(buf, FLOAT_SIZE * 3);
+    auto v1x = static_cast<TypeV>(*reinterpret_cast<float *>(buf));
+    auto v1y = static_cast<TypeV>(*reinterpret_cast<float *>(buf + FLOAT_SIZE));
+    auto v1z =
+        static_cast<TypeV>(*reinterpret_cast<float *>(buf + FLOAT_SIZE * 2));
+    assert(input.good());
+
+    // vertex 2
+    input.read(buf, FLOAT_SIZE * 3);
+    auto v2x = static_cast<TypeV>(*reinterpret_cast<float *>(buf));
+    auto v2y = static_cast<TypeV>(*reinterpret_cast<float *>(buf + FLOAT_SIZE));
+    auto v2z =
+        static_cast<TypeV>(*reinterpret_cast<float *>(buf + FLOAT_SIZE * 2));
+    assert(input.good());
+
+    // vertex 3
+    input.read(buf, FLOAT_SIZE * 3);
+    auto v3x = static_cast<TypeV>(*reinterpret_cast<float *>(buf));
+    auto v3y = static_cast<TypeV>(*reinterpret_cast<float *>(buf + FLOAT_SIZE));
+    auto v3z =
+        static_cast<TypeV>(*reinterpret_cast<float *>(buf + FLOAT_SIZE * 2));
+    assert(input.good());
+
+    // attribute (2 bytes), not sure what purpose they serve.
+    input.read(buf, 2);
+
+    N.push_back({{nx, ny, nz}});
+    V.push_back({{v1x, v1y, v1z}});
+    V.push_back({{v2x, v2y, v2z}});
+    V.push_back({{v3x, v3y, v3z}});
+
+    assert(input.good());
+    if (!input.good()) {
+      std::stringstream err_msg;
+      err_msg << "Failed to parse face " << i << " from STL file";
+      throw std::runtime_error(err_msg.str());
+    }
+  }
+  std::for_each(V.begin(), V.end(), [](const std::array<TypeV, 3> &v) {
+    for (auto x : v) {
+      if (!std::isfinite(x)) {
+        throw std::runtime_error("NaN or Inf detected in input file.");
       }
     }
-    goto close_true;
+  });
+
+  if (!V.empty()) {
+    F.resize(V.size() / 3);
+    for (size_t f = 0; f < F.size(); ++f) {
+      auto v = static_cast<TypeF>(f * 3);
+      F[f] = {{v, v + 1, v + 2}};
+    }
   }
-close_false:
-  fclose(stl_file);
-  return false;
-close_true:
-  fclose(stl_file);
+
   return true;
 }
 
+template <typename TypeV, typename TypeF, typename TypeN>
+IGL_INLINE bool readSTL(std::istream &input,
+                        std::vector<std::array<TypeV, 3>> &V,
+                        std::vector<std::array<TypeF, 3>> &F,
+                        std::vector<std::array<TypeN, 3>> &N) {
+  bool success = false;
+  if (is_stl_binary(input)) {
+    success = read_stl_binary(input, V, F, N);
+  } else {
+    success = read_stl_ascii(input, V, F, N);
+  }
+  return success;
+}
+
+template <typename DerivedV, typename DerivedF, typename DerivedN>
+IGL_INLINE bool readSTL(
+  FILE * fp,
+  Eigen::PlainObjectBase<DerivedV> & V,
+  Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedN> & N)
+{
+  std::vector<uint8_t> fileBufferBytes;
+  read_file_binary(fp,fileBufferBytes);
+  file_memory_stream stream((char*)fileBufferBytes.data(), fileBufferBytes.size());
+  return readSTL(stream, V, F, N);
+}
+
+} // namespace igl
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
-template bool igl::readSTL<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template bool igl::readSTL<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template bool igl::readSTL<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template bool igl::readSTL<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > 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<double, -1, -1, 0, -1, -1> >&);
-template bool igl::readSTL<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
-template bool igl::readSTL<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::readSTL<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(FILE*, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::readSTL<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(FILE*, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::readSTL<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(FILE*, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::readSTL<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(FILE*, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::readSTL<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(FILE*, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::readSTL<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(FILE*, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::readSTL<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(FILE*, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::readSTL<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(FILE*, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::readSTL<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(FILE*, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/readSTL.h

@@ -1,9 +1,10 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// Copyright (C) 2020 Jérémie Dumas <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_READSTL_H
 #define IGL_READSTL_H
@@ -13,10 +14,11 @@
 #  include <Eigen/Core>
 #endif
 #include <string>
-#include <cstdio>
+#include <istream>
 #include <vector>
+#include <array>
 
-namespace igl 
+namespace igl
 {
   // Read a mesh from an ascii/binary stl file.
   //
@@ -38,26 +40,24 @@ namespace igl
   //   writeOBJ("Downloads/cat.obj",V,F);
   template <typename DerivedV, typename DerivedF, typename DerivedN>
   IGL_INLINE bool readSTL(
-    const std::string & filename,
+    std::istream &input,
     Eigen::PlainObjectBase<DerivedV> & V,
     Eigen::PlainObjectBase<DerivedF> & F,
     Eigen::PlainObjectBase<DerivedN> & N);
-  // Inputs:
-  //   stl_file  pointer to already opened .stl file 
-  // Outputs:
-  //   stl_file  closed file
+
   template <typename TypeV, typename TypeF, typename TypeN>
   IGL_INLINE bool readSTL(
-    FILE * stl_file, 
-    std::vector<std::vector<TypeV> > & V,
-    std::vector<std::vector<TypeF> > & F,
-    std::vector<std::vector<TypeN> > & N);
-  template <typename TypeV, typename TypeF, typename TypeN>
+    std::istream &input,
+    std::vector<std::array<TypeV, 3> > & V,
+    std::vector<std::array<TypeF, 3> > & F,
+    std::vector<std::array<TypeN, 3> > & N);
+
+  template <typename DerivedV, typename DerivedF, typename DerivedN>
   IGL_INLINE bool readSTL(
-    const std::string & filename,
-    std::vector<std::vector<TypeV> > & V,
-    std::vector<std::vector<TypeF> > & F,
-    std::vector<std::vector<TypeN> > & N);
+    FILE * fp,
+    Eigen::PlainObjectBase<DerivedV> & V,
+    Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedN> & N);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/read_triangle_mesh.cpp

@@ -103,6 +103,7 @@ IGL_INLINE bool igl::read_triangle_mesh(
 {
   using namespace std;
   using namespace Eigen;
+  Eigen::MatrixXd N;
   vector<vector<double > > vV,vN,vTC,vC;
   vector<vector<int > > vF,vFTC,vFN;
   vector<tuple<string, int, int>> FM;
@@ -140,10 +141,10 @@ IGL_INLINE bool igl::read_triangle_mesh(
   }else if(ext == "ply")
   {
     return readPLY(fp, V, F);
-    
+
   }else if(ext == "stl")
   {
-    if(!readSTL(fp,vV,vF,vN))
+    if(!readSTL(fp,V,F,N))
     {
       return false;
     }

include/igl/repdiag.cpp

@@ -71,8 +71,7 @@ IGL_INLINE void igl::repdiag(
 {
   int m = A.rows();
   int n = A.cols();
-  B.resize(m*d,n*d);
-  B.array() *= 0;
+  B.setZero(m*d,n*d);
   for(int i = 0;i<d;i++)
   {
     B.block(i*m,i*n,m,n) = A;

include/igl/signed_distance.cpp

@@ -13,6 +13,7 @@
 #include "per_vertex_normals.h"
 #include "point_mesh_squared_distance.h"
 #include "pseudonormal_test.h"
+#include "fast_winding_number.h"
 
 
 template <
@@ -37,10 +38,17 @@ IGL_INLINE void igl::signed_distance(
 {
   using namespace Eigen;
   using namespace std;
+
   const int dim = V.cols();
+
   assert((V.cols() == 3||V.cols() == 2) && "V should have 3d or 2d positions");
   assert((P.cols() == 3||P.cols() == 2) && "P should have 3d or 2d positions");
   assert(V.cols() == P.cols() && "V should have same dimension as P");
+
+  if (sign_type == SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER){
+    assert(V.cols() == 3 && "V should be 3D for fast winding number");
+  }
+
   // Only unsigned distance is supported for non-triangles
   if(sign_type != SIGNED_DISTANCE_TYPE_UNSIGNED)
   {
@@ -67,6 +75,9 @@ IGL_INLINE void igl::signed_distance(
   Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> E;
   Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> EMAP;
   WindingNumberAABB<RowVector3S,DerivedV,DerivedF> hier3;
+  igl::FastWindingNumberBVH fwn_bvh;
+  Eigen::VectorXf W;
+
   switch(sign_type)
   {
     default:
@@ -88,6 +99,10 @@ IGL_INLINE void igl::signed_distance(
           break;
       }
       break;
+    case SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER:
+      // assert above ensures dim == 3
+      igl::fast_winding_number(V.template cast<float>().eval(), F, 2, fwn_bvh);
+
     case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
       switch(dim)
       {
@@ -187,6 +202,14 @@ IGL_INLINE void igl::signed_distance(
           }
           break;
         }
+        case SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER:
+        {
+          //assert above ensured 3D
+          Scalar w = fast_winding_number(fwn_bvh, 2, q3.template cast<float>().eval());         
+          s = 1.-2.*std::abs(w);  
+
+          break;
+        }
         case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
         {
           RowVector3S n3;
@@ -439,7 +462,6 @@ IGL_INLINE void igl::signed_distance_winding_number(
   using namespace std;
   typedef Eigen::Matrix<typename DerivedV::Scalar,1,2> RowVector2S;
   sqrd = tree.squared_distance(V,F,RowVector2S(q),i,(RowVector2S&)c);
-  Scalar w;
   // TODO: using .data() like this is very dangerous... This is assuming
   // colmajor order
   assert(!V.derived().IsRowMajor);
@@ -447,6 +469,55 @@ IGL_INLINE void igl::signed_distance_winding_number(
   s = 1.-2.*winding_number(V,F,q);
 }
 
+//Multi point by parrallel for on single point
+template <
+  typename DerivedP,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedS>
+IGL_INLINE void igl::signed_distance_fast_winding_number(
+    const Eigen::MatrixBase<DerivedP> & P,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const AABB<DerivedV,3> & tree,
+    const igl::FastWindingNumberBVH & fwn_bvh,
+    Eigen::PlainObjectBase<DerivedS> & S)
+  {
+    typedef Eigen::Matrix<typename DerivedV::Scalar,1,3> RowVector3S;
+    S.resize(P.rows(),1);
+    int min_parallel = 10000; 
+    parallel_for(P.rows(), [&](const int p)
+    {
+      RowVector3S q;
+      q.head(P.row(p).size()) = P.row(p);
+      // get sdf for single point, update result matrix
+      S(p) = signed_distance_fast_winding_number(q, V, F, tree,fwn_bvh);
+    }
+    ,min_parallel);  
+  }
+
+//Single Point
+template <
+  typename Derivedq,
+  typename DerivedV,
+  typename DerivedF>
+IGL_INLINE typename DerivedV::Scalar igl::signed_distance_fast_winding_number(
+    const Eigen::MatrixBase<Derivedq> & q,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const AABB<DerivedV,3> & tree,
+    const igl::FastWindingNumberBVH & fwn_bvh)
+  {
+    typedef typename DerivedV::Scalar Scalar;
+    Scalar s,sqrd;
+    Eigen::Matrix<Scalar,1,3> c;
+    int i = -1;
+    sqrd = tree.squared_distance(V,F,q,i,c);
+    Scalar w = fast_winding_number(fwn_bvh,2,q.template cast<float>());
+    //0.5 is on surface
+    return sqrt(sqrd)*(1.-2.*std::abs(w));
+  }
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
@@ -463,4 +534,5 @@ template void igl::signed_distance_pseudonormal<Eigen::Matrix<double, -1, -1, 0,
 template void igl::signed_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::SignedDistanceType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::signed_distance_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, double, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::WindingNumberAABB<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, double&, double&, int&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 template Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar igl::signed_distance_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::WindingNumberAABB<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&);
+template void igl::signed_distance_fast_winding_number<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(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<int, -1, -1, 0, -1, -1> > const&, igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3> const&, igl::FastWindingNumberBVH const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/signed_distance.h

@@ -11,6 +11,7 @@
 #include "igl_inline.h"
 #include "AABB.h"
 #include "WindingNumberAABB.h"
+#include "fast_winding_number.h"
 #include <Eigen/Core>
 #include <vector>
 namespace igl
@@ -18,11 +19,14 @@ namespace igl
   enum SignedDistanceType
   {
     // Use fast pseudo-normal test [Bærentzen & Aanæs 2005]
-    SIGNED_DISTANCE_TYPE_PSEUDONORMAL   = 0,
-    SIGNED_DISTANCE_TYPE_WINDING_NUMBER = 1,
-    SIGNED_DISTANCE_TYPE_DEFAULT        = 2,
-    SIGNED_DISTANCE_TYPE_UNSIGNED       = 3,
-    NUM_SIGNED_DISTANCE_TYPE            = 4
+    SIGNED_DISTANCE_TYPE_PSEUDONORMAL         = 0,
+    // Use winding number [Jacobson, Kavan Sorking-Hornug 2013]
+    SIGNED_DISTANCE_TYPE_WINDING_NUMBER       = 1,
+    SIGNED_DISTANCE_TYPE_DEFAULT              = 2,
+    SIGNED_DISTANCE_TYPE_UNSIGNED             = 3,
+    // Use Fast winding number [Barill, Dickson, Schmidt, Levin, Jacobson 2018]
+    SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER  = 4,
+    NUM_SIGNED_DISTANCE_TYPE                  = 5
   };
   // Computes signed distance to a mesh
   //
@@ -62,7 +66,22 @@ namespace igl
     Eigen::PlainObjectBase<DerivedI> & I,
     Eigen::PlainObjectBase<DerivedC> & C,
     Eigen::PlainObjectBase<DerivedN> & N);
-  // Default bounds
+  // Computes signed distance to a mesh, with default bounds
+  //
+  // Inputs:
+  //   P  #P by 3 list of query point positions
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by ss list of triangle indices, ss should be 3 unless sign_type ==
+  //     SIGNED_DISTANCE_TYPE_UNSIGNED
+  //   sign_type  method for computing distance _sign_ S
+  //   lower_bound  lower bound of distances needed {std::numeric_limits::min}
+  //   upper_bound  lower bound of distances needed {std::numeric_limits::max}
+  // Outputs:
+  //   S  #P list of smallest signed distances
+  //   I  #P list of facet indices corresponding to smallest distances
+  //   C  #P by 3 list of closest points
+  //   N  #P by 3 list of closest normals (only set if
+  //     sign_type=SIGNED_DISTANCE_TYPE_PSEUDONORMAL)
   template <
     typename DerivedP,
     typename DerivedV,
@@ -80,7 +99,7 @@ namespace igl
     Eigen::PlainObjectBase<DerivedI> & I,
     Eigen::PlainObjectBase<DerivedC> & C,
     Eigen::PlainObjectBase<DerivedN> & N);
-  // Computes signed distance to mesh
+  // Computes signed distance to mesh using pseudonormal with precomputed AABB tree and edge/vertice normals
   //
   // Inputs:
   //   tree  AABB acceleration tree (see AABB.h)
@@ -236,6 +255,63 @@ namespace igl
     Scalar & sqrd,
     int & i,
     Eigen::PlainObjectBase<Derivedc> & c);
+
+
+  // Calculates signed distance at query points P, using fast winding number
+  //   for sign.
+  //
+  // Usage:
+  //     VectorXd S;  
+  //     VectorXd V, P; //where V is mesh vertices, P are query points
+  //     VectorXi F;  
+  //     igl::FastWindingNumberBVH fwn_bvh;
+  //     igl::fast_winding_number(V.cast<float>(), F, 2, fwn_bvh);
+  //     igl::signed_distance_fast_winding_number(P,V,F,tree,fwn_bvh,S);
+  //
+  // Inputs:
+  //   P  #P by 3 list of query point positions
+  //   V  #V by 3 list of triangle indices
+  //   F  #F by 3 list of triangle normals 
+  //   tree  AABB acceleration tree (see AABB.h)
+  //   bvh fast winding precomputation (see Fast_Winding_Number.h)   
+  // Outputs:
+  //   S  #P list of signed distances of each point in P
+  template <
+    typename DerivedP,
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedS>
+  IGL_INLINE void signed_distance_fast_winding_number(
+    const Eigen::MatrixBase<DerivedP> & P,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const AABB<DerivedV,3> & tree,
+    const igl::FastWindingNumberBVH & fwn_bvh,
+    Eigen::PlainObjectBase<DerivedS> & S
+  );
+
+  // Calculates signed distance at query point q, using fast winding number
+  //   for sign.
+  //
+  // Inputs:
+  //   tree  AABB acceleration tree (see AABB.h)
+  //   V  #V by 3 list of triangle indices
+  //   F  #F by 3 list of triangle normals 
+  //   bvh fast winding precomputation (see Fast_Winding_Number.h)   
+  //   q  1 by 3 list of query point positions
+  // Outputs:
+  //   S  #P list of signed distances of each point in P
+  template <
+    typename Derivedq,
+    typename DerivedV,
+    typename DerivedF>
+  IGL_INLINE typename DerivedV::Scalar signed_distance_fast_winding_number(
+    const Eigen::MatrixBase<Derivedq> & q,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const AABB<DerivedV,3> & tree,
+    const igl::FastWindingNumberBVH & fwn_bvh
+  );
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/string_utils.cpp

@@ -0,0 +1,22 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Jérémie Dumas <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "string_utils.h"
+
+#include <cstring>
+
+namespace igl {
+
+IGL_INLINE bool starts_with(const std::string &str, const std::string &prefix) {
+  return (str.rfind(prefix, 0) == 0);
+}
+
+IGL_INLINE bool starts_with(const char *str, const char* prefix) {
+  return strncmp(prefix, str, strlen(prefix)) == 0;
+}
+
+}

include/igl/string_utils.h

@@ -0,0 +1,27 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2020 Jérémie Dumas <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_STRING_UTILS_H
+#define IGL_STRING_UTILS_H
+
+#include "igl_inline.h"
+
+#include <string>
+
+namespace igl {
+
+IGL_INLINE bool starts_with(const std::string &str, const std::string &prefix);
+
+IGL_INLINE bool starts_with(const char *str, const char* prefix);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "string_utils.cpp"
+#endif
+
+#endif

include/igl/xml/writeDAE.cpp

@@ -19,9 +19,8 @@ IGL_INLINE bool igl::xml::writeDAE(
 {
   using namespace std;
   using namespace Eigen;
-  using namespace tinyxml2;
 
-  XMLDocument* doc = new XMLDocument();
+  tinyxml2::XMLDocument* doc = new tinyxml2::XMLDocument();
 
   const auto & ele = [&doc](
     const std::string tag,
@@ -30,11 +29,11 @@ IGL_INLINE bool igl::xml::writeDAE(
     const std::map<std::string,std::string> attribs =
       std::map<std::string,std::string>(),
     const std::string text="",
-    const std::list<XMLElement *> children =
-       std::list<XMLElement *>()
-    )->XMLElement *
+    const std::list<tinyxml2::XMLElement *> children =
+       std::list<tinyxml2::XMLElement *>()
+    )->tinyxml2::XMLElement *
   {
-    XMLElement * element = doc->NewElement(tag.c_str());
+    tinyxml2::XMLElement * element = doc->NewElement(tag.c_str());
     for(const auto & key_value :  attribs)
     {
       element->SetAttribute(key_value.first.c_str(),key_value.second.c_str());
@@ -126,9 +125,9 @@ IGL_INLINE bool igl::xml::writeDAE(
     }));
   // tinyxml2 seems **not** to print the <?xml ...?> header by default, but it
   // also seems that that's OK
-  XMLError error = doc->SaveFile(filename.c_str());
+  tinyxml2::XMLError error = doc->SaveFile(filename.c_str());
   bool ret = true;
-  if(error != XML_NO_ERROR)
+  if(error != tinyxml2::XML_SUCCESS)
   {
     doc->PrintError();
     ret = false;

tests/include/igl/signed_distance.cpp

@@ -15,6 +15,7 @@ TEST_CASE("signed_distance: single_tet", "[igl]")
     0,2,1,
     0,3,2,
     1,2,3;
+
   Eigen::MatrixXd P(1,3);
   P<<0.5,0.5,0.5;
   for(const igl::SignedDistanceType type :
@@ -22,14 +23,23 @@ TEST_CASE("signed_distance: single_tet", "[igl]")
       igl::SIGNED_DISTANCE_TYPE_PSEUDONORMAL  ,
       igl::SIGNED_DISTANCE_TYPE_WINDING_NUMBER,
       igl::SIGNED_DISTANCE_TYPE_DEFAULT       ,
-      igl::SIGNED_DISTANCE_TYPE_UNSIGNED      })
+      igl::SIGNED_DISTANCE_TYPE_UNSIGNED      ,
+      igl::SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER 
+      })
   {
     Eigen::VectorXd S;
     Eigen::VectorXi I;
     Eigen::MatrixXd C,N;
     igl::signed_distance( P,V,F,type,S,I,C,N);
     Eigen::VectorXd Sexact (1,1);Sexact<<sqrt(1./12.);
-    test_common::assert_near(S,Sexact,1e-15);
+
+    if (type == igl::SIGNED_DISTANCE_TYPE_FAST_WINDING_NUMBER) {
+      // loosen tolerance on fwn. 
+      test_common::assert_near(S,Sexact,1e-7);
+    } else {
+      test_common::assert_near(S,Sexact,1e-15);
+    }
+    
   }
 }
 
@@ -52,7 +62,7 @@ TEST_CASE("signed_distance: single_triangle", "[igl]")
       igl::SIGNED_DISTANCE_TYPE_PSEUDONORMAL  ,
       igl::SIGNED_DISTANCE_TYPE_WINDING_NUMBER,
       igl::SIGNED_DISTANCE_TYPE_DEFAULT       ,
-      igl::SIGNED_DISTANCE_TYPE_UNSIGNED      
+      igl::SIGNED_DISTANCE_TYPE_UNSIGNED
       })
   {
     Eigen::VectorXd S;
@@ -63,3 +73,4 @@ TEST_CASE("signed_distance: single_triangle", "[igl]")
     test_common::assert_near(S,Sexact,1e-15);
   }
 }
+

tutorial/704_SignedDistance/main.cpp

@@ -19,7 +19,10 @@
 
 Eigen::MatrixXd V;
 Eigen::MatrixXi T,F;
+
 igl::AABB<Eigen::MatrixXd,3> tree;
+igl::FastWindingNumberBVH fwn_bvh;
+
 Eigen::MatrixXd FN,VN,EN;
 Eigen::MatrixXi E;
 Eigen::VectorXi EMAP;
@@ -28,6 +31,8 @@ double max_distance = 1;
 double slice_z = 0.5;
 bool overlay = false;
 
+bool useFastWindingNumber = false;
+
 void update_visualization(igl::opengl::glfw::Viewer & viewer)
 {
   using namespace Eigen;
@@ -72,12 +77,16 @@ void update_visualization(igl::opengl::glfw::Viewer & viewer)
 
   // Compute signed distance
   VectorXd S_vis;
+
+  if (!useFastWindingNumber)
   {
     VectorXi I;
     MatrixXd N,C;
     // Bunny is a watertight mesh so use pseudonormal for signing
     signed_distance_pseudonormal(V_vis,V,F,tree,FN,VN,EN,EMAP,S_vis,I,C,N);
-  }
+  } else {
+    signed_distance_fast_winding_number(V_vis, V, F, tree, fwn_bvh, S_vis);
+  }    
 
   const auto & append_mesh = [&F_vis,&V_vis](
     const Eigen::MatrixXd & V,
@@ -114,6 +123,12 @@ bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int mod)
     case ',':
       slice_z = std::max(slice_z-0.01,0.01);
       break;
+    case '1':
+      useFastWindingNumber = true;
+      break;
+    case '2':
+      useFastWindingNumber = false;
+      break;
   }
   update_visualization(viewer);
   return true;
@@ -127,6 +142,7 @@ int main(int argc, char *argv[])
   cout<<"Usage:"<<endl;
   cout<<"[space]  toggle showing surface."<<endl;
   cout<<"'.'/','  push back/pull forward slicing plane."<<endl;
+  cout<< "1/2 toggle between fast winding number (1) and pseudonormal (2) signing. \n";
   cout<<endl;
 
   // Load mesh: (V,T) tet-mesh of convex hull, F contains original surface
@@ -143,6 +159,9 @@ int main(int argc, char *argv[])
     max_distance = sqrt(sqrD.maxCoeff());
   }
 
+  // Fast winding and Pseudo normal depend on differnt AABB trees... We initialize both here.
+
+  // Pseudonormal setup...
   // Precompute signed distance AABB tree
   tree.init(V,F);
   // Precompute vertex,edge and face normals
@@ -152,6 +171,9 @@ int main(int argc, char *argv[])
   igl::per_edge_normals(
     V,F,igl::PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,FN,EN,E,EMAP);
 
+  // fast winding number setup (just init fwn bvh)
+  igl::fast_winding_number(V, F, 2, fwn_bvh);
+
   // Plot the generated mesh
   igl::opengl::glfw::Viewer viewer;
   update_visualization(viewer);