diff --git a/CHANGELOG.rst b/CHANGELOG.rst index eba7e9bd..e9068fd5 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -11,6 +11,12 @@ adheres to `Semantic Versioning `_. `Unreleased`_ ------------- +Added +~~~~~ + +- Integration of `libigl `_ to efficiently + handle ray-polygon intersection detection in the C++ engine. + Changed ~~~~~~~ diff --git a/CMakeLists.txt b/CMakeLists.txt index 56c2d0ad..7c714ba2 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/external/CMakeLists.txt b/external/CMakeLists.txt index 4df7d600..7cb76de7 100644 --- a/external/CMakeLists.txt +++ b/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) diff --git a/pyroomacoustics/libroom_src/geometry_utils.cpp b/pyroomacoustics/libroom_src/geometry_utils.cpp new file mode 100644 index 00000000..459bc1ad --- /dev/null +++ b/pyroomacoustics/libroom_src/geometry_utils.cpp @@ -0,0 +1,66 @@ +// Included by geometry_utils.hpp +#include +#include +#include + +std::tuple> +triangulate_walls(const std::vector> &walls) { + std::vector all_vertices; + std::vector all_faces; + std::vector 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()); + } + + 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); +} diff --git a/pyroomacoustics/libroom_src/geometry_utils.hpp b/pyroomacoustics/libroom_src/geometry_utils.hpp new file mode 100644 index 00000000..820c7426 --- /dev/null +++ b/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 +#include +#include + +/** + * @brief Triangulates a set of walls into a single mesh. + * + * @param walls Vector of Wall objects (3D). + * @return std::tuple> + * - 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> +triangulate_walls(const std::vector> &walls); + +#include "geometry_utils.cpp" + +#endif // __GEOMETRY_UTILS_H__ diff --git a/pyroomacoustics/libroom_src/libroom.cpp b/pyroomacoustics/libroom_src/libroom.cpp index 83298faf..688147bd 100644 --- a/pyroomacoustics/libroom_src/libroom.cpp +++ b/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 &, const std::vector> &, float, int, float, float, float, float, bool>()) + .def( + py::init> &, const std::vector &, + const Eigen::MatrixXf &, const Eigen::MatrixXi &, + const std::vector &, const std::vector> &, + float, int, float, float, float, float, bool>()) .def(py::init &, const Eigen::Array &, const Eigen::Array &, @@ -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"); } diff --git a/pyroomacoustics/libroom_src/room.cpp b/pyroomacoustics/libroom_src/room.cpp index 37298929..94edde88 100644 --- a/pyroomacoustics/libroom_src/room.cpp +++ b/pyroomacoustics/libroom_src/room.cpp @@ -24,11 +24,10 @@ * If not, see . */ - -#include -#include -#include "constants.hpp" #include "room.hpp" +#include "constants.hpp" +#include +#include size_t number_image_sources_2(size_t max_order) { /* @@ -104,6 +103,29 @@ Room::Room(const Vectorf &_room_size, init(); } +template +Room::Room(const std::vector> &_walls, + const std::vector &_obstructing_walls, + const Eigen::MatrixXf &_V, const Eigen::MatrixXi &_F, + const std::vector &_face_to_wall, + const std::vector> &_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, int, float> Room::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 hit; + Vectorf<3> _start = start.template cast(); + Vectorf<3> _end = end.template cast(); + 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(); + 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, int, float> Room::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::contains(const Vectorf point) { // then the point is in the room => return true return ((n_intersections % 2) == 1); } - diff --git a/pyroomacoustics/libroom_src/room.hpp b/pyroomacoustics/libroom_src/room.hpp index e1596c8b..d601334b 100644 --- a/pyroomacoustics/libroom_src/room.hpp +++ b/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 -#include -#include +#include #include #include #include +#include +#include +#include +#include #include "common.hpp" +#include "microphone.hpp" #include "wall.hpp" template @@ -93,6 +96,12 @@ class Room std::vector> 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 face_to_wall; + igl::AABB 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> &_walls, + const std::vector &_obstructing_walls, + const Eigen::MatrixXf &_V, + const Eigen::MatrixXi &_F, + const std::vector &_face_to_wall, + const std::vector> &_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> &_walls, @@ -288,7 +314,6 @@ class Room bool is_visible_dfs(const Vectorf &p, ImageSource &is, Vectorf &source_direction); bool is_obstructed_dfs(const Vectorf &p, ImageSource &is); int fill_sources(); - }; #include "room.cpp" diff --git a/pyroomacoustics/room.py b/pyroomacoustics/room.py index e0076581..72fc61aa 100644 --- a/pyroomacoustics/room.py +++ b/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