Speckle contrast from the Doppler spectrum¶

Laser speckle contrast imaging (LSCI) images blood flow. Coherent light scattered by tissue forms a speckle pattern; moving scatterers make it fluctuate, which blurs it over the camera’s exposure and lowers its contrast $K = \sigma_I / \langle I \rangle$. Slow flow keeps the speckle sharp, fast flow washes it out.

The LSCI camera records intensity, without a local oscillator: it is homodyne (self-beating) detection, integrated over the exposure. This page computes the speckle contrast from the Doppler spectrum that DopplerHistogramPathTracerInline renders, the same spectrum as in the Doppler spectrum tutorial, with no other render pass.

It uses the scene described in Doppler rendering.

$Static light fraction, speckle contrast at three exposures and contrast against exposure for four pixels$

Speckle contrast from the Doppler spectrum¶

The pass renders $S(f)$, the light that reaches a pixel at each Doppler shift $f$: the power spectrum of the optical field, relative to the laser frequency. Without a local oscillator, the photocurrent follows the intensity, and for fully developed speckle (a Gaussian field) its power spectrum is the autocorrelation of the optical one, as derived for homodyne detection in Kim et al. (2025):

\[ S_i(f) = \int S(f')\, S(f' + f)\, df' . \]

The camera integrates the intensity over the exposure $T$, which filters the fluctuations by $\operatorname{sinc}^2(fT)$, with $\operatorname{sinc}(x) = \sin(\pi x) / (\pi x)$. The speckle contrast is then

\[ K^2(T) = \beta\, \frac{\int S_i(f)\, \operatorname{sinc}^2(fT)\, df}{I^2}, \qquad I = \int S(f)\, df , \]

where $\beta$ is the coherence factor, 1 when a pixel is smaller than a speckle. This is the frequency form of the usual model, $K^2 = \frac{2\beta}{T} \int_0^T \left(1 - \frac{\tau}{T}\right) |g_1(\tau)|^2\, d\tau$ (Bandyopadhyay et al. 2005), with $g_1$ the Fourier transform of $S$, normalized. With the spectrum in bins of width $\Delta f$, $S_k$ the light in bin $k$:

\[ K^2(T) = \beta\, \frac{\sum_m A[m]\, \operatorname{sinc}^2(m\, \Delta f\, T)}{I^2}, \qquad A[m] = \sum_k S_k\, S_{k+m} , \]

which holds for exposures well below $1 / \Delta f$.

Light from static surfaces has no shift and never decorrelates: for long exposures, $K$ tends to $\sqrt{\beta}$ times the static fraction of the light. Its beat with the moving light is the term that multi-exposure speckle models add for static scattering (Parthasarathy et al. 2008), where static light acts as a local oscillator; the autocorrelation contains it with nothing added.

1. Load the scene¶

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)
# Keep every material separate, so that each object keeps its name (identical materials are merged by default).
testbed.load_scene("cornell-box/scene-v4-nolight.pbrt", falcor.SceneBuilderFlags.DontMergeMaterials)
testbed.resize_frame_buffer(128, 128)
testbed.scene.camera.aspectRatio = 1.0
testbed.clock.pause()

2. Build the render graph¶

As in the Doppler spectrum tutorial, with a point light at the camera, which lights the whole scene as LSCI does. The motion is slow, as blood flow is: at 785 nm, the tall box coming towards the camera at 1 mm/s shifts light that it reflects once by $2v/\lambda \approx 2.5$ kHz. The spectrum covers -8 to 8 kHz in 1024 bins of 15.6 Hz, fine enough for exposures up to 10 ms ($1 / \Delta f = 64$ ms).

computeDirect is off. LSCI looks into tissue, where light scatters many times, and the speckle decorrelates because those paths have different shifts. A single reflection off a rigidly moving surface has a single shift and alone would not decorrelate (see Limitations). One wavelength needs one channel: useSingleChannel.

F_MAX, BINS = 0.008, 1024  # Doppler shifts from -8 to 8 kHz (given in MHz), in bins of 15.6 Hz

graph = testbed.create_render_graph("SpeckleContrast")
graph.create_pass("VBuffer", "VBufferRT", {"samplePattern": "Center", "sampleCount": 1})
# A point light at the camera: wide-field illumination.
graph.create_pass("Light", "LaserLight", {
    "isLightSourceLaser": False, "laserCollocated": True, "laserPower": [10.0, 10.0, 10.0],
})
tracer = graph.create_pass("Tracer", "DopplerHistogramPathTracerInline", {
    "samplesPerPixel": 64, "maxBounces": 4, "computeDirect": False,
    "wavelength": 785.0,  # nm
    "frequencyMin": -F_MAX, "frequencyMax": F_MAX, "frequencyBin": BINS,
    # Slow motion, as of blood flow: the tall box comes towards the camera at 1 mm/s, and the short box turns about
    # its vertical axis at 3 mrad/s (about 1 mm/s at its edges).
    "velocities": {
        "TallBox": {"linear": [0.0, 0.0, 1e-3]},
        "ShortBox": {"angular": [0.0, 3e-3, 0.0], "center": [0.329, 0.3, 0.375]},
    },
    "accumulate": True, "useSingleChannel": True,  # one wavelength: one channel
})
graph.add_edge("VBuffer.vbuffer", "Tracer.vbuffer")
graph.add_edge("VBuffer.viewW", "Tracer.viewW")
graph.add_edge("Light", "Tracer")
graph.mark_output("Tracer.spectrum")
testbed.render_graph = graph

3. Render the spectrum twice¶

$A[m]$ multiplies the spectrum by itself, and a spectrum multiplied by itself also squares its Monte Carlo noise, which would raise the contrast. Two spectra with independent samples, one for each factor, avoid that bias. reset() starts a new accumulation, and the new frames draw new samples.

FRAMES = 128  # 128 frames x 64 spp per spectrum
BIN_MHZ = 2 * F_MAX / BINS


def render_spectrum():
    """The light in each bin, (bins, pixels): a spectrum with new samples on every call."""
    tracer.reset()
    for _ in range(FRAMES):
        testbed.frame()
    spectrum = graph.get_output("Tracer.spectrum").to_numpy()  # per MHz, summed over the frames
    return (spectrum / FRAMES * BIN_MHZ).reshape(BINS, -1)


light_a = render_spectrum()
light_b = render_spectrum()

4. Compute the speckle contrast¶

$A[m]$ is computed with FFTs, a block of pixels at a time, and the sum over $m$ is a product with the matrix of $\operatorname{sinc}^2(m\, \Delta f\, T)$ for 41 exposures from 1 µs to 10 ms. The static fraction is the light in the bin of zero shift.

BETA = 1.0  # coherence factor: 1 when a pixel is smaller than a speckle
EXPOSURES = np.logspace(-6, -2, 41)  # 1 us to 10 ms
lags = np.concatenate([np.arange(BINS), np.arange(-BINS, 0)])  # m, in the order of a length-2B FFT
# sinc^2(m df T) for every exposure T and lag m (np.sinc(x) is sin(pi x) / (pi x)).
window = np.sinc(np.outer(EXPOSURES, lags * BIN_MHZ * 1e6)) ** 2
intensity2 = light_a.sum(0, dtype=np.float64) * light_b.sum(0, dtype=np.float64)  # I^2, from the two renders
lit = intensity2 > 0

k2 = np.empty((len(EXPOSURES), intensity2.size))
for start in range(0, intensity2.size, 2048):  # 2048 pixels at a time, to bound the memory
    pixels = slice(start, start + 2048)
    # Homodyne PSD: the autocorrelation of the spectrum, A[m] = sum_k S_k S_(k+m), with S from one render and
    # S_(k+m) from the other, so that the Monte Carlo noise of each does not add up.
    spectrum_a = np.fft.rfft(light_a[:, pixels], 2 * BINS, axis=0)
    spectrum_b = np.fft.rfft(light_b[:, pixels], 2 * BINS, axis=0)
    psd = np.fft.irfft(np.conj(spectrum_a) * spectrum_b, 2 * BINS, axis=0)
    psd = 0.5 * (psd + psd[(-lags) % (2 * BINS)])  # A[m] = A[-m]
    k2[:, pixels] = window @ psd
contrast = np.sqrt(np.clip(BETA * k2 / np.where(lit, intensity2, 1.0), 0.0, None))
contrast[:, ~lit] = np.nan
contrast = contrast.reshape(len(EXPOSURES), 128, 128)
# Static light: the bin of zero shift (0 is the lower edge of bin BINS / 2).
static = np.sqrt(np.clip(light_a[BINS // 2] * light_b[BINS // 2] / np.where(lit, intensity2, 1.0), 0.0, None))
static = np.where(lit, static, np.nan).reshape(128, 128)
np.savez("lsci_basics.npz", exposures=EXPOSURES, contrast=contrast, static_fraction=static)

5. Show the contrast¶

For short exposures $K = 1$: the speckle is fully developed. It starts to decorrelate around 0.1 ms, the inverse of the kHz spread of the shifts. All the light of the boxes is shifted, and their contrast keeps falling (0.2 at 10 ms). The walls level off at their static fraction: 0.80 on the back wall and 0.64 on the right wall, whose light partly came by way of the boxes.

pixels = {"tall box": (70, 45), "short box": (100, 85), "back wall": (30, 64), "right wall": (60, 120)}
shown_exposures = [2e-4, 1e-3, 1e-2]  # s

figure, axes = plt.subplots(1, 5, figsize=(20, 3.9), gridspec_kw={"width_ratios": [1, 1, 1, 1, 1.5]})
image = axes[0].imshow(static, cmap="viridis", vmin=0, vmax=1)
axes[0].set_title("static light fraction")
figure.colorbar(image, ax=axes[0], fraction=0.046)
for axis, exposure in zip(axes[1:4], shown_exposures):
    image = axis.imshow(contrast[np.argmin(np.abs(EXPOSURES - exposure))], cmap="magma", vmin=0, vmax=1)
    axis.set_title(f"contrast K, T = {exposure * 1e3:g} ms")
    figure.colorbar(image, ax=axis, fraction=0.046)
for (name, (y, x)), color in zip(pixels.items(), ["C0", "C1", "C2", "C3"]):
    for axis in axes[:4]:
        axis.plot(x, y, "o", mfc="none", mec=color, mew=2)
    axes[4].semilogx(EXPOSURES * 1e3, contrast[:, y, x], color=color, label=name)
    axes[4].axhline(np.sqrt(BETA) * static[y, x], color=color, ls="--", lw=0.8)
axes[4].set_xlabel("exposure T (ms)")
axes[4].set_ylim(0, 1.05)
axes[4].set_title("contrast K; dashed: static light fraction")
axes[4].legend(loc="lower left")
for axis in axes[:4]:
    axis.set_xticks([])
    axis.set_yticks([])
figure.tight_layout()
figure.savefig("lsci_basics.png", dpi=75)

The full script: lsci_basics.py. It takes about 10 seconds and 1.5 GB of memory.

Limitations¶

It gives the expected contrast, not speckle images: there is no speckle pattern, so $\beta$, which depends on the speckle size against the pixel size, is an input.
Fully developed (Gaussian) speckle is assumed, as in the usual LSCI model.
Motion is rigid and deterministic: each path has one Doppler shift, and the contrast falls only through the spread of shifts between paths. Brownian (unordered) motion applies to the scatterers of media (see blood flow in tissue), not to surfaces.

References¶

R. Bandyopadhyay, A. S. Gittings, S. S. Suh, P. K. Dixon and D. J. Durian. Speckle-visibility spectroscopy: a tool to study time-varying dynamics. Review of Scientific Instruments 76, 093110, 2005.
D. A. Boas and A. K. Dunn. Laser speckle contrast imaging in biomedical optics. Journal of Biomedical Optics 15(1), 011109, 2010.
A. B. Parthasarathy, W. J. Tom, A. Gopal, X. Zhang and A. K. Dunn. Robust flow measurement with multi-exposure speckle imaging. Optics Express 16(3), 1975–1989, 2008.
J. Kim et al. A Monte Carlo rendering framework for simulating optical heterodyne detection. ACM Transactions on Graphics (SIGGRAPH), 2025. Project page