bad types in lex (#2544)

* swept volume overload

* swept volume: take transform list, templatize, expose SignedDistanceType

Replace the transform(t)+steps interface of swept_volume,
swept_volume_signed_distance and swept_volume_bounding_box with a list of
rigid transformations passed directly. The transform(t)+steps overloads are
removed rather than kept.

Templatize on Eigen::MatrixBase inputs / Eigen::PlainObjectBase outputs per
libigl style. The transform list is templated on both scalar and allocator so
std::vector<Eigen::Affine3d> and the aligned_allocator spelling both bind.

Expose SignedDistanceType on swept_volume and swept_volume_signed_distance,
dispatching like signed_distance_3 (pseudonormal / winding number / fast
winding number / unsigned) with the precomputation hoisted out of the
per-time-step loop.

swept_volume's isolevel is now a distance (typename DerivedV::Scalar) rather
than a count of grid cells. The padding cancelled out of the grid spacing
already (h == diag/(grid_res-1) regardless of pad), so pad is now derived from
the requested distance instead. Passing isolevel = k*h for integer k
reproduces the old isolevel_grid = k grid exactly. The doc had claimed
"distance level to be contoured" all along while the size_t type and
isolevel_grid implementation meant cells.

707_SweptVolume preloads the motion as a list of transforms sampled uniformly
over t in [0,1] and passes that to swept_volume. Its isolevel of 0.1 had been
truncating to 0 by the size_t conversion, so the tutorial was dilating by
nothing; it is now 10% of the bunny's largest side.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 5 <noreply@anthropic.com>

.gitattributes

@@ -0,0 +1,20 @@
+# Set the default behavior, in case people don't have core.autocrlf set.
+* text=auto
+
+# Explicitly declare text files you want to always be normalized and converted
+# to native line endings on checkout.
+*.tex text
+*.bib text
+*.svg text
+*.py text
+*.vbs text
+*.cpp text
+*.hpp text
+Makefile text
+
+# Declare files that will always have CRLF line endings on checkout.
+*.sln text eol=crlf
+
+# Denote all files that are truly binary and should not be modified.
+*.png binary
+*.jpg binary

.github/ISSUE_TEMPLATE/compilation.md

@@ -21,5 +21,6 @@ assignees: ''
 <!-- Check all that apply (change to `[x]`) -->
 
 - [ ] Windows
-- [ ] macOS
+- [ ] macOS Intel
+- [ ] macOS Arm (e.g., M1, M2)
 - [ ] Linux

.github/workflows/continuous.yml

@@ -20,24 +20,19 @@ jobs:
   ####################
 
   Unix:
-    name: ${{ matrix.name }} (${{ matrix.config }}, ${{ fromJSON('["HeaderOnly", "Static"]')[matrix.static == 'ON'] }})
+    name: ${{ matrix.os }} ${{ fromJSON('["Header-Only", "Static"]')[matrix.build-params.static == 'ON'] }} ${{ matrix.build-params.tutorials == 'ON' && 'tutorial' || ''}} ${{ matrix.build-params.tests == 'ON' && 'tests' || ''}} ${{ matrix.config }}
     runs-on: ${{ matrix.os }}
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-20.04, macos-latest]
+        os: [ubuntu-24.04, macos-15]
         config: [Release]
-        static: [ON, OFF]
-        include:
-          - os: macos-latest
-            name: macOS
-          - os: ubuntu-20.04
-            name: Linux
+        build-params: [ {static:  ON, tutorials:  ON, tests:  ON }, {static: OFF, tutorials: OFF, tests:  ON }, {static: OFF, tutorials:  ON, tests: OFF }]
     env:
       IGL_NUM_THREADS: 1  # See https://github.com/libigl/libigl/pull/996
     steps:
       - name: Checkout repository
-        uses: actions/checkout@v1
+        uses: actions/checkout@v4
         with:
           fetch-depth: 10
 
@@ -55,15 +50,20 @@ jobs:
       - name: Dependencies (macOS)
         if: runner.os == 'macOS'
         run: |
-          brew update
-          brew install ccache
+          HOMEBREW_NO_AUTO_UPDATE=1 brew install ccache
+
+      - name: Setup Xcode version
+        if: runner.os == 'macOS'
+        uses: maxim-lobanov/setup-xcode@v1
+        with:
+          xcode-version: latest-stable
 
       - name: Cache Build
         id: cache-build
-        uses: actions/cache@v2
+        uses: actions/cache@v4
         with:
           path: ~/.ccache
-          key: ${{ runner.os }}-${{ matrix.config }}-${{ matrix.static }}-cache
+          key: ${{ runner.os }}-${{ matrix.config }}-${{ matrix.build-params.static }}-cache
 
       - name: Prepare ccache
         run: |
@@ -74,76 +74,87 @@ jobs:
         run: |
           mkdir -p build
           cd build
+          # https://github.com/eclipse-ecal/ecal/issues/2041
           cmake .. \
+            -DCMAKE_POLICY_VERSION_MINIMUM=3.5 \
             -DCMAKE_CXX_COMPILER_LAUNCHER=ccache \
             -DCMAKE_BUILD_TYPE=${{ matrix.config }} \
-            -DLIBIGL_USE_STATIC_LIBRARY=${{ matrix.static }} \
-            -DLIBIGL_COPYLEFT_CGAL=ON 
+            -DLIBIGL_USE_STATIC_LIBRARY=${{ matrix.build-params.static }} \
+            -DLIBIGL_BUILD_TUTORIALS=${{ matrix.build-params.tutorials }} \
+            -DLIBIGL_GLFW_TESTS=OFF \
+            -DLIBIGL_BUILD_TESTS=${{ matrix.build-params.tests }} \
+            -DLIBIGL_COPYLEFT_CGAL=ON
 
       - name: Build
         run: cd build; make -j2; ccache --show-stats
 
       - name: Tests
-        run: cd build; ctest --verbose
+        run: cd build; ctest --show-only; ctest --verbose
 
   ####################
   # Windows
   ####################
 
   Windows:
-    name: Windows (${{ matrix.config }}, ${{ fromJSON('["HeaderOnly", "Static"]')[matrix.static == 'ON'] }})
+    name: Windows ${{ fromJSON('["Header-Only", "Static"]')[matrix.build-params.static == 'ON'] }} ${{ matrix.build-params.tutorials == 'ON' && 'tutorial' || ''}} ${{ matrix.build-params.selected_tutorial != 'NONE' && matrix.build-params.selected_tutorial || '' }} ${{ matrix.build-params.tests == 'ON' && 'tests' || ''}} ${{ matrix.config }}
     runs-on: windows-2022
-    env:
-      CC: cl.exe
-      CXX: cl.exe
     strategy:
       fail-fast: false
       matrix:
         config: [Release]
-        static: [ON, OFF]
+        build-params: [ 
+        {static:  ON, tutorials:  ON, tests:  ON, selected_tutorial: NONE}, 
+        {static: OFF, tutorials: OFF, tests:  ON, selected_tutorial: NONE}, 
+        {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 1},
+        {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 2},
+        {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 3},
+        {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 4},
+        {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 5},
+        {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 6},
+        {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 7},
+        {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 8},
+        {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 9},
+        {static: OFF, tutorials:  ON, tests: OFF, selected_tutorial: 10},
+      ]
     steps:
       - name: Checkout repository
-        uses: actions/checkout@v1
+        uses: actions/checkout@v4
         with:
           fetch-depth: 10
-      - uses: seanmiddleditch/gha-setup-ninja@master
 
-      - name: Set env
-        run: |
-          echo "appdata=$env:LOCALAPPDATA" >> ${env:GITHUB_ENV}
+      - name: Install Ninja
+        uses: seanmiddleditch/gha-setup-ninja@master
 
-      - name: Cache build
-        id: cache-build
-        uses: actions/cache@v2
-        with:
-          path: ${{ env.appdata }}\Mozilla\sccache
-          key: ${{ runner.os }}-${{ matrix.config }}-${{ matrix.static }}-cache
+      - name: Install sccache
+        uses: mozilla-actions/sccache-action@v0.0.9
 
-
-      - name: Prepare sccache
-        run: |
-          Invoke-Expression (New-Object System.Net.WebClient).DownloadString('https://get.scoop.sh')
-          scoop install sccache --global
-          # Scoop modifies the PATH so we make it available for the next steps of the job
-          echo "${env:PATH}" >> ${env:GITHUB_PATH}
-
-        # We run configure + build in the same step, since they both need to call VsDevCmd
-        # Also, cmd uses ^ to break commands into multiple lines (in powershell this is `)
       - name: Configure and build
         shell: cmd
         run: |
           call "C:\Program Files\Microsoft Visual Studio\2022\Enterprise\Common7\Tools\VsDevCmd.bat" -arch=x64
+          # https://github.com/eclipse-ecal/ecal/issues/2041
           cmake -G Ninja ^
+            -DCMAKE_POLICY_VERSION_MINIMUM=3.5 ^
             -DCMAKE_CXX_COMPILER_LAUNCHER=sccache ^
             -DCMAKE_BUILD_TYPE=${{ matrix.config }} ^
-            -DLIBIGL_USE_STATIC_LIBRARY=${{ matrix.static }} ^
+            -DLIBIGL_USE_STATIC_LIBRARY=${{ matrix.build-params.static }} ^
             -DLIBIGL_COPYLEFT_CGAL=ON ^
-            -DCMAKE_JOB_POOLS=pool-linking=1;pool-compilation=1 ^
-            -DCMAKE_JOB_POOL_COMPILE:STRING=pool-compilation ^
-            -DCMAKE_JOB_POOL_LINK:STRING=pool-linking ^
+            -DLIBIGL_BUILD_TUTORIALS=${{ matrix.build-params.tutorials }} ^
+            -DLIBIGL_BUILD_TESTS=${{ matrix.build-params.tests }} ^
+            -DLIBIGL_GLFW_TESTS=OFF ^
+            -DLIBIGL_TUTORIALS_CHAPTER1=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '1') && 'ON' || 'OFF' }} ^
+            -DLIBIGL_TUTORIALS_CHAPTER2=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '2') && 'ON' || 'OFF' }} ^
+            -DLIBIGL_TUTORIALS_CHAPTER3=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '3') && 'ON' || 'OFF' }} ^
+            -DLIBIGL_TUTORIALS_CHAPTER4=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '4') && 'ON' || 'OFF' }} ^
+            -DLIBIGL_TUTORIALS_CHAPTER5=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '5') && 'ON' || 'OFF' }} ^
+            -DLIBIGL_TUTORIALS_CHAPTER6=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '6') && 'ON' || 'OFF' }} ^
+            -DLIBIGL_TUTORIALS_CHAPTER7=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '7') && 'ON' || 'OFF' }} ^
+            -DLIBIGL_TUTORIALS_CHAPTER8=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '8') && 'ON' || 'OFF' }} ^
+            -DLIBIGL_TUTORIALS_CHAPTER9=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '9') && 'ON' || 'OFF' }} ^
+            -DLIBIGL_TUTORIALS_CHAPTER10=${{ (matrix.build-params.selected_tutorial == 'NONE' || matrix.build-params.selected_tutorial == '10') && 'ON' || 'OFF' }} ^
             -B build ^
             -S .
-          cmake --build build
+          cmake --build build -j2
 
       - name: Tests
-        run: cd build; ctest --verbose
+        run: cd build; ctest --show-only; ctest --verbose -j2

.gitignore

@@ -47,3 +47,8 @@ LibiglOptions.cmake
 
 # macos debris
 .DS_Store
+*~
+dox/
+latex/
+scripts/
+CLAUDE.md

CMakeLists.txt

@@ -6,6 +6,8 @@ else()
     set(LIBIGL_TOPLEVEL_PROJECT OFF)
 endif()
 
+
+
 # Check required CMake version
 set(REQUIRED_CMAKE_VERSION "3.16.0")
 if(LIBIGL_TOPLEVEL_PROJECT)
@@ -29,13 +31,13 @@ option(HUNTER_ENABLED "Enable Hunter package manager support" OFF)
 if(HUNTER_ENABLED)
     include("cmake/misc/HunterGate.cmake")
     HunterGate(
-        URL "https://github.com/cpp-pm/hunter/archive/v0.23.300.tar.gz"
-        SHA1 "1151d539465d9cdbc880ee30f794864aec11c448"
+        URL "https://github.com/cpp-pm/hunter/archive/v0.24.8.tar.gz"
+        SHA1 "ca7838dded9a1811b04ffd56175f629e0af82d3d"
     )
 endif()
 
 ################################################################################
-project(libigl VERSION 2.4.0)
+project(libigl VERSION 2.5.0)
 
 # CMake module path
 list(PREPEND CMAKE_MODULE_PATH
@@ -55,6 +57,8 @@ set_property(GLOBAL PROPERTY __igl_module_path ${CMAKE_MODULE_PATH})
 
 set(LIBIGL_DEFAULT_CGAL ${LIBIGL_TOPLEVEL_PROJECT})
 set(MATLAB_ADDITIONAL_VERSIONS
+    "R2023b=10.4"
+    "R2023a=10.4"
     "R2022b=10.3"
     "R2022a=10.2"
     "R2021b=10.1"
@@ -74,6 +78,14 @@ if(LIBIGL_TOPLEVEL_PROJECT)
         message(WARNING "Mosek not found, disabling igl_restricted::mosek module.")
     endif()
 endif()
+set(LIBIGL_DEFAULT_COMISO ${LIBIGL_TOPLEVEL_PROJECT})
+if(LIBIGL_TOPLEVEL_PROJECT AND (NOT APPLE) AND UNIX)
+  find_package(BLAS QUIET)
+    if(NOT BLAS_FOUND)
+        set(LIBIGL_DEFAULT_COMISO OFF)
+        message(WARNING "BLAS not found, disabling igl_copyleft::comiso module.")
+    endif()
+endif()
 
 # Build tests and tutorials
 option(LIBIGL_BUILD_TESTS      "Build libigl unit test"                ${LIBIGL_TOPLEVEL_PROJECT})
@@ -87,12 +99,14 @@ option(LIBIGL_USE_STATIC_LIBRARY "Use libigl as static library" ${LIBIGL_TOPLEVE
 
 # Permissive modules. These modules are available under MPL2 license, and their dependencies are available
 # under a permissive or public domain license.
+option(LIBIGL_CYCODEBASE       "Build target igl::cycodebase"       ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_EMBREE           "Build target igl::embree"           ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_GLFW             "Build target igl::glfw"             ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_IMGUI            "Build target igl::imgui"            ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_OPENGL           "Build target igl::opengl"           ${LIBIGL_TOPLEVEL_PROJECT})
-option(LIBIGL_PNG              "Build target igl::png"              ${LIBIGL_TOPLEVEL_PROJECT})
+option(LIBIGL_STB              "Build target igl::stb"              ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_PREDICATES       "Build target igl::predicates"       ${LIBIGL_TOPLEVEL_PROJECT})
+option(LIBIGL_SPECTRA          "Build target igl::spectra"          ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_XML              "Build target igl::xml"              ${LIBIGL_TOPLEVEL_PROJECT})
 
 # Copyleft modules. These modules are available under GPL license, and their dependencies are
@@ -108,6 +122,12 @@ option(LIBIGL_RESTRICTED_MATLAB   "Build target igl_restricted::matlab"   ${LIBI
 option(LIBIGL_RESTRICTED_MOSEK    "Build target igl_restricted::mosek"    ${LIBIGL_DEFAULT_MOSEK})
 option(LIBIGL_RESTRICTED_TRIANGLE "Build target igl_restricted::triangle" ${LIBIGL_TOPLEVEL_PROJECT})
 
+# GLFW doesn't run on headless CI machines so don't run (or build them).
+# Unfortunately on headless mac machines glfw seems to hang rather than crash
+# making it hard to catch at runtime.
+option(LIBIGL_GLFW_TESTS       "Build igl::glfw tests"             ${LIBIGL_TOPLEVEL_PROJECT})
+
+option(LIBIGL_WARNINGS_AS_ERRORS "Turn on many warnings and treat as errors" OFF)
 
 set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin")
 set(CMAKE_LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib")
@@ -167,3 +187,4 @@ if(LIBIGL_INSTALL)
     write_basic_package_version_file("${version_config_file}" COMPATIBILITY SameMajorVersion)
     install(FILES "${project_config_out}" "${version_config_file}" DESTINATION "${export_dest_dir}")
 endif()
+

LibiglOptions.cmake.sample

@@ -56,7 +56,7 @@
 # option(LIBIGL_GLFW                "Build target igl::glfw"                ON)
 # option(LIBIGL_IMGUI               "Build target igl::imgui"               ON)
 # option(LIBIGL_OPENGL              "Build target igl::opengl"              ON)
-# option(LIBIGL_PNG                 "Build target igl::png"                 ON)
+# option(LIBIGL_STB                 "Build target igl::stb"                 ON)
 # option(LIBIGL_PREDICATES          "Build target igl::predicates"          ON)
 # option(LIBIGL_XML                 "Build target igl::xml"                 ON)
 # option(LIBIGL_COPYLEFT_CGAL       "Build target igl_copyleft::cgal"       ON)

README.md

@@ -8,6 +8,6 @@
 Documentation, tutorial, and instructions at <https://libigl.github.io>.
 
 
-| 🚨 Important |
+| 🆕 Doxygen Documentation |
 |:---|
-| The latest version of libigl (v2.4.0) introduces some **breaking changes** to its CMake build system. Please read our [changelog](https://libigl.github.io/changelog/) page for instructions on how to update your project accordingly. |
+| The latest version of libigl (v2.5.0) introduces [doxygen generated detailed documentation](https://libigl.github.io/dox/index.html) |

cmake/find/FindGMP.cmake

@@ -1,95 +0,0 @@
-# Try to find the GNU Multiple Precision Arithmetic Library (GMP)
-# See http://gmplib.org/
-
-if(${CMAKE_VERSION} VERSION_LESS "3.18.0")
-    set(REQUIRED_FLAG "")
-else()
-    set(REQUIRED_FLAG REQUIRED)
-endif()
-
-# On Windows, we must use the pre-compiled versions downloaded with libigl
-if(WIN32)
-    set(NO_DEFAULT_FLAG NO_DEFAULT_PATH)
-else()
-    set(NO_DEFAULT_FLAG "")
-endif()
-
-find_path(GMP_INCLUDES
-    NAMES
-        gmp.h
-    PATHS
-        ENV GMP_DIR
-        ${INCLUDE_INSTALL_DIR}
-    PATH_SUFFIXES
-        include
-    ${REQUIRED_FLAG}
-    ${NO_DEFAULT_FLAG}
-)
-
-find_library(GMP_LIBRARIES
-    NAMES
-        gmp
-        libgmp-10
-    PATHS
-        ENV GMP_DIR
-        ${LIB_INSTALL_DIR}
-    PATH_SUFFIXES
-        lib
-    ${REQUIRED_FLAG}
-    ${NO_DEFAULT_FLAG}
-)
-
-set(GMP_EXTRA_VARS "")
-if(WIN32)
-    # Find dll file and set IMPORTED_LOCATION to the .dll file
-    find_file(GMP_RUNTIME_LIB
-        NAMES
-            gmp.dll
-            libgmp-10.dll
-        PATHS
-            ENV GMP_DIR
-            ${LIB_INSTALL_DIR}
-        PATH_SUFFIXES
-            lib
-        ${REQUIRED_FLAG}
-        ${NO_DEFAULT_FLAG}
-    )
-    list(APPEND GMP_EXTRA_VARS GMP_RUNTIME_LIB)
-endif()
-
-include(FindPackageHandleStandardArgs)
-find_package_handle_standard_args(GMP
-    REQUIRED_VARS
-        GMP_INCLUDES
-        GMP_LIBRARIES
-        ${GMP_EXTRA_VARS}
-    REASON_FAILURE_MESSAGE
-        "GMP is not installed on your system. Either install GMP using your preferred package manager, or disable libigl modules that depend on GMP, such as CGAL. See LibiglOptions.cmake.sample for configuration options. Do not forget to delete your <build>/CMakeCache.txt for the changes to take effect."
-)
-mark_as_advanced(GMP_INCLUDES GMP_LIBRARIES)
-
-if(GMP_INCLUDES AND GMP_LIBRARIES AND NOT TARGET gmp::gmp)
-    if(GMP_RUNTIME_LIB)
-        add_library(gmp::gmp SHARED IMPORTED)
-    else()
-        add_library(gmp::gmp UNKNOWN IMPORTED)
-    endif()
-
-    # Set public header location and link language
-    set_target_properties(gmp::gmp PROPERTIES
-        IMPORTED_LINK_INTERFACE_LANGUAGES "C"
-        INTERFACE_INCLUDE_DIRECTORIES "${GMP_INCLUDES}"
-    )
-
-    # Set lib location. On Windows we specify both the .lib and the .dll paths
-    if(GMP_RUNTIME_LIB)
-        set_target_properties(gmp::gmp PROPERTIES
-            IMPORTED_IMPLIB "${GMP_LIBRARIES}"
-            IMPORTED_LOCATION "${GMP_RUNTIME_LIB}"
-        )
-    else()
-        set_target_properties(gmp::gmp PROPERTIES
-            IMPORTED_LOCATION "${GMP_LIBRARIES}"
-        )
-    endif()
-endif()

cmake/find/FindMPFR.cmake

@@ -1,95 +0,0 @@
-# Try to find the MPFR library
-# See http://www.mpfr.org/
-
-if(${CMAKE_VERSION} VERSION_LESS "3.18.0")
-    set(REQUIRED_FLAG "")
-else()
-    set(REQUIRED_FLAG REQUIRED)
-endif()
-
-# On Windows, we must use the pre-compiled versions downloaded with libigl
-if(WIN32)
-    set(NO_DEFAULT_FLAG NO_DEFAULT_PATH)
-else()
-    set(NO_DEFAULT_FLAG "")
-endif()
-
-find_path(MPFR_INCLUDES
-    NAMES
-        mpfr.h
-    PATHS
-        ENV MPFR_DIR
-        ${INCLUDE_INSTALL_DIR}
-    PATH_SUFFIXES
-        include
-    ${REQUIRED_FLAG}
-    ${NO_DEFAULT_FLAG}
-)
-
-find_library(MPFR_LIBRARIES
-    NAMES
-        mpfr
-        libmpfr-4
-    PATHS
-        ENV MPFR_DIR
-        ${LIB_INSTALL_DIR}
-    PATH_SUFFIXES
-        lib
-    ${REQUIRED_FLAG}
-    ${NO_DEFAULT_FLAG}
-)
-
-set(MPFR_EXTRA_VARS "")
-if(WIN32)
-    # Find dll file and set IMPORTED_LOCATION to the .dll file
-    find_file(MPFR_RUNTIME_LIB
-        NAMES
-            mpfr.dll
-            libmpfr-4.dll
-        PATHS
-            ENV MPFR_DIR
-            ${LIB_INSTALL_DIR}
-        PATH_SUFFIXES
-            lib
-        ${REQUIRED_FLAG}
-        ${NO_DEFAULT_FLAG}
-    )
-    list(APPEND MPFR_EXTRA_VARS MPFR_RUNTIME_LIB)
-endif()
-
-include(FindPackageHandleStandardArgs)
-find_package_handle_standard_args(MPFR
-    REQUIRED_VARS
-        MPFR_INCLUDES
-        MPFR_LIBRARIES
-        ${MPFR_EXTRA_VARS}
-    REASON_FAILURE_MESSAGE
-        "MPFR is not installed on your system. Either install MPFR using your preferred package manager, or disable libigl modules that depend on MPFR, such as CGAL. See LibiglOptions.cmake.sample for configuration options. Do not forget to delete your <build>/CMakeCache.txt for the changes to take effect."
-)
-mark_as_advanced(MPFR_INCLUDES MPFR_LIBRARIES)
-
-if(MPFR_INCLUDES AND MPFR_LIBRARIES AND NOT TARGET mpfr::mpfr)
-    if(MPFR_RUNTIME_LIB)
-        add_library(mpfr::mpfr SHARED IMPORTED)
-    else()
-        add_library(mpfr::mpfr UNKNOWN IMPORTED)
-    endif()
-
-    # Set public header location and link language
-    set_target_properties(mpfr::mpfr PROPERTIES
-        IMPORTED_LINK_INTERFACE_LANGUAGES "C"
-        INTERFACE_INCLUDE_DIRECTORIES "${MPFR_INCLUDES}"
-    )
-
-    # Set lib location. On Windows we specify both the .lib and the .dll paths
-    if(MPFR_RUNTIME_LIB)
-        set_target_properties(mpfr::mpfr PROPERTIES
-            IMPORTED_IMPLIB "${MPFR_LIBRARIES}"
-            IMPORTED_LOCATION "${MPFR_RUNTIME_LIB}"
-        )
-    else()
-        set_target_properties(mpfr::mpfr PROPERTIES
-            IMPORTED_LOCATION "${MPFR_LIBRARIES}"
-        )
-    endif()
-endif()

cmake/igl/igl_add_library.cmake

@@ -33,8 +33,12 @@ function(igl_add_library module_name)
         target_compile_definitions(${module_name} ${IGL_SCOPE} -DIGL_STATIC_LIBRARY)
     endif()
 
-    # C++11 features
-    target_compile_features(${module_name} ${IGL_SCOPE} cxx_std_11)
+    # C++17 features
+    target_compile_features(${module_name} ${IGL_SCOPE} cxx_std_17)
+
+    if(LIBIGL_WARNINGS_AS_ERRORS)
+      target_compile_options(${module_name} PRIVATE -Wall -Wextra -Wpedantic -Wno-sign-compare -Werror -Wno-gnu -Wno-unknown-pragmas)
+    endif()
 
     # Other compilation flags
     if(MSVC)
@@ -42,6 +46,7 @@ function(igl_add_library module_name)
         target_compile_options(${module_name} ${IGL_SCOPE} $<$<COMPILE_LANGUAGE:CXX>:/MP> $<$<COMPILE_LANGUAGE:CXX>:/bigobj>)
         target_compile_definitions(${module_name} ${IGL_SCOPE} -DNOMINMAX)
 
+
         # Silencing some compilation warnings
         if(LIBIGL_USE_STATIC_LIBRARY)
             target_compile_options(${module_name} PRIVATE

cmake/igl/igl_add_tutorial.cmake

@@ -7,7 +7,9 @@ function(igl_add_tutorial name)
     endforeach()
 
     message(STATUS "Creating libigl tutorial: ${name}")
-    add_executable(${name} ${CMAKE_CURRENT_SOURCE_DIR}/${name}/main.cpp)
+    # get all cpp files in ${CMAKE_CURRENT_SOURCE_DIR}/${name}/
+    file(GLOB SRCFILES ${CMAKE_CURRENT_SOURCE_DIR}/${name}/*.cpp)
+    add_executable(${name} ${SRCFILES})
     target_link_libraries(${name} PRIVATE
         igl::core
         igl::tutorial_data
@@ -15,4 +17,13 @@ function(igl_add_tutorial name)
     )
 
     set_target_properties(${name} PROPERTIES FOLDER Libigl_Tutorials)
+
+    # Do this codesign only on macOS
+    # add_custom_command(TARGET your_target POST_BUILD COMMAND codesign -s - $<TARGET_FILE:your_target>
+    if(APPLE)
+      add_custom_command(TARGET ${name} POST_BUILD
+        COMMAND codesign -f -s - $<TARGET_FILE:${name}>
+        COMMENT "Codesigning ${name}"
+      )
+    endif()
 endfunction()

cmake/igl/igl_copy_dll.cmake

@@ -50,6 +50,10 @@ function(igl_copy_dll target)
     if(NOT WIN32)
         return()
     endif()
+    if(NOT TARGET ${target})
+      message(STATUS "igl_copy_dll() was called with a non-target: ${target}")
+      return()
+    endif()
 
     # Sanity checks
     get_target_property(TYPE ${target} TYPE)

cmake/igl/libigl-config.cmake.in

@@ -1,3 +1,7 @@
 @PACKAGE_INIT@
 
+include(CMakeFindDependencyMacro)
+find_dependency(Eigen3 REQUIRED)
+find_dependency(Threads REQUIRED)
+include("${CMAKE_CURRENT_LIST_DIR}/LibiglConfigTargets.cmake")
 check_required_components(Libigl)

cmake/igl/modules/copyleft/comiso.cmake

@@ -1,27 +0,0 @@
-# 1. Define module
-igl_add_library(igl_copyleft_comiso)
-
-# 2. Include headers
-include(GNUInstallDirs)
-target_include_directories(igl_copyleft_comiso ${IGL_SCOPE}
-    $<BUILD_INTERFACE:${libigl_SOURCE_DIR}/include>
-    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
-)
-
-# 3. Target sources
-file(GLOB INC_FILES "${libigl_SOURCE_DIR}/include/igl/copyleft/comiso/*.h")
-file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/include/igl/copyleft/comiso/*.cpp")
-igl_target_sources(igl_copyleft_comiso ${INC_FILES} ${SRC_FILES})
-
-# 4. Dependencies
-include(comiso)
-igl_include(copyleft core)
-target_link_libraries(igl_copyleft_comiso ${IGL_SCOPE}
-    igl::core
-    igl_copyleft::core
-    CoMISo::CoMISo
-)
-
-# 5. Unit tests
-file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/tests/include/igl/copyleft/comiso/*.cpp")
-igl_add_test(igl_copyleft_comiso ${SRC_FILES})

cmake/igl/modules/cycodebase.cmake

@@ -0,0 +1,26 @@
+# 1. Define module
+igl_add_library(igl_cycodebase)
+
+# 2. Include headers
+include(GNUInstallDirs)
+target_include_directories(igl_cycodebase ${IGL_SCOPE}
+    $<BUILD_INTERFACE:${libigl_SOURCE_DIR}/include>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+
+# 3. Target sources
+file(GLOB INC_FILES "${libigl_SOURCE_DIR}/include/igl/cycodebase/*.h")
+file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/include/igl/cycodebase/*.cpp")
+igl_target_sources(igl_cycodebase ${INC_FILES} ${SRC_FILES})
+
+# 4. Dependencies
+include(cycodebase)
+target_link_libraries(igl_cycodebase ${IGL_SCOPE}
+    igl::core
+    cyCodeBase::cyCodeBase
+)
+
+# 5. Unit tests
+file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/tests/include/igl/cycodebase/*.cpp")
+igl_add_test(igl_cycodebase ${SRC_FILES})
+

cmake/igl/modules/glfw.cmake

@@ -21,3 +21,9 @@ target_link_libraries(igl_glfw ${IGL_SCOPE}
     igl::opengl
     glfw::glfw
 )
+
+# 5. Unit tests
+if(LIBIGL_GLFW_TESTS)
+  file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/tests/include/igl/opengl/glfw/*.cpp")
+  igl_add_test(igl_glfw ${SRC_FILES})
+endif()

cmake/igl/modules/png.cmake

@@ -1,23 +0,0 @@
-# 1. Define module
-igl_add_library(igl_png)
-
-# 2. Include headers
-include(GNUInstallDirs)
-target_include_directories(igl_png ${IGL_SCOPE}
-    $<BUILD_INTERFACE:${libigl_SOURCE_DIR}/include>
-    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
-)
-
-# 3. Target sources
-file(GLOB INC_FILES "${libigl_SOURCE_DIR}/include/igl/png/*.h")
-file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/include/igl/png/*.cpp")
-igl_target_sources(igl_png ${INC_FILES} ${SRC_FILES})
-
-# 4. Dependencies
-include(stb)
-igl_include(opengl)
-target_link_libraries(igl_png ${IGL_SCOPE}
-    igl::core
-    igl::opengl
-    stb::stb
-)

cmake/igl/modules/spectra.cmake

@@ -0,0 +1,25 @@
+# 1. Define module
+igl_add_library(igl_spectra)
+
+# 2. Include headers
+include(GNUInstallDirs)
+target_include_directories(igl_spectra ${IGL_SCOPE}
+    $<BUILD_INTERFACE:${libigl_SOURCE_DIR}/include>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+
+# 3. Target sources
+file(GLOB INC_FILES "${libigl_SOURCE_DIR}/include/igl/spectra/*.h")
+file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/include/igl/spectra/*.cpp")
+igl_target_sources(igl_spectra ${INC_FILES} ${SRC_FILES})
+
+# 4. Dependencies
+include(spectra)
+target_link_libraries(igl_spectra ${IGL_SCOPE}
+    igl::core
+    spectra::spectra
+)
+
+# 5. Unit tests
+file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/tests/include/igl/spectra/*.cpp")
+igl_add_test(igl_spectra ${SRC_FILES})

cmake/igl/modules/stb.cmake

@@ -0,0 +1,35 @@
+# 1. Define module
+igl_add_library(igl_stb)
+
+# 2. Include headers
+include(GNUInstallDirs)
+target_include_directories(igl_stb ${IGL_SCOPE}
+    $<BUILD_INTERFACE:${libigl_SOURCE_DIR}/include>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+
+# 3. Target sources
+file(GLOB INC_FILES "${libigl_SOURCE_DIR}/include/igl/stb/*.h")
+file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/include/igl/stb/*.cpp")
+if(LIBIGL_OPENGL)
+  message(STATUS "Including igl/opengl/stb support")
+  file(GLOB OPENGL_INC_FILES "${libigl_SOURCE_DIR}/include/igl/opengl/stb/*.h")
+  file(GLOB OPENGL_SRC_FILES "${libigl_SOURCE_DIR}/include/igl/opengl/stb/*.cpp")
+  list(APPEND INC_FILES ${OPENGL_INC_FILES})
+  list(APPEND SRC_FILES ${OPENGL_SRC_FILES})
+endif()
+igl_target_sources(igl_stb ${INC_FILES} ${SRC_FILES})
+
+# 4. Dependencies
+include(stb)
+target_link_libraries(igl_stb ${IGL_SCOPE}
+    igl::core
+    stb::stb
+)
+
+if(LIBIGL_OPENGL)
+  igl_include(opengl)
+  target_link_libraries(igl_stb ${IGL_SCOPE}
+      igl::opengl
+  )
+endif()

cmake/libigl.cmake

@@ -13,18 +13,19 @@ include(igl_windows)
 
 # Libigl permissive modules
 igl_include(core)
+igl_include_optional(cycodebase)
 igl_include_optional(embree)
 igl_include_optional(opengl)
 igl_include_optional(glfw)
 igl_include_optional(imgui)
 igl_include_optional(predicates)
-igl_include_optional(png)
+igl_include_optional(stb)
+igl_include_optional(spectra)
 igl_include_optional(xml)
 
 # Libigl copyleft modules
 igl_include_optional(copyleft core)
 igl_include_optional(copyleft cgal)
-igl_include_optional(copyleft comiso)
 igl_include_optional(copyleft tetgen)
 
 # Libigl restricted modules

cmake/recipes/external/boost.cmake

@@ -4,32 +4,27 @@ endif()
 
 message(STATUS "Third-party: creating targets 'Boost::boost'...")
 
+cmake_minimum_required(VERSION 3.24) # Ensure modern FetchContent features
+project(BoostFetchExample)
+
 include(FetchContent)
+
+# Define the Boost library to fetch
 FetchContent_Declare(
-    boost-cmake
-    GIT_REPOSITORY https://github.com/Orphis/boost-cmake.git
-    GIT_TAG 7f97a08b64bd5d2e53e932ddf80c40544cf45edf
+    Boost
+    URL https://archives.boost.io/release/1.86.0/source/boost_1_86_0.tar.gz
+    URL_HASH MD5=ac857d73bb754b718a039830b07b9624
 )
+# Fetch Boost
+FetchContent_MakeAvailable(Boost)
 
-set(PREVIOUS_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS})
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fPIC")
-set(OLD_CMAKE_POSITION_INDEPENDENT_CODE ${CMAKE_POSITION_INDEPENDENT_CODE})
-set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 
-# This guy will download boost using FetchContent
-FetchContent_GetProperties(boost-cmake)
-if(NOT boost-cmake_POPULATED)
-    FetchContent_Populate(boost-cmake)
-    # File lcid.cpp from Boost_locale.cpp doesn't compile on MSVC, so we exclude them from the default
-    # targets being built by the project (only targets explicitly used by other targets will be built).
-    add_subdirectory(${boost-cmake_SOURCE_DIR} ${boost-cmake_BINARY_DIR} EXCLUDE_FROM_ALL)
-endif()
+# Ensure Boost paths are set before CGAL
+set(Boost_INCLUDE_DIR ${boost_SOURCE_DIR})
+set(Boost_LIBRARY_DIR ${boost_BINARY_DIR})
 
-set(CMAKE_POSITION_INDEPENDENT_CODE ${OLD_CMAKE_POSITION_INDEPENDENT_CODE})
-set(CMAKE_CXX_FLAGS "${PREVIOUS_CMAKE_CXX_FLAGS}")
-
-# Set VS target folders
-set(boost_modules
+# Add Boost libraries needed for your project
+set(BOOST_LIBRARIES
     container
     regex
     atomic
@@ -48,6 +43,7 @@ set(boost_modules
     log_setup
     unit_test_framework
     math
+    multiprecision
     program_options
     timer
     random
@@ -55,10 +51,10 @@ set(boost_modules
     system
     thread
     type_erasure
-)
-foreach(module IN ITEMS ${boost_modules})
-    if(TARGET Boost_${module})
-        set_target_properties(Boost_${module} PROPERTIES FOLDER ThirdParty/Boost)
-    endif()
-endforeach()
+  )
 
+foreach(lib IN LISTS BOOST_LIBRARIES)
+    add_library(boost_${lib} INTERFACE)
+    target_include_directories(boost_${lib} INTERFACE ${Boost_SOURCE_DIR})
+    target_link_libraries(boost_${lib} INTERFACE Boost::${lib})
+endforeach()

cmake/recipes/external/cgal.cmake

@@ -7,12 +7,8 @@ message(STATUS "Third-party: creating target 'CGAL::CGAL'")
 include(FetchContent)
 FetchContent_Declare(
     cgal
-    #GIT_REPOSITORY https://github.com/CGAL/cgal.git
-    #GIT_TAG        f7c3c8212b56c0d6dae63787efc99093f4383415
-    URL https://github.com/CGAL/cgal/releases/download/releases%2FCGAL-4.12.2/CGAL-4.12.2.tar.xz
-    URL_MD5 c94a0081c3836fd01ccb4d1e8bdd5d4f
-    # URL https://github.com/CGAL/cgal/releases/download/v5.2.1/CGAL-5.2.1-library.tar.xz
-    # URL_MD5 c1c3a9abe9106b5f3ff8dccaf2ddc0b7
+    URL https://github.com/CGAL/cgal/releases/download/v6.0.1/CGAL-6.0.1-library.tar.xz
+    URL_MD5 ea827f6778063e00554ae41f4c845492
 )
 FetchContent_GetProperties(cgal)
 if(cgal_POPULATED)
@@ -33,25 +29,21 @@ function(cgal_import_target)
         set(${NAME}_ROOT ${CMAKE_CURRENT_BINARY_DIR}/${NAME} CACHE PATH "")
     endmacro()
 
-    include(gmp)
-    include(mpfr)
     include(boost)
 
-    ignore_package(GMP 5.0.1)
-    set(GMP_INCLUDE_DIR "")
-    set(GMP_LIBRARIES gmp::gmp)
-
-    ignore_package(MPFR 3.0.0)
-    set(MPFR_INCLUDE_DIR "")
-    set(MPFR_LIBRARIES mpfr::mpfr)
 
     ignore_package(Boost 1.71.0)
     set(Boost_INCLUDE_DIRS "")
-    set(Boost_LIBRARIES Boost::thread Boost::system)
+    set(Boost_LIBRARIES Boost::thread Boost::system Boost::multiprecision)
 
     # Prefer Config mode before Module mode to prevent CGAL from loading its own FindXXX.cmake
     set(CMAKE_FIND_PACKAGE_PREFER_CONFIG TRUE)
 
+    # https://stackoverflow.com/a/71714947/148668
+    set(CGAL_DATA_DIR "unspecified")
+
+    set(CGAL_CMAKE_EXACT_NT_BACKEND "BOOST_BACKEND" CACHE STRING "CGAL exact NT backend")
+    set(CGAL_DISABLE_GMP ON CACHE BOOL "Disable GMP")
     find_package(CGAL CONFIG COMPONENTS Core PATHS ${cgal_SOURCE_DIR} NO_DEFAULT_PATH)
 endfunction()
 

cmake/recipes/external/comiso.cmake

@@ -1,31 +0,0 @@
-if(TARGET CoMISo::CoMISo)
-    return()
-endif()
-
-message(STATUS "Third-party: creating target 'CoMISo::CoMISo'")
-
-include(FetchContent)
-FetchContent_Declare(
-    comiso
-    GIT_REPOSITORY https://github.com/libigl/CoMISo.git
-    GIT_TAG 536440e714f412e7ef6c0b96b90ba37b1531bb39
-)
-
-include(eigen)
-
-FetchContent_MakeAvailable(comiso)
-
-add_library(CoMISo::CoMISo ALIAS CoMISo)
-
-# Copy .hh headers into a subfolder `CoMISo/`
-file(GLOB_RECURSE INC_FILES "${comiso_SOURCE_DIR}/*.hh" "${comiso_SOURCE_DIR}/*.cc")
-set(output_folder "${CMAKE_CURRENT_BINARY_DIR}/CoMISo/include/CoMISo")
-message(VERBOSE "Copying CoMISo headers to '${output_folder}'")
-foreach(filepath IN ITEMS ${INC_FILES})
-    file(RELATIVE_PATH filename "${comiso_SOURCE_DIR}" ${filepath})
-    configure_file(${filepath} "${output_folder}/${filename}" COPYONLY)
-endforeach()
-
-target_include_directories(CoMISo PUBLIC ${CMAKE_CURRENT_BINARY_DIR}/CoMISo/include)
-
-set_target_properties(CoMISo PROPERTIES FOLDER ThirdParty)

cmake/recipes/external/cycodebase.cmake

@@ -0,0 +1,18 @@
+if(TARGET cycodebase::cycodebase)
+    return()
+endif()
+
+FetchContent_Declare(
+  cyCodeBase
+  GIT_REPOSITORY https://github.com/cemyuksel/cyCodeBase/
+  GIT_TAG e36f3cffca65eb12a8a071f0443128b7de6ed75d
+)
+FetchContent_Populate(cyCodeBase)
+add_library(cyCodeBase_interface INTERFACE)
+target_include_directories(cyCodeBase_interface INTERFACE ${cycodebase_SOURCE_DIR})
+
+if(NOT (CMAKE_SYSTEM_PROCESSOR MATCHES "x86_64|AMD64|i[3-6]86"))
+  target_compile_definitions(cyCodeBase_interface INTERFACE CY_NO_INTRIN_H)
+endif()
+
+add_library(cyCodeBase::cyCodeBase ALIAS cyCodeBase_interface)

cmake/recipes/external/eigen.cmake

@@ -8,7 +8,7 @@ include(FetchContent)
 FetchContent_Declare(
     eigen
     GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
-    GIT_TAG tags/3.3.7
+    GIT_TAG tags/5.0.1
     GIT_SHALLOW TRUE
 )
 FetchContent_GetProperties(eigen)

cmake/recipes/external/embree.cmake

@@ -8,7 +8,7 @@ include(FetchContent)
 FetchContent_Declare(
     embree
     GIT_REPOSITORY https://github.com/embree/embree.git
-    GIT_TAG        v3.13.3
+    GIT_TAG        v4.4.0
     GIT_SHALLOW    TRUE
 )
 

cmake/recipes/external/glad.cmake

@@ -8,7 +8,7 @@ include(FetchContent)
 FetchContent_Declare(
     glad
     GIT_REPOSITORY https://github.com/libigl/libigl-glad.git
-    GIT_TAG        09b4969c56779f7ddf8e6176ec1873184aec890f
+    GIT_TAG        651a425101365aa6e8504988ef9bb363d066c5ee
 )
 
 FetchContent_MakeAvailable(glad)

cmake/recipes/external/gmp.cmake

@@ -1,54 +0,0 @@
-if(TARGET gmp::gmp)
-    return()
-endif()
-
-# Download precompiled .dll on Windows
-if(WIN32)
-  include(gmp_mpfr)
-  # Find_package will look for our downloaded lib on Windows, and system-wide on Linux/macOS
-  find_package(GMP REQUIRED)
-else()
-  message(STATUS "Third-party: creating target 'gmp::gmp'")
-
-  include(FetchContent)
-  include(ProcessorCount)
-  ProcessorCount(Ncpu)
-  include(ExternalProject)
-  set(prefix ${FETCHCONTENT_BASE_DIR}/gmp)
-  set(gmp_INSTALL ${prefix}/install)
-  set(gmp_LIB_DIR ${gmp_INSTALL}/lib)
-  set(gmp_LIBRARY ${gmp_LIB_DIR}/${CMAKE_STATIC_LIBRARY_PREFIX}gmp${CMAKE_STATIC_LIBRARY_SUFFIX})
-  set(gmp_INCLUDE_DIR ${gmp_INSTALL}/include)
-
-  ExternalProject_Add(gmp
-    PREFIX ${prefix}
-    URL  https://gmplib.org/download/gmp/gmp-6.2.1.tar.xz
-    URL_MD5 0b82665c4a92fd2ade7440c13fcaa42b
-    UPDATE_DISCONNECTED true  # need this to avoid constant rebuild
-    PATCH_COMMAND 
-      curl "https://gmplib.org/repo/gmp/raw-rev/5f32dbc41afc" "|" git apply -v
-    CONFIGURE_HANDLED_BY_BUILD ON  # avoid constant reconfigure
-    CONFIGURE_COMMAND 
-      ${prefix}/src/gmp/configure 
-      --disable-debug --disable-dependency-tracking --enable-cxx --with-pic
-      --prefix=${gmp_INSTALL}
-      --disable-shared
-    BUILD_COMMAND make -j${Ncpu}
-    INSTALL_COMMAND make -j${Ncpu} install
-    INSTALL_DIR ${gmp_INSTALL}
-    TEST_COMMAND ""
-    BUILD_BYPRODUCTS ${gmp_LIBRARY}
-  )
-  ExternalProject_Get_Property(gmp SOURCE_DIR)
-  set(gmp_LIBRARIES ${gmp_LIBRARY})
-  add_library(gmp::gmp INTERFACE IMPORTED GLOBAL)
-  file(MAKE_DIRECTORY ${gmp_INCLUDE_DIR})  # avoid race condition
-  target_include_directories(gmp::gmp INTERFACE ${gmp_INCLUDE_DIR})
-  target_link_libraries(gmp::gmp INTERFACE "${gmp_LIBRARIES}")  # need the quotes to expand list
-  add_dependencies(gmp::gmp gmp)
-
-endif()
-
-if(NOT TARGET gmp::gmp)
-    message(FATAL_ERROR "Creation of target 'gmp::gmp' failed")
-endif()

cmake/recipes/external/gmp_mpfr.cmake

@@ -1,34 +0,0 @@
-if(WIN32)
-    message(STATUS "Third-party: downloading gmp + mpfr")
-
-    include(FetchContent)
-
-    # CGAL 5+ ships with a single .zip combining GMP + MPFR's precompiled dlls.
-    # For now we still download them separately.
-
-    FetchContent_Declare(
-        gmp
-        URL     https://cgal.geometryfactory.com/CGAL/precompiled_libs/auxiliary/x64/GMP/5.0.1/gmp-all-CGAL-3.9.zip
-        URL_MD5 508c1292319c832609329116a8234c9f
-    )
-    FetchContent_MakeAvailable(gmp)
-
-    FetchContent_Declare(
-        mpfr
-        URL https://cgal.geometryfactory.com/CGAL/precompiled_libs/auxiliary/x64/MPFR/3.0.0/mpfr-all-CGAL-3.9.zip
-        URL_MD5 48840454eef0ff18730050c05028734b
-    )
-    FetchContent_MakeAvailable(mpfr)
-
-    # FetchContent_Declare(
-    #     gmp_mpfr
-    #     URL https://github.com/CGAL/cgal/releases/download/v5.2.1/CGAL-5.2.1-win64-auxiliary-libraries-gmp-mpfr.zip
-    #     URL_MD5 247f4dca741c6b9a9be76286414070fa
-    # )
-
-    # For CGAL
-    set(ENV{GMP_DIR} "${gmp_SOURCE_DIR}")
-    set(ENV{MPFR_DIR} "${mpfr_SOURCE_DIR}")
-else()
-    # On Linux/macOS, gmp+mpfr will be fetched and compiled
-endif()

cmake/recipes/external/imgui.cmake

@@ -9,6 +9,7 @@ FetchContent_Declare(
     imgui
     GIT_REPOSITORY https://github.com/ocornut/imgui.git
     GIT_TAG v1.85
+    GIT_SHALLOW TRUE
 )
 FetchContent_MakeAvailable(imgui)
 

cmake/recipes/external/libigl_tests_data.cmake

@@ -6,13 +6,13 @@ message(STATUS "Third-party: creating target 'igl::tests_data'")
 
 include(FetchContent)
 FetchContent_Declare(
-    libigl_tests_tata
+    libigl_tests_data
     GIT_REPOSITORY https://github.com/libigl/libigl-tests-data
     GIT_TAG        19cedf96d70702d8b3a83eb27934780c542356fe
 )
-FetchContent_MakeAvailable(libigl_tests_tata)
+FetchContent_MakeAvailable(libigl_tests_data)
 
 add_library(igl_tests_data INTERFACE)
 add_library(igl::tests_data ALIAS igl_tests_data)
 
-target_compile_definitions(igl_tests_data INTERFACE LIBIGL_DATA_DIR=\"${libigl_tests_tata_SOURCE_DIR}\")
+target_compile_definitions(igl_tests_data INTERFACE LIBIGL_DATA_DIR=\"${libigl_tests_data_SOURCE_DIR}\")

cmake/recipes/external/libigl_tutorial_data.cmake

@@ -6,13 +6,13 @@ message(STATUS "Third-party: creating target 'igl::tutorial_data'")
 
 include(FetchContent)
 FetchContent_Declare(
-    libigl_tutorial_tata
+    libigl_tutorial_data
     GIT_REPOSITORY https://github.com/libigl/libigl-tutorial-data
-    GIT_TAG        c1f9ede366d02e3531ecbaec5e3769312f31cccd
+    GIT_TAG        644dd4104843b6d736745d9dafbd70bf8d175648
 )
-FetchContent_MakeAvailable(libigl_tutorial_tata)
+FetchContent_MakeAvailable(libigl_tutorial_data)
 
 add_library(igl_tutorial_data INTERFACE)
 add_library(igl::tutorial_data ALIAS igl_tutorial_data)
 
-target_compile_definitions(igl_tutorial_data INTERFACE "-DTUTORIAL_SHARED_PATH=\"${libigl_tutorial_tata_SOURCE_DIR}\"")
+target_compile_definitions(igl_tutorial_data INTERFACE "-DTUTORIAL_SHARED_PATH=\"${libigl_tutorial_data_SOURCE_DIR}\"")

cmake/recipes/external/mpfr.cmake

@@ -1,61 +0,0 @@
-# Expects
-#   gmp_INCLUDE_DIR
-#   gmp_LIB_DIR
-#   gmp_LIBRARIES
-if(TARGET mpfr::mpfr)
-    return()
-endif()
-
-# Download precompiled .dll on Windows
-if(WIN32)
-  include(gmp_mpfr)
-  # Find_package will look for our downloaded lib on Windows, and system-wide on Linux/macOS
-  find_package(MPFR REQUIRED)
-else()
-  message(STATUS "Third-party: creating target 'mpfr::mpfr'")
-
-  include(FetchContent)
-  include(ProcessorCount)
-  ProcessorCount(Ncpu)
-  include(ExternalProject)
-  set(prefix ${FETCHCONTENT_BASE_DIR}/mpfr)
-  set(mpfr_INSTALL ${prefix}/install)
-  set(mpfr_LIBRARY ${mpfr_INSTALL}/lib/${CMAKE_STATIC_LIBRARY_PREFIX}mpfr${CMAKE_STATIC_LIBRARY_SUFFIX})
-  set(mpfr_INCLUDE_DIR ${mpfr_INSTALL}/include)
-
-  ExternalProject_Add(mpfr
-    PREFIX ${prefix}
-    DEPENDS gmp
-    URL  https://ftp.gnu.org/gnu/mpfr/mpfr-4.1.0.tar.xz
-    URL_MD5 bdd3d5efba9c17da8d83a35ec552baef
-    UPDATE_DISCONNECTED true  # need this to avoid constant rebuild
-    CONFIGURE_HANDLED_BY_BUILD ON  # avoid constant reconfigure
-    CONFIGURE_COMMAND 
-      ${prefix}/src/mpfr/configure 
-      --disable-debug --disable-dependency-tracking  --disable-silent-rules --enable-cxx --with-pic
-      --with-gmp-include=${gmp_INCLUDE_DIR} --with-gmp-lib=${gmp_LIB_DIR}
-      --disable-shared
-      --prefix=${mpfr_INSTALL}
-      --disable-shared
-    BUILD_COMMAND make -j${Ncpu}
-    INSTALL_COMMAND make -j${Ncpu} install
-    INSTALL_DIR ${mpfr_INSTALL}
-    TEST_COMMAND ""
-    BUILD_BYPRODUCTS ${mpfr_LIBRARY}
-  )
-  #PATCH_COMMAND  curl "https://raw.githubusercontent.com/Homebrew/formula-patches/03cf8088210822aa2c1ab544ed58ea04c897d9c4/libtool/configure-big_sur.diff" "|" sed -e "s/configure.orig/configure/g" "|" git apply -v
-  ExternalProject_Get_Property(mpfr SOURCE_DIR)
-  set(mpfr_LIBRARIES ${mpfr_LIBRARY})
-  add_library(mpfr::mpfr INTERFACE IMPORTED GLOBAL)
-  file(MAKE_DIRECTORY ${mpfr_INCLUDE_DIR})  # avoid race condition
-  target_include_directories(mpfr::mpfr INTERFACE ${mpfr_INCLUDE_DIR})
-  target_link_libraries(mpfr::mpfr INTERFACE "${mpfr_LIBRARIES}")  # need the quotes to expand list
-  # This is necessary to ensure that mpfr appears before gmp in link order.
-  # Otherwise undefined reference errors occur at link time on Linux with gcc
-  target_link_libraries(mpfr::mpfr INTERFACE "${gmp_LIBRARIES}") 
-  add_dependencies(mpfr::mpfr mpfr)
-endif()
-
-if(NOT TARGET mpfr::mpfr)
-    message(FATAL_ERROR "Creation of target 'mpfr::mpfr' failed")
-endif()

cmake/recipes/external/predicates.cmake

@@ -8,7 +8,7 @@ include(FetchContent)
 FetchContent_Declare(
     predicates
     GIT_REPOSITORY https://github.com/libigl/libigl-predicates.git
-    GIT_TAG        488242fa2b1f98a9c5bd1441297fb4a99a6a9ae4
+    GIT_TAG        decb7bc1260e689cbe008109e3cc5d3a5a433aea
 )
 
 FetchContent_MakeAvailable(predicates)

cmake/recipes/external/spectra.cmake

@@ -0,0 +1,16 @@
+if(TARGET spectra::spectra)
+  return()
+endif()
+include(FetchContent)
+
+message(STATUS "Third-party: creating target 'spectra::spectra'")
+
+# Use fork because yixuan/spectra struggles to find Eigen3
+FetchContent_Declare(
+  Spectra
+  GIT_REPOSITORY https://github.com/alecjacobson/spectra/
+  GIT_TAG bbdc521b70a733c52ebfc0ac1484c82e13c3d140
+)
+FetchContent_MakeAvailable(Spectra)
+
+add_library(spectra::spectra ALIAS Spectra)

cmake/recipes/external/tetgen.cmake

@@ -8,7 +8,7 @@ include(FetchContent)
 FetchContent_Declare(
     tetgen
     GIT_REPOSITORY https://github.com/libigl/tetgen.git
-    GIT_TAG        4f3bfba3997f20aa1f96cfaff604313a8c2c85b6
+    GIT_TAG        e05aca7df74e3f531bc35733ed87d36d437266c5
 )
 
 FetchContent_MakeAvailable(tetgen)

cmake/recipes/external/triangle.cmake

@@ -8,7 +8,7 @@ include(FetchContent)
 FetchContent_Declare(
     triangle
     GIT_REPOSITORY https://github.com/libigl/triangle.git
-    GIT_TAG        3ee6cac2230f0fe1413879574f741c7b6da11221
+    GIT_TAG        62f02db9ab4ff4b62d5ff82a77c8ea458c84c23a
 )
 
 FetchContent_MakeAvailable(triangle)

docs/index.md

@@ -0,0 +1,54 @@
+# libigl - A simple C++ geometry processing library
+
+This detailed documentation browser is automatically generated from the comments
+in libigl header (.h) files.
+
+In general, each [libigl function](./namespaceigl.html#func-members) (e.g., `igl::func`) will be defined in a
+correspondingly named header file (e.g., `#include <igl/func.h>`).
+
+The _core_ library only depends on the standard template library (`std::`) and 
+Eigen. These functions reside directly the [`igl::` namespace](./namespaceigl.html)
+
+Functions with further dependencies reside in a corresonding sub-namespace. For
+example, the function `igl::spectra::lscm` depends on the Spectra library so it
+resides in the [`igl::spectra::` namespace](./namespaceigl_1_1spectra.html).
+
+Functions which depend on external code under a copyleft license reside in the
+[`igl::copyleft::` namepsace](file:///Users/alecjacobson/Repos/libigl/dox/namespaceigl_1_1copyleft.html).
+
+
+Most libigl functions are templated over the Eigen matrix inputs and outputs.
+Callers can choose their own scalar types (e.g., `double`/`float`) and storage
+orders (`Eigen::ColMajor`/`Eigen::RowMajor`). Libigl can be used as a:
+
+ - **header only library** (via CMake, make sure
+   `LIBIGL_USE_STATIC_LIBRARY=OFF`) and insure that `IGL_STATIC_LIBRARY` is
+   _not_ defined_ when compiling --- easiest if you're new to libigl, or 
+ - **static library** (`LIBIGL_USE_STATIC_LIBRARY=ON` → `IGL_STATIC_LIBRARY` is
+  defined) --- speeds up repeated compilation.
+
+The libigl static library is filled with _explicit template instantiations_ for
+common Eigen inputs and outputs. If the library doesn't contain your types, you
+may get some form of linker error (e.g., `Undefined symbols for architecture`,
+`undefined reference to` or `unresolved external symbol`).
+
+You can fix this by:
+
+1. Switching to header only mode for your project,
+2. Making  a file in your project to compile the missing templates. E.g., `my_templates.cpp` 
+```cpp
+#ifdef IGL_STATIC_LIBRARY
+#undef IGL_STATIC_LIBRARY
+#endif
+#include <igl/per_vertex_normals.h>
+template void igl::per_vertex_normals<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&); 
+```
+3. [Submit a PR](https://github.com/libigl/libigl/pulls) containing your missing template to the development branch of libigl
+4. Change your input/output types to match existing templates (e.g., `Eigen::MatrixXd`).
+
+
+
+
+https://libigl.github.io/
+
+https://github.com/libigl/libigl/

include/igl/ARAPEnergyType.h

@@ -9,27 +9,24 @@
 #define IGL_ARAPENERGYTYPE_H
 namespace igl
 {
-  //     ARAP_ENERGY_TYPE_SPOKES  "As-rigid-as-possible Surface Modeling" by [Sorkine and
-  //       Alexa 2007], rotations defined at vertices affecting incident edges,
-  //       default
-  //     ARAP_ENERGY_TYPE_SPOKES-AND-RIMS  Adapted version of "As-rigid-as-possible Surface
-  //       Modeling" by [Sorkine and Alexa 2007] presented in section 4.2 of or
-  //       "A simple geometric model for elastic deformation" by [Chao et al.
-  //       2010], rotations defined at vertices affecting incident edges and
-  //       opposite edges
-  //     ARAP_ENERGY_TYPE_ELEMENTS  "A local-global approach to mesh parameterization" by
-  //       [Liu et al.  2010] or "A simple geometric model for elastic
-  //       deformation" by [Chao et al.  2010], rotations defined at elements
-  //       (triangles or tets) 
-  //     ARAP_ENERGY_TYPE_DEFAULT  Choose one automatically: spokes and rims
-  //       for surfaces, elements for planar meshes and tets (not fully
-  //       supported)
+  /// Enum for choosing ARAP energy type
   enum ARAPEnergyType
   {
+    /// "As-rigid-as-possible Surface Modeling" by [Sorkine and Alexa 2007],
+    /// rotations defined at vertices affecting incident edges, default
     ARAP_ENERGY_TYPE_SPOKES = 0,
+    /// Adapted version of "As-rigid-as-possible Surface Modeling" by [Sorkine
+    /// and Alexa 2007] presented in section 4.2 of or "A simple geometric model
+    /// for elastic deformation" by [Chao et al.\ 2010], rotations defined at
+    /// vertices affecting incident edges and opposite edges
     ARAP_ENERGY_TYPE_SPOKES_AND_RIMS = 1,
+    /// "A local-global approach to mesh parameterization" by [Liu et al.\ 2010]
+    /// or "A simple geometric model for elastic deformation" by [Chao et al.\ 2010], rotations defined at elements (triangles or tets) 
     ARAP_ENERGY_TYPE_ELEMENTS = 2,
+    /// Choose one automatically: spokes and rims for surfaces, elements for
+    /// planar meshes and tets (not fully supported)
     ARAP_ENERGY_TYPE_DEFAULT = 3,
+    /// Total number of types
     NUM_ARAP_ENERGY_TYPES = 4
   };
 }

include/igl/AtA_cached.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 "AtA_cached.h"
+#include "IGL_ASSERT.h"
 
 #include <iostream>
 #include <vector>
@@ -34,12 +35,12 @@ IGL_INLINE void igl::AtA_cached_precompute(
       int col = k;
       int row = *(A.innerIndexPtr()+l);
       int value_index = l;
-      assert(col < A.cols());
-      assert(col >= 0);
-      assert(row < A.rows());
-      assert(row >= 0);
-      assert(value_index >= 0);
-      assert(value_index < A.nonZeros());
+      IGL_ASSERT(col < A.cols());
+      IGL_ASSERT(col >= 0);
+      IGL_ASSERT(row < A.rows());
+      IGL_ASSERT(row >= 0);
+      IGL_ASSERT(value_index >= 0);
+      IGL_ASSERT(value_index < A.nonZeros());
 
       Col_RowPtr[col].push_back(row);
       Col_IndexPtr[col].push_back(value_index);
@@ -74,12 +75,12 @@ IGL_INLINE void igl::AtA_cached_precompute(
       int col = k;
       int row = *(AtA.innerIndexPtr()+l);
       int value_index = l;
-      assert(col < AtA.cols());
-      assert(col >= 0);
-      assert(row < AtA.rows());
-      assert(row >= 0);
-      assert(value_index >= 0);
-      assert(value_index < AtA.nonZeros());
+      IGL_ASSERT(col < AtA.cols());
+      IGL_ASSERT(col >= 0);
+      IGL_ASSERT(row < AtA.rows());
+      IGL_ASSERT(row >= 0);
+      IGL_ASSERT(value_index >= 0);
+      IGL_ASSERT(value_index < AtA.nonZeros());
 
       data.I_outer.push_back(data.I_row.size());
 

include/igl/AtA_cached.h

@@ -13,40 +13,47 @@
 #include <Eigen/Sparse>
 namespace igl
 {  
+  /// Hold precomputed data for AtA_cached
   struct AtA_cached_data
   {
-    // Weights
+    /// Weights (diagonal of W)
     Eigen::VectorXd W;
 
     // Flatten composition rules
+    /// @private
     std::vector<int> I_row;
+    /// @private
     std::vector<int> I_col;
+    /// @private
     std::vector<int> I_w;
 
     // For each entry of AtA, points to the beginning
     // of the composition rules
+    /// @private
     std::vector<int> I_outer;
   };
 
-  // Computes At * W * A, where A is sparse and W is diagonal. Divides the 
-  // construction in two phases, one
-  // for fixing the sparsity pattern, and one to populate it with values. Compared to
-  // evaluating it directly, this version is slower for the first time (since it requires a
-  // precomputation), but faster to the subsequent evaluations.
-  //
-  // Input:
-  //   A m x n sparse matrix
-  //   data stores the precomputed sparsity pattern, data.W contains the optional diagonal weights (stored as a dense vector). If W is not provided, it is replaced by the identity.
-  // Outputs:
-  //   AtA  m by m matrix computed as AtA * W * A
-  //
-  // Example:
-  // AtA_data = igl::AtA_cached_data();
-  // AtA_data.W = W;
-  // if (s.AtA.rows() == 0)
-  //   igl::AtA_cached_precompute(s.A,s.AtA_data,s.AtA);
-  // else
-  //   igl::AtA_cached(s.A,s.AtA_data,s.AtA);
+  /// Computes At * W * A, where A is sparse and W is diagonal.
+  ///
+  /// Divides the construction in two phases, one for fixing the sparsity
+  /// pattern, and one to populate it with values. Compared to evaluating it
+  /// directly, this version is slower for the first time (since it requires a
+  /// precomputation), but faster to the subsequent evaluations.
+  ///
+  /// @param[in] A m x n sparse matrix
+  /// @param[in,out] data stores the precomputed sparsity pattern, data.W contains the optional diagonal weights (stored as a dense vector). If W is not provided, it is replaced by the identity.
+  /// @param[out] AtA  m by m matrix computed as AtA * W * A
+  ///
+  /// #### Example:
+  ///
+  /// \code{cpp}
+  /// AtA_data = igl::AtA_cached_data();
+  /// AtA_data.W = W;
+  /// if (s.AtA.rows() == 0)
+  ///   igl::AtA_cached_precompute(s.A,s.AtA_data,s.AtA);
+  /// else
+  ///   igl::AtA_cached(s.A,s.AtA_data,s.AtA);
+  /// \endcode
   template <typename Scalar>
   IGL_INLINE void AtA_cached_precompute(
     const Eigen::SparseMatrix<Scalar>& A,
@@ -54,6 +61,11 @@ namespace igl
     Eigen::SparseMatrix<Scalar>& AtA
     );
 
+  /// Computes At * W * A, where A is sparse and W is diagonal precomputed into data.
+  ///
+  /// @param[in] A m x n sparse matrix
+  /// @param[in]  data stores the precomputed sparsity pattern, data.W contains the optional diagonal weights (stored as a dense vector). If W is not provided, it is replaced by the identity.
+  /// @param[out] AtA  m by m matrix computed as AtA * W * A
   template <typename Scalar>
   IGL_INLINE void AtA_cached(
     const Eigen::SparseMatrix<Scalar>& A,

include/igl/COLLAPSE_EDGE_NULL.h

@@ -0,0 +1,19 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2025 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_COLLAPSE_EDGE_NULL_H
+#define IGL_COLLAPSE_EDGE_NULL_H
+namespace igl
+{
+#ifndef IGL_COLLAPSE_EDGE_NULL
+  /// Special value for indicating a null vertex index as the result of a
+  /// collapsed edge.
+  #define IGL_COLLAPSE_EDGE_NULL 0
+#endif
+}
+#endif

include/igl/C_STR.h

@@ -7,12 +7,23 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_C_STR_H
 #define IGL_C_STR_H
-// http://stackoverflow.com/a/2433143/148668
-// Suppose you have a function:
-//   void func(const char * c);
-// Then you can write:
-//   func(C_STR("foo"<<1<<"bar"));
 #include <sstream>
 #include <string>
+/// Convert a stream of things to a const char *.
+///
+/// Suppose you have a function:
+/// \code{cpp}
+/// void func(const char * c);
+/// \endcode
+/// Then you can write:
+/// \code{cpp}
+///   func(C_STR("foo"<<1<<"bar"));
+/// \endcode
+/// which is equivalent to:
+/// \code{cpp}
+///   func("foo1bar");
+/// \endcode
+///
+// http://stackoverflow.com/a/2433143/148668
 #define C_STR(X) static_cast<std::ostringstream&>(std::ostringstream().flush() << X).str().c_str()
 #endif

include/igl/Camera.h

@@ -1,359 +0,0 @@
-// 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_CAMERA_H
-#define IGL_CAMERA_H
-
-// you're idiot, M$!
-#if defined(_WIN32)
-#undef far
-#undef near
-#endif
-
-#include <Eigen/Geometry>
-#include <Eigen/Core>
-#include "PI.h"
-
-#define IGL_CAMERA_MIN_ANGLE 5.0
-namespace igl
-{
-
-  // A simple camera class. The camera stores projection parameters (field of
-  // view angle, aspect ratio, near and far clips) as well as a rigid
-  // transformation *of the camera as if it were also a scene object*. Thus, the
-  // **inverse** of this rigid transformation is the modelview transformation.
-  class Camera
-  {
-    public:
-      // On windows you might need: -fno-delayed-template-parsing
-      //static constexpr double IGL_CAMERA_MIN_ANGLE = 5.;
-      //  m_angle  Field of view angle in degrees {45}
-      //  m_aspect  Aspect ratio {1}
-      //  m_near  near clipping plane {1e-2}
-      //  m_far  far clipping plane {100}
-      //  m_at_dist  distance of looking at point {1}
-      //  m_orthographic  whether to use othrographic projection {false}
-      //  m_rotation_conj  Conjugate of rotation part of rigid transformation of
-      //    camera {identity}. Note: we purposefully store the conjugate because
-      //    this is what TW_TYPE_QUAT4D is expecting.
-      //  m_translation  Translation part of rigid transformation of camera
-      //    {(0,0,1)}
-      double m_angle, m_aspect, m_near, m_far, m_at_dist;
-      bool m_orthographic;
-      Eigen::Quaterniond m_rotation_conj;
-      Eigen::Vector3d m_translation;
-    public:
-      inline Camera();
-      inline virtual ~Camera(){}
-      // Return projection matrix that takes relative camera coordinates and
-      // transforms it to viewport coordinates
-      //
-      // Note:
-      //
-      //     if(m_angle > 0)
-      //     {
-      //       gluPerspective(m_angle,m_aspect,m_near,m_at_dist+m_far);
-      //     }else
-      //     {
-      //       gluOrtho(-0.5*aspect,0.5*aspect,-0.5,0.5,m_at_dist+m_near,m_far);
-      //     }
-      //
-      // Is equivalent to
-      //
-      //     glMultMatrixd(projection().data());
-      //
-      inline Eigen::Matrix4d projection() const;
-      // Return an Affine transformation (rigid actually) that 
-      // takes relative coordinates and tramsforms them into world 3d
-      // coordinates: moves the camera into the scene.
-      inline Eigen::Affine3d affine() const;
-      // Return an Affine transformation (rigid actually) that puts the takes a
-      // world 3d coordinate and transforms it into the relative camera
-      // coordinates: moves the scene in front of the camera.
-      //
-      // Note:
-      //
-      //     gluLookAt(
-      //       eye()(0), eye()(1), eye()(2),
-      //       at()(0), at()(1), at()(2),
-      //       up()(0), up()(1), up()(2));
-      //
-      // Is equivalent to
-      //
-      //     glMultMatrixd(camera.inverse().matrix().data());
-      //
-      // See also: affine, eye, at, up
-      inline Eigen::Affine3d inverse() const;
-      // Returns world coordinates position of center or "eye" of camera.
-      inline Eigen::Vector3d eye() const;
-      // Returns world coordinate position of a point "eye" is looking at.
-      inline Eigen::Vector3d at() const;
-      // Returns world coordinate unit vector of "up" vector
-      inline Eigen::Vector3d up() const;
-      // Return top right corner of unit plane in relative coordinates, that is
-      // (w/2,h/2,1)
-      inline Eigen::Vector3d unit_plane() const;
-      // Move dv in the relative coordinate frame of the camera (move the FPS)
-      //
-      // Inputs:
-      //   dv  (x,y,z) displacement vector
-      //
-      inline void dolly(const Eigen::Vector3d & dv);
-      // "Scale zoom": Move `eye`, but leave `at`
-      //
-      // Input:
-      //   s  amount to scale distance to at
-      inline void push_away(const double s);
-      // Aka "Hitchcock", "Vertigo", "Spielberg" or "Trombone" zoom:
-      // simultaneously dolly while changing angle so that `at` not only stays
-      // put in relative coordinates but also projected coordinates. That is
-      //
-      // Inputs:
-      //   da  change in angle in degrees
-      inline void dolly_zoom(const double da);
-      // Turn around eye so that rotation is now q
-      //
-      // Inputs:
-      //   q  new rotation as quaternion
-      inline void turn_eye(const Eigen::Quaterniond & q);
-      // Orbit around at so that rotation is now q
-      //
-      // Inputs:
-      //   q  new rotation as quaternion
-      inline void orbit(const Eigen::Quaterniond & q);
-      // Rotate and translate so that camera is situated at "eye" looking at "at"
-      // with "up" pointing up.
-      //
-      // Inputs:
-      //   eye  (x,y,z) coordinates of eye position
-      //   at   (x,y,z) coordinates of at position
-      //   up   (x,y,z) coordinates of up vector
-      inline void look_at(
-        const Eigen::Vector3d & eye,
-        const Eigen::Vector3d & at,
-        const Eigen::Vector3d & up);
-    // Needed any time Eigen Structures are used as class members
-    // http://eigen.tuxfamily.org/dox-devel/group__TopicStructHavingEigenMembers.html
-    public:
-        EIGEN_MAKE_ALIGNED_OPERATOR_NEW
-  };
-}
-
-// Implementation
-#include "PI.h"
-#include "EPS.h"
-#include <cmath>
-#include <iostream>
-#include <cassert>
-
-inline igl::Camera::Camera():
-  m_angle(45.0),m_aspect(1),m_near(1e-2),m_far(100),m_at_dist(1),
-  m_orthographic(false),
-  m_rotation_conj(1,0,0,0),
-  m_translation(0,0,1)
-{
-}
-
-inline Eigen::Matrix4d igl::Camera::projection() const
-{
-  Eigen::Matrix4d P;
-  using namespace std;
-  const double far = m_at_dist + m_far;
-  const double near = m_near;
-  // http://stackoverflow.com/a/3738696/148668
-  if(m_orthographic)
-  {
-    const double f = 0.5;
-    const double left = -f*m_aspect;
-    const double right = f*m_aspect;
-    const double bottom = -f;
-    const double top = f;
-    const double tx = (right+left)/(right-left);
-    const double ty = (top+bottom)/(top-bottom);
-    const double tz = (far+near)/(far-near);
-    const double z_fix = 0.5 /m_at_dist / tan(m_angle*0.5 * (igl::PI/180.) );
-    P<<
-      z_fix*2./(right-left), 0, 0, -tx,
-      0, z_fix*2./(top-bottom), 0, -ty,
-      0, 0, -z_fix*2./(far-near),  -tz,
-      0, 0, 0, 1;
-  }else
-  {
-    const double yScale = tan(PI*0.5 - 0.5*m_angle*PI/180.);
-    // http://stackoverflow.com/a/14975139/148668
-    const double xScale = yScale/m_aspect;
-    P<< 
-      xScale, 0, 0, 0,
-      0, yScale, 0, 0,
-      0, 0, -(far+near)/(far-near), -1,
-      0, 0, -2.*near*far/(far-near), 0;
-    P = P.transpose().eval();
-  }
-  return P;
-}
-
-inline Eigen::Affine3d igl::Camera::affine() const
-{
-  using namespace Eigen;
-  Affine3d t = Affine3d::Identity();
-  t.rotate(m_rotation_conj.conjugate());
-  t.translate(m_translation);
-  return t;
-}
-
-inline Eigen::Affine3d igl::Camera::inverse() const
-{
-  using namespace Eigen;
-  Affine3d t = Affine3d::Identity();
-  t.translate(-m_translation);
-  t.rotate(m_rotation_conj);
-  return t;
-}
-
-inline Eigen::Vector3d igl::Camera::eye() const
-{
-  using namespace Eigen;
-  return affine() * Vector3d(0,0,0);
-}
-
-inline Eigen::Vector3d igl::Camera::at() const
-{
-  using namespace Eigen;
-  return affine() * (Vector3d(0,0,-1)*m_at_dist);
-}
-
-inline Eigen::Vector3d igl::Camera::up() const
-{
-  using namespace Eigen;
-  Affine3d t = Affine3d::Identity();
-  t.rotate(m_rotation_conj.conjugate());
-  return t * Vector3d(0,1,0);
-}
-
-inline Eigen::Vector3d igl::Camera::unit_plane() const
-{
-  // Distance of center pixel to eye
-  const double d = 1.0;
-  const double a = m_aspect;
-  const double theta = m_angle*PI/180.;
-  const double w =
-    2.*sqrt(-d*d/(a*a*pow(tan(0.5*theta),2.)-1.))*a*tan(0.5*theta);
-  const double h = w/a;
-  return Eigen::Vector3d(w*0.5,h*0.5,-d);
-}
-
-inline void igl::Camera::dolly(const Eigen::Vector3d & dv)
-{
-  m_translation += dv;
-}
-
-inline void igl::Camera::push_away(const double s)
-{
-  using namespace Eigen;
-#ifndef NDEBUG
-  Vector3d old_at = at();
-#endif
-  const double old_at_dist = m_at_dist;
-  m_at_dist = old_at_dist * s;
-  dolly(Vector3d(0,0,1)*(m_at_dist - old_at_dist));
-  assert((old_at-at()).squaredNorm() < DOUBLE_EPS);
-}
-
-inline void igl::Camera::dolly_zoom(const double da)
-{
-  using namespace std;
-  using namespace Eigen;
-#ifndef NDEBUG
-  Vector3d old_at = at();
-#endif
-  const double old_angle = m_angle;
-  if(old_angle + da < IGL_CAMERA_MIN_ANGLE)
-  {
-    m_orthographic = true;
-  }else if(old_angle + da > IGL_CAMERA_MIN_ANGLE)
-  {
-    m_orthographic = false;
-  }
-  if(!m_orthographic)
-  {
-    m_angle += da;
-    m_angle = min(89.,max(IGL_CAMERA_MIN_ANGLE,m_angle));
-    // change in distance
-    const double s = 
-      (2.*tan(old_angle/2./180.*igl::PI)) /
-      (2.*tan(m_angle/2./180.*igl::PI)) ;
-    const double old_at_dist = m_at_dist;
-    m_at_dist = old_at_dist * s;
-    dolly(Vector3d(0,0,1)*(m_at_dist - old_at_dist));
-    assert((old_at-at()).squaredNorm() < DOUBLE_EPS);
-  }
-}
-
-inline void igl::Camera::turn_eye(const Eigen::Quaterniond & q)
-{
-  using namespace Eigen;
-  Vector3d old_eye = eye();
-  // eye should be fixed
-  //
-  // eye_1 = R_1 * t_1 = eye_0
-  // t_1 = R_1' * eye_0
-  m_rotation_conj = q.conjugate();
-  m_translation = m_rotation_conj * old_eye;
-  assert((old_eye - eye()).squaredNorm() < DOUBLE_EPS);
-}
-
-inline void igl::Camera::orbit(const Eigen::Quaterniond & q)
-{
-  using namespace Eigen;
-  Vector3d old_at = at();
-  // at should be fixed
-  //
-  // at_1 = R_1 * t_1 - R_1 * z = at_0
-  // t_1 = R_1' * (at_0 + R_1 * z)
-  m_rotation_conj = q.conjugate();
-  m_translation = 
-    m_rotation_conj * 
-      (old_at + 
-         m_rotation_conj.conjugate() * Vector3d(0,0,1) * m_at_dist);
-  assert((old_at - at()).squaredNorm() < DOUBLE_EPS);
-}
-
-inline void igl::Camera::look_at(
-  const Eigen::Vector3d & eye,
-  const Eigen::Vector3d & at,
-  const Eigen::Vector3d & up)
-{
-  using namespace Eigen;
-  using namespace std;
-  // http://www.opengl.org/sdk/docs/man2/xhtml/gluLookAt.xml
-  // Normalize vector from at to eye
-  Vector3d F = eye-at;
-  m_at_dist = F.norm();
-  F.normalize();
-  // Project up onto plane orthogonal to F and normalize
-  assert(up.cross(F).norm() > DOUBLE_EPS && "(eye-at) x up ≈ 0");
-  const Vector3d proj_up = (up-(up.dot(F))*F).normalized();
-  Quaterniond a,b;
-  a.setFromTwoVectors(Vector3d(0,0,-1),-F);
-  b.setFromTwoVectors(a*Vector3d(0,1,0),proj_up);
-  m_rotation_conj = (b*a).conjugate();
-  m_translation = m_rotation_conj * eye;
-  //cout<<"m_at_dist: "<<m_at_dist<<endl;
-  //cout<<"proj_up: "<<proj_up.transpose()<<endl;
-  //cout<<"F: "<<F.transpose()<<endl;
-  //cout<<"eye(): "<<this->eye().transpose()<<endl;
-  //cout<<"at(): "<<this->at().transpose()<<endl;
-  //cout<<"eye()-at(): "<<(this->eye()-this->at()).normalized().transpose()<<endl;
-  //cout<<"eye-this->eye(): "<<(eye-this->eye()).squaredNorm()<<endl;
-  assert(           (eye-this->eye()).squaredNorm() < DOUBLE_EPS);
-  //assert((F-(this->eye()-this->at()).normalized()).squaredNorm() < 
-  //  DOUBLE_EPS);
-  assert(           (at-this->at()).squaredNorm() < DOUBLE_EPS);
-  //assert(        (proj_up-this->up()).squaredNorm() < DOUBLE_EPS);
-}
-
-#endif

include/igl/EPS.h

@@ -10,13 +10,17 @@
 #include "igl_inline.h"
 namespace igl
 {
-  // Define a standard value for double epsilon
+  /// Standard value for double epsilon
   const double DOUBLE_EPS    = 1.0e-14;
+  /// Standard value for double epsilon²
   const double DOUBLE_EPS_SQ = 1.0e-28;
+  /// Standard value for single epsilon
   const float FLOAT_EPS    = 1.0e-7f;
+  /// Standard value for single epsilon²
   const float FLOAT_EPS_SQ = 1.0e-14f;
-  // Function returning EPS for corresponding type
+  /// Function returning EPS for corresponding type
   template <typename S_type> IGL_INLINE S_type EPS();
+  /// Function returning EPS_SQ for corresponding type
   template <typename S_type> IGL_INLINE S_type EPS_SQ();
   // Template specializations for float and double
   template <> IGL_INLINE float EPS<float>();

include/igl/FileEncoding.h

@@ -10,7 +10,7 @@
 
 namespace igl
 {
-
+  /// File encoding types for writing files.
 enum class FileEncoding {
   Binary,
   Ascii

include/igl/FileMemoryStream.h

@@ -32,7 +32,7 @@ namespace igl {
     pos_type seekoff(
       off_type off, 
       std::ios_base::seekdir dir,
-      std::ios_base::openmode which) override 
+      std::ios_base::openmode /*which*/) override 
     {
       if (dir == std::ios_base::cur) 
       {
@@ -50,6 +50,7 @@ namespace igl {
     }
   };
   
+  /// Class to convert a FILE * to an std::istream
   struct FileMemoryStream : virtual FileMemoryBuffer, public std::istream 
   {
     FileMemoryStream( char const *first_elem, size_t size)

include/igl/HalfEdgeIterator.cpp

@@ -8,6 +8,8 @@
 
 #include "HalfEdgeIterator.h"
 
+#include <cassert>
+
 template <typename DerivedF, typename DerivedFF, typename DerivedFFi>
 IGL_INLINE igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::HalfEdgeIterator(
     const Eigen::MatrixBase<DerivedF>& _F,
@@ -63,10 +65,10 @@ IGL_INLINE bool igl::HalfEdgeIterator<DerivedF,DerivedFF,DerivedFFi>::isBorder()
 /*!
  * Returns the next edge skipping the border
  *      _________
- *     /\ c | b /\
- *    /  \  |  /  \
- *   / d  \ | / a  \
- *  /______\|/______\
+ *     ╱╲ c | b ╱╲
+ *    ╱  ╲  |  ╱  ╲
+ *   ╱ d  ╲ | ╱ a  ╲
+ *  ╱______╲|╱______╲
  *          v
  * In this example, if a and d are of-border and the pos is iterating counterclockwise, this method iterate through the faces incident on vertex v,
  * producing the sequence a, b, c, d, a, b, c, ...

include/igl/HalfEdgeIterator.h

@@ -11,35 +11,24 @@
 #include <Eigen/Core>
 
 #include <vector>
-#include <igl/igl_inline.h>
+#include "igl_inline.h"
 
-// This file violates many of the libigl style guidelines.
 
 namespace igl
 {
-  // HalfEdgeIterator - Fake halfedge for fast and easy navigation
-  // on triangle meshes with vertex_triangle_adjacency and
-  // triangle_triangle adjacency
-  //
-  // Note: this is different to classical Half Edge data structure.
-  //    Instead, it follows cell-tuple in [Brisson, 1989]
-  //    "Representing geometric structures in d dimensions: topology and order."
-  //    This class can achieve local navigation similar to half edge in OpenMesh
-  //    But the logic behind each atom operation is different.
-  //    So this should be more properly called TriangleTupleIterator.
-  //
-  // Each tuple contains information on (face, edge, vertex)
-  //    and encoded by (face, edge \in {0,1,2}, bool reverse)
-  //
-  // Inputs:
-  //    F #F by 3 list of "faces"
-  //    FF #F by 3 list of triangle-triangle adjacency.
-  //    FFi #F by 3 list of FF inverse. For FF and FFi, refer to
-  //        "triangle_triangle_adjacency.h"
-  // Usages:
-  //    FlipF/E/V changes solely one actual face/edge/vertex resp.
-  //    NextFE iterates through one-ring of a vertex robustly.
-  //
+  /// Fake halfedge for fast and easy navigation
+  /// on triangle meshes with vertex_triangle_adjacency and
+  /// triangle_triangle adjacency
+  ///
+  /// Note: this is different to classical Half Edge data structure.
+  ///    Instead, it follows cell-tuple in [Brisson, 1989]
+  ///    "Representing geometric structures in d dimensions: topology and order."
+  ///    This class can achieve local navigation similar to half edge in OpenMesh
+  ///    But the logic behind each atom operation is different.
+  ///    So this should be more properly called TriangleTupleIterator.
+  ///
+  /// Each tuple contains information on (face, edge, vertex)
+  ///    and encoded by (face, edge \in {0,1,2}, bool reverse)
   template <
     typename DerivedF,
     typename DerivedFF,
@@ -47,7 +36,15 @@ namespace igl
   class HalfEdgeIterator
   {
   public:
-    // Init the HalfEdgeIterator by specifying Face,Edge Index and Orientation
+    /// Init the HalfEdgeIterator by specifying Face,Edge Index and Orientation
+    ///
+    /// @param[in] F #F by 3 list of "faces"
+    /// @param[in] FF #F by 3 list of triangle-triangle adjacency.
+    /// @param[in] FFi #F by 3 list of FF inverse. For FF and FFi, refer to
+    ///        "triangle_triangle_adjacency.h"
+    /// @param[in] _fi index of the selected face
+    /// @param[in] _ii index of the selected face
+    /// @param[in] _reverse orientation of the selected face
     IGL_INLINE HalfEdgeIterator(
         const Eigen::MatrixBase<DerivedF>& _F,
         const Eigen::MatrixBase<DerivedFF>& _FF,
@@ -57,41 +54,48 @@ namespace igl
         bool _reverse = false
         );
 
-    // Change Face
+    /// Change Face
     IGL_INLINE void flipF();
 
-    // Change Edge
+    /// Change Edge
     IGL_INLINE void flipE();
 
-    // Change Vertex
+    /// Change Vertex
     IGL_INLINE void flipV();
 
+    /// Determine if on border.
+    /// @returns true if the current edge is on the border
     IGL_INLINE bool isBorder();
 
-    /*!
-     * Returns the next edge skipping the border
-     *      _________
-     *     /\ c | b /\
-     *    /  \  |  /  \
-     *   / d  \ | / a  \
-     *  /______\|/______\
-     *          v
-     * In this example, if a and d are of-border and the pos is iterating
-     counterclockwise, this method iterate through the faces incident on vertex
-     v,
-     * producing the sequence a, b, c, d, a, b, c, ...
-     */
+    /// Change to next edge skipping the border
+    ///      _________
+    ///     ╱╲ c | b ╱╲
+    ///    ╱  ╲  |  ╱  ╲
+    ///   ╱ d  ╲ | ╱ a  ╲
+    ///  ╱______╲|╱______╲
+    ///          v
+    /// In this example, if a and d are of-border and the pos is iterating
+    /// counterclockwise, this method iterate through the faces incident on vertex
+    /// v,
+    /// producing the sequence a, b, c, d, a, b, c, ...
+    ///
+    /// @returns true if the next edge is not on the border
     IGL_INLINE bool NextFE();
 
-    // Get vertex index
+    /// Get vertex index
+    /// @return vertex index
     IGL_INLINE int Vi();
 
-    // Get face index
+    /// Get face index
+    /// @return face index
     IGL_INLINE int Fi();
 
-    // Get edge index
+    /// Get edge index
+    /// @return edge index
     IGL_INLINE int Ei();
 
+    /// Check if two HalfEdgeIterator are the same
+    /// @return true if two HalfEdgeIterator are the same
     IGL_INLINE bool operator==(HalfEdgeIterator& p2);
 
   private:

include/igl/Hit.h

@@ -11,19 +11,21 @@
 
 namespace igl
 {
-  // Reimplementation of the embree::Hit struct from embree1.0
-  // 
-  // TODO: template on floating point type
+  /// Reimplementation of the embree::Hit struct from embree1.0
+  /// 
+  template <typename Scalar>
   struct Hit
   {
-    int id; // primitive id
-    int gid; // geometry id (not used)
-    // barycentric coordinates so that 
-    //   pos = V.row(F(id,0))*(1-u-v)+V.row(F(id,1))*u+V.row(F(id,2))*v;
-    float u,v; 
-    // parametric distance so that
-    //   pos = origin + t * dir
-    float t; 
+    /// primitive id
+    int id; 
+    /// geometry id (not used)
+    int gid; 
+    /// barycentric coordinates so that 
+    ///   pos = V.row(F(id,0))*(1-u-v)+V.row(F(id,1))*u+V.row(F(id,2))*v;
+    Scalar u,v; 
+    /// parametric distance so that
+    ///   pos = origin + t * dir
+    Scalar t; 
   };
 }
 #endif 

include/igl/IGL_ASSERT.h

@@ -0,0 +1,9 @@
+// https://stackoverflow.com/a/985807/148668
+#include <cassert>
+#ifndef IGL_ASSERT
+#ifdef NDEBUG
+#define IGL_ASSERT(x) do { (void)sizeof(x);} while (0)
+#else
+#define IGL_ASSERT(x) assert(x)
+#endif
+#endif

include/igl/IndexComparison.h

@@ -8,10 +8,8 @@
 #ifndef IGL_INDEXCOMPARISON_H
 #define IGL_INDEXCOMPARISON_H
 namespace igl{
-  // Comparison struct used by sort
-  // http://bytes.com/topic/c/answers/132045-sort-get-index
-
-  // For use with functions like std::sort
+  /// Comparison struct used by sort
+  /// http://bytes.com/topic/c/answers/132045-sort-get-index
   template<class T> struct IndexLessThan
   {
     IndexLessThan(const T arr) : arr(arr) {}
@@ -22,7 +20,7 @@ namespace igl{
     const T arr;
   };
 
-  // For use with functions like std::unique
+  /// Comparison struct used by unique
   template<class T> struct IndexEquals
   {
     IndexEquals(const T arr) : arr(arr) {}
@@ -33,7 +31,7 @@ namespace igl{
     const T arr;
   };
 
-  // For use with functions like std::sort
+  /// Comparison struct for vectors for use with functions like std::sort
   template<class T> struct IndexVectorLessThan
   {
     IndexVectorLessThan(const T & vec) : vec ( vec) {}
@@ -44,7 +42,7 @@ namespace igl{
     const T & vec;
   };
 
-  // For use with functions like std::sort
+  /// Comparison struct for use with functions like std::sort
   template<class T> struct IndexDimLessThan
   {
     IndexDimLessThan(const T & mat,const int & dim, const int & j) : 
@@ -67,7 +65,7 @@ namespace igl{
     const int & j;
   };
 
-  // For use with functions like std::sort
+  /// Comparison struct For use with functions like std::sort
   template<class T> struct IndexRowLessThan
   {
     IndexRowLessThan(const T & mat) : mat ( mat) {}
@@ -91,7 +89,7 @@ namespace igl{
     const T & mat;
   };
 
-  // For use with functions like std::sort
+  /// Comparison struct for use with functions like std::sort
   template<class T> struct IndexRowEquals
   {
     IndexRowEquals(const T & mat) : mat ( mat) {}

include/igl/LinSpaced.h

@@ -1,33 +1,42 @@
 #ifndef IGL_LINSPACED_H
 #define IGL_LINSPACED_H
 #include <Eigen/Core>
-// This function is not intended to be a permanent function of libigl. Rather
-// it is a "drop-in" workaround for documented bug in Eigen:
-// http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1383
-//
-// Replace: 
-//
-//     Eigen::VectorXi::LinSpaced(size,low,high);
-//
-// With:
-//
-//     igl::LinSpaced<Eigen::VectorXi>(size,low,high);
-//
-// Specifcally, this version will _always_ return an empty vector if size==0,
-// regardless of the values for low and high. If size != 0, then this simply
-// returns the result of Eigen::Derived::LinSpaced.
-//
-// Until this bug is fixed, we should also avoid calls to the member function
-// `.setLinSpaced`. This means replacing:
-//
-//     a.setLinSpaced(size,low,high);
-//
-// with
-//
-//     a = igl::LinSpaced<decltype(a) >(size,low,high);
-//
+/// @file LinSpaced.h
+///
+/// This function is not intended to be a permanent function of libigl. Rather
+/// it is a "drop-in" workaround for documented bug in Eigen:
+/// http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1383
+///
+/// Replace: 
+///
+///     Eigen::VectorXi::LinSpaced(size,low,high);
+///
+/// With:
+///
+///     igl::LinSpaced<Eigen::VectorXi>(size,low,high);
+///
+/// Specifcally, this version will _always_ return an empty vector if size==0,
+/// regardless of the values for low and high. If size != 0, then this simply
+/// returns the result of Eigen::Derived::LinSpaced.
+///
+/// Until this bug is fixed, we should also avoid calls to the member function
+/// `.setLinSpaced`. This means replacing:
+///
+///     a.setLinSpaced(size,low,high);
+///
+/// with
+///
+///     a = igl::LinSpaced<decltype(a) >(size,low,high);
+///
 namespace igl
 {
+  /// Replacement for Eigen::DenseBase::LinSpaced
+  /// @param[in] size  number of elements
+  /// @param[in] low   first element
+  /// @param[in] high  last element
+  /// @return          vector of size elements linearly spaced between low and
+  ///
+  /// \fileinfo
   template <typename Derived>
   //inline typename Eigen::DenseBase< Derived >::RandomAccessLinSpacedReturnType 
   inline Derived LinSpaced(

include/igl/MappingEnergyType.h

@@ -9,10 +9,9 @@
 #define IGL_MAPPINGENERGYTYPE_H
 namespace igl
 {
-  // Energy Types used for Parameterization/Mapping. 
-  // Refer to SLIM [Rabinovich et al. 2017] for more details
+  /// Energy Types used for Parameterization/Mapping. 
+  /// Refer to SLIM [Rabinovich et al. 2017] for more details
   // Todo: Integrate with ARAPEnergyType
-  
   enum MappingEnergyType
   {
     ARAP = 0,

include/igl/MeshBooleanType.h

@@ -9,13 +9,20 @@
 #define IGL_MESH_BOOLEAN_TYPE_H
 namespace igl
 {
+  /// Boolean operation types
   enum MeshBooleanType
   {
+    /// A ∪ B
     MESH_BOOLEAN_TYPE_UNION = 0,
+    /// A ∩ B
     MESH_BOOLEAN_TYPE_INTERSECT = 1,
+    /// A \ B
     MESH_BOOLEAN_TYPE_MINUS = 2,
+    /// A ⊕ B
     MESH_BOOLEAN_TYPE_XOR = 3,
+    /// Resolve intersections without removing any non-coplanar faces
     MESH_BOOLEAN_TYPE_RESOLVE = 4,
+    /// Total number of Boolean options
     NUM_MESH_BOOLEAN_TYPES = 5
   };
 };

include/igl/MshLoader.cpp

@@ -114,10 +114,31 @@ IGL_INLINE igl::MshLoader::MshLoader(const std::string &filename) {
     fin.close();
 }
 
+IGL_INLINE int igl::MshLoader::node_dense_index(int node_tag) const {
+    const auto it = m_node_tag_to_dense.find(node_tag);
+    if (it == m_node_tag_to_dense.end()) {
+        std::stringstream err_msg;
+        err_msg << "Unknown node tag: " << node_tag;
+        throw std::runtime_error(err_msg.str());
+    }
+	return it->second;
+}
+
+IGL_INLINE int igl::MshLoader::element_dense_index(int elem_tag) const {
+	const auto it = m_element_tag_to_dense.find(elem_tag);
+    if (it == m_element_tag_to_dense.end()) {
+		std::stringstream err_msg;
+        err_msg << "Unknown element tag: " << elem_tag;
+		throw std::runtime_error(err_msg.str());
+    }
+	return it->second;
+}
+
 IGL_INLINE void igl::MshLoader::parse_nodes(std::ifstream& fin) {
     size_t num_nodes;
     fin >> num_nodes;
     m_nodes.resize(num_nodes*3);
+    m_node_tag_to_dense.clear();
 
     if (m_binary) {
 		size_t stride = (4+3*m_data_size);
@@ -127,23 +148,37 @@ IGL_INLINE void igl::MshLoader::parse_nodes(std::ifstream& fin) {
         fin.read(data, num_bytes);
 
         for (size_t i=0; i<num_nodes; i++) {
-            int node_idx;
-			memcpy(&node_idx, data+i*stride, sizeof(int));
-			node_idx-=1;
-			// directly move into vector storage
-			// this works only when m_data_size==sizeof(Float)==sizeof(double)
-			memcpy(&m_nodes[node_idx*3], data+i*stride + 4, m_data_size*3);
+            int node_tag;
+			memcpy(&node_tag, data+i*stride, sizeof(int));
+
+            if (node_tag <= 0) {
+                throw std::runtime_error("Invalid node tag");
+            }
+            if (m_node_tag_to_dense.find(node_tag) != m_node_tag_to_dense.end()) {
+                throw std::runtime_error("Duplicate node tag");
+            }
+            m_node_tag_to_dense[node_tag] = static_cast<int>(i);
+            // directly move into vector storage
+            // this works only when m_data_size==sizeof(Float)==sizeof(double)
+			memcpy(&m_nodes[i*3], data+i*stride + 4, m_data_size*3);
         }
         delete [] data;
     } else {
-        int node_idx;
+        int node_tag;
         for (size_t i=0; i<num_nodes; i++) {
-            fin >> node_idx;
-            node_idx -= 1;
+            fin >> node_tag;
+
+            if (node_tag <= 0) {
+				throw std::runtime_error("Invalid node tag");
+            }
+            if (m_node_tag_to_dense.find(node_tag) != m_node_tag_to_dense.end()) {
+                throw std::runtime_error("Duplicate node tag");
+            }
+            m_node_tag_to_dense[node_tag] = static_cast<int>(i);
             // here it's 3D node explicitly
-            fin >> m_nodes[node_idx*3]
-                >> m_nodes[node_idx*3+1]
-                >> m_nodes[node_idx*3+2];
+            fin >> m_nodes[i*3]
+                >> m_nodes[i*3+1]
+                >> m_nodes[i*3+2];
         }
     }
 }
@@ -152,6 +187,7 @@ IGL_INLINE void igl::MshLoader::parse_elements(std::ifstream& fin) {
     m_elements_tags.resize(2); //hardcoded to have 2 tags
     size_t num_elements;
     fin >> num_elements;
+    m_element_tag_to_dense.clear();
 
     size_t nodes_per_element;
 
@@ -168,15 +204,24 @@ IGL_INLINE void igl::MshLoader::parse_elements(std::ifstream& fin) {
 
             // store node info
             for (size_t i=0; i<num_elems; i++) {
-                int elem_idx;
+                int elem_tag;
 
                 // all elements in the segment share the same elem_type and number of nodes per element
                 m_elements_types.push_back(elem_type);
                 m_elements_lengths.push_back(nodes_per_element);
 
-                fin.read((char*)&elem_idx, sizeof(int));
-                elem_idx -= 1;
-                m_elements_ids.push_back(elem_idx);
+                fin.read((char*)&elem_tag, sizeof(int));
+
+                if (elem_tag <= 0) {
+                    throw std::runtime_error("Invalid element tag");
+                }
+                if (m_element_tag_to_dense.find(elem_tag) != m_element_tag_to_dense.end()) {
+                    throw std::runtime_error("Duplicate element tag");
+                }
+                m_element_tag_to_dense[elem_tag] = static_cast<int>(m_elements_ids.size());
+
+                elem_tag -= 1;
+                m_elements_ids.push_back(elem_tag);
 
                 // read first two tags
                 for (size_t j=0; j<num_tags; j++) {
@@ -191,10 +236,10 @@ IGL_INLINE void igl::MshLoader::parse_elements(std::ifstream& fin) {
                 m_elements_nodes_idx.push_back(m_elements.size());
                 // Element values.
                 for (size_t j=0; j<nodes_per_element; j++) {
-                    int idx;
-                    fin.read((char*)&idx, sizeof(int));
+                    int node_tag;
+                    fin.read((char*)&node_tag, sizeof(int));
                     
-                    m_elements.push_back(idx-1);
+                    m_elements.push_back(node_dense_index(node_tag));
                 }
             }
             elem_read += num_elems;
@@ -202,8 +247,16 @@ IGL_INLINE void igl::MshLoader::parse_elements(std::ifstream& fin) {
     } else {
         for (size_t i=0; i<num_elements; i++) {
             // Parse per element header
-            int elem_num, elem_type, num_tags;
-            fin >> elem_num >> elem_type >> num_tags;
+            int elem_tag, elem_type, num_tags;
+            fin >> elem_tag >> elem_type >> num_tags;
+
+            if (elem_tag <= 0) {
+                throw std::runtime_error("Invalid element tag");
+            }
+            if (m_element_tag_to_dense.find(elem_tag) != m_element_tag_to_dense.end()) {
+                throw std::runtime_error("Duplicate element tag");
+            }
+            m_element_tag_to_dense[elem_tag] = static_cast<int>(m_elements_ids.size());
 
             // read tags.
             for (size_t j=0; j<num_tags; j++) {
@@ -218,14 +271,14 @@ IGL_INLINE void igl::MshLoader::parse_elements(std::ifstream& fin) {
             m_elements_types.push_back(elem_type);
             m_elements_lengths.push_back(nodes_per_element);
 
-            elem_num -= 1;
-            m_elements_ids.push_back(elem_num);
+            elem_tag -= 1;
+            m_elements_ids.push_back(elem_tag);
             m_elements_nodes_idx.push_back(m_elements.size());
             // Parse node idx.
             for (size_t j=0; j<nodes_per_element; j++) {
-                int idx;
-                fin >> idx;
-                m_elements.push_back(idx-1); // msh index starts from 1.
+                int node_tag;
+                fin >> node_tag;
+                m_elements.push_back(node_dense_index(node_tag)); // msh index starts from 1.
             }
         }
     }
@@ -274,7 +327,7 @@ IGL_INLINE void igl::MshLoader::parse_node_field( std::ifstream& fin ) {
     int num_components    = int_tags[1];
     int num_entries       = int_tags[2];
 
-    std::vector<Float> field( num_entries*num_components );
+    std::vector<Float> field((m_nodes.size()/3)*num_components);
 
     if (m_binary) {
         size_t num_bytes = (num_components * m_data_size + 4) * num_entries;
@@ -282,23 +335,20 @@ IGL_INLINE void igl::MshLoader::parse_node_field( std::ifstream& fin ) {
         igl::_msh_eat_white_space(fin);
         fin.read(data, num_bytes);
         for (size_t i=0; i<num_entries; i++) {
-			int node_idx;
-			memcpy(&node_idx,&data[i*(4+num_components*m_data_size)],4);
-			
-            if(node_idx<1) throw std::runtime_error("Negative or zero index");
-            node_idx -= 1;
-			
-            if(node_idx>=num_entries) throw std::runtime_error("Index too big");
+			int node_tag;
+			memcpy(&node_tag,&data[i*(4+num_components*m_data_size)],4);
+
+            const int node_idx = node_dense_index(node_tag);
             size_t base_idx = i*(4+num_components*m_data_size) + 4;
             // TODO: make this work when m_data_size != sizeof(double) ?
 			memcpy(&field[node_idx*num_components], &data[base_idx], num_components*m_data_size);
         }
         delete [] data;
     } else {
-        int node_idx;
+        int node_tag;
         for (size_t i=0; i<num_entries; i++) {
-            fin >> node_idx;
-            node_idx -= 1;
+            fin >> node_tag;
+            const int node_idx = node_dense_index(node_tag);
             for (size_t j=0; j<num_components; j++) {
                 fin >> field[node_idx*num_components+j];
             }
@@ -346,7 +396,7 @@ IGL_INLINE void igl::MshLoader::parse_element_field(std::ifstream& fin) {
     std::string fieldname = str_tags[0];
     int num_components = int_tags[1];
     int num_entries = int_tags[2];
-    std::vector<Float> field(num_entries*num_components);
+    std::vector<Float> field(m_elements_ids.size()*num_components);
 
     if (m_binary) {
         size_t num_bytes = (num_components * m_data_size + 4) * num_entries;
@@ -354,20 +404,20 @@ IGL_INLINE void igl::MshLoader::parse_element_field(std::ifstream& fin) {
         igl::_msh_eat_white_space(fin);
         fin.read(data, num_bytes);
         for (int i=0; i<num_entries; i++) {
-			int elem_idx;
+			int elem_tag;
 			// works with sizeof(int)==4
-			memcpy(&elem_idx, &data[i*(4+num_components*m_data_size)],4);
-            elem_idx -= 1;
+			memcpy(&elem_tag, &data[i*(4+num_components*m_data_size)],4);
+            const int elem_idx = element_dense_index(elem_tag);
 			
 			// directly copy data into vector storage space
 			memcpy(&field[elem_idx*num_components], &data[i*(4+num_components*m_data_size) + 4], m_data_size*num_components);
         }
         delete [] data;
     } else {
-        int elem_idx;
+        int elem_tag;
         for (size_t i=0; i<num_entries; i++) {
-            fin >> elem_idx;
-            elem_idx -= 1;
+            fin >> elem_tag;
+            const int elem_idx = element_dense_index(elem_tag);
             for (size_t j=0; j<num_components; j++) {
                 fin >> field[elem_idx*num_components+j];
             }

include/igl/MshLoader.h

@@ -15,11 +15,12 @@
 #include <string>
 #include <vector>
 #include <algorithm>
+#include <unordered_map>
 
 namespace igl {
 
-// Class for loading information from .msh file
-// depends only on c++stl library
+/// Class for loading information from .msh file
+/// depends only on c++stl library
 class MshLoader {
     public:
 
@@ -60,6 +61,8 @@ class MshLoader {
               // other elements
               ELEMENT_POINT=15 };
     public:
+        /// Load a .msh file from a given path
+        /// @param[in] filename  path to .msh
         MshLoader(const std::string &filename);
 
     public:
@@ -153,11 +156,15 @@ class MshLoader {
         void parse_element_field(std::ifstream& fin);
         void parse_unknown_field(std::ifstream& fin,
                 const std::string& fieldname);
+        int node_dense_index(int node_tag) const;
+        int element_dense_index(int elem_tag) const;
 
     private:
         bool   m_binary;
         size_t m_data_size;
-        
+        std::unordered_map<int, int> m_node_tag_to_dense;
+        std::unordered_map<int, int> m_element_tag_to_dense;
+
         FloatVector m_nodes;    // len x 3 vector 
 
         IndexVector m_elements; // linear array for nodes corresponding to each element 
@@ -187,4 +194,4 @@ class MshLoader {
 #  include "MshLoader.cpp"
 #endif
 
-#endif //IGL_MSH_LOADER_H
+#endif //IGL_MSH_LOADER_H

include/igl/MshSaver.cpp

@@ -105,7 +105,6 @@ IGL_INLINE void igl::MshSaver::save_elements(const IndexVector& elements,
 
     if (m_num_elements > 0) {
         //int elem_type = el_type;
-        int num_elems = m_num_elements;
         //int tags = 0;
         if (!m_binary) {
             size_t el_ptr=0;
@@ -213,7 +212,6 @@ IGL_INLINE void igl::MshSaver::save_vector_field(const std::string& fieldname, c
     fout << "3" << std::endl; // 3-component vector field.
     fout << m_num_nodes << std::endl; // number of nodes
 
-    const Float zero = 0.0;
     if (m_binary) {
         for (size_t i=0; i<m_num_nodes; i++) {
             int node_idx = i+1;
@@ -275,7 +273,6 @@ IGL_INLINE void igl::MshSaver::save_elem_vector_field(const std::string& fieldna
     fout << "3" << std::endl; // 3-component vector field.
     fout << m_num_elements << std::endl; // number of elements
 
-    const Float zero = 0.0;
     if (m_binary) {
         for (size_t i=0; i<m_num_elements; ++i) {
             int elem_idx = i+1;
@@ -310,7 +307,6 @@ IGL_INLINE void igl::MshSaver::save_elem_tensor_field(const std::string& fieldna
     fout << "9" << std::endl; // 9-component tensor field.
     fout << m_num_elements << std::endl; // number of elements
 
-    const Float zero = 0.0;
     
     if (m_binary) {
         for (size_t i=0; i<m_num_elements; i++) {

include/igl/MshSaver.h

@@ -16,9 +16,9 @@
 
 namespace igl {
 
-// Class for dumping information to .msh file
-// depends only on c++stl library
-// current implementation works only with 3D information
+/// Class for dumping information to .msh file
+/// depends only on c++stl library
+/// current implementation works only with 3D information
 class MshSaver {
     public:
         typedef double Float;
@@ -30,6 +30,9 @@ class MshSaver {
         typedef std::vector<IntVector>   IntField;
         typedef std::vector<std::string> FieldNames;
 
+        /// Write a .msh to a given path
+        /// @param[in] filename  path to output file
+        /// @param[in] binary    whether to write in binary format
         MshSaver(const std::string& filename, bool binary=true);
         ~MshSaver();
 

include/igl/NormalType.h

@@ -10,14 +10,15 @@
 
 namespace igl
 {
-  // PER_VERTEX_NORMALS  Normals computed per vertex based on incident faces
-  // PER_FACE_NORMALS  Normals computed per face
-  // PER_CORNER_NORMALS  Normals computed per corner (aka wedge) based on
-  //   incident faces without sharp edge
+  /// Type of mesh normal computation method
   enum NormalType
   {
+    /// Normals computed per vertex based on incident faces
     PER_VERTEX_NORMALS,
+    /// Normals computed per face
     PER_FACE_NORMALS,
+    /// Normals computed per corner (aka wedge) based on incident faces without
+    /// sharp edge
     PER_CORNER_NORMALS
   };
 #  define NUM_NORMAL_TYPE 3

include/igl/ONE.h

@@ -9,9 +9,9 @@
 #define IGL_ONE_H
 namespace igl
 {
-  // Often one needs a reference to a dummy variable containing one as its
-  // value, for example when using AntTweakBar's
-  // TwSetParam( "3D View", "opened", TW_PARAM_INT32, 1, &INT_ONE);
+  /// Often one needs a reference to a dummy variable containing one as its
+  /// value, for example when using AntTweakBar's
+  /// TwSetParam( "3D View", "opened", TW_PARAM_INT32, 1, &INT_ONE);
   const char CHAR_ONE = 1;
   const int INT_ONE = 1;
   const unsigned int UNSIGNED_INT_ONE = 1;

include/igl/Orientation.h

@@ -0,0 +1,28 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Qingnan Zhou <qnzhou@gmail.com>
+// Copyright (C) 2025 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/.
+#pragma once
+#ifndef IGL_ORIENTATION_H
+#define IGL_ORIENTATION_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+  /// Types of orientations and other predicate results.
+  ///
+  /// \fileinfo
+  enum class Orientation {
+    POSITIVE=1, INSIDE=1,
+    NEGATIVE=-1, OUTSIDE=-1,
+    COLLINEAR=0, COPLANAR=0, COCIRCULAR=0, COSPHERICAL=0, DEGENERATE=0
+  };
+}
+
+
+#endif

include/igl/PI.h

@@ -11,8 +11,10 @@ namespace igl
 {
   // Use standard mathematical constants' M_PI if available
 #ifdef M_PI
+  /// π
   constexpr double PI = M_PI;
 #else
+  /// π
   constexpr double PI = 3.1415926535897932384626433832795;
 #endif
 }

include/igl/PlainMatrix.h

@@ -0,0 +1,95 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2024 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_PLAINMATRIX_H
+#define IGL_PLAINMATRIX_H
+#include <Eigen/Core>
+
+#include <type_traits>
+#include <Eigen/Dense>
+
+// Define void_t for compatibility if it's not in the standard library (C++11 and later)
+#if __cplusplus < 201703L
+namespace std {
+  template <typename... Ts>
+  using void_t = void;
+}
+#endif
+
+
+#ifndef IGL_DEFAULT_MAJORING
+#define IGL_DEFAULT_MAJORING Eigen::ColMajor
+#endif
+
+namespace igl
+{
+  template <typename Derived, int Rows, int Cols, int Options>
+  struct PlainMatrixHelper {
+    using Type = Eigen::Matrix<typename Derived::Scalar,Rows,Cols,((Rows == 1 && Cols != 1) ? Eigen::RowMajor : ((Cols == 1 && Rows != 1) ? Eigen::ColMajor : Options))>;
+  };
+  template <typename Derived, typename = void>
+  struct get_options {
+    static constexpr int value = IGL_DEFAULT_MAJORING;
+  };
+  
+  template <typename Derived>
+  struct get_options<Derived, std::void_t<decltype(Derived::Options)>> {
+    static constexpr int value = Derived::Options;
+  };
+  /// Some libigl implementations would (still do?) use a pattern like:
+  ///
+  /// template <typename DerivedA>
+  /// void foo(const Eigen::MatrixBase<DerivedA>& A)
+  /// {
+  ///    DerivedA B;
+  ///    igl::unique_rows(A,true,B);
+  /// }
+  ///
+  /// If `DerivedA` is `Eigen::Matrix`, then this may compile, but `DerivedA` might be
+  /// from a Eigen::Map or Eigen::Ref and fail to compile due to missing
+  /// construtor.
+  ///
+  /// Even worse, the code above will work if `DerivedA` has dynamic rows, but will
+  /// throw a runtime error if `DerivedA` has fixed number of rows.
+  ///
+  /// Instead it's better to declare `B` as a `Eigen::Matrix`
+  ///
+  /// Eigen::Matrix<typename DerivedA::Scalar,Eigen::Dynamic,DerivedA::ColsAtCompileTime,DerivedA::Options> B;
+  ///
+  /// Using `Eigen::Dynamic` for dimensions that may not be known at compile
+  /// time (or may be different from A).
+  ///
+  /// `igl::PlainMatrix` is just a helper to make this easier. So in this case
+  /// we could write:
+  ///
+  /// igl::PlainMatrix<DerivedA,Eigen::Dynamic> B;
+  ///
+  /// IIUC, if the code in question looks like:
+  ///
+  /// template <typename DerivedC>
+  /// void foo(Eigen::PlainObjectBase<DerivedC>& C)
+  /// {
+  ///    DerivedC B;
+  ///    …
+  ///    C.resize(not_known_at_compile_time,also_not_known_at_compile_time);
+  /// }
+  ///
+  /// Then it's probably fine. If C can be resized to different sizes, then
+  /// `DerivedC` should be `Eigen::Matrix`-like .
+  // Helper to check if `Options` exists in Derived
+
+  // Modify PlainMatrix to use get_options
+  template <typename Derived, 
+            int Rows = Derived::RowsAtCompileTime, 
+            int Cols = Derived::ColsAtCompileTime, 
+            int Options = get_options<Derived>::value>
+  using PlainMatrix = typename PlainMatrixHelper<Derived, Rows, Cols, Options>::Type;
+
+}
+
+#endif

include/igl/PlainVector.h

@@ -0,0 +1,36 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2024 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_PLAINVECTOR_H
+#define IGL_PLAINVECTOR_H
+#include <Eigen/Core>
+#include "PlainMatrix.h"
+
+namespace igl
+{
+  // PlainVectorHelper to determine correct matrix type based on Derived and Size
+  template <typename Derived, int Size, int Options>
+    struct PlainVectorHelper {
+      // Conditional Type: Column vector if is_column_vector is true, otherwise row vector
+      using Type = Eigen::Matrix<
+        typename Derived::Scalar,
+        (Derived::ColsAtCompileTime == 1 && Derived::RowsAtCompileTime != 1) ? Size : 1,
+        (Derived::ColsAtCompileTime == 1 && Derived::RowsAtCompileTime != 1) ? 1 : Size,
+        Options>;
+    };
+
+  /// \see PlainMatrix
+  template <
+    typename Derived, 
+    int Size = (Derived::ColsAtCompileTime == 1 && Derived::RowsAtCompileTime != 1) ? Derived::RowsAtCompileTime : Derived::ColsAtCompileTime,
+    int Options = get_options<Derived>::value>
+  using PlainVector = typename PlainVectorHelper<Derived, Size, Options>::Type;
+
+}
+#endif
+

include/igl/REDRUM.h

@@ -35,9 +35,17 @@
 
 #else
 
+/// Bold red colored text
+/// @param[in] X  text to color
+/// @returns colored text as "stream"
+/// #### Example:
+///
+/// \code{cpp}
+/// std::cout<<REDRUM("File "<<filename<<" not found.")<<std::endl;
+/// \endcode
+#define REDRUM(X)      "\e[1m\e[31m"<<X<<"\e[m"
 // Bold Red, etc.
 #define NORUM(X)       ""<<X<<""
-#define REDRUM(X)      "\e[1m\e[31m"<<X<<"\e[m"
 #define GREENRUM(X)    "\e[1m\e[32m"<<X<<"\e[m"
 #define YELLOWRUM(X)   "\e[1m\e[33m"<<X<<"\e[m"
 #define BLUERUM(X)     "\e[1m\e[34m"<<X<<"\e[m"

include/igl/STR.h

@@ -7,12 +7,23 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_STR_H
 #define IGL_STR_H
-// http://stackoverflow.com/a/2433143/148668
 #include <string>
 #include <sstream>
-// Suppose you have a function:
-//   void func(std::string c);
-// Then you can write:
-//   func(STR("foo"<<1<<"bar"));
+/// Convert a stream of things to std:;string
+///
+/// Suppose you have a function:
+/// \code{cpp}
+/// void func(std::string s);
+/// \endcode
+/// Then you can write:
+/// \code{cpp}
+///   func(C_STR("foo"<<1<<"bar"));
+/// \endcode
+/// which is equivalent to:
+/// \code{cpp}
+///   func("foo1bar");
+/// \endcode
+///
+// http://stackoverflow.com/a/2433143/148668
 #define STR(X) static_cast<std::ostringstream&>(std::ostringstream().flush() << X).str()
 #endif 

include/igl/Singular_Value_Decomposition_Main_Kernel_Body.hpp

@@ -18,6 +18,9 @@
 #pragma warning( disable : 592 )
 #endif
 
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wuninitialized"
+
 // #define USE_ACCURATE_RSQRT_IN_JACOBI_CONJUGATION
 // #define PERFORM_STRICT_QUATERNION_RENORMALIZATION
 
@@ -1272,6 +1275,8 @@
 #endif
 #endif
 
+#pragma clang diagnostic pop
+
 #ifdef __INTEL_COMPILER
 #pragma warning( default : 592 )
 #endif

include/igl/SolverStatus.h

@@ -9,14 +9,16 @@
 #define IGL_SOLVER_STATUS_H
 namespace igl
 {
+  /// Solver status type used by min_quad_with_fixed
   enum SolverStatus
   {
-    // Good
+    // Good. Solver declared convergence 
     SOLVER_STATUS_CONVERGED = 0,
-    // OK
+    // OK. Solver reached max iterations
     SOLVER_STATUS_MAX_ITER = 1,
-    // Bad
+    // Bad. Solver reported failure
     SOLVER_STATUS_ERROR = 2,
+    // Total number of solver types
     NUM_SOLVER_STATUSES = 3,
   };
 };

include/igl/SortableRow.h

@@ -14,16 +14,23 @@
 
 namespace igl
 {
-  // Templates:
-  //   T  should be a matrix that implements .size(), and operator(int i)
+  /// A row of things that can be sorted against other rows
+  /// @tparam T  should be a vector/matrix/array that implements .size(), and operator(int i)
   template <typename T>
   class SortableRow
   {
     public:
+      /// The data
       T data;
     public:
+      /// Default constructor
       SortableRow():data(){};
+      /// Constructor
+      /// @param[in] data  the data
       SortableRow(const T & data):data(data){};
+      /// Less than comparison
+      /// @param[in] that  the other row
+      /// @returns true if this row is less than that row
       bool operator<(const SortableRow & that) const
       {
         // Lexicographical
@@ -41,6 +48,9 @@ namespace igl
         // All characters the same, comes done to length
         return this->data.size()<that.data.size();
       };
+      /// Equality comparison
+      /// @param[in] that  the other row
+      /// @returns true if this row is equal to that row
       bool operator==(const SortableRow & that) const
       {
         if(this->data.size() != that.data.size())
@@ -56,6 +66,9 @@ namespace igl
         }
         return true;
       };
+      /// Inequality comparison
+      /// @param[in] that  the other row
+      /// @returns true if this row is not equal to that row
       bool operator!=(const SortableRow & that) const
       {
         return !(*this == that);

include/igl/SphereMeshWedge.cpp

@@ -0,0 +1,232 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2025 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 "SphereMeshWedge.h"
+#include "round_cone_signed_distance.h"
+#include "sign.h"
+#include <cassert>
+#include <Eigen/QR>
+#include <Eigen/Geometry>
+
+template <typename Scalar>
+IGL_INLINE igl::SphereMeshWedge<Scalar>::SphereMeshWedge(
+  const RowVector3S & V0,
+  const RowVector3S & V1,
+  const RowVector3S & V2,
+  const Scalar r0,
+  const Scalar r1,
+  const Scalar r2)
+{
+  // Internal copy
+  V.row(0) = V0;
+  V.row(1) = V1;
+  V.row(2) = V2;
+  r(0) = r0;
+  r(1) = r1;
+  r(2) = r2;
+
+  flavor = FULL;
+  // By default use full
+
+  EV.row(0) = V.row(2) - V.row(1);
+  EV.row(1) = V.row(0) - V.row(2);
+  EV.row(2) = V.row(1) - V.row(0);
+  l = EV.rowwise().norm();
+  l2 = l.array().square();
+  rr << r(1) - r(2), r(2) - r(0), r(0) - r(1);
+  a2 = l2.array() - rr.array().square();
+  il2 = 1.0/l2.array();
+
+  /////////////////////////////////////////////
+  /// BIG_VERTEX ?
+  /////////////////////////////////////////////
+  {
+    r.maxCoeff(&max_i);
+    int j = (max_i+1)%3;
+    int k = (max_i+2)%3;
+    if((l(k) + r(j) < r(max_i)) && (l(j) + r(k) < r(max_i)))
+    {
+      flavor = BIG_VERTEX;
+    }
+  }
+
+  /////////////////////////////////////////////
+  /// BIG_EDGE ?
+  /////////////////////////////////////////////
+  if(flavor == FULL)
+  {
+    // Case where one edge's roundCone containes the others
+    for(int e = 0;e<3;e++)
+    {
+      const int i = (e+1)%3;
+      const int j = (e+2)%3;
+      const int k = (e+3)%3;
+      const Scalar s = 
+        igl::round_cone_signed_distance(V.row(k),V.row(i),V.row(j),r(i),r(j));
+      if(-s > r(k))
+      {
+        flavor = BIG_EDGE;
+        max_i = i;
+        break;
+      }
+    }
+  }
+
+  if(flavor == FULL && !compute_planes())
+  {
+    flavor = NO_TRIANGLE;
+  }
+}
+
+template <typename Scalar>
+IGL_INLINE Scalar igl::SphereMeshWedge<Scalar>::operator()(const RowVector3S & p) const
+{
+  if(flavor == BIG_VERTEX)
+  {
+    // Case 0: Vertex i
+    return (p - V.row(max_i)).norm() - r(max_i);
+  }
+
+  if(flavor == BIG_EDGE)
+  {
+    const int i = max_i;
+    const int j = (i+1)%3;
+    // Case 1: Edge e
+    return this->round_cone_signed_distance(p,i,j);
+  }
+
+  Scalar s = std::numeric_limits<Scalar>::infinity();
+  if(flavor == FULL)
+  {
+    // This is possibly the bottleneck and could be turned into precomputed
+    // plane equations.
+
+    // signed distance to triangle plane (this is immediately recomputed later in
+    // sdSkewedExtrudedTriangle...)
+    const auto plane_sdf = [](
+        const RowVector3S & p,
+        const Eigen::RowVector4d & plane)
+    {
+      return plane.head<3>().dot(p) + plane(3);
+    };
+    Scalar d0 = plane_sdf(p, planes.row(0));
+    Scalar planes_s = -std::abs(d0);
+    // Reflect if necessary so that q is always on negative side of plane
+    RowVector3S q = p - (d0 - planes_s) * planes.row(0).template head<3>();
+    // Other planes (for negative side slab, by symmetry)
+    for(int i = 1;i<planes.rows();i++)
+    {
+      planes_s = std::max(planes_s,plane_sdf(q, planes.row(i)));
+    }
+    // This produces correct interior distance
+    if(planes_s <= 0)
+    {
+      s = std::min(s,planes_s);
+    }else
+    {
+
+      const auto & nor = planes.row(1).template head<3>();
+      const RowVector3S q0 = q - T.row(0);
+      const RowVector3S q1 = q - T.row(1);
+      const RowVector3S q2 = q - T.row(2);
+
+      if(!(sign(C.row(0).dot(q0)) + 
+           sign(C.row(1).dot(q1)) + 
+           sign(C.row(2).dot(q2))<2.0))
+      {
+        s = std::min(s,planes_s);
+      }
+    }
+
+    //s = std::min(s,sdSkewedExtrudedTriangle(q,V,T));
+    //s = std::min(s,sdSkewedExtrudedTriangle(p,B,V));
+  }
+  assert(flavor == FULL || flavor == NO_TRIANGLE);
+
+  for(int e = 0;e<3;e++)
+  {
+    const int i = (e+1)%3;
+    const int j = (e+2)%3;
+    s = std::min(s,this->round_cone_signed_distance(p,i,j));
+  }
+  return s;
+}
+
+template <typename Scalar>
+IGL_INLINE bool igl::SphereMeshWedge<Scalar>::compute_planes()
+{
+  // Non-degenerate case
+  const RowVector3S & a = V.row(0);
+  const RowVector3S & b = V.row(1);
+  const RowVector3S & c = V.row(2);
+  const Scalar & ra = r(0);
+  const Scalar & rb = r(1);
+  const Scalar & rc = r(2);
+  Eigen::Matrix<Scalar,2,3,Eigen::RowMajor> A;
+  A<<
+    b-a,
+    c-a;
+  const Eigen::Vector2d d(rb-ra,rc-ra);
+  const RowVector3S N = (A.row(0).cross(A.row(1))).normalized();
+  //const Eigen::CompleteOrthogonalDecomposition<decltype(A)> cod(A);
+  const RowVector3S n0 = A.completeOrthogonalDecomposition().solve(d);
+  const Scalar qA = N.squaredNorm();
+  // qB is zeros by construction. We could delete all terms involving qB
+  // It's not even clear if keeping them would lead to more accurate results.
+  const Scalar qB = 2 * N.dot(n0);
+  const Scalar qC = n0.squaredNorm() - 1;
+  const Scalar qD = qB*qB - 4*qA*qC;
+
+  if(qD<0) { return false; }
+
+  Scalar t_sol_1 = (-qB + std::sqrt(qD)) / (2*qA);
+  RowVector3S n1 = -(t_sol_1 * N + n0);
+  T  = V + r * n1;
+
+  const auto plane_equation = [](
+      const RowVector3S & a,
+      const RowVector3S & b,
+      const RowVector3S & c)->Eigen::RowVector4d
+  {
+    RowVector3S n = (b-a).cross(c-a).normalized();
+    n.normalize();
+    Scalar d = -n.dot(a);
+    return Eigen::RowVector4d(n(0),n(1),n(2),d);
+  };
+
+  planes.row(0) = plane_equation(V.row(0),V.row(1),V.row(2));
+  planes.row(1) = plane_equation(T.row(2),T.row(1),T.row(0));
+  planes.row(2) = plane_equation(V.row(1),V.row(0),T.row(0));
+  planes.row(3) = plane_equation(V.row(2),V.row(1),T.row(1));
+  planes.row(4) = plane_equation(V.row(0),V.row(2),T.row(2));
+
+  // Determine if the closest point is on the face.
+  const RowVector3S v10 = T.row(1) - T.row(0); 
+  const RowVector3S v21 = T.row(2) - T.row(1); 
+  const RowVector3S v02 = T.row(0) - T.row(2); 
+  const auto & nor = planes.row(1).template head<3>();
+  const RowVector3S c10 = v10.cross(nor);
+  const RowVector3S c21 = v21.cross(nor);
+  const RowVector3S c02 = v02.cross(nor);
+  C<<c10,c21,c02;
+  return true;
+}
+
+template <typename Scalar>
+IGL_INLINE Scalar igl::SphereMeshWedge<Scalar>::round_cone_signed_distance(const RowVector3S & p, const int i, const int j) const
+{
+  const int e = (j+1)%3;
+  return igl::round_cone_signed_distance(
+    p, V.row(i), r(i), r(j), EV.row(e), l2(e), rr(e), a2(e), il2(e));
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+/// Explicit template instantiation
+template class igl::SphereMeshWedge<double>;
+#endif
+

include/igl/SphereMeshWedge.h

@@ -0,0 +1,87 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2025 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_SPHERE_MESH_WEDGE_H
+#define IGL_SPHERE_MESH_WEDGE_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// A class to compute the signed distance to a "Sphere-Mesh Wedge" as seen in
+  /// variable radius offset surfaces or Sphere-Meshes. Each wedge is defined
+  /// by three vertices and three radii, one at each vertex. The wedge is
+  /// the union of all spheres at points on the triangle with radius linearly
+  /// interpolated. See, e.g., "Sphere-Meshes for Real-Time Hand Modeling and
+  /// Tracking" or "A Multilinear Model for Bidirectional Craniofacial
+  /// Reconstruction" or "Sphere-Meshes: Shape Approximation using Spherical
+  /// Quadric Error Metrics" or "Variable-Radius Offset Surface Approximation on
+  /// the GPU".
+  ///
+  template <typename Scalar>
+  class SphereMeshWedge
+  {
+    public: 
+      using RowVector3S = Eigen::Matrix<Scalar, 1, 3>;
+      // Fields
+      enum 
+      {
+        BIG_VERTEX = 0,
+        BIG_EDGE = 1,
+        NO_TRIANGLE = 2,
+        FULL = 3
+      } flavor;
+      Eigen::Matrix<Scalar,3,3,Eigen::RowMajor> V;
+      Eigen::Matrix<Scalar,3,1> r;
+      Eigen::Matrix<Scalar,3,3,Eigen::RowMajor> EV;
+      Eigen::Matrix<Scalar,3,1> l,l2,rr,a2,il2;
+      int max_i;
+      Eigen::Matrix<Scalar,5,4,Eigen::RowMajor> planes;
+      Eigen::Matrix<Scalar,3,3,Eigen::RowMajor> T;
+      Eigen::Matrix<Scalar,3,3,Eigen::RowMajor> C;
+
+      SphereMeshWedge(){}
+      /// Constructor that takes three vertices and three radii
+      ///
+      /// @param V0  first vertex position
+      /// @param V1  second vertex position
+      /// @param V2  third vertex position
+      /// @param r0  radius at first vertex
+      /// @param r1  radius at second vertex
+      /// @param r2  radius at third vertex
+      IGL_INLINE SphereMeshWedge(
+        const RowVector3S & V0,
+        const RowVector3S & V1,
+        const RowVector3S & V2,
+        const Scalar r0,
+        const Scalar r1,
+        const Scalar r2);
+      /// @param[in] p  3-vector query point
+      /// @return signed distance to the wedge at point p
+      IGL_INLINE Scalar operator()(const RowVector3S & p) const;
+    private:
+      /// Precompute planes used for determining bounded signed to the skewed
+      /// triangular slab portion.
+      ///
+      /// @return true if planes are well defined (false implies this slab has
+      /// no contribution).
+      IGL_INLINE bool compute_planes();
+      /// Compute the signed distance to the wedge at a point p for the edge
+      /// (i,j) 
+      ///
+      /// @param[in] p  3-vector query point
+      /// @param[in] i  index of first vertex (0,1,2)
+      /// @param[in] j  index of second vertex (0,1,2)
+      IGL_INLINE Scalar round_cone_signed_distance(const RowVector3S & p, const int i, const int j) const;
+  };
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "SphereMeshWedge.cpp"
+#endif
+#endif 

include/igl/Timer.h

@@ -22,13 +22,15 @@
 #include <sys/time.h>
 #endif
 #include <cstddef>
+#include <cstdint>
 
 namespace igl
 {
+  /// Simple timer class
   class Timer
   {
   public:
-    // default constructor
+    /// default constructor
     Timer():
       stopped(0),
 #ifdef WIN32
@@ -64,10 +66,13 @@ namespace igl
     }
 
 #ifdef __APPLE__
-    //Raw mach_absolute_times going in, difference in seconds out
-    double subtractTimes( uint64_t endTime, uint64_t startTime )
+    /// Raw mach_absolute_times going in, difference in seconds out
+    /// @param[in] endTime   end time
+    /// @param[in] startTime start time
+    /// @return time
+    double subtractTimes( std::uint64_t endTime, std::uint64_t startTime )
     {
-      uint64_t difference = endTime - startTime;
+      std::uint64_t difference = endTime - startTime;
       static double conversion = 0.0;
 
       if( conversion == 0.0 )
@@ -84,7 +89,7 @@ namespace igl
     }
 #endif
 
-    // start timer
+    /// start timer
     void   start()               
     {
       stopped = 0; // reset stop flag
@@ -98,7 +103,7 @@ namespace igl
 
     }
 
-    // stop the timer
+    /// stop the timer
     void   stop()                
     {
       stopped = 1; // set timer stopped flag
@@ -112,23 +117,27 @@ namespace igl
 #endif
 
     }
-    // get elapsed time in second
+    /// get elapsed time in second
+    /// @return time in seconds
     double getElapsedTime()      
     {
       return this->getElapsedTimeInSec();
     }
-    // get elapsed time in second (same as getElapsedTime)
+    /// get elapsed time in second (same as getElapsedTime)
+    /// @return time
     double getElapsedTimeInSec() 
     {
       return this->getElapsedTimeInMicroSec() * 0.000001;
     }
 
-    // get elapsed time in milli-second
+    /// get elapsed time in milli-second
+    /// @return time
     double getElapsedTimeInMilliSec()
     {
       return this->getElapsedTimeInMicroSec() * 0.001;
     }
-    // get elapsed time in micro-second
+    /// get elapsed time in micro-second
+    /// @return time
     double getElapsedTimeInMicroSec()          
     {
       double startTimeInMicroSec = 0;
@@ -167,8 +176,8 @@ namespace igl
     LARGE_INTEGER startCount;     
     LARGE_INTEGER endCount;       
 #elif __APPLE__
-    uint64_t startCount;           
-    uint64_t endCount;             
+    std::uint64_t startCount;           
+    std::uint64_t endCount;             
 #else
     timeval startCount;           
     timeval endCount;             

include/igl/Viewport.h

@@ -10,6 +10,7 @@
 
 namespace igl
 {
+  /// @private
   // Simple Viewport class for an opengl context. Handles reshaping and mouse.
   struct Viewport
   {

include/igl/WindingNumberAABB.h

@@ -13,19 +13,24 @@
 #ifndef IGL_WINDINGNUMBERAABB_H
 #define IGL_WINDINGNUMBERAABB_H
 #include "WindingNumberTree.h"
+#include "PlainMatrix.h"
 
 namespace igl
 {
+  /// Class for building an AABB tree to implement the divide and conquer
+  /// algorithm described in [Jacobson et al. 2013]. 
   template <
-    typename Point,
-    typename DerivedV, 
-    typename DerivedF >
-  class WindingNumberAABB : public WindingNumberTree<Point,DerivedV,DerivedF>
+    typename Scalar, 
+    typename Index>
+  class WindingNumberAABB : public WindingNumberTree<Scalar,Index>
   {
     protected:
+      // WindingNumberTree defines Point
+      using Point = typename WindingNumberTree<Scalar,Index>::Point;
+      using MatrixXF = typename WindingNumberTree<Scalar,Index>::MatrixXF;
       Point min_corner;
       Point max_corner;
-      typename DerivedV::Scalar total_positive_area;
+      Scalar total_positive_area;
     public: 
       enum SplitMethod
       {
@@ -35,16 +40,25 @@ namespace igl
       } split_method;
     public:
       inline WindingNumberAABB():
-        total_positive_area(std::numeric_limits<typename DerivedV::Scalar>::infinity()),
+        total_positive_area(std::numeric_limits<Scalar>::infinity()),
         split_method(MEDIAN_ON_LONGEST_AXIS)
       {}
+      /// Constructor
+      ///
+      /// @param[in] V  #V by 3 list of vertex positions
+      /// @param[in] F  #F by 3 list of triangle indices into V
+      template <typename DerivedV, typename DerivedF>
       inline WindingNumberAABB(
         const Eigen::MatrixBase<DerivedV> & V,
         const Eigen::MatrixBase<DerivedF> & F);
       inline WindingNumberAABB(
-        const WindingNumberTree<Point,DerivedV,DerivedF> & parent,
-        const Eigen::MatrixBase<DerivedF> & F);
-      // Initialize some things
+        const WindingNumberTree<Scalar,Index> & parent,
+        const typename WindingNumberTree<Scalar,Index>::MatrixXF & F);
+      /// Initialize the hierarchy to a given mesh
+      ///
+      /// @param[in] V  #V by 3 list of vertex positions
+      /// @param[in] F  #F by 3 list of triangle indices into V
+      template <typename DerivedV, typename DerivedF>
       inline void set_mesh(
         const Eigen::MatrixBase<DerivedV> & V,
         const Eigen::MatrixBase<DerivedF> & F);
@@ -53,8 +67,8 @@ namespace igl
       inline virtual void grow();
       // Compute min and max corners
       inline void compute_min_max_corners();
-      inline typename DerivedV::Scalar max_abs_winding_number(const Point & p) const;
-      inline typename DerivedV::Scalar max_simple_abs_winding_number(const Point & p) const;
+      inline Scalar max_abs_winding_number(const Point & p) const;
+      inline Scalar max_simple_abs_winding_number(const Point & p) const;
   };
 }
 
@@ -77,70 +91,73 @@ namespace igl
 #  define WindingNumberAABB_MIN_F 100
 #endif
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline void igl::WindingNumberAABB<Point,DerivedV,DerivedF>::set_mesh(
+template <typename Scalar, typename Index>
+  template <typename DerivedV, typename DerivedF>
+inline void igl::WindingNumberAABB<Scalar,Index>::set_mesh(
     const Eigen::MatrixBase<DerivedV> & V,
     const Eigen::MatrixBase<DerivedF> & F)
 {
-  igl::WindingNumberTree<Point,DerivedV,DerivedF>::set_mesh(V,F);
+  // static assert that DerivedF::ColsAtCompileTime == 3 or Eigen::Dynamic
+  static_assert(
+    DerivedF::ColsAtCompileTime == 3 || DerivedF::ColsAtCompileTime == Eigen::Dynamic,
+    "F should have 3 or Dynamic columns");
+  igl::WindingNumberTree<Scalar,Index>::set_mesh(V,F);
   init();
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline void igl::WindingNumberAABB<Point,DerivedV,DerivedF>::init()
+template <typename Scalar, typename Index>
+inline void igl::WindingNumberAABB<Scalar,Index>::init()
 {
-  using namespace Eigen;
   assert(max_corner.size() == 3);
   assert(min_corner.size() == 3);
   compute_min_max_corners();
-  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,1> dblA;
-  doublearea(this->getV(),this->getF(),dblA);
+  Eigen::Matrix<Scalar,Eigen::Dynamic,1> dblA;
+  doublearea((*this->Vptr),(this->F),dblA);
   total_positive_area = dblA.sum()/2.0;
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline igl::WindingNumberAABB<Point,DerivedV,DerivedF>::WindingNumberAABB(
+template <typename Scalar, typename Index>
+template <typename DerivedV, typename DerivedF>
+inline igl::WindingNumberAABB<Scalar,Index>::WindingNumberAABB(
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedF> & F):
-  WindingNumberTree<Point,DerivedV,DerivedF>(V,F),
+  WindingNumberTree<Scalar,Index>(V,F),
   min_corner(),
   max_corner(),
   total_positive_area(
-    std::numeric_limits<typename DerivedV::Scalar>::infinity()),
+    std::numeric_limits<Scalar>::infinity()),
   split_method(MEDIAN_ON_LONGEST_AXIS)
 {
   init();
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline igl::WindingNumberAABB<Point,DerivedV,DerivedF>::WindingNumberAABB(
-  const WindingNumberTree<Point,DerivedV,DerivedF> & parent,
-  const Eigen::MatrixBase<DerivedF> & F):
-  WindingNumberTree<Point,DerivedV,DerivedF>(parent,F),
+template <typename Scalar, typename Index>
+inline igl::WindingNumberAABB<Scalar,Index>::WindingNumberAABB(
+  const WindingNumberTree<Scalar,Index> & parent,
+  const typename WindingNumberTree<Scalar,Index>::MatrixXF & F):
+  WindingNumberTree<Scalar,Index>(parent,F),
   min_corner(),
   max_corner(),
   total_positive_area(
-    std::numeric_limits<typename DerivedV::Scalar>::infinity()),
+    std::numeric_limits<Scalar>::infinity()),
   split_method(MEDIAN_ON_LONGEST_AXIS)
 {
   init();
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline void igl::WindingNumberAABB<Point,DerivedV,DerivedF>::grow()
+template <typename Scalar, typename Index>
+inline void igl::WindingNumberAABB<Scalar,Index>::grow()
 {
-  using namespace std;
-  using namespace Eigen;
   // Clear anything that already exists
   this->delete_children();
 
-  //cout<<"cap.rows(): "<<this->getcap().rows()<<endl;
-  //cout<<"F.rows(): "<<this->getF().rows()<<endl;
+  //cout<<"cap.rows(): "<<(this->cap).rows()<<endl;
+  //cout<<"F.rows(): "<<(this->F).rows()<<endl;
 
   // Base cases
   if(
-    this->getF().rows() <= (WindingNumberAABB_MIN_F>0?WindingNumberAABB_MIN_F:0) ||
-    (this->getcap().rows() - 2) >= this->getF().rows())
+    (this->F).rows() <= (WindingNumberAABB_MIN_F>0?WindingNumberAABB_MIN_F:0) ||
+    ((this->cap).rows() - 2) >= (this->F).rows())
   {
     // Don't grow
     return;
@@ -148,8 +165,8 @@ inline void igl::WindingNumberAABB<Point,DerivedV,DerivedF>::grow()
 
   // Compute longest direction
   int max_d = -1;
-  typename DerivedV::Scalar max_len = 
-    -numeric_limits<typename DerivedV::Scalar>::infinity();
+  Scalar max_len = 
+    -std::numeric_limits<Scalar>::infinity();
   for(int d = 0;d<min_corner.size();d++)
   {
     if( (max_corner[d] - min_corner[d]) > max_len )
@@ -159,13 +176,13 @@ inline void igl::WindingNumberAABB<Point,DerivedV,DerivedF>::grow()
     }
   }
   // Compute facet barycenters
-  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,Eigen::Dynamic> BC;
-  barycenter(this->getV(),this->getF(),BC);
+  Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic> BC;
+  barycenter((*this->Vptr),(this->F),BC);
 
 
   // Blerg, why is selecting rows so difficult
 
-  typename DerivedV::Scalar split_value;
+  Scalar split_value;
   // Split in longest direction
   switch(split_method)
   {
@@ -182,8 +199,8 @@ inline void igl::WindingNumberAABB<Point,DerivedV,DerivedF>::grow()
   //cout<<"c: "<<0.5*(max_corner[max_d] + min_corner[max_d])<<" "<<
   //  "m: "<<split_value<<endl;;
 
-  vector<int> id( this->getF().rows());
-  for(int i = 0;i<this->getF().rows();i++)
+  std::vector<int> id( (this->F).rows());
+  for(int i = 0;i<(this->F).rows();i++)
   {
     if(BC(i,max_d) <= split_value)
     {
@@ -201,19 +218,19 @@ inline void igl::WindingNumberAABB<Point,DerivedV,DerivedF>::grow()
     // badly balanced base case (could try to recut)
     return;
   }
-  assert(lefts+rights == this->getF().rows());
-  DerivedF leftF(lefts,  this->getF().cols());
-  DerivedF rightF(rights,this->getF().cols());
+  assert(lefts+rights == (this->F).rows());
+  MatrixXF leftF(lefts,  (this->F).cols());
+  MatrixXF rightF(rights,(this->F).cols());
   int left_i = 0;
   int right_i = 0;
-  for(int i = 0;i<this->getF().rows();i++)
+  for(int i = 0;i<(this->F).rows();i++)
   {
     if(id[i] == 0)
     {
-      leftF.row(left_i++) = this->getF().row(i);
+      leftF.row(left_i++) = (this->F).row(i);
     }else if(id[i] == 1)
     {
-      rightF.row(right_i++) = this->getF().row(i);
+      rightF.row(right_i++) = (this->F).row(i);
     }else
     {
       assert(false);
@@ -222,18 +239,18 @@ inline void igl::WindingNumberAABB<Point,DerivedV,DerivedF>::grow()
   assert(right_i == rightF.rows());
   assert(left_i == leftF.rows());
   // Finally actually grow children and Recursively grow
-  WindingNumberAABB<Point,DerivedV,DerivedF> * leftWindingNumberAABB = 
-    new WindingNumberAABB<Point,DerivedV,DerivedF>(*this,leftF);
+  WindingNumberAABB<Scalar,Index> * leftWindingNumberAABB = 
+    new WindingNumberAABB<Scalar,Index>(*this,leftF);
   leftWindingNumberAABB->grow();
   this->children.push_back(leftWindingNumberAABB);
-  WindingNumberAABB<Point,DerivedV,DerivedF> * rightWindingNumberAABB = 
-    new WindingNumberAABB<Point,DerivedV,DerivedF>(*this,rightF);
+  WindingNumberAABB<Scalar,Index> * rightWindingNumberAABB = 
+    new WindingNumberAABB<Scalar,Index>(*this,rightF);
   rightWindingNumberAABB->grow();
   this->children.push_back(rightWindingNumberAABB);
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline bool igl::WindingNumberAABB<Point,DerivedV,DerivedF>::inside(const Point & p) const
+template <typename Scalar, typename Index>
+inline bool igl::WindingNumberAABB<Scalar,Index>::inside(const Point & p) const
 {
   assert(p.size() == max_corner.size());
   assert(p.size() == min_corner.size());
@@ -250,39 +267,38 @@ inline bool igl::WindingNumberAABB<Point,DerivedV,DerivedF>::inside(const Point
   return true;
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline void igl::WindingNumberAABB<Point,DerivedV,DerivedF>::compute_min_max_corners()
+template <typename Scalar, typename Index>
+inline void igl::WindingNumberAABB<Scalar,Index>::compute_min_max_corners()
 {
-  using namespace std;
   // initialize corners
   for(int d = 0;d<min_corner.size();d++)
   {
-    min_corner[d] =  numeric_limits<typename Point::Scalar>::infinity();
-    max_corner[d] = -numeric_limits<typename Point::Scalar>::infinity();
+    min_corner[d] =  std::numeric_limits<typename Point::Scalar>::infinity();
+    max_corner[d] = -std::numeric_limits<typename Point::Scalar>::infinity();
   }
 
   this->center = Point(0,0,0);
   // Loop over facets
-  for(int i = 0;i<this->getF().rows();i++)
+  for(int i = 0;i<(this->F).rows();i++)
   {
-    for(int j = 0;j<this->getF().cols();j++)
+    for(int j = 0;j<(this->F).cols();j++)
     {
       for(int d = 0;d<min_corner.size();d++)
       {
         min_corner[d] = 
-          this->getV()(this->getF()(i,j),d) < min_corner[d] ?  
-            this->getV()(this->getF()(i,j),d) : min_corner[d];
+          (*this->Vptr)((this->F)(i,j),d) < min_corner[d] ?  
+            (*this->Vptr)((this->F)(i,j),d) : min_corner[d];
         max_corner[d] = 
-          this->getV()(this->getF()(i,j),d) > max_corner[d] ?  
-            this->getV()(this->getF()(i,j),d) : max_corner[d];
+          (*this->Vptr)((this->F)(i,j),d) > max_corner[d] ?  
+            (*this->Vptr)((this->F)(i,j),d) : max_corner[d];
       }
       // This is biased toward vertices incident on more than one face, but
       // perhaps that's good
-      this->center += this->getV().row(this->getF()(i,j));
+      this->center += (*this->Vptr).row((this->F)(i,j));
     }
   }
   // Average
-  this->center.array() /= this->getF().size();
+  this->center.array() /= (this->F).size();
 
   //cout<<"min_corner: "<<this->min_corner.transpose()<<endl;
   //cout<<"Center: "<<this->center.transpose()<<endl;
@@ -293,32 +309,29 @@ inline void igl::WindingNumberAABB<Point,DerivedV,DerivedF>::compute_min_max_cor
   this->radius = (max_corner-min_corner).norm()/2.0;
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline typename DerivedV::Scalar
-igl::WindingNumberAABB<Point,DerivedV,DerivedF>::max_abs_winding_number(const Point & p) const
+template <typename Scalar, typename Index>
+inline Scalar
+igl::WindingNumberAABB<Scalar,Index>::max_abs_winding_number(const Point & p) const
 {
-  using namespace std;
   // Only valid if not inside
   if(inside(p))
   {
-    return numeric_limits<typename DerivedV::Scalar>::infinity();
+    return std::numeric_limits<Scalar>::infinity();
   }
   // Q: we know the total positive area so what's the most this could project
   // to? Remember it could be layered in the same direction.
-  return numeric_limits<typename DerivedV::Scalar>::infinity();
+  return std::numeric_limits<Scalar>::infinity();
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline typename DerivedV::Scalar 
-  igl::WindingNumberAABB<Point,DerivedV,DerivedF>::max_simple_abs_winding_number(
+template <typename Scalar, typename Index>
+inline Scalar 
+  igl::WindingNumberAABB<Scalar,Index>::max_simple_abs_winding_number(
   const Point & p) const
 {
-  using namespace std;
-  using namespace Eigen;
   // Only valid if not inside
   if(inside(p))
   {
-    return numeric_limits<typename DerivedV::Scalar>::infinity();
+    return std::numeric_limits<Scalar>::infinity();
   }
   // Max simple is the same as sum of positive winding number contributions of
   // bounding box
@@ -326,10 +339,10 @@ inline typename DerivedV::Scalar
   // begin precomputation
   //MatrixXd BV((int)pow(2,3),3);
   typedef
-    Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,Eigen::Dynamic>
+    Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>
     MatrixXS;
   typedef
-    Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,Eigen::Dynamic>
+    Eigen::Matrix<Index,Eigen::Dynamic,Eigen::Dynamic>
     MatrixXF;
   MatrixXS BV((int)(1<<3),3);
   BV <<
@@ -374,16 +387,4 @@ inline typename DerivedV::Scalar
   return igl::winding_number(BV,PBF,p);
 }
 
-// This is a bullshit template because AABB annoyingly needs templates for bad
-// combinations of 3D V with DIM=2 AABB
-//
-// _Define_ as a no-op rather than monkeying around with the proper code above
-namespace igl
-{
-  template <> inline igl::WindingNumberAABB<Eigen::Matrix<double, 1, 3, 1, 1, 3>,Eigen::Matrix<double, -1, 2, 0, -1, 2>,Eigen::Matrix<int, -1, 2, 0, -1, 2>>::WindingNumberAABB(const Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 0, -1, 2>> & V, const Eigen::MatrixBase<Eigen::Matrix<int, -1, 2, 0, -1, 2>> & F){};
-  template <> inline void igl::WindingNumberAABB<Eigen::Matrix<double, 1, 3, 1, 1, 3>,Eigen::Matrix<double, -1, 2, 0, -1, 2>,Eigen::Matrix<int, -1, 2, 0, -1, 2>>::grow(){};
-  template <> inline void igl::WindingNumberAABB<Eigen::Matrix<double, 1, 3, 1, 1, 3>,Eigen::Matrix<double, -1, 2, 0, -1, 2>,Eigen::Matrix<int, -1, 2, 0, -1, 2>>::init(){};
-
-}
-
 #endif

include/igl/WindingNumberAntipodalScene.h

@@ -0,0 +1,259 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2026 Philip Trettner <trettner@shapedcode.com>, Cedric Martens <cedric.martens@umontreal.ca>
+//
+// 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_WINDINGNUMBERANTIPODALSCENE_H
+#define IGL_WINDINGNUMBERANTIPODALSCENE_H
+
+#include "PI.h"
+#include "parallel_for.h"
+
+#include <Eigen/Core>
+
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+#include <unordered_map>
+#include <vector>
+
+namespace igl
+{
+  /// Precomputed scene for the Antipodal Method generalized winding number
+  ///
+  /// The scene stores only the weighted open-boundary edges of the input
+  /// triangle mesh and a fixed antipodal reference direction `x0`. Closed
+  /// (manifold) meshes have an empty boundary and the fractional term is
+  /// exactly zero.
+  ///
+  /// Querying the winding number additionally requires an `Intersector` that
+  /// returns the signed ray-mesh crossing count along `(p, x0)` over the
+  /// original mesh; the scene itself is intersector-agnostic. See
+  /// `igl::embree::EmbreeIntersector::signedIntersectionsRay` for an
+  /// optimized concrete implementation.
+  ///
+  /// ### Intersector concept
+  /// A type `I` satisfies the concept when it exposes:
+  ///   - `using OriginType    = ...;`     (3D row vector type)
+  ///   - `using DirectionType = ...;`     (3D row vector type)
+  ///   - `int signedIntersectionsRay(
+  ///         OriginType origin, DirectionType direction,
+  ///         /* defaulted tnear, tfar, mask */) const;`
+  ///
+  /// `winding_number` casts query point/direction to the intersector's types
+  /// at the call site, so a `double`-precision scene against a float-only
+  /// intersector works without an adaptor.
+  template <typename Scalar>
+  class WindingNumberAntipodalScene
+  {
+  public:
+    using Point     = Eigen::Matrix<Scalar, 1, 3>;
+    using Direction = Eigen::Matrix<Scalar, 1, 3>;
+
+  private:
+    struct WeightedSeg
+    {
+      Point a;
+      Point b;
+      Scalar w;
+    };
+
+  public:
+    /// Build a scene from a 3D triangle mesh.
+    ///
+    /// @param[in] V   #V by 3 list of vertex positions
+    /// @param[in] F   #F by 3 list of triangle indices
+    /// @param[in] x0  unit reference direction (defaults to a fixed non-axis-aligned vector)
+    template <typename DerivedV, typename DerivedF>
+    WindingNumberAntipodalScene(
+      const Eigen::MatrixBase<DerivedV> & V,
+      const Eigen::MatrixBase<DerivedF> & F,
+      const Direction & x0 = default_x0())
+      : m_x0(x0), m_face_count(static_cast<size_t>(F.rows()))
+    {
+      assert(V.cols() == 3 && "WindingNumberAntipodalScene: only 3D vertex positions are supported");
+      assert(F.cols() == 3 && "WindingNumberAntipodalScene: only triangle meshes are supported");
+      build_boundary_segments(V, F, m_boundary);
+    }
+
+    /// Single-point query: full generalized winding number at `p`
+    /// (fractional + signed integer crossings).
+    template <typename Intersector, typename Derivedp>
+    Scalar winding_number(
+      const Intersector & intersector,
+      const Eigen::MatrixBase<Derivedp> & p) const
+    {
+      const Point pp(static_cast<Scalar>(p(0)),
+                     static_cast<Scalar>(p(1)),
+                     static_cast<Scalar>(p(2)));
+      const Direction x1 = -m_x0;
+      Scalar area = Scalar(0);
+      for (const auto & ws : m_boundary)
+      {
+        const Point v0 = ws.a - pp;
+        const Point v1 = ws.b - pp;
+        area += ws.w * half_solid_angle_unorm(x1, v0, v1);
+      }
+      const Scalar frac = area / (Scalar(2) * Scalar(igl::PI));
+
+      using IO = typename Intersector::OriginType;
+      using ID = typename Intersector::DirectionType;
+      const IO io(static_cast<typename IO::Scalar>(pp(0)),
+                  static_cast<typename IO::Scalar>(pp(1)),
+                  static_cast<typename IO::Scalar>(pp(2)));
+      const ID id(static_cast<typename ID::Scalar>(m_x0(0)),
+                  static_cast<typename ID::Scalar>(m_x0(1)),
+                  static_cast<typename ID::Scalar>(m_x0(2)));
+      const int c = intersector.signedIntersectionsRay(io, id);
+
+      return frac + Scalar(c);
+    }
+
+    /// Batch query, parallelized via `igl::parallel_for`.
+    ///
+    /// @param[in]  intersector  Concept-compatible intersector built over the
+    ///                          same mesh used to construct the scene. Must
+    ///                          be safe to query concurrently.
+    /// @param[in]  O            #O by 3 list of query points
+    /// @param[out] W            #O by 1 list of winding numbers
+    template <typename Intersector, typename DerivedO, typename DerivedW>
+    void winding_number(
+      const Intersector & intersector,
+      const Eigen::MatrixBase<DerivedO> & O,
+      Eigen::PlainObjectBase<DerivedW> & W) const
+    {
+      W.resize(O.rows(), 1);
+
+      // Adaptive parallel-for threshold.
+      //
+      // The libigl thread pool has a roughly fixed ~1 ms TOTAL overhead per
+      // parallel_for invocation (not per iteration). So we only spawn the
+      // pool when the WHOLE batch is expected to take ≥ 1 ms.
+      //
+      // Per-query work heuristic:
+      //   t_q ≈ 50 * B + 100 * sqrt(F)  ns
+      // with B = boundary segment count, F = triangle count. Pool wins once
+      //   t_q * O > 10^6 ns  ⇒  O > 10^6 / t_q
+      // which is exactly parallel_for's `min_parallel` semantics.
+      //
+      // (This is a rough heuristic and should be revisited once parallel_for becomes lower-overhead)
+      const double t_q_ns =
+          50.0  * static_cast<double>(m_boundary.size()) +
+          100.0 * std::sqrt(static_cast<double>(m_face_count));
+      const size_t min_parallel = static_cast<size_t>(
+          std::ceil(1.0e6 / std::max(t_q_ns, 1.0)));
+
+      igl::parallel_for(O.rows(), [&](const int o)
+      {
+        W(o) = winding_number(intersector, O.row(o));
+      }, min_parallel);
+    }
+
+    /// Reference direction `x0` used to evaluate this scene.
+    const Direction & x0() const { return m_x0; }
+    /// Number of weighted boundary edge segments.
+    size_t num_boundary_segments() const { return m_boundary.size(); }
+    /// Number of triangles in the original mesh.
+    size_t num_faces() const { return m_face_count; }
+
+    /// Default reference direction: normalize(1, sqrt(2), sqrt(3)). A fixed
+    /// non-axis-aligned unit vector; any unit vector works per the paper;
+    /// this choice avoids accidental alignment with axis-aligned geometry.
+    static Direction default_x0()
+    {
+      Direction d(Scalar(1),
+                  Scalar(std::sqrt(2.0)),
+                  Scalar(std::sqrt(3.0)));
+      return d / d.norm();
+    }
+
+    /// Half the signed solid angle subtended by the spherical triangle
+    /// (x1, v0, v1) at the origin, via the unnormalized
+    /// Van Oosterom-Strackee formula. `x1` is expected to be a unit vector
+    /// (the antipodal "south pole" `-x0`); `v0`, `v1` need not be normalized.
+    /// Callers accumulate per-edge contributions and divide by 2π.
+    static Scalar half_solid_angle_unorm(
+      const Direction & x1, const Point & v0, const Point & v1)
+    {
+      const Scalar l0 = v0.norm();
+      const Scalar l1 = v1.norm();
+      const Scalar num   = x1.dot(v0.cross(v1));
+      const Scalar denom = l0 * l1
+                         + l1 * x1.dot(v0)
+                         + l0 * x1.dot(v1)
+                         + v0.dot(v1);
+      return std::atan2(num, denom);
+    }
+
+    /// Extract the open boundary as oriented, weighted edge segments.
+    ///
+    /// Each undirected edge {min, max} accumulates +1 for every triangle
+    /// that traverses it as min→max and -1 for every traversal max→min.
+    /// Interior edges of an oriented manifold cancel to zero and are
+    /// dropped. Surviving edges are emitted with positive weight; the
+    /// segment direction is flipped when the net count is negative so the
+    /// weight is always > 0. Non-manifold (≥3 incident triangles per edge)
+    /// is handled the same way. The surviving net count becomes the weight.
+    template <typename DerivedV, typename DerivedF>
+    static void build_boundary_segments(
+      const Eigen::MatrixBase<DerivedV> & V,
+      const Eigen::MatrixBase<DerivedF> & F,
+      std::vector<WeightedSeg> & out)
+    {
+      out.clear();
+      std::unordered_map<std::uint64_t, int> counts;
+      counts.reserve(static_cast<size_t>(F.rows()) * 3);
+
+      const auto pack = [](int lo, int hi) -> std::uint64_t
+      {
+        return (static_cast<std::uint64_t>(static_cast<std::uint32_t>(lo)) << 32)
+             |  static_cast<std::uint64_t>(static_cast<std::uint32_t>(hi));
+      };
+
+      for (Eigen::Index t = 0; t < F.rows(); ++t)
+      {
+        const int tri[3] = {
+          static_cast<int>(F(t, 0)),
+          static_cast<int>(F(t, 1)),
+          static_cast<int>(F(t, 2))
+        };
+        for (int e = 0; e < 3; ++e)
+        {
+          const int u = tri[e];
+          const int v = tri[(e + 1) % 3];
+          const int lo = u < v ? u : v;
+          const int hi = u < v ? v : u;
+          const std::uint64_t k = pack(lo, hi);
+          counts[k] += (u < v) ? +1 : -1;
+        }
+      }
+
+      out.reserve(counts.size());
+      for (const auto & kv : counts)
+      {
+        const int sum = kv.second;
+        if (sum == 0) continue;
+        const std::uint32_t lo = static_cast<std::uint32_t>(kv.first >> 32);
+        const std::uint32_t hi = static_cast<std::uint32_t>(kv.first & 0xFFFFFFFFu);
+        const Point a(static_cast<Scalar>(V(lo, 0)),
+                      static_cast<Scalar>(V(lo, 1)),
+                      static_cast<Scalar>(V(lo, 2)));
+        const Point b(static_cast<Scalar>(V(hi, 0)),
+                      static_cast<Scalar>(V(hi, 1)),
+                      static_cast<Scalar>(V(hi, 2)));
+        WeightedSeg ws;
+        if (sum > 0) { ws.a = a; ws.b = b; ws.w = static_cast<Scalar>(sum); }
+        else         { ws.a = b; ws.b = a; ws.w = static_cast<Scalar>(-sum); }
+        out.push_back(ws);
+      }
+    }
+
+    std::vector<WeightedSeg> m_boundary;
+    Direction m_x0;
+    size_t m_face_count = 0;
+  };
+}
+
+#endif

include/igl/WindingNumberMethod.h

@@ -9,14 +9,15 @@
 #define IGL_WINDINGNUMBERMETHOD_H
 namespace igl
 {
-  // EXACT_WINDING_NUMBER_METHOD  exact hierarchical evaluation
-  // APPROX_SIMPLE_WINDING_NUMBER_METHOD  poor approximation
-  // APPROX_CACHE_WINDING_NUMBER_METHOD  another poor approximation
   enum WindingNumberMethod
   {
+    // exact hierarchical evaluation
     EXACT_WINDING_NUMBER_METHOD = 0,
+    // poor approximation
     APPROX_SIMPLE_WINDING_NUMBER_METHOD = 1,
+    // another poor approximation
     APPROX_CACHE_WINDING_NUMBER_METHOD = 2,
+    /// Number of winding number methods
     NUM_WINDING_NUMBER_METHODS = 3
   };
 }

include/igl/WindingNumberTree.h

@@ -11,74 +11,68 @@
 #include <map>
 #include <Eigen/Dense>
 #include "WindingNumberMethod.h"
+#include <cassert>
+#include <memory>
 
 namespace igl
 {
-  // Space partitioning tree for computing winding number hierarchically.
-  //
-  // Templates:
-  //   Point  type for points in space, e.g. Eigen::Vector3d
+  /// Space partitioning tree for computing winding number hierarchically.
   template <
-    typename Point,
-    typename DerivedV, 
-    typename DerivedF >
+    typename Scalar, 
+    typename Index>
   class WindingNumberTree
   {
     public:
+      using Point = Eigen::Matrix<Scalar,1,3>;
       // Method to use (see enum above)
       //static double min_max_w;
       static std::map< 
         std::pair<const WindingNumberTree*,const WindingNumberTree*>, 
-        typename DerivedV::Scalar>
+        Scalar>
           cached;
-      // This is only need to fill in references, it should never actually be touched
-      // and shouldn't cause race conditions. (This is a hack, but I think it's "safe")
-      static DerivedV dummyV;
     protected:
       WindingNumberMethod method;
       const WindingNumberTree * parent;
       std::list<WindingNumberTree * > children;
       typedef 
-        Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,Eigen::Dynamic>
+        Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>
         MatrixXS;
       typedef 
-        Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,Eigen::Dynamic>
+        Eigen::Matrix<Index,Eigen::Dynamic,Eigen::Dynamic>
         MatrixXF;
-      //// List of boundary edges (recall edges are vertices in 2d)
-      //const Eigen::MatrixXi boundary;
-      // Base mesh vertices
-      DerivedV & V;
-      // Base mesh vertices with duplicates removed
+      // Base mesh vertices with duplicates removed (root will fill this in and
+      // then everyone's Vptr will point to it.
       MatrixXS SV;
+      // Shared pointer to base mesh vertices
+      std::shared_ptr<MatrixXS> Vptr;
       // Facets in this bounding volume
       MatrixXF F;
       // Tessellated boundary curve
       MatrixXF cap;
       // Upper Bound on radius of enclosing ball
-      typename DerivedV::Scalar radius;
+      Scalar radius;
       // (Approximate) center (of mass)
       Point center;
     public:
       inline WindingNumberTree();
       // For root
+      template <typename DerivedV, typename DerivedF>
       inline WindingNumberTree(
         const Eigen::MatrixBase<DerivedV> & V,
         const Eigen::MatrixBase<DerivedF> & F);
       // For chilluns 
       inline WindingNumberTree(
-        const WindingNumberTree<Point,DerivedV,DerivedF> & parent,
-        const Eigen::MatrixBase<DerivedF> & F);
+        const WindingNumberTree<Scalar,Index> & parent,
+        const typename igl::WindingNumberTree<Scalar,Index>::MatrixXF & F);
       inline virtual ~WindingNumberTree();
       inline void delete_children();
-      inline virtual void set_mesh(
+      template <typename DerivedV, typename DerivedF>
+      inline void set_mesh(
         const Eigen::MatrixBase<DerivedV> & V,
         const Eigen::MatrixBase<DerivedF> & F);
       // Set method
       inline void set_method( const WindingNumberMethod & m);
     public:
-      inline const DerivedV & getV() const;
-      inline const MatrixXF & getF() const;
-      inline const MatrixXF & getcap() const;
       // Grow the Tree recursively
       inline virtual void grow();
       // Determine whether a given point is inside the bounding 
@@ -93,12 +87,12 @@ namespace igl
       // Inputs:
       //   p  query point 
       // Returns winding number 
-      inline typename DerivedV::Scalar winding_number(const Point & p) const;
+      inline Scalar winding_number(const Point & p) const;
       // Same as above, but always computes winding number using exact method
       // (sum over every facet)
-      inline typename DerivedV::Scalar winding_number_all(const Point & p) const;
+      inline Scalar winding_number_all(const Point & p) const;
       // Same as above, but always computes using sum over tessllated boundary
-      inline typename DerivedV::Scalar winding_number_boundary(const Point & p) const;
+      inline Scalar winding_number_boundary(const Point & p) const;
       //// Same as winding_number above, but if max_simple_abs_winding_number is
       //// less than some threshold min_max_w just return 0 (colloquially the "fast
       //// multipole method)
@@ -121,10 +115,10 @@ namespace igl
       // Inputs:
       //   p  query point 
       // Returns max winding number of 
-      inline virtual typename DerivedV::Scalar max_abs_winding_number(const Point & p) const; 
+      inline virtual Scalar max_abs_winding_number(const Point & p) const; 
       // Same as above, but stronger assumptions on (V,F). Assumes (V,F) is a
       // simple polyhedron
-      inline virtual typename DerivedV::Scalar max_simple_abs_winding_number(const Point & p) const;
+      inline virtual Scalar max_simple_abs_winding_number(const Point & p) const;
       // Compute or read cached winding number for point p with respect to mesh
       // in bounding box, recursing according to approximation criteria
       //
@@ -132,7 +126,7 @@ namespace igl
       //   p  query point 
       //   that  WindingNumberTree containing mesh w.r.t. which we're computing w.n.
       // Returns cached winding number
-      inline virtual typename DerivedV::Scalar cached_winding_number(const WindingNumberTree & that, const Point & p) const;
+      inline virtual Scalar cached_winding_number(const WindingNumberTree & that, const Point & p) const;
   };
 }
 
@@ -143,89 +137,95 @@ namespace igl
 #include "triangle_fan.h"
 #include "exterior_edges.h"
 
-#include <igl/PI.h>
-#include <igl/remove_duplicate_vertices.h>
+#include "PI.h"
+#include "remove_duplicate_vertices.h"
 
 #include <iostream>
 #include <limits>
 
-//template <typename Point, typename DerivedV, typename DerivedF>
-//WindingNumberMethod WindingNumberTree<Point,DerivedV,DerivedF>::method = EXACT_WINDING_NUMBER_METHOD;
-//template <typename Point, typename DerivedV, typename DerivedF>
-//double WindingNumberTree<Point,DerivedV,DerivedF>::min_max_w = 0;
-template <typename Point, typename DerivedV, typename DerivedF>
-std::map< std::pair<const igl::WindingNumberTree<Point,DerivedV,DerivedF>*,const igl::WindingNumberTree<Point,DerivedV,DerivedF>*>, typename DerivedV::Scalar>
-  igl::WindingNumberTree<Point,DerivedV,DerivedF>::cached;
+//template <typename Scalar, typename Index>
+//WindingNumberMethod WindingNumberTree<Scalar,Index>::method = EXACT_WINDING_NUMBER_METHOD;
+//template <typename Scalar, typename Index>
+//double WindingNumberTree<Scalar,Index>::min_max_w = 0;
+template <typename Scalar, typename Index>
+std::map< std::pair<const igl::WindingNumberTree<Scalar,Index>*,const igl::WindingNumberTree<Scalar,Index>*>, Scalar>
+  igl::WindingNumberTree<Scalar,Index>::cached;
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline igl::WindingNumberTree<Point,DerivedV,DerivedF>::WindingNumberTree():
+template <typename Scalar, typename Index>
+inline igl::WindingNumberTree<Scalar,Index>::WindingNumberTree():
   method(EXACT_WINDING_NUMBER_METHOD),
   parent(NULL),
-  V(dummyV),
   SV(),
   F(),
   cap(),
-  radius(std::numeric_limits<typename DerivedV::Scalar>::infinity()),
+  radius(std::numeric_limits<Scalar>::infinity()),
   center(0,0,0)
 {
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline igl::WindingNumberTree<Point,DerivedV,DerivedF>::WindingNumberTree(
+template <typename Scalar, typename Index>
+template <typename DerivedV, typename DerivedF>
+inline igl::WindingNumberTree<Scalar,Index>::WindingNumberTree(
   const Eigen::MatrixBase<DerivedV> & _V,
   const Eigen::MatrixBase<DerivedF> & _F):
   method(EXACT_WINDING_NUMBER_METHOD),
   parent(NULL),
-  V(dummyV),
   SV(),
   F(),
   cap(),
-  radius(std::numeric_limits<typename DerivedV::Scalar>::infinity()),
+  radius(std::numeric_limits<Scalar>::infinity()),
   center(0,0,0)
 {
   set_mesh(_V,_F);
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline void igl::WindingNumberTree<Point,DerivedV,DerivedF>::set_mesh(
+template <typename Scalar, typename Index>
+template <typename DerivedV, typename DerivedF>
+inline void igl::WindingNumberTree<Scalar,Index>::set_mesh(
     const Eigen::MatrixBase<DerivedV> & _V,
     const Eigen::MatrixBase<DerivedF> & _F)
 {
-  using namespace std;
   // Remove any exactly duplicate vertices
   // Q: Can this ever increase the complexity of the boundary?
   // Q: Would we gain even more by remove almost exactly duplicate vertices?
-  MatrixXF SF,SVI,SVJ;
+  Eigen::Matrix<typename MatrixXF::Scalar,Eigen::Dynamic,1> SVI,SVJ;
   igl::remove_duplicate_vertices(_V,_F,0.0,SV,SVI,SVJ,F);
-  triangle_fan(igl::exterior_edges(F),cap);
-  V = SV;
+  {
+    Eigen::Matrix<typename MatrixXF::Scalar,Eigen::Dynamic,2> EE;
+    igl::exterior_edges(F,EE);
+    triangle_fan(EE,cap);
+  }
+  // point Vptr to SV
+  Vptr = std::make_shared<MatrixXS>(SV);
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline igl::WindingNumberTree<Point,DerivedV,DerivedF>::WindingNumberTree(
-  const igl::WindingNumberTree<Point,DerivedV,DerivedF> & parent,
-  const Eigen::MatrixBase<DerivedF> & _F):
+template <typename Scalar, typename Index>
+inline igl::WindingNumberTree<Scalar,Index>::WindingNumberTree(
+  const igl::WindingNumberTree<Scalar,Index> & parent,
+  const typename igl::WindingNumberTree<Scalar,Index>::MatrixXF & _F):
   method(parent.method),
   parent(&parent),
-  V(parent.V),
+  Vptr(parent.Vptr),
   SV(),
   F(_F),
-  cap(triangle_fan(igl::exterior_edges(_F)))
+  cap()
 {
+  Eigen::Matrix<typename MatrixXF::Scalar,Eigen::Dynamic,2> EE;
+  igl::exterior_edges(F,EE);
+  triangle_fan(EE,cap);
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline igl::WindingNumberTree<Point,DerivedV,DerivedF>::~WindingNumberTree()
+template <typename Scalar, typename Index>
+inline igl::WindingNumberTree<Scalar,Index>::~WindingNumberTree()
 {
   delete_children();
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline void igl::WindingNumberTree<Point,DerivedV,DerivedF>::delete_children()
+template <typename Scalar, typename Index>
+inline void igl::WindingNumberTree<Scalar,Index>::delete_children()
 {
-  using namespace std;
   // Delete children
-  typename list<WindingNumberTree<Point,DerivedV,DerivedF>* >::iterator cit = children.begin();
+  typename std::list<WindingNumberTree<Scalar,Index>* >::iterator cit = children.begin();
   while(cit != children.end())
   {
     // clear the memory of this item
@@ -235,8 +235,8 @@ inline void igl::WindingNumberTree<Point,DerivedV,DerivedF>::delete_children()
   }
 }
       
-template <typename Point, typename DerivedV, typename DerivedF>
-inline void igl::WindingNumberTree<Point,DerivedV,DerivedF>::set_method(const WindingNumberMethod & m)
+template <typename Scalar, typename Index>
+inline void igl::WindingNumberTree<Scalar,Index>::set_method(const WindingNumberMethod & m)
 {
   this->method = m;
   for(auto child : children)
@@ -245,44 +245,23 @@ inline void igl::WindingNumberTree<Point,DerivedV,DerivedF>::set_method(const Wi
   }
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline const DerivedV & igl::WindingNumberTree<Point,DerivedV,DerivedF>::getV() const
-{
-  return V;
-}
-
-template <typename Point, typename DerivedV, typename DerivedF>
-inline const typename igl::WindingNumberTree<Point,DerivedV,DerivedF>::MatrixXF& 
-  igl::WindingNumberTree<Point,DerivedV,DerivedF>::getF() const
-{
-  return F;
-}
-
-template <typename Point, typename DerivedV, typename DerivedF>
-inline const typename igl::WindingNumberTree<Point,DerivedV,DerivedF>::MatrixXF& 
-  igl::WindingNumberTree<Point,DerivedV,DerivedF>::getcap() const
-{
-  return cap;
-}
-
-template <typename Point, typename DerivedV, typename DerivedF>
-inline void igl::WindingNumberTree<Point,DerivedV,DerivedF>::grow()
+template <typename Scalar, typename Index>
+inline void igl::WindingNumberTree<Scalar,Index>::grow()
 {
   // Don't grow
   return;
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline bool igl::WindingNumberTree<Point,DerivedV,DerivedF>::inside(const Point & /*p*/) const
+template <typename Scalar, typename Index>
+inline bool igl::WindingNumberTree<Scalar,Index>::inside(const Point & /*p*/) const
 {
   return true;
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline typename DerivedV::Scalar 
-igl::WindingNumberTree<Point,DerivedV,DerivedF>::winding_number(const Point & p) const
+template <typename Scalar, typename Index>
+inline Scalar 
+igl::WindingNumberTree<Scalar,Index>::winding_number(const Point & p) const
 {
-  using namespace std;
   //cout<<"+"<<boundary.rows();
   // If inside then we need to be careful
   if(inside(p))
@@ -291,9 +270,9 @@ igl::WindingNumberTree<Point,DerivedV,DerivedF>::winding_number(const Point & p)
     if(children.size()>0)
     {
       // Recurse on each child and accumulate
-      typename DerivedV::Scalar sum = 0;
+      Scalar sum = 0;
       for(
-        typename list<WindingNumberTree<Point,DerivedV,DerivedF>* >::const_iterator cit = children.begin();
+        typename std::list<WindingNumberTree<Scalar,Index>* >::const_iterator cit = children.begin();
         cit != children.end();
         cit++)
       {
@@ -331,7 +310,7 @@ igl::WindingNumberTree<Point,DerivedV,DerivedF>::winding_number(const Point & p)
           return winding_number_boundary(p);
         case APPROX_SIMPLE_WINDING_NUMBER_METHOD:
         {
-          typename DerivedV::Scalar dist = (p-center).norm();
+          Scalar dist = (p-center).norm();
           // Radius is already an overestimate of inside
           if(dist>1.0*radius)
           {
@@ -356,24 +335,22 @@ igl::WindingNumberTree<Point,DerivedV,DerivedF>::winding_number(const Point & p)
   return 0;
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline typename DerivedV::Scalar 
-  igl::WindingNumberTree<Point,DerivedV,DerivedF>::winding_number_all(const Point & p) const
+template <typename Scalar, typename Index>
+inline Scalar 
+  igl::WindingNumberTree<Scalar,Index>::winding_number_all(const Point & p) const
 {
-  return igl::winding_number(V,F,p);
+  return igl::winding_number(*Vptr,F,p);
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline typename DerivedV::Scalar 
-igl::WindingNumberTree<Point,DerivedV,DerivedF>::winding_number_boundary(const Point & p) const
+template <typename Scalar, typename Index>
+inline Scalar 
+igl::WindingNumberTree<Scalar,Index>::winding_number_boundary(const Point & p) const
 {
-  using namespace Eigen;
-  using namespace std;
-  return igl::winding_number(V,cap,p);
+  return igl::winding_number(*Vptr,cap,p);
 }
 
-//template <typename Point, typename DerivedV, typename DerivedF>
-//inline double igl::WindingNumberTree<Point,DerivedV,DerivedF>::winding_number_approx_simple(
+//template <typename Scalar, typename Index>
+//inline double igl::WindingNumberTree<Scalar,Index>::winding_number_approx_simple(
 //  const Point & p, 
 //  const double min_max_w)
 //{
@@ -388,46 +365,43 @@ igl::WindingNumberTree<Point,DerivedV,DerivedF>::winding_number_boundary(const P
 //  }
 //}
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline void igl::WindingNumberTree<Point,DerivedV,DerivedF>::print(const char * tab)
+template <typename Scalar, typename Index>
+inline void igl::WindingNumberTree<Scalar,Index>::print(const char * tab)
 {
-  using namespace std;
   // Print all facets
-  cout<<tab<<"["<<endl<<F<<endl<<"]";
+  std::cout<<tab<<"["<<std::endl<<F<<std::endl<<"]";
   // Print children
   for(
-      typename list<WindingNumberTree<Point,DerivedV,DerivedF>* >::iterator cit = children.begin();
+      typename std::list<WindingNumberTree<Scalar,Index>* >::iterator cit = children.begin();
       cit != children.end();
       cit++)
   {
-    cout<<","<<endl;
-    (*cit)->print((string(tab)+"").c_str());
+    std::cout<<","<<std::endl;
+    (*cit)->print((std::string(tab)+"").c_str());
   }
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline typename DerivedV::Scalar 
-igl::WindingNumberTree<Point,DerivedV,DerivedF>::max_abs_winding_number(const Point & /*p*/) const
+template <typename Scalar, typename Index>
+inline Scalar 
+igl::WindingNumberTree<Scalar,Index>::max_abs_winding_number(const Point & /*p*/) const
 {
-  return std::numeric_limits<typename DerivedV::Scalar>::infinity();
+  return std::numeric_limits<Scalar>::infinity();
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline typename DerivedV::Scalar 
-igl::WindingNumberTree<Point,DerivedV,DerivedF>::max_simple_abs_winding_number(
+template <typename Scalar, typename Index>
+inline Scalar 
+igl::WindingNumberTree<Scalar,Index>::max_simple_abs_winding_number(
   const Point & /*p*/) const
 {
-  using namespace std;
-  return numeric_limits<typename DerivedV::Scalar>::infinity();
+  return std::numeric_limits<Scalar>::infinity();
 }
 
-template <typename Point, typename DerivedV, typename DerivedF>
-inline typename DerivedV::Scalar 
-igl::WindingNumberTree<Point,DerivedV,DerivedF>::cached_winding_number(
-  const igl::WindingNumberTree<Point,DerivedV,DerivedF> & that,
+template <typename Scalar, typename Index>
+inline Scalar 
+igl::WindingNumberTree<Scalar,Index>::cached_winding_number(
+  const igl::WindingNumberTree<Scalar,Index> & that,
   const Point & p) const
 {
-  using namespace std;
   // Simple metric for `is_far`
   //
   //   this             that
@@ -448,7 +422,7 @@ igl::WindingNumberTree<Point,DerivedV,DerivedF>::cached_winding_number(
   bool is_far = this->radius<that.radius;
   if(is_far)
   {
-    typename DerivedV::Scalar a = atan2(
+    Scalar a = atan2(
       that.radius - this->radius,
       (that.center - this->center).norm());
     assert(a>0);
@@ -458,7 +432,7 @@ igl::WindingNumberTree<Point,DerivedV,DerivedF>::cached_winding_number(
   if(is_far)
   {
     // Not implemented yet
-    pair<const WindingNumberTree*,const WindingNumberTree*> this_that(this,&that);
+    std::pair<const WindingNumberTree*,const WindingNumberTree*> this_that(this,&that);
     // Need to compute it for first time?
     if(cached.count(this_that)==0)
     {
@@ -473,7 +447,7 @@ igl::WindingNumberTree<Point,DerivedV,DerivedF>::cached_winding_number(
   }else
   {
     for(
-      typename list<WindingNumberTree<Point,DerivedV,DerivedF>* >::const_iterator cit = children.begin();
+      typename std::list<WindingNumberTree<Scalar,Index>* >::const_iterator cit = children.begin();
       cit != children.end();
       cit++)
     {
@@ -491,11 +465,4 @@ igl::WindingNumberTree<Point,DerivedV,DerivedF>::cached_winding_number(
   return 0;
 }
 
-// Explicit instantiation of static variable
-template <
-  typename Point,
-  typename DerivedV, 
-  typename DerivedF >
-DerivedV igl::WindingNumberTree<Point,DerivedV,DerivedF>::dummyV;
-
 #endif

include/igl/accumarray.h

@@ -11,14 +11,11 @@
 #include <Eigen/Core>
 namespace igl
 {
-  // ACCUMARRY Like Matlab's accumarray. Accumulate values in V using subscripts
-  // in S.
-  //
-  // Inputs:
-  //   S  #S list of subscripts
-  //   V  #V list of values
-  // Outputs:
-  //   A  max(subs)+1 list of accumulated values
+  /// Accumulate values in V using subscripts in S. Like Matlab's accumarray. 
+  ///
+  /// @param[in] S  #S list of subscripts
+  /// @param[in] V  #V list of values
+  /// @param[out] A  max(subs)+1 list of accumulated values
   template <
     typename DerivedS,
     typename DerivedV,
@@ -28,11 +25,11 @@ namespace igl
     const Eigen::MatrixBase<DerivedS> & S,
     const Eigen::MatrixBase<DerivedV> & V,
     Eigen::PlainObjectBase<DerivedA> & A);
-  // Inputs:
-  //   S  #S list of subscripts
-  //   V  single value used for all
-  // Outputs:
-  //   A  max(subs)+1 list of accumulated values
+  /// Accumulate constant value `V` using subscripts in S. Like Matlab's accumarray. 
+  ///
+  /// @param[in] S  #S list of subscripts
+  /// @param[in] V  single value used for all
+  /// @param[out] A  max(subs)+1 list of accumulated values
   template <
     typename DerivedS,
     typename DerivedA

include/igl/active_set.cpp

@@ -11,7 +11,10 @@
 #include "slice_into.h"
 #include "cat.h"
 //#include "matlab_format.h"
+#include "placeholders.h"
+#include "PlainMatrix.h"
 
+#include <cassert>
 #include <iostream>
 #include <limits>
 #include <algorithm>
@@ -31,25 +34,24 @@ template <
   >
 IGL_INLINE igl::SolverStatus igl::active_set(
   const Eigen::SparseMatrix<AT>& A,
-  const Eigen::PlainObjectBase<DerivedB> & B,
-  const Eigen::PlainObjectBase<Derivedknown> & known,
-  const Eigen::PlainObjectBase<DerivedY> & Y,
+  const Eigen::MatrixBase<DerivedB> & B,
+  const Eigen::MatrixBase<Derivedknown> & known,
+  const Eigen::MatrixBase<DerivedY> & Y,
   const Eigen::SparseMatrix<AeqT>& Aeq,
-  const Eigen::PlainObjectBase<DerivedBeq> & Beq,
+  const Eigen::MatrixBase<DerivedBeq> & Beq,
   const Eigen::SparseMatrix<AieqT>& Aieq,
-  const Eigen::PlainObjectBase<DerivedBieq> & Bieq,
-  const Eigen::PlainObjectBase<Derivedlx> & p_lx,
-  const Eigen::PlainObjectBase<Derivedux> & p_ux,
+  const Eigen::MatrixBase<DerivedBieq> & Bieq,
+  const Eigen::MatrixBase<Derivedlx> & p_lx,
+  const Eigen::MatrixBase<Derivedux> & p_ux,
   const igl::active_set_params & params,
   Eigen::PlainObjectBase<DerivedZ> & Z
   )
 {
+
 //#define ACTIVE_SET_CPP_DEBUG
 #if defined(ACTIVE_SET_CPP_DEBUG) && !defined(_MSC_VER)
 #  warning "ACTIVE_SET_CPP_DEBUG"
 #endif
-  using namespace Eigen;
-  using namespace std;
   SolverStatus ret = SOLVER_STATUS_ERROR;
   const int n = A.rows();
   assert(n == A.cols() && "A must be square");
@@ -72,7 +74,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   if(p_lx.size() == 0)
   {
     lx = Derivedlx::Constant(
-      n,1,-numeric_limits<typename Derivedlx::Scalar>::max());
+      n,1,-std::numeric_limits<typename Derivedlx::Scalar>::max());
   }else
   {
     lx = p_lx;
@@ -80,7 +82,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   if(p_ux.size() == 0)
   {
     ux = Derivedux::Constant(
-      n,1,numeric_limits<typename Derivedux::Scalar>::max());
+      n,1,std::numeric_limits<typename Derivedux::Scalar>::max());
   }else
   {
     ux = p_ux;
@@ -104,14 +106,12 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   typedef int BOOL;
 #define TRUE 1
 #define FALSE 0
-  Matrix<BOOL,Dynamic,1> as_lx = Matrix<BOOL,Dynamic,1>::Constant(n,1,FALSE);
-  Matrix<BOOL,Dynamic,1> as_ux = Matrix<BOOL,Dynamic,1>::Constant(n,1,FALSE);
-  Matrix<BOOL,Dynamic,1> as_ieq = Matrix<BOOL,Dynamic,1>::Constant(Aieq.rows(),1,FALSE);
+  Eigen::Matrix<BOOL,Eigen::Dynamic,1> as_lx = Eigen::Matrix<BOOL,Eigen::Dynamic,1>::Constant(n,1,FALSE);
+  Eigen::Matrix<BOOL,Eigen::Dynamic,1> as_ux = Eigen::Matrix<BOOL,Eigen::Dynamic,1>::Constant(n,1,FALSE);
+  Eigen::Matrix<BOOL,Eigen::Dynamic,1> as_ieq = Eigen::Matrix<BOOL,Eigen::Dynamic,1>::Constant(Aieq.rows(),1,FALSE);
 
   // Keep track of previous Z for comparison
-  DerivedZ old_Z;
-  old_Z = DerivedZ::Constant(
-      n,1,numeric_limits<typename DerivedZ::Scalar>::max());
+  PlainMatrix<DerivedZ> old_Z;
 
   int iter = 0;
   while(true)
@@ -121,35 +121,43 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     cout<<"  pre"<<endl;
 #endif
     // FIND BREACHES OF CONSTRAINTS
+#ifdef ACTIVE_SET_CPP_DEBUG
     int new_as_lx = 0;
     int new_as_ux = 0;
     int new_as_ieq = 0;
+#endif
     if(Z.size() > 0)
     {
       for(int z = 0;z < n;z++)
       {
         if(Z(z) < lx(z))
         {
+#ifdef ACTIVE_SET_CPP_DEBUG
           new_as_lx += (as_lx(z)?0:1);
+#endif
           //new_as_lx++;
           as_lx(z) = TRUE;
         }
         if(Z(z) > ux(z))
         {
+#ifdef ACTIVE_SET_CPP_DEBUG
           new_as_ux += (as_ux(z)?0:1);
+#endif
           //new_as_ux++;
           as_ux(z) = TRUE;
         }
       }
       if(Aieq.rows() > 0)
       {
-        DerivedZ AieqZ;
+        PlainMatrix<DerivedZ,Eigen::Dynamic> AieqZ;
         AieqZ = Aieq*Z;
         for(int a = 0;a<Aieq.rows();a++)
         {
           if(AieqZ(a) > Bieq(a))
           {
+#ifdef ACTIVE_SET_CPP_DEBUG
             new_as_ieq += (as_ieq(a)?0:1);
+#endif
             as_ieq(a) = TRUE;
           }
         }
@@ -158,14 +166,17 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       cout<<"  new_as_lx: "<<new_as_lx<<endl;
       cout<<"  new_as_ux: "<<new_as_ux<<endl;
 #endif
-      const double diff = (Z-old_Z).squaredNorm();
-#ifdef ACTIVE_SET_CPP_DEBUG
-      cout<<"diff: "<<diff<<endl;
-#endif
-      if(diff < params.solution_diff_threshold)
+      if(iter > 0)
       {
-        ret = SOLVER_STATUS_CONVERGED;
-        break;
+        const double diff = (Z-old_Z).squaredNorm();
+#ifdef ACTIVE_SET_CPP_DEBUG
+        cout<<"diff: "<<diff<<endl;
+#endif
+        if(diff < params.solution_diff_threshold)
+        {
+          ret = SOLVER_STATUS_CONVERGED;
+          break;
+        }
       }
       old_Z = Z;
     }
@@ -190,9 +201,9 @@ IGL_INLINE igl::SolverStatus igl::active_set(
 #endif
 
     // PREPARE FIXED VALUES
-    Derivedknown known_i;
+    Eigen::Matrix<typename Derivedknown::Scalar,Eigen::Dynamic,1> known_i;
     known_i.resize(nk + as_lx_count + as_ux_count,1);
-    DerivedY Y_i;
+    PlainMatrix<DerivedY,Eigen::Dynamic,1> Y_i;
     Y_i.resize(nk + as_lx_count + as_ux_count,1);
     {
       known_i.block(0,0,known.rows(),known.cols()) = known;
@@ -225,7 +236,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     // PREPARE EQUALITY CONSTRAINTS
     Eigen::Matrix<typename DerivedY::Scalar, Eigen::Dynamic, 1> as_ieq_list(as_ieq_count,1);
     // Gather active constraints and resp. rhss
-    DerivedBeq Beq_i;
+    PlainMatrix<DerivedBeq,Eigen::Dynamic,1> Beq_i;
     Beq_i.resize(Beq.rows()+as_ieq_count,1);
     Beq_i.head(Beq.rows()) = Beq;
     {
@@ -243,7 +254,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       assert(k == as_ieq_count);
     }
     // extract active constraint rows
-    SparseMatrix<AeqT> Aeq_i,Aieq_i;
+    Eigen::SparseMatrix<AeqT> Aeq_i,Aieq_i;
     slice(Aieq,as_ieq_list,1,Aieq_i);
     // Append to equality constraints
     cat(1,Aeq,Aieq_i,Aeq_i);
@@ -253,7 +264,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
 #ifndef NDEBUG
     {
       // NO DUPES!
-      Matrix<BOOL,Dynamic,1> fixed = Matrix<BOOL,Dynamic,1>::Constant(n,1,FALSE);
+      Eigen::Matrix<BOOL ,Eigen::Dynamic,1> fixed = Eigen::Matrix<BOOL ,Eigen::Dynamic,1>::Constant(n,1,FALSE);
       for(int k = 0;k<known_i.size();k++)
       {
         assert(!fixed[known_i(k)]);
@@ -262,7 +273,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     }
 #endif
 
-    DerivedZ sol;
+    PlainMatrix<DerivedZ,Eigen::Dynamic,Eigen::Dynamic> sol;
     if(known_i.size() == A.rows())
     {
       // Everything's fixed?
@@ -270,7 +281,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       cout<<"  everything's fixed."<<endl;
 #endif
       Z.resize(A.rows(),Y_i.cols());
-      slice_into(Y_i,known_i,1,Z);
+      Z(known_i,igl::placeholders::all) = Y_i;
       sol.resize(0,Y_i.cols());
       assert(Aeq_i.rows() == 0 && "All fixed but linearly constrained");
     }else
@@ -280,11 +291,15 @@ IGL_INLINE igl::SolverStatus igl::active_set(
 #endif
       if(!min_quad_with_fixed_precompute(A,known_i,Aeq_i,params.Auu_pd,data))
       {
+#ifdef ACTIVE_SET_CPP_DEBUG
         cerr<<"Error: min_quad_with_fixed precomputation failed."<<endl;
+#endif
         if(iter > 0 && Aeq_i.rows() > Aeq.rows())
         {
+#ifdef ACTIVE_SET_CPP_DEBUG
           cerr<<"  *Are you sure rows of [Aeq;Aieq] are linearly independent?*"<<
             endl;
+#endif
         }
         ret = SOLVER_STATUS_ERROR;
         break;
@@ -294,7 +309,9 @@ IGL_INLINE igl::SolverStatus igl::active_set(
 #endif
       if(!min_quad_with_fixed_solve(data,B,Y_i,Beq_i,Z,sol))
       {
+#ifdef ACTIVE_SET_CPP_DEBUG
         cerr<<"Error: min_quad_with_fixed solve failed."<<endl;
+#endif
         ret = SOLVER_STATUS_ERROR;
         break;
       }
@@ -308,18 +325,18 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     }
 
     // Compute Lagrange multiplier values for known_i
-    SparseMatrix<AT> Ak;
+    Eigen::SparseMatrix<AT> Ak;
     // Slow
     slice(A,known_i,1,Ak);
-    DerivedB Bk;
-    slice(B,known_i,Bk);
-    MatrixXd Lambda_known_i = -(0.5*Ak*Z + 0.5*Bk);
+    //slice(B,known_i,Bk);
+    PlainMatrix<DerivedB,Eigen::Dynamic> Bk = B(known_i,igl::placeholders::all);
+    Eigen::MatrixXd Lambda_known_i = -(0.5*Ak*Z + 0.5*Bk);
     // reverse the lambda values for lx
     Lambda_known_i.block(nk,0,as_lx_count,1) =
       (-1*Lambda_known_i.block(nk,0,as_lx_count,1)).eval();
 
     // Extract Lagrange multipliers for Aieq_i (always at back of sol)
-    VectorXd Lambda_Aieq_i(Aieq_i.rows(),1);
+    Eigen::VectorXd Lambda_Aieq_i(Aieq_i.rows(),1);
     for(int l = 0;l<Aieq_i.rows();l++)
     {
       Lambda_Aieq_i(Aieq_i.rows()-1-l) = sol(sol.rows()-1-l);
@@ -365,6 +382,6 @@ IGL_INLINE igl::SolverStatus igl::active_set(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template igl::SolverStatus igl::active_set<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 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> >(Eigen::SparseMatrix<double, 0, int> const&, 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> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, igl::active_set_params const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
-template igl::SolverStatus igl::active_set<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 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> >(Eigen::SparseMatrix<double, 0, int> const&, 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> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, igl::active_set_params const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template igl::SolverStatus igl::active_set<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 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> >(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&, Eigen::SparseMatrix<double, 0, int> 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&, igl::active_set_params const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template igl::SolverStatus igl::active_set<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, 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> >(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&, Eigen::SparseMatrix<double, 0, int> 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&, igl::active_set_params const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/active_set.h

@@ -16,39 +16,45 @@
 namespace igl
 {
   struct active_set_params;
-  // Known Bugs: rows of [Aeq;Aieq] **must** be linearly independent. Should be
-  // using QR decomposition otherwise:
-  // https://v8doc.sas.com/sashtml/ormp/chap5/sect32.htm
-  //
-  // ACTIVE_SET Minimize quadratic energy 
-  //
-  // 0.5*Z'*A*Z + Z'*B + C with constraints
-  //
-  // that Z(known) = Y, optionally also subject to the constraints Aeq*Z = Beq,
-  // and further optionally subject to the linear inequality constraints that
-  // Aieq*Z <= Bieq and constant inequality constraints lx <= x <= ux
-  //
-  // Inputs:
-  //   A  n by n matrix of quadratic coefficients
-  //   B  n by 1 column of linear coefficients
-  //   known  list of indices to known rows in Z
-  //   Y  list of fixed values corresponding to known rows in Z
-  //   Aeq  meq by n list of linear equality constraint coefficients
-  //   Beq  meq by 1 list of linear equality constraint constant values
-  //   Aieq  mieq by n list of linear inequality constraint coefficients
-  //   Bieq  mieq by 1 list of linear inequality constraint constant values
-  //   lx  n by 1 list of lower bounds [] implies -Inf
-  //   ux  n by 1 list of upper bounds [] implies Inf
-  //   params  struct of additional parameters (see below)
-  //   Z  if not empty, is taken to be an n by 1 list of initial guess values
-  //     (see output)
-  // Outputs:
-  //   Z  n by 1 list of solution values
-  // Returns true on success, false on error
-  //
-  // Benchmark: For a harmonic solve on a mesh with 325K facets, matlab 2.2
-  // secs, igl/min_quad_with_fixed.h 7.1 secs
-  //
+  ///
+  /// Minimize convex quadratic energy subject to linear inequality constraints
+  ///
+  ///     min ½ Zᵀ A Z + Zᵀ B + constant
+  ///      Z
+  ///     subject to
+  ///            Aeq Z = Beq
+  ///            Aieq Z <= Bieq
+  ///            lx <= Z <= ux
+  ///            Z(known) = Y
+  ///
+  /// that Z(known) = Y, optionally also subject to the constraints Aeq*Z = Beq,
+  /// and further optionally subject to the linear inequality constraints that
+  /// Aieq*Z <= Bieq and constant inequality constraints lx <= x <= ux
+  ///
+  /// @param[in] A  n by n matrix of quadratic coefficients
+  /// @param[in] B  n by 1 column of linear coefficients
+  /// @param[in] known  list of indices to known rows in Z
+  /// @param[in] Y  list of fixed values corresponding to known rows in Z
+  /// @param[in] Aeq  meq by n list of linear equality constraint coefficients
+  /// @param[in] Beq  meq by 1 list of linear equality constraint constant values
+  /// @param[in] Aieq  mieq by n list of linear inequality constraint coefficients
+  /// @param[in] Bieq  mieq by 1 list of linear inequality constraint constant values
+  /// @param[in] lx  n by 1 list of lower bounds [] implies -Inf
+  /// @param[in] ux  n by 1 list of upper bounds [] implies Inf
+  /// @param[in] params  struct of additional parameters (see below)
+  /// @param[in,out] Z  if not empty, is taken to be an n by 1 list of initial guess values. Set to solution on output.
+  /// @return true on success, false on error
+  ///
+  /// \note Benchmark: For a harmonic solve on a mesh with 325K facets, matlab 2.2
+  /// secs, igl/min_quad_with_fixed.h 7.1 secs
+  ///
+  /// \pre rows of [Aeq;Aieq] **must** be linearly independent. Should be
+  /// using QR decomposition otherwise:
+  /// https://v8doc.sas.com/sashtml/ormp/chap5/sect32.htm
+  ///
+  /// \warning This solver is fairly experimental. It works reasonably well for
+  /// bbw problems but doesn't generalize well to other problems.  NASOQ and
+  /// OSQP are better general purpose solvers.
   template <
     typename AT, 
     typename DerivedB,
@@ -64,37 +70,40 @@ namespace igl
     >
   IGL_INLINE igl::SolverStatus active_set(
     const Eigen::SparseMatrix<AT>& A,
-    const Eigen::PlainObjectBase<DerivedB> & B,
-    const Eigen::PlainObjectBase<Derivedknown> & known,
-    const Eigen::PlainObjectBase<DerivedY> & Y,
+    const Eigen::MatrixBase<DerivedB> & B,
+    const Eigen::MatrixBase<Derivedknown> & known,
+    const Eigen::MatrixBase<DerivedY> & Y,
     const Eigen::SparseMatrix<AeqT>& Aeq,
-    const Eigen::PlainObjectBase<DerivedBeq> & Beq,
+    const Eigen::MatrixBase<DerivedBeq> & Beq,
     const Eigen::SparseMatrix<AieqT>& Aieq,
-    const Eigen::PlainObjectBase<DerivedBieq> & Bieq,
-    const Eigen::PlainObjectBase<Derivedlx> & lx,
-    const Eigen::PlainObjectBase<Derivedux> & ux,
+    const Eigen::MatrixBase<DerivedBieq> & Bieq,
+    const Eigen::MatrixBase<Derivedlx> & lx,
+    const Eigen::MatrixBase<Derivedux> & ux,
     const igl::active_set_params & params,
     Eigen::PlainObjectBase<DerivedZ> & Z
     );
 };
 
 #include "EPS.h"
+/// Input parameters controling active_set
+///
+/// \fileinfo
 struct igl::active_set_params
 {
-  // Input parameters for active_set:
-  //   Auu_pd  whether Auu is positive definite {false}
-  //   max_iter  Maximum number of iterations (0 = Infinity, {100})
-  //   inactive_threshold  Threshold on Lagrange multiplier values to determine
-  //     whether to keep constraints active {EPS}
-  //   constraint_threshold  Threshold on whether constraints are violated (0
-  //     is perfect) {EPS}
-  //   solution_diff_threshold  Threshold on the squared norm of the difference
-  //     between two consecutive solutions {EPS}
+///  Auu_pd  whether Auu is positive definite {false}
   bool Auu_pd;
+///  max_iter  Maximum number of iterations (0 = Infinity, {100})
   int max_iter;
+///  inactive_threshold  Threshold on Lagrange multiplier values to determine
+///   whether to keep constraints active {EPS}
   double inactive_threshold;
+///  constraint_threshold  Threshold on whether constraints are violated (0
+///   is perfect) {EPS}
   double constraint_threshold;
+///  solution_diff_threshold  Threshold on the squared norm of the difference
+///    between two consecutive solutions {EPS}
   double solution_diff_threshold;
+  /// @private
   active_set_params():
     Auu_pd(false),
     max_iter(100),

include/igl/adjacency_list.cpp

@@ -70,6 +70,7 @@ IGL_INLINE void igl::adjacency_list(
     for(int v=0; v<(int)SR.size();++v)
     {
       std::vector<IndexVector>& vv = A.at(v);
+      if(vv.size() == 0){ continue; }
       std::vector<std::vector<int> >& sr = SR[v];
       
       std::vector<std::vector<int> > pn = sr;

include/igl/adjacency_list.h

@@ -14,29 +14,35 @@
 #include <vector>
 namespace igl 
 {
-  // Constructs the graph adjacency list of a given mesh (V,F)
-  // Templates:
-  //   T  should be a eigen sparse matrix primitive type like int or double
-  // Inputs:
-  //   F       #F by dim list of mesh faces (must be triangles)
-  //   sorted  flag that indicates if the list should be sorted counter-clockwise
-  // Outputs: 
-  //   A  vector<vector<T> > containing at row i the adjacent vertices of vertex i
-  //
-  // Example:
-  //   // Mesh in (V,F)
-  //   vector<vector<double> > A;
-  //   adjacency_list(F,A);
-  //
-  // See also: edges, cotmatrix, diag
+  /// Constructs the graph adjacency list of a given mesh (V,F)
+  ///
+  /// @tparam T  should be a eigen sparse matrix primitive type like int or double
+  /// @param[in] F       #F by dim list of mesh faces (must be triangles)
+  /// @param[out] A  vector<vector<T> > containing at row i the adjacent vertices of vertex i
+  /// @param[in] sorted  flag that indicates if the list should be sorted counter-clockwise. Input assumed to be manifold.
+  ///
+  /// Example:
+  /// \code{.cpp}
+  ///   // Mesh in (V,F)
+  ///   vector<vector<double> > A;
+  ///   adjacency_list(F,A);
+  /// \endcode
+  ///
+  /// \see 
+  ///   adjacency_matrix
+  ///   edges, 
+  ///   cotmatrix, 
+  ///   diag
   template <typename Index, typename IndexVector>
   IGL_INLINE void adjacency_list(
     const Eigen::MatrixBase<Index>  & F,
     std::vector<std::vector<IndexVector> >& A,
     bool sorted = false);
-
-  // Variant that accepts polygonal faces. 
-  // Each element of F is a set of indices of a polygonal face.
+  /// Constructs the graph adjacency list of a given _polygon_ mesh (V,F)
+  ///
+  /// @tparam T  should be a eigen sparse matrix primitive type like int or double
+  /// @param[in] F       #F list of polygon face index lists
+  /// @param[out] A  vector<vector<T> > containing at row i the adjacent vertices of vertex i
   template <typename Index>
   IGL_INLINE void adjacency_list(
     const std::vector<std::vector<Index> > & F,

include/igl/adjacency_matrix.cpp

@@ -9,6 +9,7 @@
 
 #include "verbose.h"
 
+#include <cassert>
 #include <vector>
 
 template <typename DerivedF, typename T>
@@ -16,12 +17,10 @@ IGL_INLINE void igl::adjacency_matrix(
   const Eigen::MatrixBase<DerivedF> & F,
   Eigen::SparseMatrix<T>& A)
 {
-  using namespace std;
-  using namespace Eigen;
   typedef typename DerivedF::Scalar Index;
 
-  typedef Triplet<T> IJV;
-  vector<IJV > ijv;
+  typedef Eigen::Triplet<T> IJV;
+  std::vector<IJV > ijv;
   ijv.reserve(F.size()*2);
   // Loop over **simplex** (i.e., **not quad**)
   for(int i = 0;i<F.rows();i++)
@@ -71,11 +70,8 @@ IGL_INLINE void igl::adjacency_matrix(
   const Eigen::MatrixBase<DerivedC> & C,
   Eigen::SparseMatrix<T>& A)
 {
-  using namespace std;
-  using namespace Eigen;
-
-  typedef Triplet<T> IJV;
-  vector<IJV > ijv;
+  typedef Eigen::Triplet<T> IJV;
+  std::vector<IJV > ijv;
   ijv.reserve(C(C.size()-1)*2);
   typedef typename DerivedI::Scalar Index;
   const Index n = I.maxCoeff()+1;
@@ -116,6 +112,8 @@ IGL_INLINE void igl::adjacency_matrix(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::adjacency_matrix<Eigen::Matrix<int, -1, 3, 1, -1, 3>, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>> const&, Eigen::SparseMatrix<int, 0, int>&);
 template void igl::adjacency_matrix<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::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<int, 0, int>& );
 // generated by autoexplicit.sh
 template void igl::adjacency_matrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, bool>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<bool, 0, int>&);

include/igl/adjacency_matrix.h

@@ -15,43 +15,44 @@
 
 namespace igl 
 {
-  // Constructs the graph adjacency matrix  of a given mesh (V,F)
-  // Templates:
-  //   T  should be a eigen sparse matrix primitive type like int or double
-  // Inputs:
-  //   F  #F by dim list of mesh simplices
-  // Outputs: 
-  //   A  max(F)+1 by max(F)+1 adjacency matrix, each row i corresponding to V(i,:)
-  //
-  // Example:
-  //   // Mesh in (V,F)
-  //   Eigen::SparseMatrix<double> A;
-  //   adjacency_matrix(F,A);
-  //   // sum each row 
-  //   SparseVector<double> Asum;
-  //   sum(A,1,Asum);
-  //   // Convert row sums into diagonal of sparse matrix
-  //   SparseMatrix<double> Adiag;
-  //   diag(Asum,Adiag);
-  //   // Build uniform laplacian
-  //   SparseMatrix<double> U;
-  //   U = A-Adiag;
-  //
-  // See also: edges, cotmatrix, diag
+  /// Constructs the graph adjacency matrix  of a given mesh (V,F)
+  ///
+  /// @tparam T  should be a eigen sparse matrix primitive type like `int` or `double`
+  /// @param[in] F  #F by dim list of mesh simplices
+  /// @param[out] A  max(F)+1 by max(F)+1 adjacency matrix, each row i corresponding to V(i,:)
+  ///
+  /// #### Example
+  /// \code{.cpp}
+  ///   // Mesh in (V,F)
+  ///   Eigen::SparseMatrix<double> A;
+  ///   adjacency_matrix(F,A);
+  ///   // sum each row 
+  ///   SparseVector<double> Asum;
+  ///   sum(A,1,Asum);
+  ///   // Convert row sums into diagonal of sparse matrix
+  ///   Eigen::SparseMatrix<double> Adiag;
+  ///   diag(Asum,Adiag);
+  ///   // Build uniform laplacian
+  ///   Eigen::SparseMatrix<double> U;
+  ///   U = A-Adiag;
+  /// \endcode
+  ///
+  /// \see 
+  ///   edges,
+  ///   cotmatrix,
+  ///   diag
   template <typename DerivedF, typename T>
   IGL_INLINE void adjacency_matrix(
     const Eigen::MatrixBase<DerivedF> & F, 
     Eigen::SparseMatrix<T>& A);
-  // Constructs an vertex adjacency for a polygon mesh.
-  //
-  // Inputs:
-  //   I  #I vectorized list of polygon corner indices into rows of some matrix V
-  //   C  #polygons+1 list of cumulative polygon sizes so that C(i+1)-C(i) =
-  //     size of the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the
-  //     indices of the ith polygon
-  // Outputs:
-  //   A  max(I)+1 by max(I)+1 adjacency matrix, each row i corresponding to V(i,:)
-  //
+  /// Constructs an vertex adjacency for a polygon mesh.
+  ///
+  /// @param[in] I  #I vectorized list of polygon corner indices into rows of some matrix V
+  /// @param[in] C  #polygons+1 list of cumulative polygon sizes so that C(i+1)-C(i) =
+  ///     size of the ith polygon, and so I(C(i)) through I(C(i+1)-1) are the
+  ///     indices of the ith polygon
+  /// @param[out] A  max(I)+1 by max(I)+1 adjacency matrix, each row i corresponding to V(i,:)
+  ///
   template <typename DerivedI, typename DerivedC, typename T>
   IGL_INLINE void adjacency_matrix(
     const Eigen::MatrixBase<DerivedI> & I,

include/igl/all.h

@@ -12,16 +12,15 @@
 #include <Eigen/Sparse>
 namespace igl
 {
-  // For Dense matrices use: A.rowwise().all() or A.colwise().all()
-  //
-  // Inputs:
-  //   A  m by n sparse matrix
-  //   dim  dimension along which to check for all (1 or 2)
-  // Output:
-  //   B  n-long vector (if dim == 1) 
-  //   or
-  //   B  m-long vector (if dim == 2)
-  //
+  /// Check whether all values are logically true along a dimension.
+  ///
+  /// \note For Dense matrices use: A.rowwise().all() or A.colwise().all()
+  ///
+  /// @param[in]  A  m by n sparse matrix
+  /// @param[in] dim  dimension along which to check for all (1 or 2)
+  /// @param[out] B  n-long vector (if dim == 1) 
+  ///   or m-long vector (if dim == 2)
+  ///
   template <typename AType, typename DerivedB>
   IGL_INLINE void all(
     const Eigen::SparseMatrix<AType> & A, 

include/igl/all_pairs_distances.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "all_pairs_distances.h"
 #include <Eigen/Dense>
+#include <cassert>
 
 template <typename Mat>
 IGL_INLINE void igl::all_pairs_distances(

include/igl/all_pairs_distances.h

@@ -11,21 +11,16 @@
 
 namespace igl
 {
-  // ALL_PAIRS_DISTANCES compute distances between each point i in V and point j
-  // in U
-  // 
-  // D = all_pairs_distances(V,U)
-  // 
-  // Templates:
-  //   Mat  matrix class like MatrixXd
-  // Inputs:
-  //   V  #V by dim list of points
-  //   U  #U by dim list of points
-  //   squared  whether to return squared distances
-  // Outputs:
-  //   D  #V by #U matrix of distances, where D(i,j) gives the distance or
-  //     squareed distance between V(i,:) and U(j,:)
-  // 
+  /// Compute distances between each point i in V and point j in U
+  /// 
+  ///     D = all_pairs_distances(V,U)
+  /// 
+  /// @tparam matrix class like Eigen::MatrixXd
+  /// @param[in] V  #V by dim list of points
+  /// @param[in] U  #U by dim list of points
+  /// @param[in] squared  whether to return squared distances
+  /// @param[out] D  #V by #U matrix of distances, where D(i,j) gives the distance or
+  ///     squareed distance between V(i,:) and U(j,:)
   template <typename Mat>
   IGL_INLINE void all_pairs_distances(
     const Mat & V,

include/igl/ambient_occlusion.cpp

@@ -22,29 +22,31 @@ template <
 IGL_INLINE void igl::ambient_occlusion(
   const std::function<
     bool(
-      const Eigen::Vector3f&,
-      const Eigen::Vector3f&)
+      const Eigen::Matrix<typename DerivedP::Scalar,3,1> &,
+      const Eigen::Matrix<typename DerivedP::Scalar,3,1> &)
       > & shoot_ray,
   const Eigen::MatrixBase<DerivedP> & P,
   const Eigen::MatrixBase<DerivedN> & N,
   const int num_samples,
   Eigen::PlainObjectBase<DerivedS> & S)
 {
-  using namespace Eigen;
   const int n = P.rows();
   // Resize output
   S.resize(n,1);
   // Embree seems to be parallel when constructing but not when tracing rays
-  const MatrixXf D = random_dir_stratified(num_samples).cast<float>();
+  typedef typename DerivedP::Scalar Scalar;
+  typedef Eigen::Matrix<Scalar,3,1> Vector3N;
+
+  const Eigen::Matrix<Scalar,Eigen::Dynamic,3> D = random_dir_stratified(num_samples).cast<Scalar>();
 
   const auto & inner = [&P,&N,&num_samples,&D,&S,&shoot_ray](const int p)
   {
-    const Vector3f origin = P.row(p).template cast<float>();
-    const Vector3f normal = N.row(p).template cast<float>();
+    const Vector3N origin = P.row(p);
+    const Vector3N normal = N.row(p);
     int num_hits = 0;
     for(int s = 0;s<num_samples;s++)
     {
-      Vector3f d = D.row(s);
+      Vector3N d = D.row(s);
       if(d.dot(normal) < 0)
       {
         // reverse ray
@@ -76,17 +78,19 @@ IGL_INLINE void igl::ambient_occlusion(
   const int num_samples,
   Eigen::PlainObjectBase<DerivedS> & S)
 {
+  typedef typename DerivedV::Scalar Scalar;
+  using Vector3S = Eigen::Matrix<Scalar,3,1>;
   const auto & shoot_ray = [&aabb,&V,&F](
-    const Eigen::Vector3f& _s,
-    const Eigen::Vector3f& dir)->bool
+    const Eigen::Matrix<Scalar,3,1> & _s,
+    const Eigen::Matrix<Scalar,3,1> & dir)->bool
   {
-    Eigen::Vector3f s = _s+1e-4*dir;
-    igl::Hit hit;
+    Vector3S s = _s+1e-4*dir;
+    igl::Hit<Scalar> hit;
     return aabb.intersect_ray(
       V,
       F,
-      s  .cast<typename DerivedV::Scalar>().eval(),
-      dir.cast<typename DerivedV::Scalar>().eval(),
+      s,
+      dir,
       hit);
   };
   return ambient_occlusion(shoot_ray,P,N,num_samples,S);
@@ -107,15 +111,17 @@ IGL_INLINE void igl::ambient_occlusion(
   const int num_samples,
   Eigen::PlainObjectBase<DerivedS> & S)
 {
+  typedef typename DerivedV::Scalar Scalar;
+  using Vector3S = Eigen::Matrix<Scalar,3,1>;
   if(F.rows() < 100)
   {
     // Super naive
     const auto & shoot_ray = [&V,&F](
-      const Eigen::Vector3f& _s,
-      const Eigen::Vector3f& dir)->bool
+      const Eigen::Matrix<Scalar,3,1> & _s,
+      const Eigen::Matrix<Scalar,3,1> & dir)->bool
     {
-      Eigen::Vector3f s = _s+1e-4*dir;
-      igl::Hit hit;
+      Vector3S s = _s+1e-4*dir;
+      igl::Hit<Scalar> hit;
       return ray_mesh_intersect(s,dir,V,F,hit);
     };
     return ambient_occlusion(shoot_ray,P,N,num_samples,S);
@@ -127,13 +133,12 @@ IGL_INLINE void igl::ambient_occlusion(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
 template void igl::ambient_occlusion<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::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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
-template void igl::ambient_occlusion<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<bool (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
-template void igl::ambient_occlusion<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<bool (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
-// generated by autoexplicit.sh
-template void igl::ambient_occlusion<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<bool (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
-template void igl::ambient_occlusion<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::function<bool (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::ambient_occlusion<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<bool (Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::ambient_occlusion<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<bool (Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::ambient_occlusion<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<bool (Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&)> const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::ambient_occlusion<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::function<bool (Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 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&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::ambient_occlusion<Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, -1, 1, 0, -1, 1>>(std::function<bool (Eigen::Matrix<Eigen::Matrix<float, 1, 3, 1, 1, 3>::Scalar, 3, 1, 0, 3, 1> const&, Eigen::Matrix<Eigen::Matrix<float, 1, 3, 1, 1, 3>::Scalar, 3, 1, 0, 3, 1> const&)> const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3>> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1>>&);
+template void igl::ambient_occlusion<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 1, 0, -1, 1>>(std::function<bool (Eigen::Matrix<Eigen::Matrix<float, -1, 3, 0, -1, 3>::Scalar, 3, 1, 0, 3, 1> const&, Eigen::Matrix<Eigen::Matrix<float, -1, 3, 0, -1, 3>::Scalar, 3, 1, 0, 3, 1> const&)> const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3>> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1>>&);
+template void igl::ambient_occlusion<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 1, 0, -1, 1>>(std::function<bool (Eigen::Matrix<Eigen::Matrix<float, -1, -1, 0, -1, -1>::Scalar, 3, 1, 0, 3, 1> const&, Eigen::Matrix<Eigen::Matrix<float, -1, -1, 0, -1, -1>::Scalar, 3, 1, 0, 3, 1> const&)> const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1>> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1>>&);
 #endif

include/igl/ambient_occlusion.h

@@ -13,17 +13,17 @@
 #include <functional>
 namespace igl
 {
-  // Compute ambient occlusion per given point
-  //
-  // Inputs:
-  //    shoot_ray  function handle that outputs hits of a given ray against a
-  //      mesh (embedded in function handles as captured variable/data)
-  //    P  #P by 3 list of origin points
-  //    N  #P by 3 list of origin normals
-  // Outputs:
-  //    S  #P list of ambient occlusion values between 1 (fully occluded) and
-  //      0 (not occluded)
-  //
+  /// Compute ambient occlusion per given point using ray-mesh intersection
+  /// function handle.
+  ///
+  /// @param[in]  shoot_ray  function handle that outputs hits of a given ray against a
+  ///               mesh (embedded in function handles as captured variable/data)
+  /// @param[in]  P  #P by 3 list of origin points
+  /// @param[in]  N  #P by 3 list of origin normals
+  /// @param[in] num_samples  number of samples to use (e.g., 1000)
+  /// @param[out]  S  #P list of ambient occlusion values between 1 (fully occluded) and
+  ///      0 (not occluded)
+  ///
   template <
     typename DerivedP,
     typename DerivedN,
@@ -31,15 +31,25 @@ namespace igl
   IGL_INLINE void ambient_occlusion(
     const std::function<
       bool(
-        const Eigen::Vector3f&,
-        const Eigen::Vector3f&)
+        const Eigen::Matrix<typename DerivedP::Scalar,3,1>&,
+        const Eigen::Matrix<typename DerivedP::Scalar,3,1>&)
         > & shoot_ray,
     const Eigen::MatrixBase<DerivedP> & P,
     const Eigen::MatrixBase<DerivedN> & N,
     const int num_samples,
     Eigen::PlainObjectBase<DerivedS> & S);
-  // Inputs:
-  //   AABB  axis-aligned bounding box hierarchy around (V,F)
+  /// Compute ambient occlusion per given point for mesh (V,F) with precomputed
+  /// AABB tree.
+  ///
+  //  @param[in] AABB  axis-aligned bounding box hierarchy around (V,F)
+  /// @param[in] V  #V by 3 list of mesh vertex positions
+  /// @param[in] F  #F by 3 list of mesh face indices into V
+  /// @param[in]  P  #P by 3 list of origin points
+  /// @param[in]  N  #P by 3 list of origin normals
+  /// @param[in] num_samples  number of samples to use (e.g., 1000)
+  /// @param[out]  S  #P list of ambient occlusion values between 1 (fully occluded) and
+  ///      0 (not occluded)
+  ///
   template <
     typename DerivedV,
     int DIM,
@@ -55,9 +65,15 @@ namespace igl
     const Eigen::MatrixBase<DerivedN> & N,
     const int num_samples,
     Eigen::PlainObjectBase<DerivedS> & S);
-  // Inputs:
-  //    V  #V by 3 list of mesh vertex positions
-  //    F  #F by 3 list of mesh face indices into V
+  /// Compute ambient occlusion per given point for mesh (V,F)
+  ///
+  /// @param[in] V  #V by 3 list of mesh vertex positions
+  /// @param[in] F  #F by 3 list of mesh face indices into V
+  /// @param[in]  P  #P by 3 list of origin points
+  /// @param[in]  N  #P by 3 list of origin normals
+  /// @param[in] num_samples  number of samples to use (e.g., 1000)
+  /// @param[out]  S  #P list of ambient occlusion values between 1 (fully occluded) and
+  ///      0 (not occluded)
   template <
     typename DerivedV,
     typename DerivedF,

include/igl/angular_distance.cpp

@@ -6,8 +6,11 @@
 // 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 "angular_distance.h"
-#include <igl/EPS.h>
-#include <igl/PI.h>
+#include "EPS.h"
+#include "PI.h"
+
+#include <cassert>
+
 IGL_INLINE double igl::angular_distance(
   const Eigen::Quaterniond & A,
   const Eigen::Quaterniond & B)

include/igl/angular_distance.h

@@ -11,13 +11,12 @@
 #include <Eigen/Geometry>
 namespace igl
 {
-  // The "angular distance" between two unit quaternions is the angle of the
-  // smallest rotation (treated as an Axis and Angle) that takes A to B.
-  //
-  // Inputs:
-  //   A  unit quaternion
-  //   B  unit quaternion
-  // Returns angular distance
+  /// The "angular distance" between two unit quaternions is the angle of the
+  /// smallest rotation (treated as an Axis and Angle) that takes A to B.
+  ///
+  /// @param[in] A  unit quaternion
+  /// @param[in] B  unit quaternion
+  /// @return angular distance
   IGL_INLINE double angular_distance(
     const Eigen::Quaterniond & A,
     const Eigen::Quaterniond & B);

include/igl/any.h

@@ -12,16 +12,15 @@
 #include <Eigen/Sparse>
 namespace igl
 {
-  // For Dense matrices use: A.rowwise().any() or A.colwise().any()
-  //
-  // Inputs:
-  //   A  m by n sparse matrix
-  //   dim  dimension along which to check for any (1 or 2)
-  // Output:
-  //   B  n-long vector (if dim == 1) 
-  //   or
-  //   B  m-long vector (if dim == 2)
-  //
+  /// Check whether any values are logically true along a dimension.
+  ///
+  /// \note Dense matrices use: A.rowwise().any() or A.colwise().any()
+  ///
+  /// @param[in]  A  m by n sparse matrix
+  /// @param[in] dim  dimension along which to check for any (1 or 2)
+  /// @param[out] B  n-long vector (if dim == 1) 
+  ///   or m-long vector (if dim == 2)
+  ///
   template <typename AType, typename DerivedB>
   IGL_INLINE void any(
     const Eigen::SparseMatrix<AType> & A, 

include/igl/any_of.cpp

@@ -1,20 +0,0 @@
-// 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 "any_of.h"
-#include <Eigen/Core>
-template <typename Mat>
-IGL_INLINE bool igl::any_of(const Mat & S)
-{
-  return std::any_of(S.data(),S.data()+S.size(),[](bool s){return s;});
-}
-
-#ifdef IGL_STATIC_LIBRARY
-// Explicit template instantiation
-template bool igl::any_of<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
-#endif
-

include/igl/any_of.h

@@ -1,26 +0,0 @@
-// 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_ANY_OF_H
-#define IGL_ANY_OF_H
-#include "igl_inline.h"
-namespace igl
-{
-  // Wrapper for STL `any_of` for matrix types
-  //
-  // Inputs:
-  //   S  matrix
-  // Returns whether any entries are true
-  //
-  // Seems that Eigen (now) implements this for `Eigen::Array` 
-  template <typename Mat>
-  IGL_INLINE bool any_of(const Mat & S);
-}
-#ifndef IGL_STATIC_LIBRARY
-#  include "any_of.cpp"
-#endif
-#endif