Doppler ToF velocity (offline)

This tutorial simulates a Doppler time-of-flight camera with DopplerToFPathTracerInline and estimates the velocity of the Cornell box’s boxes from its measurements. The tall box comes towards the camera and the short box moves away, both at 2 m/s; the light is modulated at 30 MHz, and the exposure lasts 15 ms, over which the boxes really move by 3 cm. The estimate is compared with three ground-truth velocities from VelocityGroundTruthInline.

It uses the movable Cornell box described in Doppler rendering.

Estimated D-ToF velocity next to the ground-truth D-ToF, OHD and z velocities

1. Load the scene

The boxes are built as animated in scene.pyscene (addTriangleMesh(..., True)), which keeps each in its own node and acceleration structure so that the tracer can move it. DontOptimizeGraph keeps the scene graph’s nodes apart, so that the two boxes can move differently.

import falcorcomp as falcor
import numpy as np
import matplotlib
matplotlib.use("Agg")  # no window
import matplotlib.pyplot as plt

testbed = falcor.Testbed(create_window=False)
# The boxes are built as animated (see scene.pyscene); keep every object separate and named.
testbed.load_scene("cornell-box-moving/scene.pyscene",
                   falcor.SceneBuilderFlags.DontMergeMaterials | falcor.SceneBuilderFlags.DontOptimizeGraph)
testbed.resize_frame_buffer(256, 256)
testbed.scene.camera.aspectRatio = 1.0
testbed.clock.pause()

# The tall box comes towards the camera (+z), the short box moves away, at 2 m/s.
VELOCITIES = {"TallBox": {"linear": [0.0, 0.0, 2.0]}, "ShortBox": {"linear": [0.0, 0.0, -2.0]}}
F_G, T = 30.0, 0.015  # light modulation (MHz) and exposure (s)
LIGHT = {"isLightSourceLaser": False, "laserCollocated": True, "laserPower": [10.0, 10.0, 10.0]}

2. Render a Doppler ToF measurement

Every frame, the tracer draws a time in the exposure, moves the boxes there, and traces the paths. The measurement is the mean over the exposure, so AccumulatePass averages the frames; it is created with autoReset off, because the scene changes every frame.

  • Heterodyne (\(\Delta f = 1/T\)): light from static objects cancels over the exposure. Every frame also renders the time half a heterodyne period later, where the static light has the opposite sign, with the same random numbers (antithetic = half_period, randomReplay), so that it cancels within each frame (see Antithetic time pairs).

  • Homodyne (\(\Delta f = 0\)): an ordinary CW-ToF measurement over the exposure. It has no half period to pair with, so each frame renders one time.

timeSampling = stratified spreads the frames’ times evenly over the exposure. All pixels of a frame share its time, so with random times the error of the time integral is the same over the whole image.

def render(heterodyne_frequency, phase, renders):
    graph = testbed.create_render_graph("DopplerToF")
    graph.create_pass("Light", "LaserLight", LIGHT)
    properties = {
        "samplesPerPixel": 16, "maxBounces": 3, "computeDirect": True,
        "modulationFrequency": F_G, "heterodyneFrequency": heterodyne_frequency,
        "exposureTime": T, "phase": phase, "timeSampling": "stratified",
        "velocities": VELOCITIES,
    }
    if heterodyne_frequency != 0.0:
        # Pair every time t with t + 1 / (2 df), rendered with the same random numbers.
        properties.update({"antithetic": "half_period", "randomReplay": True})
    else:
        properties["antithetic"] = "none"  # without a heterodyne frequency there is no half period
    graph.create_pass("Tracer", "DopplerToFPathTracerInline", properties)
    # The tracer moves the scene every frame: keep AccumulatePass from restarting on scene changes.
    graph.create_pass("Accumulate", "AccumulatePass", {"precisionMode": "SingleCompensated", "autoReset": False})
    graph.add_edge("Light", "Tracer")
    graph.add_edge("Tracer.color", "Accumulate.input")
    graph.mark_output("Accumulate.output")
    testbed.render_graph = graph
    frames = renders // (2 if heterodyne_frequency != 0.0 else 1)  # a paired frame is two renders
    for _ in range(frames):
        testbed.frame()
    return graph.get_output("Accumulate.output").to_numpy()[..., :3].mean(-1)

3. Heterodyne and homodyne, at two sensor phases

Each measurement is rendered at sensor phases 0 and 0.25 and stored as a complex number, \(C = I_0 + i\, I_{0.25}\).

RENDERS = 4096
heterodyne = render(1.0 / T, 0.0, RENDERS) + 1j * render(1.0 / T, 0.25, RENDERS)
homodyne = render(0.0, 0.0, RENDERS) + 1j * render(0.0, 0.25, RENDERS)

4. Velocity from the ratio

For a path whose length changes at the rate \(\mathrm{d}\ell/\mathrm{d}t\), the ratio of the two measurements is \(\rho = \varepsilon T / (1 + \varepsilon T)\) with \(\varepsilon = -f_g (\mathrm{d}\ell/\mathrm{d}t)/c\), which gives the velocity \(v = -\frac{1}{2}\,\mathrm{d}\ell/\mathrm{d}t\) along the pixel ray. With one phase only, the ratio fails where the homodyne image crosses zero; the second phase avoids that.

# A path whose length changes at rate dl/dt has the frequency shift e = -f_g (dl/dt) / c, and
# heterodyne / homodyne = eT / (1 + eT). So eT = rho / (1 - rho) and v = -(1/2) dl/dt = c e / (2 f_g).
C = 299792458.0
valid = np.abs(homodyne) > 0.05 * np.percentile(np.abs(homodyne), 90)
rho = (heterodyne / np.where(valid, homodyne, np.nan)).real
velocity = C * (rho / (1 - rho) / T) / (2 * F_G * 1e6)
np.save("doppler_tof_velocity.npy", velocity)

5. Compare with the ground truth

VelocityGroundTruthInline renders three velocities of the direct path, which differ:

  • D-ToF velocity (path_length): how fast the path length along the fixed pixel ray changes, between the scene now and 1 ms later. This is what the measurement estimates. The pixel ray meets a moving face at a different point over time, so on the boxes’ sides, seen at a grazing angle, it is much larger.

  • OHD velocity (doppler): the velocity of the hit surface point towards the camera, as in the Doppler spectrum tutorial. It falls off as the cosine of the angle between the pixel ray and the motion.

  • z velocity (projection): the boxes’ velocity along the camera axis, \(\pm 2\) m/s.

def ground_truth(properties):
    graph = testbed.create_render_graph("Truth")
    graph.create_pass("Light", "LaserLight", LIGHT)
    graph.create_pass("Truth", "VelocityGroundTruthInline", {"velocities": VELOCITIES, **properties})
    graph.add_edge("Light", "Truth")
    graph.mark_output("Truth.velocity")
    testbed.render_graph = graph
    testbed.frame()
    return graph.get_output("Truth.velocity").to_numpy().reshape(velocity.shape)


panels = [
    (np.where(valid, velocity, np.nan), "estimated: heterodyne / homodyne"),
    (ground_truth({"mode": "path_length", "dt": 1e-3}), "D-ToF velocity: -(1/2) dl/dt"),
    (ground_truth({"mode": "doppler"}), "OHD velocity: u / 2"),
    (ground_truth({"mode": "projection", "direction": [0.0, 0.0, 1.0]}), "z velocity"),
]
figure, axes = plt.subplots(1, 4, figsize=(17, 4.3))
for axis, (image, title) in zip(axes, panels):
    shown = axis.imshow(image, cmap="bwr", vmin=-3, vmax=3)
    axis.set_title(title)
    axis.set_xticks([])
    axis.set_yticks([])
figure.colorbar(shown, ax=axes, fraction=0.012, label="m/s")
figure.savefig("doppler_tof_velocity.png", dpi=90, bbox_inches="tight")

On the boxes’ front faces the D-ToF velocity is +1.98 m/s and -1.99 m/s (medians). The estimate reads +2.03 m/s and -2.07 m/s, 2.6 % and 4.3 % too large, and it does not change with more renders. The ratio assumes each path’s intensity stays constant during the exposure, but the boxes’ brightness changes as they move towards or away from the light. Relative to the Doppler term, that change shrinks as the modulation frequency grows: at 300 MHz, the estimate is within 0.5 % of the D-ToF velocity.

The full script: doppler_tof_offline.py. It takes about 20 seconds.