Adds distribution objects for measured and analytic directivities for…

… sampling rays. Adds tests for these distributions as well as for the random sub-package.

CHANGELOG.rst

@@ -17,13 +17,25 @@ tracing simulation engine.
 Added
 ~~~~~
 
-- New ``pyroomacoustics.random`` module that provides some primitives for sampling
-  at random from arbitrary distributions on the sphere. This is used for source directivities
-  in the ray tracing simulator.
+- Support for source directivities in non-shoebox rooms using the images source
+  model.
 
-- New octave filter bank with energy conservation and perfect reconstruction described
-  in Antoni, "Orthogonal-like fractional-octave-band filters," 2009.
-  The filter bank is implemented in ``pyroomacoustics.acoustics.AntoniOctaveFilterBank``.
+- New ``pyroomacoustics.random`` module that provides some primitives for
+  sampling at random from arbitrary distributions on the sphere. This is used
+  for source directivities in the ray tracing simulator.
+
+- New octave filter bank with energy conservation and perfect reconstruction
+  described in Antoni, "Orthogonal-like fractional-octave-band filters," 2009.
+  The filter bank is implemented in
+  ``pyroomacoustics.acoustics.AntoniOctaveFilterBank``.
+
+- A method ``sample_rays`` is added to the ``Directivity`` objects to provide a
+  unified interface to sample rays of sources used for ray tracing.
+
+Changed
+~~~~~~~
+
+- Bumped the numpy requirement to v1.17.0 to use the ``numpy.random.Generator`` objects.
 
 
 `0.8.3`_ - 2024-12-08

pyroomacoustics/directivities/__init__.py

@@ -84,13 +84,20 @@
 from .analytic import (
     Cardioid,
     CardioidFamily,
-    CardioidFamilySampler,
     FigureEight,
     HyperCardioid,
     Omnidirectional,
     SubCardioid,
     cardioid_func,
+    cardioid_energy,
+    CardioidEnergyDistribution,
 )
+from .histogram import SphericalHistogram
+from .integration import spherical_integral
 from .base import Directivity
 from .direction import DirectionVector, Rotation3D
-from .measured import MeasuredDirectivity, MeasuredDirectivityFile
+from .measured import (
+    MeasuredDirectivity,
+    MeasuredDirectivityFile,
+    MeasuredDirectivityEnergyDistribution,
+)

pyroomacoustics/directivities/analytic.py

@@ -62,6 +62,7 @@
 from .. import random
 from ..doa import spher2cart
 from ..utilities import all_combinations, requires_matplotlib
+from .. import random
 from .base import Directivity
 from .direction import DirectionVector
 
@@ -124,8 +125,8 @@ def __init__(self, orientation, p, gain=1.0):
         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(
+        # distribution corresponding to the source energy directivity
+        self.energy_distribution = CardioidEnergyDistribution(
             loc=self._orientation.unit_vector, p=self._p
         )
 
@@ -207,8 +208,8 @@ def get_response(
             else:
                 return resp
 
-    def sample_rays(self, n_rays):
-        return self._ray_sampler(n_rays).T
+    def sample_rays(self, n_rays, rng=None):
+        return self.energy_distribution.sample(n_rays, rng=rng).T
 
     @requires_matplotlib
     def plot_response(
@@ -393,39 +394,6 @@ 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.
@@ -486,3 +454,63 @@ def cardioid_energy(p, gain=1.0):
         The gain of the cardioid function
     """
     return gain**2 * (4.0 * np.pi / 3.0) * (4 * p**2 - 2 * p + 1)
+
+
+class CardioidEnergyDistribution(random.distributions.Distribution):
+    """
+    This object draws samples from a cardioid shaped distribution on the sphere
+
+    Parameters
+    ----------
+    loc: array_like or None
+        The unit vector pointing in the main direction of the cardioid.
+        If None, the location defaults to [1, 0, 0].
+    p: float or None
+        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. When None, p is set to 0.5.
+        The parameter must be between 0 and 1
+    """
+
+    def __init__(self, loc=None, p=None):
+        super().__init__(dim=3)
+        if loc is None:
+            loc = np.array([1.0, 0.0, 0.0])
+
+        loc_norm = np.linalg.norm(loc)
+        if abs(loc_norm - 1.0) > 1e-5:
+            raise ValueError(
+                "The location parameter should be a unit "
+                f"norm vector (got norm={loc_norm})."
+            )
+
+        self._loc = loc
+
+        if not (0.0 <= p <= 1.0):
+            raise ValueError(f"The value of p should be between 0.0 and 1.0 (got {p=})")
+
+        self._coeff = p
+
+        self._sampler = random.sampler.RejectionSampler(
+            desired_func=self._unnormalized_pdf,
+            proposal_dist=random.distributions.UnnormalizedUniformSpherical(dim=3),
+            scale=1.0,
+        )
+
+    def _unnormalized_pdf(self, x):
+        x = np.moveaxis(x, -1, 0)
+        response = cardioid_func(
+            x,
+            direction=self._loc,
+            p=self._coeff,
+            gain=1.0,
+            normalize=False,
+            magnitude=True,
+        )
+        return response**2
+
+    def pdf(self, x):
+        return self._unnormalized_pdf(x) / cardioid_energy(p=self._coeff)
+
+    def sample(self, size=None, rng=None):
+        return self._sampler(size=size, rng=rng)

pyroomacoustics/directivities/base.py

@@ -94,7 +94,7 @@ def get_response(
         raise NotImplementedError
 
     @abc.abstractmethod
-    def sample_rays(self, n_rays):
+    def sample_rays(self, n_rays, rng=None):
         """
         This method samples unit vectors from the sphere according to
         the distribution of the source
@@ -103,6 +103,10 @@ def sample_rays(self, n_rays):
         ----------
         n_rays: int
             The number of rays to sample
+        rng: numpy.random.Generator or None, optional
+            A random number generator object from numpy or None.
+            If None is passed numpy.random.default_rng is used to create
+            a Generator object.
 
         Returns
         -------

pyroomacoustics/doa/histogram.py → pyroomacoustics/directivities/histogram.py

@@ -27,7 +27,7 @@
 import numpy as np
 from scipy.spatial import SphericalVoronoi, cKDTree
 
-from .doa import GridSphere
+from ..doa import GridSphere
 
 
 class SphericalHistogram:
@@ -60,6 +60,14 @@ def __init__(self, n_bins, dim=3, enable_peak_finding=False):
         self._cache_dirty = False
         self._cache_histogram = np.zeros(self.n_bins)
 
+    @property
+    def grid(self):
+        return self._grid
+
+    @property
+    def bin_areas(self):
+        return self._areas
+
     @property
     def n_dim(self):
         return self._n_dim

pyroomacoustics/directivities/measured.py

@@ -73,12 +73,13 @@
 
 import numpy as np
 from scipy.interpolate import griddata
-from scipy.spatial import cKDTree
+from scipy.spatial import cKDTree, SphericalVoronoi
 
 from .. import random
 from ..acoustics import OctaveBandsFactory
 from ..datasets import SOFADatabase
 from ..doa import Grid, GridSphere, cart2spher, fibonacci_spherical_sampling, spher2cart
+from ..directivities.integration import spherical_integral
 from ..parameters import constants
 from ..utilities import requires_matplotlib
 from .base import Directivity
@@ -87,7 +88,7 @@
 from .sofa import open_sofa_file
 
 
-class MeasuredDirectivitySampler(random.sampler.DirectionalSampler):
+class MeasuredDirectivityEnergyDistribution(random.distributions.Distribution):
     """
     This object draws samples from the distribution defined by the energy
     of the measured directional response object.
@@ -101,16 +102,31 @@ class MeasuredDirectivitySampler(random.sampler.DirectionalSampler):
     """
 
     def __init__(self, kdtree, energy):
-        super().__init__()
+        super().__init__(kdtree.m)
         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()
 
+        self._sampler = random.sampler.RejectionSampler(
+            desired_func=self._pattern,
+            proposal_dist=random.distributions.UnnormalizedUniformSpherical(self.dim),
+            scale=1.0,
+        )
+
+        # The normalizing constant of the pdf is computed by spherical integration.
+        sv = SphericalVoronoi(kdtree.data)
+        w_ = sv.calculate_areas()
+        self._normalizing_constant = np.sum(w_ * self._energy)
+
     def _pattern(self, x):
         _, index = self._kdtree.query(x)
         return self._energy[index]
 
+    def pdf(self, x):
+        return self._pattern(x) / self._normalizing_constant
+
+    def sample(self, size=None, rng=None):
+        return self._sampler(size=size, rng=rng)
+
 
 class MeasuredDirectivity(Directivity):
     """
@@ -184,7 +200,9 @@ def set_orientation(self, orientation):
         # create the kd-tree
         self._kdtree = cKDTree(self._grid.cartesian.T)
         ir_energy = np.square(self._irs).sum(axis=-1)
-        self._ray_sampler = MeasuredDirectivitySampler(self._kdtree, ir_energy)
+        self.energy_distribution = MeasuredDirectivityEnergyDistribution(
+            self._kdtree, ir_energy
+        )
 
     def get_response(
         self, azimuth, colatitude=None, magnitude=False, frequency=None, degrees=True
@@ -214,8 +232,8 @@ 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
+    def sample_rays(self, n_rays, rng=None):
+        return self.energy_distribution.sample(size=n_rays, rng=rng).T
 
     @requires_matplotlib
     def plot(self, freq_bin=0, n_grid=100, ax=None, depth=False, offset=None):

pyroomacoustics/directivities/tests/test_cardioid_energy_distribution.py

@@ -0,0 +1,65 @@
+import numpy as np
+import pyroomacoustics as pra
+
+import pytest
+
+
+@pytest.mark.parametrize("p", (0.0, 0.25, 0.5, 0.75, 1.0))
+def test_pdf_integral(p):
+
+    loc = np.ones(3) / np.sqrt(3.0)
+
+    distribution = pra.directivities.CardioidEnergyDistribution(loc=loc, p=p)
+
+    def pdf(x):
+        return distribution.pdf(x.T)
+
+    area = pra.directivities.spherical_integral(pdf, 1000)
+
+    assert abs(area - 1.0) < 1e-5
+
+
+@pytest.mark.parametrize("p", (0.0, 0.25, 0.5, 0.75, 1.0))
+def test_rejection_sampler_cardioid(p):
+
+    loc = np.ones(3) / np.sqrt(3.0)
+
+    cardioid_energy_distribution = pra.directivities.CardioidEnergyDistribution(
+        loc=loc, p=p
+    )
+
+    rng = np.random.default_rng(94877675)
+
+    hist = pra.directivities.SphericalHistogram(n_bins=36)
+
+    random_points = cardioid_energy_distribution.sample(size=(100000,), rng=rng)
+    hist.push(random_points.T)
+
+    values_expected = cardioid_energy_distribution.pdf(hist.grid.cartesian.T)
+    values_obtained = hist.histogram
+
+    idx = np.argsort(values_expected)[::-1]
+    cdf = np.cumsum(values_expected[idx] * hist._areas[idx])
+    select = idx[cdf <= 1.0]
+
+    np.testing.assert_allclose(
+        values_expected[select], values_obtained[select], atol=0.02, rtol=0.1
+    )
+
+
+@pytest.mark.parametrize("p", (0.0, 0.25, 0.5, 0.75, 1.0))
+def test_rejection_sampler_cardioid_repeatability(p):
+
+    loc = np.ones(3) / np.sqrt(3.0)
+
+    cardioid_energy_distribution = pra.directivities.CardioidEnergyDistribution(
+        loc=loc, p=p
+    )
+
+    rng = np.random.default_rng(94877675)
+    random_points1 = cardioid_energy_distribution.sample(size=(100000,), rng=rng)
+
+    rng = np.random.default_rng(94877675)
+    random_points2 = cardioid_energy_distribution.sample(size=(100000,), rng=rng)
+
+    np.testing.assert_allclose(random_points1, random_points2)