Correctness of particle tracer

[Dream4Leo]

In mitsuba, particle tracer samples an emitter ray, intersects the scene at a point SurfaceInteraction3f si and connects it to the sensor. The contribution of this light path is evaluated by connect_sensor() and splat to the sensor.

Inside connect_sensor(), the surface interaction is computed in bsdf->eval(), where si.wi would indicate direction to emitter, and wo to sensor, which is the opposite way to a path tracer.

This is where I start to get confused: bsdf->eval() multiplies the cosine forshortening term, and in particle tracer it means to multiply dot(n,wo), the cosine angle between emitter direction and normal. At first glance it seems to be incorrect.

I do know that for reciprocal BSDFs, light transport is reversible, so we can treat as if light emits from the sensor. And a simple test on cornell box indeed shows matching renderings between path tracer and particle tracer. But it is difficult to convince myself on the following example:

Consider a perspective camera at [0,0,1] above a diffuse surface z=0, and a directional emitter at grazing angle (say, dir=[0,-1,-0.01]), one-bounce. Path tracer will produce near-zero contributions (roughly 0.01/pi) for each sample, while particle tracer will have 1/pi. Therefore, I would expect particle tracer to render a much brighter image. But again both renderings match (not strictly but close mean disregarding the variance).

Did I miss something?

[ziyi-zhang]

Hi @Dream4Leo,
Great question! The following symmetry allows us to think about light tracing as the same as path tracing but with the two endpoints reversed -- actually in Mitsuba both emitter and sensor inherits from endpoint that defines functions like sample_rays/pdf_direction.
Hope this helps.

[Dream4Leo]

Thanks for the quick reply and just an additional check:
for non-reciprocal BSDFs in light tracing, should we swap wi and wo for physical correctness, and keep using cosine factor at tracing direction (bs.wo/sensor_ray.d)?

[shinyoung-yi]

Hi,
@Dream4Leo 's question might be related to the implementation convention of BSDFs according to mi.TransportMode.

Consider the following figure:

The equations in both formulations can be derived from the same logic as @ziyi-zhang 's equations.

If ctx.mode == mi.TransportMode.Radiance (default, path tracing), bsdf.eval(ctx, si, wo) indicates
(See also here. )

$$ \left(\text{BSDF value from wo to si.wi}\right) \times\cos\left(\text{wo}\right). $$

If instead, ctx.mode == mi.TransportMode.Importance (light tracing), bsdf.eval(ctx, si, wo) indicates
(See also here. )

$$ \left(\text{BSDF value from si.wi to wo}\right) \times\cos\left(\text{wo}\right). $$

This convention allows for nearly symmetric implementations for both path and light tracing.

Script for verification

You can check that Mitsuba 3's BSDF implementation follows the above conventions by running the following script.

import numpy as np
import matplotlib.pyplot as plt
import drjit as dr
import mitsuba as mi
mi.set_variant('scalar_rgb')
print(f"{dr.__version__ = }") # dr.__version__ = '0.4.6'
print(f"{mi.__version__ = }") # mi.__version__ = '3.5.2'

bsdf_diffuse = mi.load_dict({'type': 'diffuse',
                             'reflectance': {'type': 'rgb', 'value': [1, 1, 1]}})
bsdf_rd = mi.load_dict({'type': 'roughdielectric', 'alpha': 1.0})
bsdf_list = [bsdf_diffuse, bsdf_rd]

ang_list = [0, 60, -120]
w_list = [mi.Vector3f([dr.sin(dr.deg2rad(deg)), 0, dr.cos(dr.deg2rad(deg))]) for deg in ang_list]

mode_list = [mi.TransportMode.Radiance, mi.TransportMode.Importance]
for bsdf in bsdf_list:
    print("\n---------- BSDF ----------\n", bsdf)
    for i_m, mode in enumerate(mode_list):
        print(f"\n# mi.{mode}")
        ctx = mi.BSDFContext(mode)
        for i_wi, (wi, wi_label) in enumerate(zip(w_list, ang_list)):
            for i_wo, (wo, wo_label) in enumerate(zip(w_list, ang_list)):
                si = mi.SurfaceInteraction3f()
                si.wi = wi
                val = bsdf.eval(ctx, si, wo)
                print(f"[wi:{wi_label:4d}, wo:{wo_label:4d}] bsdf.eval(si.wi, wo) == {val}")

                # ---------- Test `bsdf.eval()` for diffuse ---------- 
                if bsdf is bsdf_diffuse:
                    bsdf_GT = 1 / dr.pi # Relationship between albedo and BRDF
                    if wo.z >= 0 and wi.z >= 0: # BRDF only
                        assert dr.allclose(val / dr.abs(wo.z), bsdf_GT)

                # ---------- Test `mi.TransportMode` ---------- 
                ctx_rev = mi.BSDFContext(mode_list[1-i_m])
                si_rev = mi.SurfaceInteraction3f()
                si_rev.wi = wo
                assert dr.allclose(val * dr.abs(wi.z), bsdf.eval(ctx_rev, si_rev, wi) * dr.abs(wo.z))

Running this should pass all assertions.

Instead of reporting all the printed message, I would like to highlight several parts which help our understanding about Mitsuba 3's implementation.

bsdf.eval() contains the cosine term for wo, regardless of ctx.mode

---------- BSDF ----------
 SmoothDiffuse[
  reflectance = SRGBReflectanceSpectrum[
    value = [1, 1, 1]
  ]
]
# mi.TransportMode.Radiance
[wi:   0, wo:   0] bsdf.eval(si.wi, wo) == [0.31830987334251404, 0.31830987334251404, 0.31830987334251404]
[wi:   0, wo:  60] bsdf.eval(si.wi, wo) == [0.15915493667125702, 0.15915493667125702, 0.15915493667125702]
[wi:  60, wo:   0] bsdf.eval(si.wi, wo) == [0.31830987334251404, 0.31830987334251404, 0.31830987334251404]
[wi:  60, wo:  60] bsdf.eval(si.wi, wo) == [0.15915493667125702, 0.15915493667125702, 0.15915493667125702]
(omitted...)

# mi.TransportMode.Importance
[wi:   0, wo:   0] bsdf.eval(si.wi, wo) == [0.31830987334251404, 0.31830987334251404, 0.31830987334251404]
[wi:   0, wo:  60] bsdf.eval(si.wi, wo) == [0.15915493667125702, 0.15915493667125702, 0.15915493667125702]
[wi:  60, wo:   0] bsdf.eval(si.wi, wo) == [0.31830987334251404, 0.31830987334251404, 0.31830987334251404]
[wi:  60, wo:  60] bsdf.eval(si.wi, wo) == [0.15915493667125702, 0.15915493667125702, 0.15915493667125702]
(omitted...)

In the codes, bsdf.eval() always includes the cosine term of wo,
though its meaning ($\hat\omega_i$ or $\hat\omega_o$) depends on ctx.mode.

Roll of ctx.mode, also tested for non-reciprocal BSDF

---------- BSDF ----------
 RoughDielectric[
  distribution = beckmann,
  sample_visible = 1,
  alpha = UniformSpectrum[value=1.000000],
  eta = 1.50418
]

# mi.TransportMode.Radiance
(omitted)
[wi:  60, wo:-120] bsdf.eval(si.wi, wo) == [1.4565513134002686, 1.4565513134002686, 1.4565513134002686]
[wi:-120, wo:  60] bsdf.eval(si.wi, wo) == [3.2955477237701416, 3.2955477237701416, 3.2955477237701416]
(omitted)

# mi.TransportMode.Importance
(omitted)
[wi:  60, wo:-120] bsdf.eval(si.wi, wo) == [3.2955451011657715, 3.2955451011657715, 3.2955451011657715]
[wi:-120, wo:  60] bsdf.eval(si.wi, wo) == [1.4565526247024536, 1.4565526247024536, 1.4565526247024536]
(omitted)

Note that zenith angles 60 and -120 have the same $\left|\cos\theta\right|$, so that in the above results we can compare BSDF values
themselves without any change of cosine terms

[Dream4Leo]

Thanks for the explanation on BSDF conventions. I wrote a minimum non-reciprocal BSDF based on my understanding and passed the script test you provided.