Efficient ray tracing via libigl for geometry processing

CHANGELOG.rst

@@ -11,6 +11,12 @@ adheres to `Semantic Versioning <http://semver.org/spec/v2.0.0.html>`_.
 `Unreleased`_
 -------------
 
+Added
+~~~~~
+
+- Integration of `libigl <https://github.com/libigl/libigl>`_ to efficiently
+  handle ray-polygon intersection detection in the C++ engine.
+
 Changed
 ~~~~~~~
 

CMakeLists.txt

@@ -2,6 +2,8 @@ cmake_minimum_required(VERSION 3.15)
 set(CMAKE_POLICY_VERSION_MINIMUM 3.5)
 
 project(pyroomacoustics)
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 # Find pybind11 (handled by external/CMakeLists.txt via FetchContent)
 set(PYBIND11_FINDPYTHON ON)
@@ -30,6 +32,7 @@ target_include_directories(libroom PRIVATE
 target_link_libraries(libroom PRIVATE
     Eigen3::Eigen
     nanoflann::nanoflann
+    igl::core
 )
 
 # Compile options

external/CMakeLists.txt

@@ -21,6 +21,13 @@ FetchContent_Declare(
   GIT_TAG v3.0.1
 )
 
+# libigl
+FetchContent_Declare(
+  libigl
+  GIT_REPOSITORY https://github.com/libigl/libigl
+  GIT_TAG v2.6.0
+)
+
 # eigen
 set(EIGEN_BUILD_PKGCONFIG OFF CACHE BOOL "Disable Eigen pkg-config")
 set(EIGEN_BUILD_TESTING OFF CACHE BOOL "Disable Eigen testing")
@@ -34,4 +41,4 @@ set(BUILD_BENCHMARKS OFF CACHE BOOL "Disable nanoflann benchmarks" FORCE)
 # pybind11
 set(PYBIND11_TEST OFF CACHE BOOL "Disable pybind11 tests")
 
-FetchContent_MakeAvailable(eigen nanoflann pybind11)
+FetchContent_MakeAvailable(eigen nanoflann pybind11 libigl)

pyroomacoustics/libroom_src/geometry_utils.cpp

@@ -0,0 +1,66 @@
+// Included by geometry_utils.hpp
+#include <igl/polygons_to_triangles.h>
+#include <igl/remove_duplicate_vertices.h>
+#include <iostream>
+
+std::tuple<Eigen::MatrixXf, Eigen::MatrixXi, std::vector<int>>
+triangulate_walls(const std::vector<Wall<3>> &walls) {
+  std::vector<Eigen::Vector3f> all_vertices;
+  std::vector<Eigen::Vector3i> all_faces;
+  std::vector<int> face_to_wall;
+
+  for (size_t i = 0; i < walls.size(); ++i) {
+    const auto &wall = walls[i];
+    int n_corners = wall.flat_corners.cols();
+
+    // Local 2D coordinates for triangulation
+    Eigen::MatrixXd V2(n_corners, 2);
+    for (int j = 0; j < n_corners; ++j) {
+      V2(j, 0) = (double)wall.flat_corners(0, j);
+      V2(j, 1) = (double)wall.flat_corners(1, j);
+    }
+
+    // Cumulative indices for polygons_to_triangles
+    // For a single polygon, C = [0, n_corners]
+    Eigen::VectorXi C(2);
+    C << 0, n_corners;
+    Eigen::VectorXi I(n_corners);
+    for (int j = 0; j < n_corners; ++j)
+      I(j) = j;
+
+    Eigen::MatrixXi F_wall;
+    Eigen::VectorXi J_wall;
+    igl::polygons_to_triangles(I, C, F_wall, J_wall);
+
+    // Map wall local vertex indices to global vertex indices
+    int base_v_idx = all_vertices.size();
+    for (int j = 0; j < n_corners; ++j) {
+      all_vertices.push_back(wall.corners.col(j).cast<float>());
+    }
+
+    for (int j = 0; j < F_wall.rows(); ++j) {
+      all_faces.push_back(Eigen::Vector3i(base_v_idx + F_wall(j, 0),
+                                          base_v_idx + F_wall(j, 1),
+                                          base_v_idx + F_wall(j, 2)));
+      face_to_wall.push_back(i);
+    }
+  }
+
+  // Convert to Eigen matrices
+  Eigen::MatrixXf V(all_vertices.size(), 3);
+  for (size_t i = 0; i < all_vertices.size(); ++i)
+    V.row(i) = all_vertices[i];
+
+  Eigen::MatrixXi F(all_faces.size(), 3);
+  for (size_t i = 0; i < all_faces.size(); ++i)
+    F.row(i) = all_faces[i];
+
+  // Remove duplicate vertices to ensure a watertight mesh
+  // This is important for AABB tree reliability
+  Eigen::MatrixXf V_unique;
+  Eigen::MatrixXi F_unique;
+  Eigen::VectorXi S, _I;
+  igl::remove_duplicate_vertices(V, F, 1e-5, V_unique, S, _I, F_unique);
+
+  return std::make_tuple(V_unique, F_unique, face_to_wall);
+}

pyroomacoustics/libroom_src/geometry_utils.hpp

@@ -0,0 +1,24 @@
+#ifndef __GEOMETRY_UTILS_H__
+#define __GEOMETRY_UTILS_H__
+
+#include "common.hpp"
+#include "wall.hpp"
+#include <Eigen/Dense>
+#include <tuple>
+#include <vector>
+
+/**
+ * @brief Triangulates a set of walls into a single mesh.
+ *
+ * @param walls Vector of Wall objects (3D).
+ * @return std::tuple<Eigen::MatrixXf, Eigen::MatrixXi, std::vector<int>>
+ *         - V: Matrix of vertices (n_vertices x 3)
+ *         - F: Matrix of faces (n_faces x 3)
+ *         - face_to_wall: Vector mapping each face to its original wall index.
+ */
+std::tuple<Eigen::MatrixXf, Eigen::MatrixXi, std::vector<int>>
+triangulate_walls(const std::vector<Wall<3>> &walls);
+
+#include "geometry_utils.cpp"
+
+#endif // __GEOMETRY_UTILS_H__

pyroomacoustics/libroom_src/libroom.cpp

@@ -35,6 +35,7 @@
 
 #include "common.hpp"
 #include "geometry.hpp"
+#include "geometry_utils.hpp"
 #include "microphone.hpp"
 #include "random.hpp"
 #include "rir_builder.hpp"
@@ -56,6 +57,11 @@ PYBIND11_MODULE(libroom, m) {
       .def(py::init<const std::vector<Wall<3>> &, const std::vector<int> &,
                     const std::vector<Microphone<3>> &, float, int, float,
                     float, float, float, bool>())
+      .def(
+          py::init<const std::vector<Wall<3>> &, const std::vector<int> &,
+                   const Eigen::MatrixXf &, const Eigen::MatrixXi &,
+                   const std::vector<int> &, const std::vector<Microphone<3>> &,
+                   float, int, float, float, float, float, bool>())
       .def(py::init<const Vectorf<3> &,
                     const Eigen::Array<float, Eigen::Dynamic, 6> &,
                     const Eigen::Array<float, Eigen::Dynamic, 6> &,
@@ -330,4 +336,7 @@ PYBIND11_MODULE(libroom, m) {
   m.def("rir_builder", &rir_builder, "RIR builder");
   m.def("delay_sum", &delay_sum, "Delay and sum");
   m.def("fractional_delay", &fractional_delay, "Fractional delays");
+
+  m.def("triangulate_walls", &triangulate_walls,
+        "Triangulates a set of walls into a mesh");
 }

pyroomacoustics/libroom_src/room.cpp

@@ -24,11 +24,10 @@
  * If not, see <https://opensource.org/licenses/MIT>.
  */
 
-
-#include <cmath>
-#include <algorithm>
-#include "constants.hpp"
 #include "room.hpp"
+#include "constants.hpp"
+#include <algorithm>
+#include <cmath>
 
 size_t number_image_sources_2(size_t max_order) {
   /*
@@ -104,6 +103,29 @@ Room<D>::Room(const Vectorf<D> &_room_size,
   init();
 }
 
+template <size_t D>
+Room<D>::Room(const std::vector<Wall<D>> &_walls,
+              const std::vector<int> &_obstructing_walls,
+              const Eigen::MatrixXf &_V, const Eigen::MatrixXi &_F,
+              const std::vector<int> &_face_to_wall,
+              const std::vector<Microphone<D>> &_microphones,
+              float _sound_speed, int _ism_order, float _energy_thres,
+              float _time_thres, float _mic_radius, float _mic_hist_res,
+              bool _is_hybrid_sim)
+    : walls(_walls), obstructing_walls(_obstructing_walls),
+      microphones(_microphones), V_mesh(_V), F_mesh(_F),
+      face_to_wall(_face_to_wall), sound_speed(_sound_speed),
+      ism_order(_ism_order), energy_thres(_energy_thres),
+      time_thres(_time_thres), mic_radius(_mic_radius),
+      mic_radius_sq(_mic_radius * _mic_radius), mic_hist_res(_mic_hist_res),
+      is_hybrid_sim(_is_hybrid_sim), is_shoebox(false) {
+  init();
+  // Initialize AABB tree if mesh is provided
+  if (F_mesh.rows() > 0) {
+    mesh_aabb.init(V_mesh, F_mesh);
+  }
+}
+
 template <>
 void Room<2>::make_shoebox_walls(
     const Vectorf<2> &rs,  // room_size
@@ -669,6 +691,28 @@ std::tuple<Vectorf<D>, int, float> Room<D>::next_wall_hit(
       (void)0;  // no op
     }
   } else {
+    if constexpr (D == 3) {
+      if (F_mesh.rows() > 0 && !scattered_ray) {
+        // Use libigl AABB tree for 3D ray tracing (non-scattered)
+        igl::Hit<float> hit;
+        Vectorf<3> _start = start.template cast<float>();
+        Vectorf<3> _end = end.template cast<float>();
+        Vectorf<3> _dir = (_end - _start).normalized();
+        // Offset start to avoid hitting the same wall
+        Vectorf<3> _start_offset = _start + 1e-4f * _dir;
+        if (mesh_aabb.intersect_ray(V_mesh, F_mesh, _start_offset, _dir, hit)) {
+          // hit.t is the distance to the hit point from _start_offset
+          float actual_t = hit.t + 1e-4f;
+          if (actual_t < hit_dist) {
+            hit_dist = actual_t;
+            result = start + actual_t * _dir.template cast<float>();
+            next_wall_index = face_to_wall[hit.id];
+          }
+        }
+        return std::make_tuple(result, next_wall_index, hit_dist);
+      }
+    }
+
     // For case 1) in non-convex rooms, the segment might intersect several
     // walls. In this case, we are only interested on the closest wall to
     // 'start'. That's why we need a min_dist variable
@@ -696,7 +740,7 @@ std::tuple<Vectorf<D>, int, float> Room<D>::next_wall_hit(
         if (temp_dist > libroom_eps && temp_dist < hit_dist) {
           hit_dist = temp_dist;
           result = temp_hit;
-          next_wall_index = i;
+          next_wall_index = scattered_ray ? obstructing_walls[i] : i;
         }
       }
     }
@@ -1131,4 +1175,3 @@ bool Room<D>::contains(const Vectorf<D> point) {
   // then the point is in the room  => return true
   return ((n_intersections % 2) == 1);
 }
-

pyroomacoustics/libroom_src/room.hpp

@@ -1,4 +1,4 @@
-/* 
+/*
  * Definition of the Room class
  * Copyright (C) 2019  Robin Scheibler, Cyril Cadoux
  *
@@ -26,14 +26,17 @@
 #ifndef __ROOM_H__
 #define __ROOM_H__
 
-#include <vector>
-#include <stack>
-#include <tuple>
+#include <Eigen/Core>
 #include <Eigen/Dense>
 #include <algorithm>
 #include <ctime>
+#include <igl/AABB.h>
+#include <stack>
+#include <tuple>
+#include <vector>
 
 #include "common.hpp"
+#include "microphone.hpp"
 #include "wall.hpp"
 
 template<size_t D>
@@ -93,6 +96,12 @@ class Room
     std::vector<Microphone<D>> microphones;  // The microphones are in the room
     float sound_speed = 343.;  // the speed of sound in the room
 
+    // libigl mesh for ray tracing
+    Eigen::MatrixXf V_mesh;
+    Eigen::MatrixXi F_mesh;
+    std::vector<int> face_to_wall;
+    igl::AABB<Eigen::MatrixXf, 3> mesh_aabb;
+
     // Simulation parameters
     int ism_order = 0.;
 
@@ -127,6 +136,23 @@ class Room
     // its size is n_microphones * n_sources
     MatrixXb visible_mics;
 
+    // Constructor for general rooms with mesh
+    Room(
+        const std::vector<Wall<D>> &_walls,
+        const std::vector<int> &_obstructing_walls,
+        const Eigen::MatrixXf &_V,
+        const Eigen::MatrixXi &_F,
+        const std::vector<int> &_face_to_wall,
+        const std::vector<Microphone<D>> &_microphones,
+        float _sound_speed,
+        int _ism_order,
+        float _energy_thres,
+        float _time_thres,
+        float _mic_radius,
+        float _mic_hist_res,
+        bool _is_hybrid_sim
+        );
+
     // Constructor for general rooms
     Room(
         const std::vector<Wall<D>> &_walls,
@@ -288,7 +314,6 @@ class Room
     bool is_visible_dfs(const Vectorf<D> &p, ImageSource<D> &is, Vectorf<D> &source_direction);
     bool is_obstructed_dfs(const Vectorf<D> &p, ImageSource<D> &is);
     int fill_sources();
-
 };
 
 #include "room.cpp"

pyroomacoustics/room.py

@@ -882,11 +882,41 @@ def __init__(
         use_rand_ism=False,
         max_rand_disp=0.08,
         min_phase=False,
+        V=None,
+        F=None,
+        face_to_wall=None,
     ):
-        self.walls = walls
+        if walls is None and V is not None and F is not None:
+            # derive walls from V and F
+            # if face_to_wall is not provided, we assume each face is a wall
+            # if provided, we group faces by wall index
+            if face_to_wall is None:
+                face_to_wall = np.arange(F.shape[0])
+
+            unique_walls = np.unique(face_to_wall)
+            self.walls = []
+            for w_idx in unique_walls:
+                faces = F[face_to_wall == w_idx]
+                # For simplicity, we just use the vertices of the first face
+                # but pyroomacoustics Walls are polygons.
+                # If we have multiple triangles for a single wall, we can't
+                # easily represent them as a single polygon if it's not simple.
+                # However, for simulation, the mesh is what matters.
+                # Here we create a dummy wall for each triangle if they are not grouped.
+                for face in faces:
+                    corners = V[face].T
+                    self.walls.append(
+                        wall_factory(corners, [0.0], [0.0], name=f"wall_{w_idx}")
+                    )
+        else:
+            self.walls = walls
+
+        self.V = V
+        self.F = F
+        self.face_to_wall = face_to_wall
 
         # Get the room dimension from that of the walls
-        self.dim = walls[0].dim
+        self.dim = self.walls[0].dim
 
         # Create a mapping with friendly names for walls
         self._wall_mapping()
@@ -999,7 +1029,23 @@ def _init_room_engine(self, *args):
             # This is a polygonal room
             # find the non convex walls
             obstructing_walls = find_non_convex_walls(self.walls)
-            args += [self.walls, obstructing_walls]
+
+            if self.dim == 3:
+                # If mesh is not provided, we triangulate
+                if self.V is None or self.F is None:
+                    self.V, self.F, self.face_to_wall = libroom.triangulate_walls(
+                        self.walls
+                    )
+
+                args += [
+                    self.walls,
+                    obstructing_walls,
+                    self.V,
+                    self.F,
+                    self.face_to_wall,
+                ]
+            else:
+                args += [self.walls, obstructing_walls]
 
         # for shoebox rooms, the required arguments are passed to
         # the function