[🐛 bug report] Incorrect radiance values using the `radiancemeter` plugin

[leroyvn]

Summary

I get wrong radiance values with a simple setup consisting of a homogeneous volume above a flat surface.

System configuration

• Platform: macOS 12.1
• Compiler: Apple clang version 13.0.0 (clang-1300.0.29.30)
• Python version: 3.8.10
• Mitsuba 2 version:
• Compiled variants:
  • scalar_mono
  • scalar_mono_double
  • scalar_rgb
  • scalar_rgb_double
  • scalar_spectral
  • llvm_rgb
• Embree: OFF

Description

I'm simulating radiance values with a simple setup consisting of a homogeneous volume above a flat surface, illuminated with a directional emitter (see monochromatic illustration below).

The sensor setup consists of a set of radiancemeter plugins targeting the [0, 0, 1] point (the top of the participating medium), placed in an arc above the volume at a reasonable distance (order of magnitude 1) from the target point. I get incorrect values for some of those radiance meters. Example:

On this figure, the sensor index maps to the local zenith angle, from -75° to +75°.

Sensor configurations leading incorrect values to appear seem to depend on:

• the active variant;
• sensor origin (changing the list of directions in the script below);
• sensor target;
• volume geometry (changing the offset parameter in the script below).

Steps to reproduce

Run the script below:

# %%
# Imports

import numpy as np
import matplotlib.pyplot as plt
import mitsuba

# %%
# Utility functions


def cos_angle_to_direction(cos_theta, phi):
    cos_theta = np.atleast_1d(cos_theta).astype(float)
    phi = np.atleast_1d(phi)

    sin_theta = np.sqrt(1.0 - np.multiply(cos_theta, cos_theta))
    sin_phi, cos_phi = np.sin(phi), np.cos(phi)
    return np.vstack((sin_theta * cos_phi, sin_theta * sin_phi, cos_theta)).T


def angles_to_direction(angles):
    angles = np.atleast_1d(angles).astype(float)
    if angles.ndim < 2:
        angles = angles.reshape((angles.size // 2, 2))
    if angles.ndim > 2 or angles.shape[1] != 2:
        raise ValueError(f"array must be of shape (N, 2), got {angles.shape}")

    negative_zenith = angles[:, 0] < 0
    angles[negative_zenith, 0] *= -1
    angles[negative_zenith, 1] += np.pi

    return cos_angle_to_direction(np.cos(angles[:, 0]), angles[:, 1])


def map_cube(xmin, xmax, ymin, ymax, zmin, zmax):
    from mitsuba.core import ScalarTransform4f

    scale = ScalarTransform4f.scale(
        [
            0.5 * (xmax - xmin),
            0.5 * (ymax - ymin),
            0.5 * (zmax - zmin),
        ]
    )
    translate = ScalarTransform4f.translate(
        [
            xmin + 0.5 * (xmax - xmin),
            ymin + 0.5 * (ymax - ymin),
            zmin + 0.5 * (zmax - zmin),
        ]
    )
    return translate * scale


# %%
# Scene configuration

# Participating medium transform setup
width = 1.0  # 1 unit wide
bottom = 0.0  # surface level at z = 0
top = 1.0  # top at z = 1
offset = 0.1  # further extend the shape underneath surface


def cube_trafo():
    return map_cube(
        xmin=-0.5 * width,
        xmax=0.5 * width,
        ymin=-0.5 * width,
        ymax=0.5 * width,
        zmin=bottom - offset,
        zmax=top,
    )


# Sensor directions
angles = np.deg2rad([[theta, 0] for theta in np.arange(-75, 76, 5)])
directions = -angles_to_direction(angles)

# %%
# Scene setup


def radiancemeter_dict(direction, target=(0, 0, 0), distance=2.0):
    from mitsuba.core import ScalarTransform4f, coordinate_system

    target = np.array(target)
    origin = target - distance * direction
    _, up = coordinate_system(direction)

    return {
        "type": "radiancemeter",
        "to_world": ScalarTransform4f.look_at(
            origin=origin,
            target=target,
            up=up,
        ),
        "film": {
            "type": "hdrfilm",
            "width": 1,
            "height": 1,
            "rfilter": {"type": "box"},
            "pixel_format": "luminance",
        },
    }


def load_scene():
    from mitsuba.core import ScalarTransform4f, load_dict

    scene_dict = {
        "type": "scene",
        "surface_bsdf": {
            "type": "diffuse",
            "reflectance": 1.0,
        },
        "surface_shape": {
            "type": "rectangle",
            "bsdf": {"type": "ref", "id": "surface_bsdf"},
            "to_world": ScalarTransform4f.scale([0.5, 0.5, 1]),
        },
        "atmosphere_phase": {"type": "rayleigh"},
        "atmosphere_medium": {
            "type": "homogeneous",
            "sigma_t": 1.0,
            "albedo": 0.5,
            "phase": {"type": "ref", "id": "atmosphere_phase"},
        },
        "atmosphere_shape": {
            "type": "cube",
            "to_world": cube_trafo(),
            "interior": {"type": "ref", "id": "atmosphere_medium"},
            "bsdf": {"type": "null"},
        },
        "illumination": {
            "type": "directional",
            "direction": [0, 0, -1],
        },
        "integrator": {"type": "volpathmis"},
    }

    # Add sensors
    scene_dict["camera"] = {
        "type": "perspective",
        "to_world": ScalarTransform4f.look_at(
            origin=[2, 2, 2], target=[0, 0, 0.4], up=[0, 0, 1]
        ),
        "film": {
            "type": "hdrfilm",
            "width": 640,
            "height": 480,
            "pixel_format": "luminance",
        },
    }

    for i, direction in enumerate(directions):
        scene_dict[f"radiancemeter_{str(i).zfill(3)}"] = radiancemeter_dict(
            direction, target=[0, 0, 1], distance=1.6
        )

    return load_dict(scene_dict)


# %%
# Plot radiancemeter output

plt.figure()

for variant in [
    "scalar_mono",
    "scalar_mono_double",
    "scalar_rgb",
    "scalar_rgb_double",
    "scalar_spectral",
    # "llvm_rgb",
]:
    mitsuba.set_variant(variant)
    scene = load_scene()
    result_radiancemeter = np.empty((directions.shape[0],))

    for i in range(1, len(scene.sensors())):
        radiance = np.squeeze(scene.render(spp=10000, seed=i, sensor_index=i))
        result_radiancemeter[i - 1] = radiance

    with plt.style.context({"lines.linestyle": ":", "lines.marker": "."}):
        plt.plot(result_radiancemeter, label=variant)

plt.xlabel("Sensor index")
plt.ylabel("Radiance")
plt.legend()

plt.show()
plt.close()

# %%
# Render an image with the camera

mitsuba.set_variant("scalar_mono")
scene = load_scene()
img = np.array(scene.render(spp=32, sensor_index=0)).squeeze()
plt.imshow(img)
plt.colorbar()
plt.title("scalar_mono")
plt.show()
plt.close()

[wjakob]

Just to clarify: there is some issue with scalar_rgb_double, correct? (and the other radiancemeters give the expected result?) Is this a new bug? (i.e. something that worked previously).

[leroyvn]

I have this with all scalar variants I tried (I can also double-check with LLVM variants if you'd like). I can't tell you whether this specific example actually ever worked, I put it together specifically for this issue.

Some context. We've been having this issue since months with Eradiate. However we didn't report it because we were using the old codebase: we didn't want to waste your time with the old code.

We discussed something similar a few months ago and this led you guys to fixing intersection issues with the rectangle and sphere plugins, IIRC. I tried backporting those fixes to the Eradiate kernel version of Mitsuba, but it didn't fix it. Hoping that the quality of that port would be the problem, I decided to wait and further investigate when we'd align with the Next repo.

We originally spotted that with a modified distant sensor which basically aggregates an arbitrary number of distant radiancemeters. This is what's been buggy since, well, at least 4 months. I couldn't be sure that this wasn't on us (bad sensor implementation); now that we're moving to the new codebase, I can ensure that the MWE setup is exactly the same as the one which produces the bugs in Eradiate.

[leroyvn]

I forgot a very important point: incorrect value locations also move when I stretch the volume (varying the offset parameter in the script). I'll update the issue description.

[leroyvn]

EDIT: Answering you questions properly.

• All scalar variants are problematic somehow, with a 0 at the vertical orientation and moving incorrect values depending on the factors mentioned in the description (variant, sensor position and orientation, volume transform).
• I don't know if it's a new bug, all I can say is that we've had similar issues with the old Mitsuba code since months.

[leroyvn]

I'm coming back to this issue: I updated the sample script above with the recent changes (see this gist).

The issue can still be observed with all scalar double precision variants:

I couldn't try LLVM variants because of #92.

[Speierers]

It would be great to better understand what causes this discrepancy. Maybe could you try with a simplified setup, e.g. using a constant emitter and isotropic phase function. Also you could try removing the floor plane, to see whether this "extra energy" comes from the scattering in the volume or the light bouncing on that plane.

Also trying with LLVM would be interesting.

Let me know whether this helped pinpointing where the problem is.

🦇 🥷

[leroyvn]

Thanks to you solving #92, I can now run this with the llvm_rgb variants and I get the same result:

The outlier point in this "surface, Rayleigh, directional" has the value of the reflected radiance without the participating medium (1/π).

I ran a bunch of different configurations ranging from my original configuration to what you requested (see updated script), i.e. "no surface, isotropic, constant":

There, the outlier takes the value of the emitted radiance (1).

---

My interpretation is that the ray intersection somehow misses the cube and directly hits what's behind (either the surface or nothing depending on the configuration). The random walk consequently never enters the medium, hence the radiance values we get.

[leroyvn]

Varying the surface reflectance:

No surface, isotropic, constant (radiance without medium: 1.0)

Surface (ρ=0.5), isotropic, constant (radiance without medium: 0.5)

Surface (ρ=1.0), isotropic, constant (radiance without medium: 1.0)

(Just making sure I'm picking up the reflected radiance value when there is a surface, not the emitter.)

[leroyvn]

Other things I tried:

• I changed the integrator. volpath and volpathmis produce the same results.
• I varied the position of the lower face of the cube (i.e. extending it further down the negative Z values). It produces the same results.
• I varied the position of the upper face of the cube. It produces the same results.
• I varied the distance between the radiancemeter's location and the target point. This changes the results: the outliers move around, appear and disappear (I have the feeling that the further the radiancemeter to its target, the more frequent they get).

[Speierers]

Could you hack in the integrator and make it so that it returns 1.0 when doing the ray marching, and 0.0 otherwise? Just to validate your assumtion above. If that's the case, then you should be able to isolate the ray tracing call that misses the cube and better understand what's going on.

[leroyvn]

Alright, I did what you suggested and it's becoming more and more strange. When I run my script and chain both the single and double precision variants, I get this:

In that case, the ray cast by the problematic sensor misses the front face and hits the back face (I expect ray.o.z to be equal to 1):

2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator] ray_ = Ray3d[
2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator]   o = [-0.547232, 6.70166e-17, 2.50351],
2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator]   d = [0.34202, -4.18854e-17, -0.939693],
2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator]   maxt = 1.79769e+308,
2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator]   time = 0,
2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator] ]
2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator] ray = Ray3d[
2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator]   o = [0.400367, -2.498e-16, -0.1],
2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator]   d = [0.34202, -4.18854e-17, -0.939693],
2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator]   maxt = 1.79769e+308,
2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator]   time = 0,
2022-02-26 14:20:40 WARN  main  [VolumetricPathIntegrator] ]

However, this issue disappears when I run only the double precision variant:

I made sure that RNG seed values are reproducible independently from variant selection so I wouldn't expect to be the cause of a behaviour variation.

[Speierers]

Interesting. So it seems like the scene instance of the double-variant is interfering with the scene instance of thesingle-variant. Could you maybe try to explicitly delete and garbage collect the scene at every loop iteration in your script? E.g. using gc.collect(); gc.collect()

[leroyvn]

I added a del scene; gc.collect(); gc.collect() to my loop on variants, but this doesn't seem to do much: I still get the same results.

[Speierers]

Okay, I can take a closer look at this issue. Could you share a minimal reproducer that doesn't even use an integrator, but rather do everything in python, e.g. trace a ray and return 1.0 or 0.0. Then I will continue the investigation on my side. Thanks a lot for digging this up!

[Speierers]

Pinging @leroyvn here. Is this issue still relevant, or was that bug fixed in the latest codebase?

[leroyvn]

Hi @Speierers, I'd need some time to get back to this. I couldn't reproduce this issue without the integrator at the time, and I haven't worked on it recently so I don't know if it still happens with the latest codebase. If you wish to close the issue, I can reopen it later if necessary.

[Speierers]

No worries, let's keep it open for now