From 91029768392909e601539c4918f2a1338a66057c Mon Sep 17 00:00:00 2001 From: Sergio Aguirre Date: Sun, 1 Nov 2026 12:32:38 +0100 Subject: [PATCH] [FEATURE] add pytta.rays, stochastic ray tracing in convex rooms Co-Authored-By: Claude Fable 5.1 --- docs/pyttadoc/index.rst | 1 + docs/pyttadoc/rays.rst | 7 +++ pytta/__init__.py | 5 ++ pytta/rays.py | 113 ++++++++++++++++++++++++++++++++++++++++ tests/test_rays.py | 58 +++++++++++++++++++++ 5 files changed, 184 insertions(+) create mode 100644 docs/pyttadoc/rays.rst create mode 100644 pytta/rays.py create mode 100644 tests/test_rays.py diff --git a/docs/pyttadoc/index.rst b/docs/pyttadoc/index.rst index ca3da8d..a9142c2 100644 --- a/docs/pyttadoc/index.rst +++ b/docs/pyttadoc/index.rst @@ -64,4 +64,5 @@ Modules: functions generate rooms + rays properties diff --git a/docs/pyttadoc/rays.rst b/docs/pyttadoc/rays.rst new file mode 100644 index 0000000..bf780c8 --- /dev/null +++ b/docs/pyttadoc/rays.rst @@ -0,0 +1,7 @@ +.. currentmodule:: pytta + +Ray tracing +=========== + +.. automodule:: pytta.rays + :members: diff --git a/pytta/__init__.py b/pytta/__init__.py index d8740aa..d2be3ed 100644 --- a/pytta/__init__.py +++ b/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 diff --git a/pytta/rays.py b/pytta/rays.py new file mode 100644 index 0000000..a6bd86e --- /dev/null +++ b/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 diff --git a/tests/test_rays.py b/tests/test_rays.py new file mode 100644 index 0000000..eae18cb --- /dev/null +++ b/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()