NaN with raytracing scattering for mics close to the wall

[stiansel]

I've been running into an occasional NaN issue when simulating with the new ray tracing scattering. I'm using Pyroomacoustics 0.10.0. Managed to boil it down to the following example.

The NaN occurs when negative values reach the np.sqrt in rt.py at:

        # Impulse response for every octave band for each microphone.
        rir_bands[b, :] = seq_bp * np.sqrt(hist)

Only occurs when scattering materials are used with raytracing enabled and the receiver is closer to a wall than the receiver_radius.

Example where it occurs:

import warnings

import numpy as np
import pyroomacoustics as pra

for receiver_radius in [0.5, 0.1]:  # the 0.1 case leaves the mic farther from the wall than the receiver_radius
    for use_scattering in [False, True]:
        pra.random.seed(140)  # this seed seems to reproduce the NaN issue reliably on my machine.

        room = pra.ShoeBox(
            [4.144901752471924, 6.039876461029053, 2.7347633838653564],
            fs=24000,
            materials=pra.Material(energy_absorption=0.13035287, scattering=0.1 if use_scattering else 0.0),
            max_order=0,
            ray_tracing=True,
            air_absorption=True,
        )
        room.set_ray_tracing(receiver_radius=receiver_radius)

        # The microphone is fairly close to the wall. Within the default receiver_radius (0.5 m)
        room.add_microphone([2.304156303866167, 0.17025121891320288, 1.4909059943866196])
        room.add_source([1.9948741012396916, 3.8261834557486685, 0.8465485818395717])

        with warnings.catch_warnings(record=True) as w:
            warnings.simplefilter("always")
            room.compute_rir()

            sqrt_warns = [x for x in w if "invalid value encountered in sqrt" in str(x.message)]
            has_nan = np.any(np.isnan(room.rir[0][0]))
            print(f"RT scattering: {use_scattering}. Receiver radius: {receiver_radius}. Sqrt warn: {len(sqrt_warns)}, RIR has NaN: {has_nan}")

Script output:

RT scattering: False. Receiver radius: 0.5. Sqrt warn: 0, RIR has NaN: False
RT scattering: True. Receiver radius: 0.5. Sqrt warn: 1, RIR has NaN: True
RT scattering: False. Receiver radius: 0.1. Sqrt warn: 0, RIR has NaN: False
RT scattering: True. Receiver radius: 0.1. Sqrt warn: 0, RIR has NaN: False

[fakufaku]

Hi @stiansel , thank you for flagging this! I had seen the warning, but was not overly worried about it.
NaN's are not good. I'm wondering why we'd have negative numbers there and will investigate.

[stiansel]

Thanks for looking into it!

I noticed that scat_ray in room.cpp handles distance differently than in the simul_ray function.

pyroomacoustics/pyroomacoustics/libroom_src/room.cpp

Line 764 in c7365b8

float p_hit_equal = 1.f - sqrt(1.f - mic_radius_sq / h_sq);

vs

pyroomacoustics/pyroomacoustics/libroom_src/room.cpp

Lines 892 to 893 in c7365b8

	auto p_hit =
	(1 - sqrt(1 - mic_radius_sq / std::max(mic_radius_sq, r_sq)));

I'm not familiar with the codebase though, so might not be it at all.

[fakufaku]

@stiansel I think in this case #434 was the solution kindly provided by @jguice . Essentially, tiny negative values were introduced by the histogram interpolation. Just clamping at zeros solves the nan issues without other regressions.

Thanks for pointing out the two different computations of p_hit in the case of direct and scattered rays. I think that the name p_hit represents different quantities in these two functions. You can check the reference below for details. I think the form of the book is slightly different than what's in the code but should be equivalent.

Dirk Schröder, “Physically Based Real-Time Auralization of Interactive Virtual Environments.”

• p_hit_equal from scat_ray is described in eq. 5.15.
• p_hit from simul_ray comes from the equalization of image source model and ray tracing, as described in eq. 5.54.

But, the difference between the two may only be a factor of two. And, indeed, if there is scattering happening due to a wall reflection that falls within the detector boundary, then h_sq will be smaller than mic_radius_sq and we'll get some nan...
Good find. I will protect this with an std::max.

Edit: Just tested it and it does not seem to be the cause of the nan here. It is good to fix it regardless.