[FEATURE] add pytta.rays, stochastic ray tracing in convex rooms

docs/pyttadoc/index.rst

@@ -64,4 +64,5 @@ Modules:
    functions
    generate
    rooms
+   rays
    properties

docs/pyttadoc/rays.rst

@@ -0,0 +1,7 @@
+.. currentmodule:: pytta
+
+Ray tracing
+===========
+
+.. automodule:: pytta.rays
+   :members:

pytta/__init__.py

@@ -72,6 +72,9 @@
     * pytta.iso3741:
         calculations according to the standard;
 
+    * pytta.rays:
+        stochastic ray tracing in convex rooms;
+
 
 You can find usage examples in the toolbox's examples folder.
 
@@ -108,6 +111,7 @@
 from . import utils
 from . import rooms
 from . import iso3741
+from . import rays
 
 from .functions import read_wav, write_wav, merge, split, \
     list_devices, print_devices, get_device_from_user,\
@@ -130,6 +134,7 @@
            'generate',
            'utils',
            'iso3741',
+           'rays',
            'rooms',
 
            # Functions

pytta/rays.py

@@ -0,0 +1,113 @@
+# -*- coding: utf-8 -*-
+"""
+Stochastic ray tracing in convex rooms (Krokstad, Strom and Sorsdal,
+Calculating the acoustical room response by the use of a ray tracing
+technique, J. Sound Vib. 8(1), 1968; Vorlaender, Auralization, 2008),
+as in the ray tracing part of RAVEN, which the ITA-Toolbox drives.
+
+A room is a list of planes, each with its absorption and scattering
+coefficients. Rays leave the source in random directions with equal
+energy, lose the absorbed share at each wall and are reflected
+specularly, or, with the probability of the scattering coefficient, in a
+random direction drawn from Lambert's law. A sphere around the receiver
+counts each ray crossing it with its energy times the chord length over
+the sphere volume, so the histogram is the energy per unit area arriving
+in each time bin (E0 / (4 pi r^2) for the direct sound).
+
+>>> room = pytta.rays.shoebox(6, 5, 4, alpha=0.2, scattering=0.5)
+>>> t, h = pytta.rays.trace(room, [1, 1, 1.5], [4, 3, 1.2], duration=1.)
+
+"""
+
+import numpy as np
+
+__all__ = ['shoebox', 'convex_room', 'trace']
+
+
+def shoebox(Lx, Ly, Lz, alpha, scattering=0.):
+    """Room [0, Lx] x [0, Ly] x [0, Lz]; alpha and scattering are one value
+    for all walls or six, in the order x = 0, x = Lx, y = 0, y = Ly, z = 0,
+    z = Lz."""
+    normals = np.array([[-1, 0, 0], [1, 0, 0], [0, -1, 0], [0, 1, 0],
+                        [0, 0, -1], [0, 0, 1]], float)
+    offsets = np.array([0, Lx, 0, Ly, 0, Lz], float)
+    return _room(normals, offsets, alpha, scattering)
+
+
+def convex_room(vertices, alpha, scattering=0.):
+    """Convex room spanned by its corner points (M, 3), one plane per face of
+    their convex hull; alpha and scattering are one value for all faces or
+    one per face, in the order of room['normals']."""
+    from scipy.spatial import ConvexHull
+    eq = ConvexHull(np.asarray(vertices, float)).equations
+    normals, offsets = eq[:, :3], -eq[:, 3]
+    keep = np.unique(np.round(np.column_stack([normals, offsets]), 9), axis=0, return_index=True)[1]
+    keep.sort()
+    return _room(normals[keep], offsets[keep], alpha, scattering)
+
+
+def _room(normals, offsets, alpha, scattering):
+    n = len(offsets)
+    return {'normals': normals, 'offsets': offsets,
+            'alpha': np.broadcast_to(np.asarray(alpha, float), (n,)).copy(),
+            'scattering': np.broadcast_to(np.asarray(scattering, float), (n,)).copy()}
+
+
+def _lambert(normals, rng):
+    # Directions about the inward normals -normals with density cos(theta).
+    u, v = rng.uniform(size=(2, len(normals)))
+    z = np.sqrt(1 - u)
+    r = np.sqrt(u)
+    phi = 2 * np.pi * v
+    axis = -normals
+    helper = np.where(np.abs(axis[:, :1]) < 0.9, [[1., 0, 0]], [[0, 1., 0]])
+    e1 = np.cross(axis, helper)
+    e1 /= np.linalg.norm(e1, axis=1, keepdims=True)
+    e2 = np.cross(axis, e1)
+    return (r * np.cos(phi))[:, None] * e1 + (r * np.sin(phi))[:, None] * e2 + z[:, None] * axis
+
+
+def trace(room, source, receiver, numRays=20000, duration=1., radius=0.5,
+          dt=1e-3, c=343., m=0., seed=None):
+    """
+    Energy histogram at the receiver of a unit energy impulse at the
+    source: returns (bin start times in s, energy per unit area per bin).
+    m is the power attenuation of air in 1/m.
+    """
+    rng = np.random.default_rng(seed)
+    src, rcv = np.asarray(source, float), np.asarray(receiver, float)
+    D = rng.standard_normal((numRays, 3))
+    D /= np.linalg.norm(D, axis=1, keepdims=True)
+    P = np.tile(src, (numRays, 1))
+    E = np.full(numRays, 1 / numRays)
+    L = np.zeros(numRays)
+    hist = np.zeros(int(np.ceil(duration / dt)))
+    volume = 4 / 3 * np.pi * radius**3
+    nrm, off = room['normals'], room['offsets']
+    alive = np.ones(numRays, bool)
+    while alive.any():
+        P, D, E, L = P[alive], D[alive], E[alive], L[alive]
+        # Distance to the wall the ray leaves through.
+        speed = D @ nrm.T
+        with np.errstate(divide='ignore', invalid='ignore'):
+            t = np.where(speed > 1e-12, (off - P @ nrm.T) / speed, np.inf)
+        wall = np.argmin(t, axis=1)
+        t = t[np.arange(len(t)), wall]
+        # Crossing of the receiver sphere on the way.
+        s = np.sum((rcv - P) * D, axis=1)
+        d2 = np.sum((rcv - P)**2, axis=1) - s**2
+        hit = (s > 0) & (s < t) & (d2 < radius**2)
+        when = ((L + s)[hit] / c / dt).astype(int)
+        ok = when < hist.size
+        chord = 2 * np.sqrt(radius**2 - d2[hit])
+        np.add.at(hist, when[ok], (E[hit] * chord / volume * np.exp(-m * (L + s)[hit]))[ok])
+        # Reflection.
+        P = P + t[:, None] * D
+        L = L + t
+        E = E * (1 - room['alpha'][wall]) * np.exp(-m * t)
+        n = nrm[wall]
+        specular = D - 2 * np.sum(D * n, axis=1, keepdims=True) * n
+        diffuse = rng.uniform(size=len(D)) < room['scattering'][wall]
+        D = np.where(diffuse[:, None], _lambert(n, rng), specular)
+        alive = (L < c * duration) & (E > 0)
+    return np.arange(hist.size) * dt, hist

tests/test_rays.py

@@ -0,0 +1,58 @@
+"""
+Tests for pytta.rays, stochastic ray tracing in convex rooms. Each test
+compares the energy histogram with a result of statistical room acoustics
+(Kuttruff, Room Acoustics).
+"""
+import unittest
+
+import numpy as np
+
+from pytta import rays
+
+C, DT = 343., 1e-3
+SRC, RCV = np.array([1.5, 1.2, 1.4]), np.array([4.2, 3.1, 2.2])
+V, S = 6 * 5 * 4, 2 * (30 + 24 + 20)
+
+
+class TestRays(unittest.TestCase):
+
+    def test_direct_sound(self):
+        # The direct sound carries E0 / (4 pi r^2) per unit area.
+        t, h = rays.trace(rays.shoebox(6, 5, 4, 0.2, 1.), SRC, RCV, 200000, 0.02, seed=1)
+        r = np.linalg.norm(RCV - SRC)
+        i = int(r / C / DT)
+        ratio = h[i - 1:i + 2].sum() * 4 * np.pi * r**2
+        print(f'\n  direct sound: {ratio:.3f} of E0 / (4 pi r^2)', end='')
+        self.assertLess(abs(ratio - 1), 0.08)
+
+    def test_lossless_diffuse_room_holds_e0_over_v(self):
+        # Without absorption the energy stays in the room and spreads
+        # evenly: each bin gets the energy density E0 / V times c dt.
+        t, h = rays.trace(rays.shoebox(6, 5, 4, 0., 1.), SRC, RCV, 20000, 0.5, seed=2)
+        late = h[(t > 0.2) & (t < 0.5)].mean() / (C * DT / V)
+        print(f'\n  lossless room: {late:.3f} of E0 c dt / V', end='')
+        self.assertLess(abs(late - 1), 0.03)
+
+    def test_decay_lies_between_eyring_and_sabine(self):
+        # With alpha = 0.2 on every wall and diffuse reflections the decay
+        # falls between Eyring and Sabine, as Kuttruff's correction for
+        # the spread of the free paths predicts.
+        t, h = rays.trace(rays.shoebox(6, 5, 4, 0.2, 1.), SRC, RCV, 20000, 1.2, seed=3)
+        e = np.cumsum(h[::-1])[::-1]
+        L = 10 * np.log10(e / e[0])
+        i5, i35 = np.argmax(L < -5), np.argmax(L < -35)
+        T30 = -60 / np.polyfit(t[i5:i35], L[i5:i35], 1)[0]
+        eyring = 0.161 * V / (-S * np.log(0.8))
+        sabine = 0.161 * V / (S * 0.2)
+        print(f'\n  T30 {T30:.3f} s, Eyring {eyring:.3f} s, Sabine {sabine:.3f} s', end='')
+        self.assertTrue(eyring < T30 < sabine)
+
+    def test_convex_room_from_corners(self):
+        # The hull of the corners of the shoebox gives its six planes.
+        corners = np.array([[x, y, z] for x in (0, 6) for y in (0, 5) for z in (0, 4)])
+        room = rays.convex_room(corners, 0.1)
+        self.assertEqual(len(room['offsets']), 6)
+
+
+if __name__ == '__main__':
+    unittest.main()