Fast Directional Radiosity

This examples shows a simple diffuse simulation of a shoebox room and compare the the result with the analytical solution. First we import all our dependencies.

"""Test the radiosity.Radiosity module."""
import numpy as np
import pyfar as pf
import sparrowpy as sp
import matplotlib.pyplot as plt

%matplotlib inline
# %matplotlib ipympl

Lets define our room and source position.

# Define parameters
X = 5
Y = 6
Z = 4
patch_size = 1
etc_duration = 1
etc_time_resolution = 1/1000
max_reflection_order = 150
speed_of_sound = 343.2
absorption = 0.1

# create geometry
walls = sp.testing.shoebox_room_stub(X, Y, Z)
source = pf.Coordinates(2, 2, 2)
receiver = pf.Coordinates(2, 3, 2)

Let’s create the instance of the simulation class called DirectionalRadiosityFast from the previously defined walls.

# create object
radiosity_fast = sp.DirectionalRadiosityFast.from_polygon(walls, patch_size)

Now we define the scattering coefficient to be 1 for each wall:

# create directional scattering data (totally diffuse)
brdf_sources = pf.Coordinates(0, 0, 1, weights=1)
brdf_receivers = pf.Coordinates(0, 0, 1, weights=1)
frequencies = np.array([1000])
brdf = sp.brdf.create_from_scattering(
    brdf_sources,
    brdf_receivers,
    pf.FrequencyData(1, frequencies),
    pf.FrequencyData(absorption, frequencies))

# set directional scattering data
radiosity_fast.set_wall_brdf(
    np.arange(len(walls)), brdf, brdf_sources, brdf_receivers)

/var/folders/7q/_50y0lw50_x1j732_5yyyx7h0000gn/T/ipykernel_67300/1682180744.py:12: UserWarning: Gimbal lock detected. Setting third angle to zero since it is not possible to uniquely determine all angles.
  radiosity_fast.set_wall_brdf(

The air attenuation and absorption can be defined as well.

# set air absorption
radiosity_fast.set_air_attenuation(
    pf.FrequencyData(
        np.zeros_like(brdf.frequencies),
        brdf.frequencies))

Now the simulation start, first the geometry is baked, where all patch to patch relationships are precalculated. E.g. Form-factor including the geometrical relationship between the patches as well as the BRDFs. This is the most heavy part of the simulation and is independent of the source and receiver position.

# calculate from factors including brdfs
radiosity_fast.bake_geometry()

OMP: Info #270: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
/Users/anne/git/_pyfar/sparrowpy/sparrowpy/form_factor/universal.py:48: NumbaPerformanceWarning: np.dot() is faster on contiguous arrays, called on (Array(float64, 1, 'A', False, aligned=True), Array(float64, 1, 'C', False, aligned=True))
  form_factors[i,j] = universal_form_factor(

it is also possible to write and read the FastRadiosity object at andy state of the simulation pipeline.

radiosity_fast.write('radiosity_fast.far')
radiosity_read = sp.DirectionalRadiosityFast.from_read('radiosity_fast.far')

We can now also continue with the object, which we have just read from disk. Next, the source energy is initialized.

radiosity_read.init_source_energy(source)

Next we calculate the energy exchange between the patches, then we can collect the energy at the receiver.

radiosity_read.calculate_energy_exchange(
        speed_of_sound=speed_of_sound,
        etc_time_resolution=etc_time_resolution,
        etc_duration=etc_duration,
        max_reflection_order=max_reflection_order)

The energy is collected at the receiver, this is quite fast and can be done for as many receivers as required.

etc_radiosity = radiosity_read.collect_energy_receiver_mono(
    receivers=receiver)

print(etc_radiosity)

TimeData:
(1, 1) channels with 1000 samples

To compare our energy time curve, we can calculate the analytical solution based on the diffuse sound field in the room after Kuttruff.

S = (2*X*Y) + (2*X*Z) + (2*Y*Z)
A = S*absorption
alpha_dash = A/S
r_h = 1/4*np.sqrt(A/np.pi)
print(f'reverberation distance is {r_h:.1f}m')
V = X*Y*Z
RT = 24*np.log(10)/(speed_of_sound)*V/(-S*np.log(1-alpha_dash))
print(f'reverberation time is {RT:.2f}s')
E_reverb_analytical = 4/A
t = etc_radiosity.times
# Kuttruff Eq 4.7
w_0 = E_reverb_analytical/ V
t_0 = 0.03
# Kuttruff Eq 4.10
reverberation_analytic = w_0 * np.exp(+(
    speed_of_sound*S*np.log(1-alpha_dash)/(4*V))*(t-t_0))
reverberation_analytic = pf.TimeData(reverberation_analytic, t)

reverberation distance is 0.5m
reverberation time is 1.24s

Lets compare these to results.

plt.figure()
pf.plot.time(
    reverberation_analytic, dB=True, log_prefix=10,
    label=f'analytical E_rev={E_reverb_analytical:0.2f}')

pf.plot.time(
    etc_radiosity, dB=True, log_prefix=10,
    label='simulated',
    linestyle='--')

plt.legend()
plt.show()

If the order would be increased, this would match even after 0.5s.

%load_ext watermark
%watermark -v -m -iv

Python implementation: CPython
Python version       : 3.11.11
IPython version      : 9.1.0

Compiler    : Clang 14.0.6
OS          : Darwin
Release     : 23.5.0
Machine     : arm64
Processor   : arm
CPU cores   : 8
Architecture: 64bit

pyfar     : 0.7.2
sparrowpy : 0.1.0
matplotlib: 3.10.1
numpy     : 1.26.4