diff --git a/.gitmodules b/.gitmodules index 4c6e4cba..05df7510 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,6 @@ [submodule "pyroomacoustics/libroom_src/ext/eigen"] path = pyroomacoustics/libroom_src/ext/eigen url = https://github.com/eigenteam/eigen-git-mirror.git +[submodule "pyroomacoustics/libroom_src/ext/nanoflann"] + path = pyroomacoustics/libroom_src/ext/nanoflann + url = https://github.com/jlblancoc/nanoflann.git diff --git a/MANIFEST.in b/MANIFEST.in index 0b82c592..0e910b89 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -11,9 +11,11 @@ include pyroomacoustics/libroom_src/*.h include pyroomacoustics/libroom_src/*.hpp include pyroomacoustics/libroom_src/*.cpp -# The compiled extension code rely on Eigen, which is included +# The compiled extension code rely on Eigen and nanoflann, +# which are included include pyroomacoustics/libroom_src/ext/eigen/COPYING.* graft pyroomacoustics/libroom_src/ext/eigen/Eigen +graft pyroomacoustics/libroom_src/ext/nanoflann/include include pyproject.toml include requirements.txt diff --git a/examples/plot_directivity_2D.py b/examples/plot_directivity_2D.py new file mode 100644 index 00000000..b380003c --- /dev/null +++ b/examples/plot_directivity_2D.py @@ -0,0 +1,33 @@ +import matplotlib.pyplot as plt +import numpy as np + +from pyroomacoustics import dB +from pyroomacoustics.directivities import ( + CardioidFamily, + DirectionVector, + DirectivityPattern, +) +from pyroomacoustics.doa import spher2cart + +ORIENTATION = DirectionVector(azimuth=0, colatitude=90, degrees=True) +LOWER_GAIN = -20 + +# plot each directivity +angles = np.linspace(start=0, stop=360, num=361, endpoint=True) +angles = np.radians(angles) + +# plot each pattern +fig = plt.figure() +ax = plt.subplot(111, projection="polar") +for pattern in DirectivityPattern: + + dir_obj = CardioidFamily(orientation=ORIENTATION, pattern_enum=pattern) + resp = dir_obj.get_response(angles, magnitude=True, degrees=False) + resp_db = dB(np.array(resp)) + ax.plot(angles, resp_db, label=pattern.name) + +plt.legend(bbox_to_anchor=(1, 1)) +plt.ylim([LOWER_GAIN, 0]) +ax.yaxis.set_ticks(np.arange(start=LOWER_GAIN, stop=5, step=5)) +plt.tight_layout() +plt.show() diff --git a/examples/plot_directivity_3D.py b/examples/plot_directivity_3D.py new file mode 100644 index 00000000..7887ac57 --- /dev/null +++ b/examples/plot_directivity_3D.py @@ -0,0 +1,21 @@ +import matplotlib.pyplot as plt +import numpy as np + +from pyroomacoustics.directivities import ( + CardioidFamily, + DirectionVector, + DirectivityPattern, +) + +PATTERN = DirectivityPattern.HYPERCARDIOID +ORIENTATION = DirectionVector(azimuth=0, colatitude=45, degrees=True) + +# create cardioid object +dir_obj = CardioidFamily(orientation=ORIENTATION, pattern_enum=PATTERN) + +# plot +azimuth = np.linspace(start=0, stop=360, num=361, endpoint=True) +colatitude = np.linspace(start=0, stop=180, num=180, endpoint=True) +# colatitude = None # for 2D plot +dir_obj.plot_response(azimuth=azimuth, colatitude=colatitude, degrees=True) +plt.show() diff --git a/examples/plot_directivity_one_shot.py b/examples/plot_directivity_one_shot.py new file mode 100644 index 00000000..9a41c2c7 --- /dev/null +++ b/examples/plot_directivity_one_shot.py @@ -0,0 +1,60 @@ +import matplotlib.pyplot as plt +import numpy as np + +from pyroomacoustics import all_combinations, dB +from pyroomacoustics.directivities import ( + CardioidFamily, + DirectionVector, + DirectivityPattern, + cardioid_func, +) +from pyroomacoustics.doa import spher2cart + +ORIENTATION = DirectionVector(azimuth=0, colatitude=90, degrees=True) +azimuth = np.radians(np.linspace(start=0, stop=360, num=361, endpoint=True)) +colatitude = np.radians(np.linspace(start=0, stop=180, num=180, endpoint=True)) +LOWER_GAIN = -40 + +""" 2D """ +# get cartesian coordinates +cart = spher2cart(azimuth=azimuth) +direction = spher2cart(azimuth=225, degrees=True) + +# compute response +resp = cardioid_func(x=cart, direction=direction, coef=0.5, magnitude=True) +resp_db = dB(np.array(resp)) + +# plot +plt.figure() +plt.polar(azimuth, resp_db) +plt.ylim([LOWER_GAIN, 0]) +ax = plt.gca() +ax.yaxis.set_ticks(np.arange(start=LOWER_GAIN, stop=5, step=10)) +plt.tight_layout() + +""" 3D """ +# get cartesian coordinates +spher_coord = all_combinations(azimuth, colatitude) +cart = spher2cart(azimuth=spher_coord[:, 0], colatitude=spher_coord[:, 1]) +direction = spher2cart(azimuth=0, colatitude=45, degrees=True) + +# compute response +resp = cardioid_func(x=cart, direction=direction, coef=0.25, magnitude=True) + +# plot (surface plot) +fig = plt.figure() +RESP_2D = resp.reshape(len(azimuth), len(colatitude)) +AZI, COL = np.meshgrid(azimuth, colatitude) +X = RESP_2D.T * np.sin(COL) * np.cos(AZI) +Y = RESP_2D.T * np.sin(COL) * np.sin(AZI) +Z = RESP_2D.T * np.cos(COL) +ax = fig.add_subplot(1, 1, 1, projection="3d") +ax.plot_surface(X, Y, Z) +ax.set_xlabel("x") +ax.set_ylabel("y") +ax.set_zlabel("z") +ax.set_xlim([-1, 1]) +ax.set_ylim([-1, 1]) +ax.set_zlim([-1, 1]) + +plt.show() diff --git a/pyroomacoustics/__init__.py b/pyroomacoustics/__init__.py index dcc9125c..19b4e633 100644 --- a/pyroomacoustics/__init__.py +++ b/pyroomacoustics/__init__.py @@ -122,7 +122,7 @@ from . import adaptive, bss, datasets, denoise, doa, experimental from . import libroom as libroom -from . import phase, transform +from . import phase, random, transform from .acoustics import * from .beamforming import * from .directivities import * diff --git a/pyroomacoustics/directivities/__init__.py b/pyroomacoustics/directivities/__init__.py index 671b7a25..27c0c0f3 100644 --- a/pyroomacoustics/directivities/__init__.py +++ b/pyroomacoustics/directivities/__init__.py @@ -84,6 +84,7 @@ from .analytic import ( Cardioid, CardioidFamily, + CardioidFamilySampler, FigureEight, HyperCardioid, Omnidirectional, diff --git a/pyroomacoustics/directivities/analytic.py b/pyroomacoustics/directivities/analytic.py index 59a3b28d..73c35b21 100644 --- a/pyroomacoustics/directivities/analytic.py +++ b/pyroomacoustics/directivities/analytic.py @@ -59,6 +59,7 @@ """ import numpy as np +from .. import random from ..doa import spher2cart from ..utilities import all_combinations, requires_matplotlib from .base import Directivity @@ -122,6 +123,12 @@ def __init__(self, orientation, p, gain=1.0): self._gain = gain self._pattern_name = f"cardioid family, p={self._p}" + # this is the object that will allow to sample rays according to the + # distribution corresponding to the source directivity + self._ray_sampler = CardioidFamilySampler( + loc=self._orientation.unit_vector, p=self._p + ) + @property def is_impulse_response(self): # this is not an impulse response, do not make docstring to avoid clutter in the @@ -200,6 +207,9 @@ def get_response( else: return resp + def sample_rays(self, n_rays): + return self._ray_sampler(n_rays).T + @requires_matplotlib def plot_response( self, azimuth, colatitude=None, degrees=True, ax=None, offset=None @@ -383,6 +393,39 @@ def __init__(self, gain=1.0): self._pattern_name = "omni" +class CardioidFamilySampler(random.sampler.DirectionalSampler): + """ + This object draws samples from a cardioid shaped distribution on the sphere + + Parameters + ---------- + loc: array_like + The unit vector pointing in the main direction of the cardioid + p: float + Parameter of the cardioid pattern. A value of 0 corresponds to a + figure-eight pattern, 0.5 to a cardioid pattern, and 1 to an omni + pattern + The parameter must be between 0 and 1 + """ + + def __init__(self, loc=None, p=0.5): + super().__init__(loc=loc) + self._coeff = p + + def _pattern(self, x): + response = cardioid_func( + x.T, + direction=self._loc, + p=self._coeff, + gain=1.0, + normalize=False, + magnitude=True, + ) + # The number of rays needs to be proportional to the + # response energy. + return response**2 + + def cardioid_func(x, direction, p, gain=1.0, normalize=True, magnitude=False): """ One-shot function for computing cardioid response. diff --git a/pyroomacoustics/directivities/base.py b/pyroomacoustics/directivities/base.py index 349f25f9..4c5e8c6a 100644 --- a/pyroomacoustics/directivities/base.py +++ b/pyroomacoustics/directivities/base.py @@ -92,3 +92,21 @@ def get_response( Response at provided angles. """ raise NotImplementedError + + @abc.abstractmethod + def sample_rays(self, n_rays): + """ + This method samples unit vectors from the sphere according to + the distribution of the source + + Parameters + ---------- + n_rays: int + The number of rays to sample + + Returns + ------- + ray_directions: numpy.ndarray, shape (n_dim, n_rays) + An array containing the unit vectors in its columns + """ + raise NotImplementedError diff --git a/pyroomacoustics/directivities/measured.py b/pyroomacoustics/directivities/measured.py index 2f95f04a..d4b1f4fd 100644 --- a/pyroomacoustics/directivities/measured.py +++ b/pyroomacoustics/directivities/measured.py @@ -75,6 +75,7 @@ from scipy.interpolate import griddata from scipy.spatial import cKDTree +from .. import random from ..datasets import SOFADatabase from ..doa import Grid, GridSphere, cart2spher, fibonacci_spherical_sampling, spher2cart from ..utilities import requires_matplotlib @@ -84,6 +85,34 @@ from .sofa import open_sofa_file +class MeasuredDirectivitySampler(random.sampler.DirectionalSampler): + """ + This object draws samples from the distribution defined by the energy + of the measured directional response object. + + Parameters + ---------- + loc: array_like + The unit vector pointing in the main direction of the cardioid + p: float + Parameter of the cardioid pattern. A value of 0 corresponds to a + figure-eight pattern, 0.5 to a cardioid pattern, and 1 to an omni + pattern + The parameter must be between 0 and 1 + """ + + def __init__(self, kdtree, energy): + super().__init__() + self._kdtree = kdtree + # Normalize to maximum energy 1 because that is also the + # maximum value of the proposal unnormalized uniform distribution. + self._energy = energy / energy.max() + + def _pattern(self, x): + _, index = self._kdtree.query(x) + return self._energy[index] + + class MeasuredDirectivity(Directivity): """ A class to store directivity patterns obtained by measurements. @@ -144,6 +173,10 @@ def set_orientation(self, orientation): # create the kd-tree self._kdtree = cKDTree(self._grid.cartesian.T) + # create the ray sampler + ir_energy = np.square(self._irs).mean(axis=-1) + self._ray_sampler = MeasuredDirectivitySampler(self._kdtree, ir_energy) + def get_response( self, azimuth, colatitude=None, magnitude=False, frequency=None, degrees=True ): @@ -172,6 +205,9 @@ def get_response( _, index = self._kdtree.query(cart.T) return self._irs[index, :] + def sample_rays(self, n_rays): + return self._ray_sampler(n_rays).T + @requires_matplotlib def plot(self, freq_bin=0, n_grid=100, ax=None, depth=False, offset=None): """ @@ -238,7 +274,7 @@ def plot(self, freq_bin=0, n_grid=100, ax=None, depth=False, offset=None): Y *= V Z *= V - surf = ax.plot_surface( + _ = ax.plot_surface( X, Y, Z, facecolors=cmap.to_rgba(V), linewidth=0, antialiased=False ) diff --git a/pyroomacoustics/doa/__init__.py b/pyroomacoustics/doa/__init__.py index 057d52c6..bdf9b461 100644 --- a/pyroomacoustics/doa/__init__.py +++ b/pyroomacoustics/doa/__init__.py @@ -112,6 +112,7 @@ from .doa import * from .frida import * from .grid import * +from .histogram import SphericalHistogram from .music import * from .normmusic import * from .srp import * diff --git a/pyroomacoustics/doa/grid.py b/pyroomacoustics/doa/grid.py index c0c4ed85..36b1cec4 100644 --- a/pyroomacoustics/doa/grid.py +++ b/pyroomacoustics/doa/grid.py @@ -209,7 +209,7 @@ def __init__( # If no list was provided, samples points on the sphere # as uniformly as possible - self.x, self.y, self.z = fibonacci_spherical_sampling(n_points) + self.x[:], self.y[:], self.z[:] = fibonacci_spherical_sampling(n_points) # Create convenient arrays # to access both in cartesian and spherical coordinates @@ -389,10 +389,7 @@ def plot( def plot_old(self, plot_points=False, mark_peaks=0): """Plot the points on the sphere with their values""" - from scipy import rand - try: - import matplotlib.colors as colors import matplotlib.pyplot as plt # from mpl_toolkits.mplot3d import Axes3D diff --git a/pyroomacoustics/doa/histogram.py b/pyroomacoustics/doa/histogram.py new file mode 100644 index 00000000..5c9c4cf4 --- /dev/null +++ b/pyroomacoustics/doa/histogram.py @@ -0,0 +1,110 @@ +# This module provides a class to plot histograms of data on the sphere. +# Copyright (C) 2024 Robin Scheibler +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. +# +# You should have received a copy of the MIT License along with this program. If +# not, see . +r""" +A class to make collect histograms of data distributed on the sphere. +""" +import numpy as np +from scipy.spatial import SphericalVoronoi, cKDTree + +from .doa import GridSphere + + +class SphericalHistogram: + def __init__(self, n_bins, dim=3, enable_peak_finding=False): + + self._n_dim = 3 + self._n_bins = n_bins + + if self.n_dim == 3: + self._grid = GridSphere( + n_points=self.n_bins, precompute_neighbors=enable_peak_finding + ) + else: + raise NotImplementedError("Only 3D histogram has been implemented") + + # we need to know the area of each bin + self._voronoi = SphericalVoronoi(self._grid.cartesian.T) + self._areas = self._voronoi.calculate_areas() + + # now we also need a KD-tree to do nearest neighbor search + self._kd_tree = cKDTree(self._grid.cartesian.T) + + # the counter variables for every bin + self._bins = np.zeros(self.n_bins, dtype=int) + + # the total number of points in the histogram + self._total_count = 0 + + # we cache the histogram bins + self._cache_dirty = False + self._cache_histogram = np.zeros(self.n_bins) + + @property + def n_dim(self): + return self._n_dim + + @property + def n_bins(self): + return self._n_bins + + @property + def histogram(self): + if self._cache_dirty: + # if the cache is dirty, we need to recompute + Z = np.sum(self._areas * self._bins) # partitioning constant + self._cache_histogram[:] = self._bins / Z + self._cache_dirty = False + + return self._cache_histogram + + @property + def raw_counts(self): + return self._bins + + @property + def total_count(self): + return self._total_count + + def find_peak(self, *args, **kwargs): + return self._grid.find_peaks(self, *args, **kwargs) + + def plot(self): + self._grid.set_values(self.histogram) + self._grid.plot_old() + + def push(self, points): + """ + Add new data into the histogram + + Parameters + ---------- + points: array_like, shape (n_dim, n_points) + The points to add to the histogram + """ + self._total_count += points.shape[1] + self._cache_dirty = True + + _, matches = self._kd_tree.query(points.T) + bin_indices, counts = np.unique(matches, return_counts=True) + self._bins[bin_indices] += counts diff --git a/pyroomacoustics/libroom_src/ext/nanoflann b/pyroomacoustics/libroom_src/ext/nanoflann new file mode 160000 index 00000000..b3e81831 --- /dev/null +++ b/pyroomacoustics/libroom_src/ext/nanoflann @@ -0,0 +1 @@ +Subproject commit b3e81831b847a77473e0da2ad7ee266b4f4e0d9a diff --git a/pyroomacoustics/libroom_src/libroom.cpp b/pyroomacoustics/libroom_src/libroom.cpp index 444c6be2..6f2e55c6 100644 --- a/pyroomacoustics/libroom_src/libroom.cpp +++ b/pyroomacoustics/libroom_src/libroom.cpp @@ -65,17 +65,29 @@ PYBIND11_MODULE(libroom, m) { .def("get_max_distance", &Room<3>::get_max_distance) .def("next_wall_hit", &Room<3>::next_wall_hit) .def("scat_ray", &Room<3>::scat_ray) - .def("simul_ray", &Room<3>::simul_ray) + .def("simul_ray", + (void (Room<3>::*)(float phi, float theta, + const Vectorf<3> &source_pos, float energy_0)) & + Room<3>::simul_ray) + .def("simul_ray", + (void (Room<3>::*)(const Vectorf<3> &ray_direction, + const Vectorf<3> &source_pos, float energy_0)) & + Room<3>::simul_ray) .def("ray_tracing", - (void(Room<3>::*)( + (void (Room<3>::*)( const Eigen::Matrix &angles, - const Vectorf<3> source_pos)) & + const Vectorf<3> &source_pos)) & Room<3>::ray_tracing) - .def("ray_tracing", (void(Room<3>::*)(size_t nb_phis, size_t nb_thetas, - const Vectorf<3> source_pos)) & + .def("ray_tracing", + (void (Room<3>::*)( + const Eigen::Matrix &vectors, + const Vectorf<3> &source_pos)) & + Room<3>::ray_tracing) + .def("ray_tracing", (void (Room<3>::*)(size_t nb_phis, size_t nb_thetas, + const Vectorf<3> &source_pos)) & Room<3>::ray_tracing) .def("ray_tracing", - (void(Room<3>::*)(size_t nb_rays, const Vectorf<3> source_pos)) & + (void (Room<3>::*)(size_t nb_rays, const Vectorf<3> &source_pos)) & Room<3>::ray_tracing) .def("contains", &Room<3>::contains) .def_property("is_hybrid_sim", &Room<3>::get_is_hybrid_sim, @@ -86,6 +98,7 @@ PYBIND11_MODULE(libroom, m) { .def_readonly("sources", &Room<3>::sources) .def_readonly("orders", &Room<3>::orders) .def_readonly("orders_xyz", &Room<3>::orders_xyz) + .def_readonly("source_directions", &Room<3>::source_directions) .def_readonly("attenuations", &Room<3>::attenuations) .def_readonly("gen_walls", &Room<3>::gen_walls) .def_readonly("visible_mics", &Room<3>::visible_mics) @@ -113,17 +126,29 @@ PYBIND11_MODULE(libroom, m) { .def("get_max_distance", &Room<2>::get_max_distance) .def("next_wall_hit", &Room<2>::next_wall_hit) .def("scat_ray", &Room<2>::scat_ray) - .def("simul_ray", &Room<2>::simul_ray) + .def("simul_ray", + (void (Room<2>::*)(float phi, float theta, + const Vectorf<2> &source_pos, float energy_0)) & + Room<2>::simul_ray) + .def("simul_ray", + (void (Room<2>::*)(const Vectorf<2> &ray_direction, + const Vectorf<2> &source_pos, float energy_0)) & + Room<2>::simul_ray) .def("ray_tracing", - (void(Room<2>::*)( + (void (Room<2>::*)( const Eigen::Matrix &angles, - const Vectorf<2> source_pos)) & + const Vectorf<2> &source_pos)) & + Room<2>::ray_tracing) + .def("ray_tracing", + (void (Room<2>::*)( + const Eigen::Matrix &vectors, + const Vectorf<2> &source_pos)) & Room<2>::ray_tracing) - .def("ray_tracing", (void(Room<2>::*)(size_t nb_phis, size_t nb_thetas, - const Vectorf<2> source_pos)) & + .def("ray_tracing", (void (Room<2>::*)(size_t nb_phis, size_t nb_thetas, + const Vectorf<2> &source_pos)) & Room<2>::ray_tracing) .def("ray_tracing", - (void(Room<2>::*)(size_t n_rays, const Vectorf<2> source_pos)) & + (void (Room<2>::*)(size_t n_rays, const Vectorf<2> &source_pos)) & Room<2>::ray_tracing) .def("contains", &Room<2>::contains) .def_property_readonly_static("dim", @@ -134,6 +159,7 @@ PYBIND11_MODULE(libroom, m) { .def_readonly("sources", &Room<2>::sources) .def_readonly("orders", &Room<2>::orders) .def_readonly("orders_xyz", &Room<2>::orders_xyz) + .def_readonly("source_directions", &Room<2>::source_directions) .def_readonly("attenuations", &Room<2>::attenuations) .def_readonly("gen_walls", &Room<2>::gen_walls) .def_readonly("visible_mics", &Room<2>::visible_mics) @@ -225,16 +251,18 @@ PYBIND11_MODULE(libroom, m) { // The microphone class py::class_>(m, "Microphone") - .def(py::init &, int, float, float>()) - .def_readonly("loc", &Microphone<3>::loc) - .def_readonly("hits", &Microphone<3>::hits) - .def_readonly("histograms", &Microphone<3>::histograms); + .def(py::init &, int, float, float>()) + .def("set_directions", &Microphone<3>::set_directions) + .def("log_histogram", (void(Microphone<3>::*)(float, const Eigen::ArrayXf &, const Vectorf<3> &))&Microphone<3>::log_histogram) + .def_readonly("loc", &Microphone<3>::loc) + .def_readonly("histograms", &Microphone<3>::histograms) + ; py::class_>(m, "Microphone2D") - .def(py::init &, int, float, float>()) - .def_readonly("loc", &Microphone<2>::loc) - .def_readonly("hits", &Microphone<2>::hits) - .def_readonly("histograms", &Microphone<2>::histograms); + .def(py::init &, int, float, float>()) + .def_readonly("loc", &Microphone<2>::loc) + .def_readonly("histograms", &Microphone<2>::histograms) + ; // The 2D histogram class py::class_(m, "Histogram2D") @@ -252,8 +280,10 @@ PYBIND11_MODULE(libroom, m) { .def_readonly("distance", &Hit::distance); // getter and setter for geometric epsilon - m.def("set_eps", [](const float &eps) { libroom_eps = eps; }); - m.def("get_eps", []() { return libroom_eps; }); + m.def("set_eps", [](const float &eps) { + libroom_eps = eps; }); + m.def("get_eps", []() { + return libroom_eps; }); // Routines for the geometry packages m.def("ccw3p", &ccw3p, "Determines the orientation of three points"); @@ -280,7 +310,7 @@ PYBIND11_MODULE(libroom, m) { "Computes the angle between two 2D or 3D vectors"); m.def("dist_line_point", &dist_line_point, - "Computes the distance between a point and an infinite line"); + "Computes the distance between a point and an infinite line"); m.def("rir_builder", &rir_builder, "RIR builder"); m.def("delay_sum", &delay_sum, "Delay and sum"); diff --git a/pyroomacoustics/libroom_src/microphone.hpp b/pyroomacoustics/libroom_src/microphone.hpp index b62d71f7..89452775 100644 --- a/pyroomacoustics/libroom_src/microphone.hpp +++ b/pyroomacoustics/libroom_src/microphone.hpp @@ -1,4 +1,4 @@ -/* +/* * Microphone and receiver classes * Copyright (C) 2019 Robin Scheibler * @@ -9,8 +9,8 @@ * copies of the Software, and to permit persons to whom the Software is * furnished to do so, subject to the following conditions: * - * The above copyright notice and this permission notice shall be included in all - * copies or substantial portions of the Software. + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, @@ -20,95 +20,126 @@ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * SOFTWARE. * - * You should have received a copy of the MIT License along with this program. If - * not, see . + * You should have received a copy of the MIT License along with this program. + * If not, see . */ #ifndef __MICROPHONE_HPP__ #define __MICROPHONE_HPP__ -#include -#include +#include #include #include -#include +#include +#include +#include #include "common.hpp" -template -class Microphone -{ +template +class Microphone { /* - * This is the basic microphone class. It works as an omnidirectional microphone. + * This is the basic microphone class. It works as an omnidirectional + * microphone. */ - public: - Vectorf loc; - - int n_dirs = 1; - int n_bands = 1; // the number of frequency bands in the histogram - float hist_resolution; // the size of one bin in meters - std::vector distance_bins = { 0.f }; // a list of distances forming the boundaries of the bins in the time histogram - - // We keep a log of discrete hits - std::list hits; - - // and an Energy histogram for the tail - std::vector histograms; - - Microphone(const Vectorf &_loc, int _n_bands, float _hist_res, float max_dist_init) - : loc(_loc), n_dirs(1), n_bands(_n_bands), hist_resolution(_hist_res) - { - size_t n_dist_bins_init = size_t(max_dist_init / hist_resolution) + 1; - // Initialize the histograms - histograms.resize(n_dirs); - for (auto &hist : histograms) - hist.init(n_bands, n_dist_bins_init); - } - ~Microphone() {}; + private: + using my_kd_tree_t = + nanoflann::KDTreeEigenMatrixAdaptor; + // the kd-tree + my_kd_tree_t *kdtree = nullptr; - void reset() - { - for (auto h = histograms.begin() ; h != histograms.end() ; ++h) - h->reset(); - } + MatrixXf directions; + size_t n_dist_bins_init = 1; + + public: + Vectorf loc; + + int n_dirs = 1; + int n_bands = 1; // the number of frequency bands in the histogram + float hist_resolution; // the size of one bin in meters + std::vector distance_bins = { + 0.f}; // a list of distances forming the boundaries of the bins in the + // time histogram + + // and an Energy histogram for the tail + std::vector histograms; - const Vectorf &get_loc() const - { - return loc; - }; + // Constructor for omni microphones (kd-tree not necessary) + Microphone(const Vectorf &_loc, int _n_bands, float _hist_res, + float max_dist_init) + : loc(_loc), n_dirs(1), n_bands(_n_bands), hist_resolution(_hist_res) { + n_dist_bins_init = size_t(max_dist_init / hist_resolution) + 1; + // Initialize the histograms + histograms.resize(n_dirs); + for (auto &hist : histograms) hist.init(n_bands, n_dist_bins_init); + } - float get_dir_gain(const Vectorf &origin, int band_index) const - { - return 1.; // omnidirectional + Microphone(const Vectorf &_loc, int _n_bands, float _hist_res, + float max_dist_init, const MatrixXf &directions) + : Microphone(_loc, _n_bands, _hist_res, max_dist_init) { + set_directions(directions); + } + + ~Microphone() { delete_kdtree(); }; + + void set_directions(const MatrixXf &_directions) { + delete_kdtree(); + directions = _directions; + kdtree = new my_kd_tree_t(D, std::cref(directions), 10 /* max leaf */); + n_dirs = directions.rows(); + + // Initialize the histograms + histograms.resize(n_dirs); + for (auto &hist : histograms) hist.init(n_bands, n_dist_bins_init); + } + + void make_omni() { + delete_kdtree(); + n_dirs = 1; + directions.setZero(1, D); + } + + void delete_kdtree() { + if (kdtree) { + delete kdtree; + kdtree = nullptr; } - - float get_dir_bin(const Vectorf &origin) const - { + } + + void reset() { + for (auto h = histograms.begin(); h != histograms.end(); ++h) h->reset(); + } + + const Vectorf &get_loc() const { return loc; }; + + int get_dir_bin(const Vectorf &origin) const { + if (n_dirs == 1) { return 0; // only one direction is logged (omni) - } + } else { + // direction of incoming ray + auto dir = (origin - loc).normalized(); - void log_hit(const Hit &the_hit, const Vectorf &origin) - { - Hit copy_hit(the_hit); + // find in kd-tree + MatrixXf::Index ret_index = 0; + float out_dist_sqr = 0.0f; + float *ptr = dir.data(); + kdtree->query(ptr, 1, &ret_index, &out_dist_sqr); - // Correct transmitted amplitude with directivity - for (int f(0) ; f < n_bands ; f++) - copy_hit.transmitted[f] *= get_dir_gain(origin, f); - - hits.push_back(copy_hit); + return int(ret_index); } + } - void log_histogram(float distance, const Eigen::ArrayXf &energy, const Vectorf &origin) - { - // first find the bin index - auto dist_bin_index = size_t(distance / hist_resolution); - auto dir_index = get_dir_bin(origin); - histograms[dir_index].log_col(dist_bin_index, energy); - } + void log_histogram(float distance, const Eigen::ArrayXf &energy, + const Vectorf &origin) { + // first find the bin index + auto dist_bin_index = size_t(distance / hist_resolution); + auto dir_index = get_dir_bin(origin); + histograms[dir_index].log_col(dist_bin_index, energy); + } - void log_histogram(const Hit &the_hit, const Vectorf &origin) - { - log_histogram(the_hit.distance, the_hit.transmitted, origin); - } + void log_histogram(const Hit &the_hit, const Vectorf &origin) { + log_histogram(the_hit.distance, the_hit.transmitted, origin); + } }; -#endif // __MICROPHONE_HPP__ +#endif // __MICROPHONE_HPP__ diff --git a/pyroomacoustics/libroom_src/room.cpp b/pyroomacoustics/libroom_src/room.cpp index 7e99c7ff..2218b36e 100644 --- a/pyroomacoustics/libroom_src/room.cpp +++ b/pyroomacoustics/libroom_src/room.cpp @@ -1,4 +1,4 @@ -/* +/* * Implementation of the Room class * Copyright (C) 2019 Robin Scheibler, Cyril Cadoux * @@ -9,8 +9,8 @@ * copies of the Software, and to permit persons to whom the Software is * furnished to do so, subject to the following conditions: * - * The above copyright notice and this permission notice shall be included in all - * copies or substantial portions of the Software. + * The above copyright notice and this permission notice shall be included in + * all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, @@ -20,15 +20,16 @@ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE * SOFTWARE. * - * You should have received a copy of the MIT License along with this program. If - * not, see . + * You should have received a copy of the MIT License along with this program. + * If not, see . */ -#include -#include -#include #include "room.hpp" +#include +#include +#include + const double pi = 3.14159265358979323846; const double pi_2 = 1.57079632679489661923; @@ -47,71 +48,71 @@ size_t number_image_sources_3(size_t max_order) { return 1 + 2 * max_order * (2 * max_order_sq + 3 * max_order + 4) / 3; } -template -Room::Room( - const std::vector> &_walls, - const std::vector &_obstructing_walls, - const std::vector> &_microphones, - float _sound_speed, - // parameters for the image source model - int _ism_order, - // parameters for the ray tracing - 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), - 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) -{ +template +Room::Room(const std::vector> &_walls, + const std::vector &_obstructing_walls, + const std::vector> &_microphones, + float _sound_speed, + // parameters for the image source model + int _ism_order, + // parameters for the ray tracing + 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), + 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(); } -template -Room::Room( - const Vectorf &_room_size, - const Eigen::Array &_absorption, - const Eigen::Array &_scattering, - const std::vector> &_microphones, - float _sound_speed, - // parameters for the image source model - int _ism_order, - // parameters for the ray tracing - float _energy_thres, - float _time_thres, - float _mic_radius, - float _mic_hist_res, - bool _is_hybrid_sim - ) - : microphones(_microphones), - 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(true), shoebox_size(_room_size), shoebox_absorption(_absorption), - shoebox_scattering(_scattering) -{ - if (shoebox_absorption.rows() != _scattering.rows()) - { - throw std::runtime_error("Error: The same number of absorption and scattering coefficients are required"); +template +Room::Room(const Vectorf &_room_size, + const Eigen::Array &_absorption, + const Eigen::Array &_scattering, + const std::vector> &_microphones, + float _sound_speed, + // parameters for the image source model + int _ism_order, + // parameters for the ray tracing + float _energy_thres, float _time_thres, float _mic_radius, + float _mic_hist_res, bool _is_hybrid_sim) + : microphones(_microphones), + 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(true), + shoebox_size(_room_size), + shoebox_absorption(_absorption), + shoebox_scattering(_scattering) { + if (shoebox_absorption.rows() != _scattering.rows()) { + throw std::runtime_error( + "Error: The same number of absorption and scattering coefficients are " + "required"); } make_shoebox_walls(shoebox_size, _absorption, _scattering); init(); } - -template<> +template <> void Room<2>::make_shoebox_walls( const Vectorf<2> &rs, // room_size - const Eigen::Array &abs, - const Eigen::Array &scat - ) -{ - Eigen::Matrix corners; + const Eigen::Array &abs, + const Eigen::Array &scat) { + Eigen::Matrix corners; corners.resize(2, 2); corners << 0.f, rs[0], 0.f, 0.f; @@ -127,98 +128,76 @@ void Room<2>::make_shoebox_walls( walls.push_back(Wall<2>(corners, abs.col(0), scat.col(0), "west")); } - -template<> +template <> void Room<3>::make_shoebox_walls( const Vectorf<3> &rs, // room_size - const Eigen::Array &abs, - const Eigen::Array &scat - ) -{ - Eigen::Matrix corners; - corners.resize(3,4); - - corners << 0.f, 0.f, 0.f, 0.f, - rs[1], 0.f, 0.f, rs[1], - 0.f, 0.f, rs[2], rs[2]; + const Eigen::Array &abs, + const Eigen::Array &scat) { + Eigen::Matrix corners; + corners.resize(3, 4); + + corners << 0.f, 0.f, 0.f, 0.f, rs[1], 0.f, 0.f, rs[1], 0.f, 0.f, rs[2], rs[2]; walls.push_back(Wall<3>(corners, abs.col(0), scat.col(0), "west")); - corners << rs[0], rs[0], rs[0], rs[0], - 0.f, rs[1], rs[1], 0.f, - 0.f, 0.f, rs[2], rs[2]; + corners << rs[0], rs[0], rs[0], rs[0], 0.f, rs[1], rs[1], 0.f, 0.f, 0.f, + rs[2], rs[2]; walls.push_back(Wall<3>(corners, abs.col(1), scat.col(1), "east")); - corners << 0.f, rs[0], rs[0], 0.f, - 0.f, 0.f, 0.f, 0.f, - 0.f, 0.f, rs[2], rs[2]; + corners << 0.f, rs[0], rs[0], 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, rs[2], rs[2]; walls.push_back(Wall<3>(corners, abs.col(2), scat.col(2), "south")); - corners << rs[0], 0.f, 0.f, rs[0], - rs[1], rs[1], rs[1], rs[1], - 0.f, 0.f, rs[2], rs[2]; + corners << rs[0], 0.f, 0.f, rs[0], rs[1], rs[1], rs[1], rs[1], 0.f, 0.f, + rs[2], rs[2]; walls.push_back(Wall<3>(corners, abs.col(3), scat.col(3), "north")); - corners << rs[0], 0.f, 0.f, rs[0], - 0.f, 0.f, rs[1], rs[1], - 0.f, 0.f, 0.f, 0.f; + corners << rs[0], 0.f, 0.f, rs[0], 0.f, 0.f, rs[1], rs[1], 0.f, 0.f, 0.f, 0.f; walls.push_back(Wall<3>(corners, abs.col(4), scat.col(4), "floor")); - corners << rs[0], rs[0], 0.f, 0.f, - 0.f, rs[1], rs[1], 0.f, - rs[2], rs[2], rs[2], rs[2]; + corners << rs[0], rs[0], 0.f, 0.f, 0.f, rs[1], rs[1], 0.f, rs[2], rs[2], + rs[2], rs[2]; walls.push_back(Wall<3>(corners, abs.col(5), scat.col(5), "ceiling")); } - -template -void Room::init() -{ +template +void Room::init() { /* * Constructor for non-shoebox rooms */ - if (walls.size() > D) - { + if (walls.size() > D) { n_bands = walls[0].get_n_bands(); for (auto &wall : walls) - if (n_bands != wall.get_n_bands()) - { - throw std::runtime_error("Error: All walls should have the same number of frequency bands"); + if (n_bands != wall.get_n_bands()) { + throw std::runtime_error( + "Error: All walls should have the same number of frequency bands"); } - } - else - { + } else { if (D == 2) throw std::runtime_error("Error: The minimum number of walls is 3"); else if (D == 3) throw std::runtime_error("Error: The minimum number of walls is 4"); else - throw std::runtime_error("Rooms of dimension other than 2 and 3 not supported"); + throw std::runtime_error( + "Rooms of dimension other than 2 and 3 not supported"); } // Useful for ray tracing max_dist = get_max_distance(); } - -template -int Room::image_source_model(const Vectorf &source_location) -{ +template +int Room::image_source_model(const Vectorf &source_location) { /* * This is the top-level method to run the image source model */ // make sure the list is empty - while (visible_sources.size() > 0) - visible_sources.pop(); + while (visible_sources.size() > 0) visible_sources.pop(); - if (is_shoebox) - { + if (is_shoebox) { return image_source_shoebox(source_location); - } - else - { + } else { // add the original (real) source - ImageSource real_source(source_location, n_bands); + ImageSource real_source(source_location, n_bands, microphones.size()); // Run the image source model algorithm image_sources_dfs(real_source, ism_order); @@ -228,24 +207,21 @@ int Room::image_source_model(const Vectorf &source_location) } } - -template -int Room::fill_sources() -{ +template +int Room::fill_sources() { int n_sources = visible_sources.size(); // Create linear arrays to store the image sources - if (n_sources > 0) - { + if (n_sources > 0) { // resize all the arrays sources.resize(D, n_sources); orders.resize(n_sources); gen_walls.resize(n_sources); attenuations.resize(n_bands, n_sources); visible_mics.resize(microphones.size(), n_sources); + source_directions.resize(D, microphones.size() * n_sources); - for (int i = n_sources - 1 ; i >= 0 ; i--) - { + for (int i = n_sources - 1; i >= 0; i--) { ImageSource &top = visible_sources.top(); // sample top of stack // fill the arrays @@ -255,6 +231,11 @@ int Room::fill_sources() attenuations.col(i) = top.attenuation; visible_mics.col(i) = top.visible_mics; + int m = 0; + for (auto mic = microphones.begin(); mic != microphones.end(); ++mic, ++m) + source_directions.col(i * microphones.size() + m) = + top.source_direction[m]; + visible_sources.pop(); // unstack } } @@ -262,84 +243,75 @@ int Room::fill_sources() return n_sources; } - -template -void Room::image_sources_dfs(ImageSource &is, int max_order) -{ +template +void Room::image_sources_dfs(ImageSource &is, int max_order) { /* * This function runs a depth first search (DFS) on the tree of image sources */ - ImageSource new_is(n_bands); + ImageSource new_is(n_bands, microphones.size()); // Check the visibility of the source from the different microphones bool any_visible = false; int m = 0; - for (auto mic = microphones.begin() ; mic != microphones.end() ; ++mic, ++m) - { - bool is_visible = is_visible_dfs(mic->get_loc(), is); - if (is_visible && !any_visible) - { - any_visible = is_visible; - is.visible_mics.resize(microphones.size()); - is.visible_mics.setZero(); - } - if (any_visible) - is.visible_mics.coeffRef(m) = is_visible; + for (auto mic = microphones.begin(); mic != microphones.end(); ++mic, ++m) { + // is_visible_dfs checks that the source is visible from the microphone + // and also captures the direction vector from the image source towards + // the microphone + bool is_visible = + is_visible_dfs(mic->get_loc(), is, is.source_direction[m]); + is.visible_mics.coeffRef(m) = is_visible; + any_visible |= is_visible; } if (any_visible) visible_sources.push(is); // this should push a copy onto the stack - + // If we reached maximal depth, stop - if (max_order == 0) - return; + if (max_order == 0) return; // Then, check all the reflections across the walls - for (size_t wi=0 ; wi < walls.size() ; wi++) - { + for (size_t wi = 0; wi < walls.size(); wi++) { int dir = walls[wi].reflect(is.loc, new_is.loc); // the reflected location - // We only check valid reflections (normals should point outward from the room - if (dir <= 0) - continue; + // We only check valid reflections (normals should point outward from the + // room + if (dir <= 0) continue; // The reflection is valid, fill in the image source attributes new_is.attenuation = is.attenuation * walls[wi].get_transmission(); - if (walls[wi].scatter.maxCoeff() > 0.f && is_hybrid_sim) - { + if (walls[wi].scatter.maxCoeff() > 0.f && is_hybrid_sim) { new_is.attenuation *= (1 - walls[wi].scatter).sqrt(); } new_is.order = is.order + 1; new_is.gen_wall = wi; new_is.parent = &is; - // Run the DFS recursion (on the last element in the array, the one we just added) + // Run the DFS recursion (on the last element in the array, the one we just + // added) image_sources_dfs(new_is, max_order - 1); } } - -template -bool Room::is_visible_dfs(const Vectorf &p, ImageSource &is) -{ +template +bool Room::is_visible_dfs(const Vectorf &p, ImageSource &is, + Vectorf &source_direction) { /* - Returns true if the given sound source (with image source id) is visible from point p. + Returns true if the given sound source (with image source id) is visible + from point p. room - the structure that contains all the sources and stuff p (np.array size 2 or 3) coordinates of the point where we check visibility imageId (int) id of the image within the SoundSource object - + Returns False (0) : not visible True (1) : visible */ - if (is_obstructed_dfs(p, is)) - return false; + if (is_obstructed_dfs(p, is)) return false; - if (is.parent != NULL) - { + if (is.parent != NULL) { Vectorf intersection; // get generating wall id @@ -352,25 +324,28 @@ bool Room::is_visible_dfs(const Vectorf &p, ImageSource &is) // the generating wall if (ret >= 0) // Check visibility of intersection point from parent source - return is_visible_dfs(intersection, *(is.parent)); + return is_visible_dfs(intersection, *(is.parent), source_direction); else return false; } + // the source direction is the line between the original source and p + // p is the intersection of the ray with the first reflection wall + source_direction = (p - is.loc).normalized(); + // If we get here this is the original, unobstructed, source return true; } - -template -bool Room::is_obstructed_dfs(const Vectorf &p, ImageSource &is) -{ +template +bool Room::is_obstructed_dfs(const Vectorf &p, ImageSource &is) { /* - Checks if there is a wall obstructing the line of sight going from a source to a point. + Checks if there is a wall obstructing the line of sight going from a source + to a point. room - the structure that contains all the sources and stuff - p (np.array size 2 or 3) coordinates of the point where we check obstruction - imageId (int) id of the image within the SoundSource object + p (np.array size 2 or 3) coordinates of the point where we check + obstruction imageId (int) id of the image within the SoundSource object Returns (bool) False (0) : not obstructed @@ -379,23 +354,19 @@ bool Room::is_obstructed_dfs(const Vectorf &p, ImageSource &is) int gen_wall_id = is.gen_wall; // Check candidate walls for obstructions - for (size_t ow = 0 ; ow < obstructing_walls.size() ; ow++) - { + for (size_t ow = 0; ow < obstructing_walls.size(); ow++) { int wall_id = obstructing_walls[ow]; // generating wall can't be obstructive - if (wall_id != gen_wall_id) - { + if (wall_id != gen_wall_id) { Vectorf intersection; int ret = walls[wall_id].intersection(is.loc, p, intersection); // There is an intersection and it is distinct from segment endpoints - if (ret == Wall::Isect::VALID || ret == Wall::Isect::BNDRY) - { - if (is.parent != NULL) - { + if (ret == Wall::Isect::VALID || ret == Wall::Isect::BNDRY) { + if (is.parent != NULL) { // Test if the intersection point and the image are on - // opposite sides of the generating wall + // opposite sides of the generating wall // We ignore the obstruction if it is inside the // generating wall (it is what happens in a corner) int img_side = walls[is.gen_wall].side(is.loc); @@ -403,8 +374,7 @@ bool Room::is_obstructed_dfs(const Vectorf &p, ImageSource &is) if (img_side != intersection_side && intersection_side != 0) return true; - } - else + } else return true; } } @@ -413,91 +383,98 @@ bool Room::is_obstructed_dfs(const Vectorf &p, ImageSource &is) return false; } - -template -int Room::image_source_shoebox(const Vectorf &source) -{ +template +int Room::image_source_shoebox(const Vectorf &source) { // precompute powers of the transmission coefficients std::vector transmission_pwr; - for (int i(0) ; i <= ism_order ; ++i) - transmission_pwr.push_back(Eigen::ArrayXXf(n_bands, 2*D)); + for (int i(0); i <= ism_order; ++i) + transmission_pwr.push_back(Eigen::ArrayXXf(n_bands, 2 * D)); transmission_pwr[0].setOnes(); - if (ism_order > 0) - { + if (ism_order > 0) { transmission_pwr[1] = (1.f - shoebox_absorption).sqrt(); - if (shoebox_scattering.maxCoeff() > 0.f && is_hybrid_sim) - { + if (shoebox_scattering.maxCoeff() > 0.f && is_hybrid_sim) { transmission_pwr[1] *= (1 - shoebox_scattering).sqrt(); } } - for (int i = 2 ; i <= ism_order ; ++i) - transmission_pwr[i] = transmission_pwr[i-1] * transmission_pwr[1]; + for (int i = 2; i <= ism_order; ++i) + transmission_pwr[i] = transmission_pwr[i - 1] * transmission_pwr[1]; // make sure the list is empty int n_image_sources = number_image_sources_3(ism_order); - if (D == 2) - n_image_sources = number_image_sources_2(ism_order); + if (D == 2) n_image_sources = number_image_sources_2(ism_order); - std::vector> image_sources(n_image_sources, ImageSource(n_bands)); + std::vector> image_sources( + n_image_sources, ImageSource(n_bands, microphones.size())); int img_src_index = 0; - + // L1 ball of room images int point[3] = {0, 0, 0}; // Take 2D case into account int z_max = ism_order; - if (D == 2) - z_max = 0; + if (D == 2) z_max = 0; // Walk on all the points of the discrete L! ball of radius ism_order - for (point[2] = -z_max ; point[2] <= z_max ; point[2]++) - { + for (point[2] = -z_max; point[2] <= z_max; point[2]++) { int y_max = ism_order - abs(point[2]); - for (point[1] = -y_max ; point[1] <= y_max ; point[1]++) - { + for (point[1] = -y_max; point[1] <= y_max; point[1]++) { int x_max = y_max - abs(point[1]); if (x_max < 0) x_max = 0; - for (point[0] = -x_max ; point[0] <= x_max ; point[0]++) - { + for (point[0] = -x_max; point[0] <= x_max; point[0]++) { ImageSource &is = image_sources[img_src_index++]; - is.visible_mics.resize(microphones.size()); + // is.visible_mics.resize(microphones.size()); is.visible_mics.setOnes(); // everything is visible - // Now compute the reflection, the order, and the multiplicative constant - for (size_t d = 0 ; d < D ; d++) - { + // Now compute the reflection, the order, and the multiplicative + // constant + for (size_t d = 0; d < D; d++) { // Compute the reflected source - float step = abs(point[d]) % 2 == 1 ? shoebox_size.coeff(d) - source.coeff(d) : source.coeff(d); + int abs_point = abs(point[d]); + float step = abs_point % 2 == 1 + ? shoebox_size.coeff(d) - source.coeff(d) + : source.coeff(d); + is.loc[d] = point[d] * shoebox_size.coeff(d) + step; + + // compute the source direction is.order_xyz[d] = point[d]; + int m = 0; + for (auto mic = microphones.begin(); mic != microphones.end(); + ++mic, ++m) { + if (abs_point % 2 == 1) + is.source_direction[m].coeffRef(d) = + is.loc[d] - mic->get_loc().coeff(d); + else + is.source_direction[m].coeffRef(d) = + mic->get_loc().coeff(d) - is.loc[d]; + } - // source order is just the sum of absolute values of reflection indices - is.order += abs(point[d]); + // source order is just the sum of absolute values of reflection + // indices + is.order += abs_point; // attenuation can also be computed this way int p1 = 0, p2 = 0; - if (point[d] > 0) - { - p1 = point[d]/2; - p2 = (point[d]+1)/2; - } - else if (point[d] < 0) - { - p1 = abs((point[d]-1)/2); - p2 = abs(point[d]/2); + if (point[d] > 0) { + p1 = point[d] / 2; + p2 = (point[d] + 1) / 2; + } else if (point[d] < 0) { + p1 = abs((point[d] - 1) / 2); + p2 = abs(point[d] / 2); } - is.attenuation *= transmission_pwr[p1].col(2*d); // 'west' absorption factor - is.attenuation *= transmission_pwr[p2].col(2*d+1); // 'east' absorption factor + is.attenuation *= + transmission_pwr[p1].col(2 * d); // 'west' absorption factor + is.attenuation *= + transmission_pwr[p2].col(2 * d + 1); // 'east' absorption factor } } } } // Create linear arrays to store the image sources - if (n_image_sources > 0) - { + if (n_image_sources > 0) { // resize all the arrays sources.resize(D, n_image_sources); orders.resize(n_image_sources); @@ -505,8 +482,10 @@ int Room::image_source_shoebox(const Vectorf &source) gen_walls.resize(n_image_sources); attenuations.resize(n_bands, n_image_sources); visible_mics.resize(microphones.size(), n_image_sources); + source_directions.resize(D, microphones.size() * n_image_sources); - for (std::vector::size_type idx = 0; idx < image_sources.size(); ++idx) { + for (std::vector::size_type idx = 0; idx < image_sources.size(); + ++idx) { ImageSource &top = image_sources[idx]; // sample top of stack // fill the arrays @@ -516,6 +495,11 @@ int Room::image_source_shoebox(const Vectorf &source) orders_xyz.col(idx) = top.order_xyz; attenuations.col(idx) = top.attenuation; visible_mics.col(idx) = top.visible_mics; + + int m = 0; + for (auto mic = microphones.begin(); mic != microphones.end(); ++mic, ++m) + source_directions.col(idx * microphones.size() + m) = + top.source_direction[m].normalized(); } } @@ -523,16 +507,13 @@ int Room::image_source_shoebox(const Vectorf &source) return n_image_sources; } - -template -float Room::get_max_distance() -{ - +template +float Room::get_max_distance() { /* This function outputs a value L that is strictly larger than any distance that a ray can travel in straight line in the room without hitting anything. - In other words, this function outputs the diagonal of the bounding box + In other words, this function outputs the diagonal of the bounding box of this room. - + As we are sure that a ray reflecting from a hit_point H will hit another wall in less than L meters, we can use Wall::intersection() with the segment starting at H and of length L @@ -541,30 +522,26 @@ float Room::get_max_distance() float maxX(0), maxY(0), maxZ(0); float minX(0), minY(0), minZ(0); size_t n_walls = walls.size(); - + // The first step is to go over the corners of all the walls to extract // the min x, y and z and the max x,y and z. - for (size_t i(0); i < n_walls; ++i) - { - Wall wi = this -> get_wall(i); + for (size_t i(0); i < n_walls; ++i) { + Wall wi = this->get_wall(i); Eigen::Vector3f max_coord(0, 0, 0); Eigen::Vector3f min_coord(0, 0, 0); - if (D == 2) - { + if (D == 2) { max_coord.head(2) = wi.corners.topRows(D).rowwise().maxCoeff(); min_coord.head(2) = wi.corners.topRows(D).rowwise().minCoeff(); - } else - { + } else { max_coord = wi.corners.topRows(D).rowwise().maxCoeff(); min_coord = wi.corners.topRows(D).rowwise().minCoeff(); } - // For the first wall, we have nothing to compare to - if (i == 0) - { + // For the first wall, we have nothing to compare to + if (i == 0) { maxX = max_coord[0]; minX = min_coord[0]; maxY = max_coord[1]; @@ -574,24 +551,16 @@ float Room::get_max_distance() } // For the next walls, we compare to the stored min and max x,y and z - else - { - + else { // Update max point - if (max_coord[0] > maxX) - maxX = max_coord[0]; - if (max_coord[1] > maxY) - maxY = max_coord[1]; - if (max_coord[2] > maxZ) - maxZ = max_coord[2]; + if (max_coord[0] > maxX) maxX = max_coord[0]; + if (max_coord[1] > maxY) maxY = max_coord[1]; + if (max_coord[2] > maxZ) maxZ = max_coord[2]; // Update min point - if (min_coord[0] < minX) - minX = min_coord[0]; - if (min_coord[1] < minY) - minY = min_coord[1]; - if (min_coord[2] < minZ) - minZ = min_coord[2]; + if (min_coord[0] < minX) minX = min_coord[0]; + if (min_coord[1] < minY) minY = min_coord[1]; + if (min_coord[2] < minZ) minZ = min_coord[2]; } } @@ -602,62 +571,56 @@ float Room::get_max_distance() // Return the length of the diagonal of the bounding box, + 1 return (min_point - max_point).norm() + 1; - } - -template -std::tuple < Vectorf, int, float > Room::next_wall_hit( - const Vectorf &start, - const Vectorf &end, - bool scattered_ray) -{ - +template +std::tuple, int, float> Room::next_wall_hit( + const Vectorf &start, const Vectorf &end, bool scattered_ray) { /* This function is called in 2 different contexts : - * + * * 1) When we trace the main ray : it computes the next next wall_hit * position given a segment defined by its endpoints.It also returns the * index of the wall that contains this next hit point. - * + * * start: (array size 2 or 3) the start point of the segment. Except for - * the very first ray where this point is the source, this point + * the very first ray where this point is the source, this point * is located on a wall. - * end: (array size 2 or 3) the end point of the segment. Recall that + * end: (array size 2 or 3) the end point of the segment. Recall that * this point must be set so that a wall is intersected between * start and end (dist(start,end) = Room->max_dist). * scattered_ray: false - * + * * :return: tuple - * + * * ================================================================== - * + * * 2) When we trace a scattered ray from the previous wall_hit point * back to the microphone : it checks if an (obstructing) wall stands * between that previous wall_hit position and the center of the microphone. * In that case, the scattered ray cannot reach the microphone. - * + * * start: (array size 2 or 3) the last wall_hit position - * end: (array size 2 or 3) the end point of the segment, ie the center + * end: (array size 2 or 3) the end point of the segment, ie the center * of the microphone in this case. * scattered_ray: true - * - * :return: tuple - * If no wall is intersected, then potentially_obstructing_wall_index will be -1 - * + * + * :return: tuple If no wall is intersected, then + * potentially_obstructing_wall_index will be -1 + * * In fact here we are only interested in the second element of the tuple. * */ - const static std::vector> shoebox_orders = {{0, 1, 2}, {1, 0, 2}, {2, 0, 1}}; + const static std::vector> shoebox_orders = { + {0, 1, 2}, {1, 0, 2}, {2, 0, 1}}; Vectorf result; int next_wall_index = -1; float hit_dist(max_dist); - if (is_shoebox) - { + if (is_shoebox) { // There are no obstructing walls in shoebox rooms - if (scattered_ray) - return std::make_tuple(result, -1, 0.); + if (scattered_ray) return std::make_tuple(result, -1, 0.); // The direction vector Vectorf dir = end - start; @@ -667,8 +630,7 @@ std::tuple < Vectorf, int, float > Room::next_wall_hit( auto d = shoebox_orders[shoebox_orders_idx]; float abs_dir0 = std::abs(dir[d[0]]); - if (abs_dir0 < libroom_eps) - continue; + if (abs_dir0 < libroom_eps) continue; // distance to plane float distance = 0.; @@ -676,30 +638,26 @@ std::tuple < Vectorf, int, float > Room::next_wall_hit( // this wil tell us if the front or back plane is hit int ind_inc = 0; - if (dir[d[0]] < 0) - { + if (dir[d[0]] < 0) { result[d[0]] = 0.; distance = start[d[0]]; ind_inc = 0; - } - else - { + } else { result[d[0]] = shoebox_size[d[0]]; distance = shoebox_size[d[0]] - start[d[0]]; ind_inc = 1; } - if (distance < libroom_eps) - continue; + if (distance < libroom_eps) continue; float ratio = distance / abs_dir0; // Now compute the intersection point and verify if intersection happens - for (size_t i = 1 ; i < D ; ++i) - { + for (size_t i = 1; i < D; ++i) { result[d[i]] = start[d[i]] + ratio * dir[d[i]]; // when there is no intersection, we jump to the next plane - if (result[d[i]] <= -libroom_eps || shoebox_size[d[i]] + libroom_eps <= result[d[i]]) + if (result[d[i]] <= -libroom_eps || + shoebox_size[d[i]] + libroom_eps <= result[d[i]]) goto next_plane; } @@ -710,12 +668,10 @@ std::tuple < Vectorf, int, float > Room::next_wall_hit( break; -next_plane: + next_plane: (void)0; // no op } - } - else - { + } else { // 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 @@ -724,9 +680,8 @@ std::tuple < Vectorf, int, float > Room::next_wall_hit( // For a scattered ray, we only check the obstructing walls size_t n_walls = scattered_ray ? obstructing_walls.size() : walls.size(); - for (size_t i(0) ; i < n_walls ; ++i) - { - Wall & w = scattered_ray ? walls[obstructing_walls[i]] : walls[i]; + for (size_t i(0); i < n_walls; ++i) { + Wall &w = scattered_ray ? walls[obstructing_walls[i]] : walls[i]; // To store the result of this iteration Vectorf temp_hit = Vectorf::Zero(); @@ -734,40 +689,29 @@ std::tuple < Vectorf, int, float > Room::next_wall_hit( // As a side effect, temp_hit gets a value (VectorXf) here int ret = w.intersection(start, end, temp_hit); - if (ret > -1) - { + if (ret > -1) { float temp_dist = (temp_hit - start).norm(); // Compare to min dist to see if this wall is the closest to 'start' // Compare to libroom_eps to be sure that this wall w is not the wall // where 'start' is located ('intersects' could be true because of // rounding errors) - if (temp_dist > libroom_eps && temp_dist < hit_dist) - { + if (temp_dist > libroom_eps && temp_dist < hit_dist) { hit_dist = temp_dist; result = temp_hit; next_wall_index = i; } } } - } return std::make_tuple(result, next_wall_index, hit_dist); - } - -template -bool Room::scat_ray( - const Eigen::ArrayXf &transmitted, - const Wall &wall, - const Vectorf &prev_last_hit, - const Vectorf &hit_point, - float travel_dist - ) -{ - +template +bool Room::scat_ray(const Eigen::ArrayXf &transmitted, const Wall &wall, + const Vectorf &prev_last_hit, + const Vectorf &hit_point, float travel_dist) { /* Traces a one-hop scattered ray from the last wall hit to each microphone. In case the scattering ray can indeed reach the microphone (no wall in @@ -776,29 +720,27 @@ bool Room::scat_ray( float energy: The energy of the ray right after last_wall has absorbed a part of it wall: The wall object where last_hit is located - prev_last_hit: (array size 2 or 3) the previous last wall hit_point position (needed to check that - the wall normal is correctly oriented) - hit_point: (array size 2 or 3) defines the last wall hit position - travel_dist: The total distance travelled by the ray from source to hit_point + prev_last_hit: (array size 2 or 3) the previous last wall hit_point position + (needed to check that the wall normal is correctly oriented) hit_point: (array + size 2 or 3) defines the last wall hit position travel_dist: The total + distance travelled by the ray from source to hit_point :return : true if the scattered ray reached ALL the microphones, false otw */ - // Convert the energy threshold to transmission threshold (make this more efficient at some point) + // Convert the energy threshold to transmission threshold (make this more + // efficient at some point) float distance_thres = time_thres * sound_speed; - bool ret = true; - for(size_t k(0); k < microphones.size(); ++k) - { - + bool ret = true; + for (size_t k(0); k < microphones.size(); ++k) { Vectorf mic_pos = microphones[k].get_loc(); - /* + /* * We also need to check that both the microphone and the * previous hit point are on the same side of the wall */ - if (wall.side(mic_pos) != wall.side(prev_last_hit)) - { + if (wall.side(mic_pos) != wall.side(prev_last_hit)) { ret = false; continue; } @@ -809,11 +751,11 @@ bool Room::scat_ray( float hit_distance(0.); if (!is_shoebox) - std::tie(dont_care, next_wall_index, hit_distance) = next_wall_hit(hit_point, mic_pos, true); + std::tie(dont_care, next_wall_index, hit_distance) = + next_wall_hit(hit_point, mic_pos, true); // If no wall obstructs the scattered ray - if (next_wall_index == -1) - { + if (next_wall_index == -1) { // As the ray is shot towards the microphone center, // the hop dist can be easily computed Vectorf hit_point_to_mic = mic_pos - hit_point; @@ -837,15 +779,13 @@ bool Room::scat_ray( // Ref: D. Schroeder, "Physically based real-time auralization of interactive virtual environments", // section 5.4, eq. 5.54. double r_sq = double(travel_dist_at_mic) * travel_dist_at_mic; - auto p_hit = (1 - sqrt(1 - mic_radius_sq / std::max(mic_radius_sq, r_sq))); - Eigen::ArrayXf energy = scat_trans / (r_sq * p_hit) ; + auto p_hit = + (1 - sqrt(1 - mic_radius_sq / std::max(mic_radius_sq, r_sq))); + Eigen::ArrayXf energy = scat_trans / (r_sq * p_hit); microphones[k].log_histogram(travel_dist_at_mic, energy, hit_point); - } - else + } else ret = false; - } - else - { + } else { ret = false; // if a wall intersects the scattered ray, we return false } } @@ -853,34 +793,39 @@ bool Room::scat_ray( return ret; } -template -void Room::simul_ray( - float phi, - float theta, - const Vectorf source_pos, - float energy_0 - ) -{ +template +void Room::simul_ray(float phi, float theta, const Vectorf &source_pos, + float energy_0) { + // the direction of the ray (unit vector) + Vectorf dir; + if (D == 2) + dir.head(2) = Eigen::Vector2f(cos(phi), sin(phi)); + else if (D == 3) + dir.head(3) = Eigen::Vector3f(sin(theta) * cos(phi), sin(theta) * sin(phi), + cos(theta)); - /*This function simulates one ray and fills the output vectors of + simul_ray(dir, source_pos, energy_0); +} + +template +void Room::simul_ray(const Vectorf &ray_direction, + const Vectorf &source_pos, float energy_0) { + /*This function simulates one ray and fills the output vectors of every microphone with all the entries produced by this ray (any specular or scattered ray reaching a microphone) - - phi (azimuth) and theta (colatitude) : give the orientation of the ray (2D or 3D) - source_pos: (array size 2 or 3) is the location of the sound source (NOT AN IMAGE SOURCE) - energy_0: (float) the initial energy of one ray - output: is the std::vector that contains the entries for all the simulated rays */ + + phi (azimuth) and theta (colatitude) : give the orientation of the ray (2D or + 3D) source_pos: (array size 2 or 3) is the location of the sound source (NOT + AN IMAGE SOURCE) energy_0: (float) the initial energy of one ray output: is + the std::vector that contains the entries for all the simulated rays */ // ------------------ INIT -------------------- // What we need to trace the ray // the origin of the ray Vectorf start = source_pos; + // the direction of the ray (unit vector) - Vectorf dir; - if(D == 2) - dir.head(2) = Eigen::Vector2f(cos(phi), sin(phi)); - else if (D == 3) - dir.head(3) = Eigen::Vector3f(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta)); + Vectorf dir = ray_direction; // The following initializations are arbitrary and does not count since we set // the boolean to false @@ -890,7 +835,7 @@ void Room::simul_ray( Eigen::ArrayXf transmitted = Eigen::ArrayXf::Ones(n_bands) * energy_0; Eigen::ArrayXf energy = Eigen::ArrayXf::Ones(n_bands); float travel_dist = 0; - + // To count the number of times the ray bounces on the walls // For hybrid generation we add a ray to output only if specular_counter // is higher than the ism order. @@ -902,47 +847,44 @@ void Room::simul_ray( //--------------------------------------------- - //------------------ RAY TRACING -------------------- Vectorf hit_point; - while (true) - { + while (true) { // Find the next hit point float hit_distance(0); - std::tie(hit_point, next_wall_index, hit_distance) = next_wall_hit(start, start + dir * max_dist, false); + std::tie(hit_point, next_wall_index, hit_distance) = + next_wall_hit(start, start + dir * max_dist, false); // If no wall is hit (rounding errors), stop the ray - if (next_wall_index == -1) - break; + if (next_wall_index == -1) break; // Intersected wall Wall &wall = walls[next_wall_index]; // Check if the specular ray hits any of the microphone - if (!(is_hybrid_sim && specular_counter < ism_order)) - { - for(size_t k(0) ; k < microphones.size() ; k++) - { + if (!(is_hybrid_sim && specular_counter < ism_order)) { + for (size_t k(0); k < microphones.size(); k++) { // Compute the distance between the line defined by (start, hit_point) // and the center of the microphone (mic_pos) Vectorf to_mic = microphones[k].get_loc() - start; float impact_distance = to_mic.dot(dir); - bool impacts = -libroom_eps < impact_distance && impact_distance < hit_distance + libroom_eps; + bool impacts = -libroom_eps < impact_distance && + impact_distance < hit_distance + libroom_eps; // If yes, we compute the ray's transmitted amplitude at the mic // and we continue the ray - if (impacts && - (to_mic - dir * impact_distance).norm() < mic_radius + libroom_eps) - { + if (impacts && (to_mic - dir * impact_distance).norm() < + mic_radius + libroom_eps) { // The length of this last hop float distance = fabsf(impact_distance); // Updating travel_time and transmitted amplitude for this ray - // We DON'T want to modify the variables transmitted amplitude and travel_dist - // because the ray will continue its way + // We DON'T want to modify the variables transmitted amplitude and + // travel_dist + // because the ray will continue its way float travel_dist_at_mic = travel_dist + distance; // The (r_sq * p_hit) normalization factor is necessary to equalize the energy @@ -950,10 +892,9 @@ void Room::simul_ray( // Ref: D. Schroeder, "Physically based real-time auralization of interactive virtual environments", // section 5.4, eq. 5.54. double r_sq = double(travel_dist_at_mic) * travel_dist_at_mic; - auto p_hit = (1 - sqrt(1 - mic_radius_sq / std::max(mic_radius_sq, r_sq))); + auto p_hit = + (1 - sqrt(1 - mic_radius_sq / std::max(mic_radius_sq, r_sq))); energy = transmitted / (r_sq * p_hit); - // energy = transmitted / (travel_dist_at_mic - sqrtf(fmaxf(0.f, travel_dist_at_mic * travel_dist_at_mic - mic_radius_sq))); - microphones[k].log_histogram(travel_dist_at_mic, energy, start); } } @@ -964,16 +905,9 @@ void Room::simul_ray( transmitted *= wall.get_energy_reflection(); // Let's shoot the scattered ray induced by the rebound on the wall - if (wall.scatter.maxCoeff() > 0.f) - { + if (wall.scatter.maxCoeff() > 0.f) { // Shoot the scattered ray - scat_ray( - transmitted, - wall, - start, - hit_point, - travel_dist - ); + scat_ray(transmitted, wall, start, hit_point, travel_dist); // The overall ray's energy gets decreased by the total // amount of scattered energy @@ -981,8 +915,7 @@ void Room::simul_ray( } // Check if we reach the thresholds for this ray - if (travel_dist > distance_thres || transmitted.maxCoeff() < e_thres) - break; + if (travel_dist > distance_thres || transmitted.maxCoeff() < e_thres) break; // set up for next iteration specular_counter += 1; @@ -991,38 +924,36 @@ void Room::simul_ray( } } - -template +template void Room::ray_tracing( - const Eigen::Matrix &angles, - const Vectorf source_pos - ) -{ + const Eigen::Matrix &angles, + const Vectorf &source_pos) { // float energy_0 = 2.f / (mic_radius * mic_radius * angles.cols()); float energy_0 = 2.f / angles.cols(); - for (int k(0) ; k < angles.cols() ; k++) - { - float phi = angles.coeff(0,k); + for (int k(0); k < angles.cols(); k++) { + float phi = angles.coeff(0, k); float theta = pi_2; - if (D == 3) - theta = angles.coeff(1,k); + if (D == 3) theta = angles.coeff(1, k); simul_ray(phi, theta, source_pos, energy_0); } } - -template +template void Room::ray_tracing( - size_t nb_phis, - size_t nb_thetas, - const Vectorf source_pos - ) -{ + const Eigen::Matrix &unit_vectors, + const Vectorf &source_pos) { + float energy_0 = 2.f / unit_vectors.cols(); + for (int c = 0; c < unit_vectors.cols(); c++) + simul_ray(unit_vectors.col(c), source_pos, energy_0); +} +template +void Room::ray_tracing(size_t nb_phis, size_t nb_thetas, + const Vectorf &source_pos) { /*This method produced all the time/energy entries needed to compute the RIR using ray-tracing with the following parameters @@ -1031,26 +962,23 @@ void Room::ray_tracing( (NOTE: nb_phis*nb_thetas is the number of simulated rays source_pos: (array size 2 or 3) represents the position of the sound source - :returns: + :returns: a std::vector where each entry is a tuple (time, energy) reprensenting a ray (scattered or not) reaching the microphone */ - // ------------------ INIT -------------------- // initial energy of one ray float energy_0 = 2.f / (nb_phis * nb_thetas); // ------------------ RAY TRACING -------------------- - for (size_t i(0); i < nb_phis; ++i) - { - float phi = 2 * pi * (float) i / nb_phis; + for (size_t i(0); i < nb_phis; ++i) { + float phi = 2 * pi * (float)i / nb_phis; - for (size_t j(0); j < nb_thetas; ++j) - { + for (size_t j(0); j < nb_thetas; ++j) { // Having a 3D uniform sampling of the sphere surrounding the room - float theta = std::acos(2 * ((float) j / nb_thetas) - 1); + float theta = std::acos(2 * ((float)j / nb_thetas) - 1); // For 2D, this parameter means nothing, but we set it to // PI/2 to be consistent @@ -1063,21 +991,15 @@ void Room::ray_tracing( // if we work in 2D rooms, only 1 elevation angle is needed // => get out of the theta loop - if (D == 2) - { + if (D == 2) { break; } } } } - -template -void Room::ray_tracing( - size_t n_rays, - const Vectorf source_pos - ) -{ +template +void Room::ray_tracing(size_t n_rays, const Vectorf &source_pos) { /*This method produced all the time/energy entries needed to compute the RIR using ray-tracing with the following parameters @@ -1086,19 +1008,16 @@ void Room::ray_tracing( source_pos: (array size 2 or 3) represents the position of the sound source */ - // ------------------ INIT -------------------- // initial energy of one ray float energy_0 = 2.f / n_rays; // ------------------ RAY TRACING -------------------- - if (D == 3) - { + if (D == 3) { auto offset = 2.f / n_rays; auto increment = pi * (3.f - sqrt(5.f)); // phi increment - for (size_t i(0); i < n_rays ; ++i) - { + for (size_t i(0); i < n_rays; ++i) { auto z = (i * offset - 1) + offset / 2.f; auto rho = sqrt(1.f - z * z); @@ -1112,47 +1031,38 @@ void Room::ray_tracing( simul_ray(azimuth, colatitude, source_pos, energy_0); } - } - else if (D == 2) - { + } else if (D == 2) { float offset = 2. * pi / n_rays; - for (size_t i(0) ; i < n_rays ; ++i) + for (size_t i(0); i < n_rays; ++i) simul_ray(i * offset, 0.f, source_pos, energy_0); } } - -template -bool Room::contains(const Vectorf point) -{ - +template +bool Room::contains(const Vectorf point) { /*This methods checks if a point is contained in the room - + point: (array size 2 or 3) representing a point in the room - + :returs: true if the point is inside the room, false otherwise*/ - - + // ------- USING RAY CASTING ALGO ------- // First we need to build a point outside the room - // For this we take the min (x,y,z) point and subtract 1 (arbitrary) to each coordinate + // For this we take the min (x,y,z) point and subtract 1 (arbitrary) to each + // coordinate size_t n_walls = walls.size(); - Eigen::Matrix min_coord; + Eigen::Matrix min_coord; min_coord.setZero(); - for (size_t i(0); i < n_walls; ++i) - { + for (size_t i(0); i < n_walls; ++i) { Wall &wi = this->get_wall(i); - // First iteration - if (i == 0) - { + // First iteration + if (i == 0) { min_coord.col(0) = wi.corners.rowwise().minCoeff(); - } - else - { + } else { min_coord.col(1) = wi.corners.rowwise().minCoeff(); min_coord.col(0) = min_coord.rowwise().minCoeff(); } @@ -1162,39 +1072,34 @@ bool Room::contains(const Vectorf point) // ------------------------------------------ - // Now we build a segment between 'outside_point' and 'point' + // Now we build a segment between 'outside_point' and 'point' // We must look at the number of walls that intersect this segment size_t n_intersections(0); bool ambiguous_intersection = false; - do - { - // We restart the computation with a modified output_point as long as we find - // an ambiguous intersection (on the edge or on the endpoint of a segment) - // Note : ambiguous intersection means that Wall::intersects() method - // gives a result strictly above 0 - + do { + // We restart the computation with a modified output_point as long as we + // find an ambiguous intersection (on the edge or on the endpoint of a + // segment) Note : ambiguous intersection means that Wall::intersects() + // method gives a result strictly above 0 + n_intersections = 0; ambiguous_intersection = false; outside_point[0] -= (float)(rand() % 27) / 50; outside_point[1] -= (float)(rand() % 22) / 26; - if (D == 3) - { + if (D == 3) { outside_point[2] -= (float)(rand() % 24 / 47); } - for (size_t i(0); i < n_walls; ++i) - { - - Wall & w = walls[i]; + for (size_t i(0); i < n_walls; ++i) { + Wall &w = walls[i]; int result = w.intersects(outside_point, point); ambiguous_intersection = ambiguous_intersection || (result > 0); - if (result > -1) - { + if (result > -1) { n_intersections++; } } diff --git a/pyroomacoustics/libroom_src/room.hpp b/pyroomacoustics/libroom_src/room.hpp index d7c9398d..4de6add8 100644 --- a/pyroomacoustics/libroom_src/room.hpp +++ b/pyroomacoustics/libroom_src/room.hpp @@ -52,25 +52,29 @@ struct ImageSource // this is a unit vector from the center of the source pointing // in the direction of the path to the microphone - Vectorf source_impact_dir; + std::vector> source_direction; // This contains the reflection orders with respect to x/y/z axis // for the shoebox image source model Vectori order_xyz; - ImageSource(size_t n_bands) + ImageSource(size_t n_bands, size_t n_mics) : order(0), gen_wall(-1), parent(NULL) { loc.setZero(); attenuation.resize(n_bands); attenuation.setOnes(); + source_direction.resize(n_mics); + visible_mics.setZero(n_mics); } - ImageSource(const Vectorf &_loc, size_t n_bands) + ImageSource(const Vectorf &_loc, size_t n_bands, size_t n_mics) : loc(_loc), order(0), gen_wall(-1), parent(NULL) { attenuation.resize(n_bands); attenuation.setOnes(); + source_direction.resize(n_mics); + visible_mics.setZero(n_mics); } }; @@ -116,6 +120,7 @@ class Room Eigen::VectorXi gen_walls; Eigen::VectorXi orders; Eigen::Matrix orders_xyz; + Eigen::Matrix source_directions; Eigen::MatrixXf attenuations; // This array will get filled by visibility status @@ -223,24 +228,35 @@ class Room void simul_ray( float phi, float theta, - const Vectorf source_pos, + const Vectorf &source_pos, + float energy_0 + ); + + void simul_ray( + const Vectorf &ray_direction, + const Vectorf &source_pos, float energy_0 ); void ray_tracing( const Eigen::Matrix &angles, - const Vectorf source_pos + const Vectorf &source_pos + ); + + void ray_tracing( + const Eigen::Matrix &unit_vectors, + const Vectorf &source_pos ); void ray_tracing( size_t nb_phis, size_t nb_thetas, - const Vectorf source_pos + const Vectorf &source_pos ); void ray_tracing( size_t n_rays, - const Vectorf source_pos + const Vectorf &source_pos ); bool contains(const Vectorf point); @@ -254,7 +270,7 @@ class Room // Image source model internal methods void image_sources_dfs(ImageSource &is, int max_order); - bool is_visible_dfs(const Vectorf &p, ImageSource &is); + bool is_visible_dfs(const Vectorf &p, ImageSource &is, Vectorf &source_direction); bool is_obstructed_dfs(const Vectorf &p, ImageSource &is); int fill_sources(); diff --git a/pyroomacoustics/random/__init__.py b/pyroomacoustics/random/__init__.py new file mode 100644 index 00000000..2e126f17 --- /dev/null +++ b/pyroomacoustics/random/__init__.py @@ -0,0 +1,2 @@ +from . import distributions, sampler +from .spherical import power_spherical, uniform_spherical diff --git a/pyroomacoustics/random/distributions.py b/pyroomacoustics/random/distributions.py new file mode 100644 index 00000000..c64f96e9 --- /dev/null +++ b/pyroomacoustics/random/distributions.py @@ -0,0 +1,116 @@ +import abc + +import numpy as np +import scipy + +from .spherical import power_spherical, uniform_spherical + + +class Distribution(abc.ABC): + """ + Abstract base class for distributions. Derived classes + should implement a ``pdf`` method that provides the pdf value for samples + and a ``sample`` method that allows to sample from the distribution + """ + + def __init__(self, dim): + self._dim = dim + + @property + def dim(self): + return self._dim + + @abc.abstractmethod + def pdf(self): + pass + + @abc.abstractmethod + def sample(self): + pass + + def __call__(self, *args, **kwargs): + return self.sample(*args, **kwargs) + + +class AdHoc(Distribution): + """ + A helper class to construct distribution from separate pdf/sample + functions + """ + + def __init__(self, pdf, sampler, dim): + super().__init__(dim) + + self._pdf = pdf + self._sampler = sampler + + def pdf(self, *args, **kwargs): + return self._pdf(*args, **kwargs) + + def sample(self, *args, **kwargs): + return self._sampler(*args, **kwargs) + + +class UniformSpherical(Distribution): + """ + Uniform distribution on the n-sphere + """ + + def __init__(self, dim=3): + super().__init__(dim) + self._area = ( + self._dim + * np.pi ** (self._dim / 2) + / scipy.special.gamma(self._dim / 2 + 1) + ) + + @property + def _n_sphere_area(self): + return self._area + + def pdf(self, x): + return 1.0 / self._n_sphere_area + + def sample(self, size=None): + return uniform_spherical(dim=self.dim, size=size) + + +class PowerSpherical(Distribution): + def __init__(self, dim=3, loc=None, scale=None): + super().__init__(dim) + + # first canonical basis vector + if loc is None: + self._loc = np.zeros(self.dim) + self._loc[0] = 1.0 + else: + self._loc = loc + + if not scale: + self._scale = 1.0 + else: + self._scale = scale + + # computation the normalizing constant + a = (self.dim - 1) / 2 + self.scale + b = (self.dim - 1) / 2 + self._Z_inv = 1.0 / ( + 2 ** (a + b) + * np.pi**b + * scipy.special.gamma(a) + / scipy.special.gamma(a + b) + ) + + @property + def loc(self): + return self._loc + + @property + def scale(self): + return self._scale + + def pdf(self, x): + assert ( + x.shape[-1] == self.dim + ), "Input dimension does not match distribution dimension" + return self._Z_inv * (1.0 + np.matmul(x, self.loc)) ** self.scale diff --git a/pyroomacoustics/random/sampler.py b/pyroomacoustics/random/sampler.py new file mode 100644 index 00000000..12316cb6 --- /dev/null +++ b/pyroomacoustics/random/sampler.py @@ -0,0 +1,111 @@ +import abc + +import numpy as np +from numpy.random import default_rng + +from . import distributions +from .spherical import uniform_spherical + + +class RejectionSampler: + def __init__(self, desired_func, proposal_dist, scale=None): + + if not scale: + self._scale = 1.0 + else: + self._scale = scale + + self._desired_func = desired_func + self._proposal_dist = proposal_dist + + self._dim = self._proposal_dist.dim + + # get a random number generator + self._rng = default_rng() + + # keep track of the efficiency for analysis purpose + self._n_proposed = 0 + self._n_accepted = 0 + + @property + def efficiency(self): + if self._n_proposed == 0: + # before starting to sample, efficiency is maximum for consistency + # i.e., we haven't rejected anything + return 1.0 + else: + return self._n_accepted / self._n_proposed + + @property + def dim(self): + return self._dim + + def __call__(self, size=None): + + if not size: + size = 1 + + flat_size = int(np.prod(size)) # flat size + + offset = 0 + samples = np.zeros((flat_size, self.dim)) + + while offset < flat_size: + + n_propose = flat_size - offset + + proposal = self._proposal_dist.sample(size=n_propose) + proposal_pdf_value = self._proposal_dist.pdf(proposal) + desired_pdf_value = self._desired_func(proposal) + u = ( + self._rng.uniform(size=proposal_pdf_value.shape) + * proposal_pdf_value + * self._scale + ) + accept = np.where(u < desired_pdf_value)[0] + + n_accept = len(accept) + samples[offset : offset + n_accept, :] = proposal[accept, :] + + offset += n_accept + + # accounting + self._n_proposed += n_propose + self._n_accepted += n_accept + + if isinstance(size, int): + final_shape = [size, self.dim] + else: + final_shape = list(size) + [self.dim] + return samples.reshape(final_shape) + + +class DirectionalSampler(RejectionSampler): + def __init__(self, loc=None): + + if loc is None: + self._dim = 3 + self._loc = np.zeros(self._dim, dtype=float) + self._loc[0] = 1.0 + else: + self._loc = np.array(loc, dtype=float) + assert self._loc.ndim == 1 + self._dim = len(self._loc) + + self._loc /= np.linalg.norm(self._loc) + + # proposal distribution is the unnormalized uniform spherical distribution + unnormalized_uniform = distributions.AdHoc( + lambda x: np.ones_like(x.shape[-1]), uniform_spherical, self._dim + ) + super().__init__(self._pattern, unnormalized_uniform) + + @abc.abstractmethod + def _pattern(self, x): + """ + Parameters + ---------- + x: array_like, shape (..., n_dim) + Cartesian coordinates + """ + pass diff --git a/pyroomacoustics/random/spherical.py b/pyroomacoustics/random/spherical.py new file mode 100644 index 00000000..58653f13 --- /dev/null +++ b/pyroomacoustics/random/spherical.py @@ -0,0 +1,96 @@ +import numpy as np +from numpy.random import default_rng + +_eps = 1e-7 + + +def uniform_spherical(dim=3, size=None): + """ + Generates uniform samples on the n-sphere + + Parameters + ---------- + size: int or tuple of ints, optional + The number of samples to generate + dim: int, optional + The number of dimensions of the sphere, the default is dim=3 + + Returns + ------- + out: ndarray, shape (*size, dim) + The samples draw from the uniform distribution on the n-sphere + """ + if size is None: + size = [1, dim] + elif isinstance(size, int): + size = [size, dim] + else: + size = list(size) + [dim] + + rng = default_rng() + + out = rng.standard_normal(size=size) + out /= np.linalg.norm(out, axis=-1, keepdims=True) + + return out + + +def power_spherical(loc=None, scale=None, dim=3, size=None): + """ + Generates power spherical samples on the (n-1)-sphere according to + + Nicola De Cao, Wilker Aziz, "The Power Spherical distribution", arXiv, 2020. + http://arxiv.org/abs/2006.04437v1 + + Parameters + ---------- + loc: float or array_like of floats, optional + The location (i.e., direction) unit vector + scale: float or array_like of floats + The scale parameter descibing the spread of the distribution + dim: int, optional + The number of dimensions of the sphere, the default is dim=3 + size: int or tuple of ints, optional + The number of samples to generate + + Returns + ------- + out: ndarray, shape (*size, dim) + The samples draw from the uniform distribution on the n-sphere + """ + + if size is None: + size = [1, dim] + elif isinstance(size, int): + size = [size, dim] + else: + size = list(size) + [dim] + + e1 = np.zeros(dim) + e1[0] = 1.0 + + if loc is None: + loc = e1.copy() + + if scale is None: + scale = 1.0 + + rng = default_rng() + + z = rng.beta((dim - 1.0) / 2.0 + scale, (dim - 1) / 2.0, size=size[:-1]) + v = uniform_spherical(size=size[:-1], dim=dim - 1) + + t = 2 * z[..., None] - 1 + y = np.concatenate((t, np.sqrt(1 - t**2) * v), axis=-1) # shape (*size, dim) + + u_hat = e1 - loc + # here the _eps is to avoid division by zero so that it is fine that when + # e1 == loc the vector u is zero + # this was verified in the code from the original paper + # https://github.com/nicola-decao/power_spherical/blob/master/power_spherical/distributions.py + u = u_hat / (np.linalg.norm(u_hat, axis=-1, keepdims=True) + _eps) + + # out = y - u * (u[..., None, :] @ y[..., :, None])[..., 0] + out = y - 2 * u * (u * y).sum(axis=-1, keepdims=True) + + return out diff --git a/pyroomacoustics/random/tests/test_rejection.py b/pyroomacoustics/random/tests/test_rejection.py new file mode 100644 index 00000000..ee989dbe --- /dev/null +++ b/pyroomacoustics/random/tests/test_rejection.py @@ -0,0 +1,56 @@ +import pyroomacoustics as pra +from pyroomacoustics.directivities import CardioidFamilySampler + +if __name__ == "__main__": + import matplotlib.pyplot as plt + import mpl_toolkits.mplot3d.axes3d as axes3d + import numpy as np + + # let's visualize samples in the sphere + + theta, phi = np.linspace(0, 2 * np.pi, 20), np.linspace(0, np.pi, 20) + THETA, PHI = np.meshgrid(theta, phi) + X, Y, Z = np.sin(PHI) * np.cos(THETA), np.sin(PHI) * np.sin(THETA), np.cos(PHI) + + """ + fig = plt.figure(figsize=(8, 8)) + ax = fig.gca(projection="3d") + ax.plot_wireframe(X, Y, Z, linewidth=1, alpha=0.25, color="gray") + """ + + for loc, scale in zip([[1, 1, 1]], [10, 1, 0.1]): + print(loc, scale) + # Figure-of-eight + sampler = CardioidFamilySampler(loc=loc, p=0) + + # Measured eigenmike response + eigenmike = pra.MeasuredDirectivityFile("EM32_Directivity", fs=16000) + rot_54_73 = pra.Rotation3D([73, 54], "yz", degrees=True) + dir_obj_Emic = eigenmike.get_mic_directivity("EM_32_9", orientation=rot_54_73) + sampler = dir_obj_Emic._ray_sampler + + points = sampler(size=100000).T # shape (n_dim, n_points) + + # Create a spherical histogram + hist = pra.doa.SphericalHistogram(n_bins=500) + hist.push(points) + hist.plot() + + """ + ax.scatter(X, Y, Z, s=50) + ax.plot( + *np.stack((torch.zeros_like(loc), loc)).T, + linewidth=4, + label="$\kappa={}$".format(scale) + ) + """ + + print("Sampler's efficiency:", sampler.efficiency) + + """ + ax.view_init(30, 45) + ax.tick_params(axis="both") + plt.legend() + """ + + plt.show() diff --git a/pyroomacoustics/random/tests/test_spherical.py b/pyroomacoustics/random/tests/test_spherical.py new file mode 100644 index 00000000..b119244b --- /dev/null +++ b/pyroomacoustics/random/tests/test_spherical.py @@ -0,0 +1,49 @@ +import numpy as np +from numpy.random import default_rng + +import pyroomacoustics as pra + + +def test_stat_prop_power_spherical(): + + rng = default_rng() + + for d in range(2, 4): + + pass + + +if __name__ == "__main__": + + import matplotlib.pyplot as plt + import mpl_toolkits.mplot3d.axes3d as axes3d + import numpy as np + + # let's visualize samples in the sphere + + theta, phi = np.linspace(0, 2 * np.pi, 20), np.linspace(0, np.pi, 20) + THETA, PHI = np.meshgrid(theta, phi) + X, Y, Z = np.sin(PHI) * np.cos(THETA), np.sin(PHI) * np.sin(THETA), np.cos(PHI) + + fig = plt.figure(figsize=(8, 8)) + ax = fig.gca(projection="3d") + ax.plot_wireframe(X, Y, Z, linewidth=1, alpha=0.25, color="gray") + + for loc, scale in zip(np.eye(3), [10, 1, 0.1]): + print(loc, scale) + X, Y, Z = pra.random.power_spherical(loc=loc, scale=scale, size=1000).T + + ax.scatter(X, Y, Z, s=50) + """ + ax.plot( + *np.stack((torch.zeros_like(loc), loc)).T, + linewidth=4, + label="$\kappa={}$".format(scale) + ) + """ + + ax.view_init(30, 45) + ax.tick_params(axis="both") + plt.legend() + + plt.show() diff --git a/pyroomacoustics/room.py b/pyroomacoustics/room.py index 329842d2..852c99d4 100644 --- a/pyroomacoustics/room.py +++ b/pyroomacoustics/room.py @@ -2155,22 +2155,25 @@ def image_source_model(self): if n_visible_sources > 0: # Copy to python managed memory - - source.images = ( - self.room_engine.sources.copy() - ) # Positions of the image source (3,n) n: n_sources - source.orders = ( - self.room_engine.orders.copy() - ) # Reflection order for each image source shape n:n_sources - source.orders_xyz = ( - self.room_engine.orders_xyz.copy() - ) # Reflection order for each image source for each coordinate shape (3,n) n:n_sources - source.walls = ( - self.room_engine.gen_walls.copy() - ) # Something that i don't get [-1,-1,-1,-1,-1...] shape n:n_sources - source.damping = ( - self.room_engine.attenuations.copy() - ) # Octave band damping's shape (no_of_octave_bands*n_sources) damping value for each image source for each octave bands + # image source locations + source.images = self.room_engine.sources.copy() + # image source reflection orders + source.orders = self.room_engine.orders.copy() + # image source reflection orders per x/y/z axis (shoebox only) + source.orders_xyz = self.room_engine.orders_xyz.copy() + # the microphone direction from the perspective of the image source + source.directions = ( + self.room_engine.source_directions.reshape( + (self.dim, -1, self.mic_array.M) + ) + .transpose(2, 0, 1) + .copy() + ) # (n_mics, dim, n_image_sources) + # the index of the wall that produced the image source + source.walls = self.room_engine.gen_walls.copy() + # the attenuation coefficient associated with the image source + source.damping = self.room_engine.attenuations.copy() + # --deprecated-- source.generators = -np.ones(source.walls.shape) # if randomized image method is selected, add a small random diff --git a/pyroomacoustics/tests/test_source_directivity_flipping.py b/pyroomacoustics/tests/test_source_directivity_flipping.py index 0495e617..c65abb87 100644 --- a/pyroomacoustics/tests/test_source_directivity_flipping.py +++ b/pyroomacoustics/tests/test_source_directivity_flipping.py @@ -7,6 +7,88 @@ class TestSourceDirectivityFlipping(TestCase): + def test_source_directions(self): + # create room + mic_loc = [5, 12, 5] + source_loc = [5, 2, 5] + room = ( + pra.ShoeBox( + p=[10, 14, 10], + max_order=2, + ) + .add_source(source_loc) + .add_microphone(mic_loc) + ) + + # compute image sources + room.image_source_model() + + # compute azimuth_s and colatitude_s pair for images along x-axis + source_angle_array = source_angle_shoebox( + image_source_loc=room.sources[0].images, + wall_flips=abs(room.sources[0].orders_xyz), + mic_loc=mic_loc, + ) + source_angle_array = np.array(source_angle_array) + + source_dir = room.sources[0].directions[0] + azimuth = np.arctan2(source_dir[1], source_dir[0]) + colatitude = np.pi / 2 - np.arcsin(source_dir[2]) + source_angle_array_2 = np.array([azimuth, colatitude]) + + np.testing.assert_almost_equal( + source_angle_array, source_angle_array_2, decimal=4 + ) + + def test_robin_check(self): + # create room + mic_loc = [5, 12, 5] + source_loc = [5, 2, 5] + room = ( + pra.ShoeBox( + p=[10, 14, 10], + max_order=1, + ) + .add_source(source_loc) + .add_microphone(mic_loc) + ) + + # compute image sources + room.image_source_model() + + # compute azimuth_s and colatitude_s pair for images along x-axis + source_angle_array = source_angle_shoebox( + image_source_loc=room.sources[0].images, + wall_flips=abs(room.sources[0].orders_xyz), + mic_loc=mic_loc, + ) + source_angle_array = np.array(source_angle_array) + + skipped = 0 + for i in range(7): + img = room.sources[0].images[:, i] + if np.allclose(img, [5, -2, 5]): + sd = [-np.pi / 2, np.pi / 2] + elif np.allclose(img, [-5, 2, 5]): + sd = [3 * np.pi / 4, np.pi / 2] + elif np.allclose(img, [5, 26, 5]): + sd = [np.pi / 2, np.pi / 2] + elif np.allclose(img, [15, 2, 5]): + sd = [np.pi / 4, np.pi / 2] + elif np.allclose(img, [5, 2, -5]): + sd = [np.pi / 2, 3 * np.pi / 4] + elif np.allclose(img, [5, 2, 15]): + sd = [np.pi / 2, np.pi / 4] + elif np.allclose(img, [5, 2, 5]): + sd = [np.pi / 2, np.pi / 2] + else: + skipped += 1 + continue + + np.testing.assert_almost_equal(source_angle_array[:, i], sd) + + assert skipped == 0 + def test_x(self): # create room mic_loc = [5, 12, 5] diff --git a/pyroomacoustics/tests/tests_libroom/test_microphone.py b/pyroomacoustics/tests/tests_libroom/test_microphone.py new file mode 100644 index 00000000..269229a9 --- /dev/null +++ b/pyroomacoustics/tests/tests_libroom/test_microphone.py @@ -0,0 +1,86 @@ +import matplotlib.pyplot as plt +import numpy as np + +import pyroomacoustics as pra + + +def inspect_directional_receiver_pattern_anechoic(): + room_dim = [10, 10, 4] + materials = pra.Material(1.0) # fully absorbant + + materials = pra.make_materials( + north=0.0, + south=1.0, + east=1.0, + west=1.0, + floor=1.0, + ceiling=1.0, + ) + max_order = -1 # deactivate ism + + room = pra.ShoeBox(room_dim, materials=materials, max_order=-1) + room.set_ray_tracing( + n_rays=100000, receiver_radius=0.1, time_thres=0.1, hist_bin_size=1.0 + ) + room.add_microphone([8, 5, 2]) + + # add directional receiver + grid = pra.doa.GridSphere(n_points=256) + room.room_engine.microphones[0].set_directions(grid.cartesian.T) + + room.room_engine.ray_tracing(100000, [2, 5, 2]) + + bins = [] + for hist in room.room_engine.microphones[0].histograms: + bins.append(np.sum(hist.get_hist())) + + grid.set_values(bins) + print(bins) + + grid.plot_old() + + +def inspect_directional_receiver_pattern_reverb(): + room_dim = [10, 10, 4] + materials = pra.Material(0.15) # fully absorbant + max_order = -1 # deactivate ism + + room = pra.ShoeBox(room_dim, materials=materials, max_order=-1) + room.set_ray_tracing( + n_rays=100000, receiver_radius=0.1, time_thres=0.1, hist_bin_size=1.0 + ) + room.add_microphone([8, 5, 2]) + + # add directional receiver + grid = pra.doa.GridSphere(n_points=256) + room.room_engine.microphones[0].set_directions(grid.cartesian.T) + + room.room_engine.ray_tracing(100000, [2, 5, 2]) + + bins = [] + for hist in room.room_engine.microphones[0].histograms: + bins.append(np.sum(hist.get_hist())) + + grid.set_values(bins) + print(bins) + + grid.plot_old() + + +if __name__ == "__main__": + + mic = pra.libroom.Microphone([0.0, 0.0, 0.0], 1, 0.1, 1.0) + grid = pra.doa.GridSphere(n_points=16) + dirs = grid.cartesian.T.copy() + mic.set_directions(dirs) + + test_dir = dirs[0] + + mic.log_histogram(0.5, [1.0], test_dir) + + print(mic.histograms[0].get_hist()) + print(mic.histograms[1].get_hist()) + + inspect_directional_receiver_pattern_anechoic() + inspect_directional_receiver_pattern_reverb() + plt.show() diff --git a/setup.py b/setup.py index 544a6f08..340b934d 100644 --- a/setup.py +++ b/setup.py @@ -68,6 +68,7 @@ def __str__(self): str(get_pybind_include()), str(get_pybind_include(user=True)), os.path.join(libroom_src_dir, "ext/eigen"), + os.path.join(libroom_src_dir, "ext/nanoflann/include"), ], language="c++", extra_compile_args=[