Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction Alfonso López-Ruiza,∗ , Diego Royoa

arXiv:2608.28183v1 [cs.DC] 28 Aug 2026

a Universidad de Zaragoza–I3A, Zaragoza, Spain

ARTICLE INFO

ABSTRACT

Keywords: Non-line-of-sight imaging GPU computing CUDA Real-time-reconstruction Memory-efficient algorithms Multiple-frame aggregation

Non-line-of-sight (NLOS) imaging reconstructs scenes hidden around a corner from indirect light recorded by a single-photon avalanche diode (SPAD). A single reconstruction is a large inverse problem: billions of photon timestamps must be binned, moved through memory, transformed and inverted. As SPAD arrays raise acquisition throughput, reconstruction becomes the limiting stage. We rebuild the GPU execution of two established wave-based algorithms, 𝑓 –𝑘 migration and phasor-fields, for both streaming and offline processing. On the phasor-fields side we assemble the ring-and-radius kernels of previous work once and offline, using the analytic Fourier transform of a ring, so the propagation kernel never exists in dense form at runtime, reducing the memory and bandwidth. We reorganize the pipeline of both algorithms with fused kernels, warp-level photon binning, batched transforms, CUDA graph replay, and FP16 storage applied only where it reduces the actual bottleneck. Our implementations are up to 42× faster than the reference streaming pipeline and up to 14× faster than the fastest published GPU baseline, all while using a fraction of the memory (down to 2.5%), enabling vastly larger and finer reconstructions on the same hardware, or comparable ones within a much lower memory budget. We report an ablation of each implementation choice and propose three denoising strategies enabled by the resulting frame budget for next-generation NLOS video processing.

1. Introduction Non-line-of-sight (NLOS) imaging reconstructs hidden objects by analyzing indirect light scattered by a visible relay wall. In an active time-of-flight system, a pulsed laser illuminates the relay wall and a single-photon avalanche diode (SPAD) records, with picosecond resolution, the roundtrip time of flight of each photon traveling to the hidden scene and back. Repeating this measurement across many wall positions yields a three-dimensional signal indexed by the two spatial coordinates of the wall and the photon time of flight, from which NLOS algorithms recover a volumetric reconstruction. From a computing perspective, a single frame is therefore a large inverse problem whose cost grows with the sampling of the relay wall, the temporal resolution of the sensor, and the size of the reconstructed volume. In practice, it can comprise billions of photon timestamps (Liu et al., 2019), every one of which must be binned, moved through memory, transformed, and inverted. With current technology, the main limitation is the low signal-to-noise ratio (SNR) inherent to three-bounce light transport: only a tiny fraction of the emitted light ever reaches the sensor, and the limited detection rate and noise level of the SPAD force setups to repeat the same measurement many times to accumulate a usable signal. SPAD arrays relieve this acquisition bottleneck by recording photon events in parallel (Schwartz et al., 2008). Still, current technology forces a trade-off between pixel count and timing resolution. Arrays ∗ Corresponding author

alfonso.lopezr@unizar.es (A. López-Ruiz); droyo@unizar.es (D.

Royo) ORCID (s): 0000-0003-1423-9496 (A. López-Ruiz); 0000-0001-6880-322X (D. Royo)

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

with sufficient timing resolution for NLOS reconstruction (∼10ps), only reach 16×16 pixels (Conca et al., 2019). On the other hand, newer 3D-stacked architectures reach far more pixels at coarser timing, such as the 256×128, 60ps sensor of Efe et al. (2026) built for LiDAR. As 3D stacking technology advances, arrays of 128 × 128 pixels and beyond with precise timing would raise pixel-level capture parallelism by nearly two orders of magnitude over current high-timing-resolution arrays. Our work does not address the SNR problem. We observe instead that every credible remedy for low SNR produces more data, and that reconstruction has already become a limiting factor in its own right. The throughput of NLOS imaging algorithms determines the size and resolution of the volumes that can be reconstructed at all, and the SPAD arrays now under development will only push that demand further. Previous work has accelerated NLOS reconstruction with GPU rasterization (Arellano et al., 2017), CUDA (Nam et al., 2021; Sun et al., 2026), field-programmable gate array (FPGA) hardware (Liao et al., 2022), and high-level GPU frameworks such as PyTorch (Mu et al., 2025). We use two main baselines: the streaming phasor-fields pipeline of Nam et al. (2021), which can capture and reconstruct up to five frames per second, and Physics to the Rescue (PttR) (Mu et al., 2025), the fastest published GPU baseline for preprocessed data, built using PyTorch. In this work, we take two established algorithms, 𝑓 –𝑘 migration (Lindell et al., 2019) and phasor-fields (Liu et al., 2019), and rebuild their execution on the GPU. We combine existing algorithmic optimizations with new ones, most notably an offline construction of the ring-and-radius kernels of Jiang et al. (2022) that builds them from the analytic Fourier transform of a ring and so removes dense kernel reconstruction from the runtime entirely, and three strategies Page 1 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction Dense volume (720×520×512) Non-Line-of-Sight setup Ultra-fast laser / SPAD device Relay wall

Occluder

Hidden object

Coarse volume (360×260×512) Hidden object

Sun et al.’26

Our f-k migration

Nam et al. ‘21

PttR (Mu et al. ‘25)

Our phasor fields

Figure 1: Performance overview on the office scene, reconstructed at two volume sizes: a dense 720 × 520 × 512, and a coarser 360 × 260 × 512 which therefore shows fewer details in the reconstruction. For each, the left panel reports reconstruction time per frame and the right panel peak video memory, with the methods grouped into three families: backprojection (Sun et al., 2026), 𝑓 –𝑘 migration, and phasor-fields (Nam et al., 2021; Mu et al., 2025). Hatched bars are previous work and solid bars are ours. Backprojection is three to four orders of magnitude slower than the rest, so its bar runs off the top of the time panel. Against PttR, the fastest published GPU baseline, our 𝑓 –𝑘 is 2.7× faster at the coarser volume and our phasor-fields about 14×, in a fraction of the memory; at the dense volume neither phasor-fields baseline fits in the 16 GB of our GPU at all (OOM).

to reduce background noise in continuous reconstruction. Their workload is highly parallel, yet not automatically fast on a GPU: we analyze how performance and memory scale with capture resolution, reconstruction size, and float precision. We identify the dominant bottlenecks and accelerate the reconstructions with CUDA. As a result, on streaming reconstruction we are 42× faster than that reference pipeline using 8.2% of its video memory (0.27 against 3.34 GB), and on full multi-plane volumes we are roughly an order of magnitude faster than PttR (7.7× and 14.0×) using only 2.5% of its peak memory on average (0.22–0.29 against 8.9– 12.7 GB). The two gains compound: the memory savings lift the ceiling on what can be reconstructed at all, letting us handle reconstruction volumes far larger than prior work, or comparable ones in commodity GPUs. We evaluate two operating modes: streaming reconstruction, where one frame is processed while the next is captured, and offline reconstruction, where all transient data is already available. Streaming tests whether the method can keep up with future capture hardware, while offline reconstruction tests scalability for high-resolution experiments. We will release our implementation upon acceptance to make these comparisons reproducible.

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

In summary, our contributions are: • High-throughput, low-memory GPU pipeline design for 𝑓 –𝑘 migration and phasor-fields. • Algorithmic optimizations, notably an offline construction of the ring-and-radius kernel of Jiang et al. (2022) from the analytic Fourier transform of a ring, which removes the runtime rebuild. • CUDA implementations and benchmarking for both streaming and offline settings. • An analysis of performance and memory usage across reconstruction sizes, precision modes, and padding settings, and against previous accelerated NLOS imaging implementations. • Three strategies for removing noise in streamed reconstructions.

2. Related Work 2.1. Time-of-flight NLOS imaging We focus on active NLOS imaging, where a laser source illuminates the scene. The first NLOS imaging method Page 2 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

relied on expensive streak cameras to capture time-resolved indirect light (Velten et al., 2012, 2013). The introduction of SPADs (Buttafava et al., 2015; Schwartz et al., 2008), far more affordable than streak cameras, broadened access to NLOS imaging and incentivized the development of different reconstruction methods (Faccio et al., 2020). Broadly, there are two possible capture configurations: confocal, where the illuminated and measured points on the relay wall are colocated, and non-confocal, which relaxes this restriction. In our work, we accelerate two wave-based methods: the phasorfields formulation (Liu et al., 2019) and the 𝑓 –𝑘 migration approach (Lindell et al., 2019). Our implementation supports both confocal and non-confocal capture.

Algorithmic optimizations for NLOS imaging. Wavebased NLOS methods let us bring tools from conventional optics into the NLOS domain. For example, some works implement convolutional versions of plane-to-plane light propagation operators (Liu et al., 2020), use approximations targeting efficiency (Ahn et al., 2019) or reducing memory usage (Luesia-Lahoz et al., 2023), or even use the phase information of light to achieve similar quality with fewer operations (Luesia-Lahoz et al., 2025). Most relevant to our work, Jiang et al. (2022) exploit radial symmetry in propagation operators and propose efficient ring-based representations. Their formulation only keeps a single radial profile and rebuilds the full propagation kernel at reconstruction time, which recovers the memory but pays for it in bandwidth. We take this representation to the GPU and, additionally, remove the need to rebuild. The Fourier transform of a ring is analytic, so the radial kernel can be assembled once offline and the propagation kernel then samples it directly (Section 4.3). The kernel never exists in dense form at runtime, so in our implementation the memory saving and the traffic saving apply at the same time. On a different topic, Nam et al. (2021) reduce noise in consecutive frames of streamed NLOS reconstruction using depth-dependent averages. Other works have explored learning-based techniques (Ye et al., 2024), but show limited generalization capabilities. In our work, we propose two novel denoising techniques enabled by our high-throughput capabilities by combining different reconstructions of the same frame.

Hardware acceleration for NLOS imaging. Arellano et al. (2017) used GPU rasterization to accelerate backprojection by tessellating the ellipsoids used for reconstruction. Nam et al. (2021) achieved NLOS imaging at five frames per second, introducing several optimizations that accelerate capture and improve reconstruction quality, but their main speed-up comes from their CUDA implementation of a phasor-fields method. Liao et al. (2022) reported real-time phasor-fields reconstruction on FPGA hardware, but their implementation is not publicly available. Their accelerator implements the ring-and-radius formulation of Jiang et al. (2022) in hardware, including its kernel reconstruction at runtime, so the rebuild cost that our offline construction

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

removes applies to this baseline as well. More recently, Sun et al. (2026) implemented CUDA backprojection for irregular relay surfaces. Still, all these approaches target relatively small reconstruction volumes, and do not scale well. Our closest point of comparison is Physics to the Rescue (Mu et al., 2025) (PttR), which implements phasorfields in PyTorch and reports very high throughput, but for single-plane reconstruction. However, reconstruction of an entire volume requires multiple planes at different depths, especially when the depth of the hidden object is unknown. Their learned feature embedding pipeline avoids this explicit depth sweep and reports real-time performance at 11.8 FPS, but it solves the problem with a learned model. We achieve a 40% speed-up with respect to PttR on single-plane reconstruction. Still, our gains are much more apparent in high-throughput multiple-plane (volumetric) reconstruction, where we reconstruct 7.7× faster on the Z scene and 14.0× faster on the office scene while using on average 2.5% of its peak video memory (Table 3).

Baselines. In summary, we compare against the baseline NLOS imaging works listed in Table 1. We distinguish between methods designed for real-time streaming and offline references. Across all these works the reported numbers differ in reconstruction algorithm, hardware platform, and depth settings, so our evaluation separates experiments we can rerun on our hardware from values reported in the literature.

2.2. GPU pipelines for streaming and memory-bound applications The bottlenecks we observe in previous NLOS reconstruction implementations (many short kernel launches, keeping the GPU fed, volumes larger than available memory, etc.) are not specific to NLOS imaging. They appear in any system where an instrument feeds a GPU faster than software can consume the data (Magro, 2014), and most of our solutions are taken from that literature. Our kernels are short, so launch overhead is a significant fraction of the frame time. We record the per-frame sequence into a CUDA graph and replay it, following Ekelund et al. (2025), who use graphs for iterative solvers, and Li et al. (2026), who add event chaining across streams. To keep the device busy we overlap capture and compute in a producer-consumer front end (Yuan et al., 2026; Xue et al., 2024). Badawood and Aldosari (2026) argue that pipelines of this kind should be evaluated on latency and not only on throughput, so we report frame periods rather than kernel times. Memory capacity is a more serious constraint. Volumes that do not fit in device memory are usually staged in tiles (Jurado et al., 2022; Jaros et al., 2026), or the access pattern is restructured (Liu et al., 2026). We never build the propagation kernel densely, so there is nothing to stage. We treat precision the same way and lower it only where it buys performance, as in mixedprecision solvers that choose the precision of each component (Tsai et al., 2023), batching the small operations that remain (Luszczek et al., 2024). Page 3 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction Table 1 Previous NLOS imaging methods considered in our evaluation. We distinguish between baselines rerun in our hardware and methods included as reported literature values. Partially modified denotes code changes required to run the method on our data, mainly for loading HDF5 capture datasets. Every method is run with its own published reconstruction operator. Compared work

Model

Comparison use

Main limitation for comparison

Fast backprojection (Arellano et al., 2017) CUDA irregular-relay NLOS (Sun et al., 2026) NLOS at 5 FPS (Nam et al., 2021) FPGA-accelerated NLOS (Liao et al., 2022) Physics to the Rescue (Mu et al., 2025)

Backprojection Backprojection Phasor-fields Phasor-fields Phasor-fields / learned

Rerun (partially modified) Rerun Rerun Reported only Rerun

Different reconstruction model; costly for many ellipsoids Different reconstruction model Real-time but low-resolution setting Specialized hardware; code unavailable Highest FPS reported for single-depth reconstruction

𝑓 –𝑘 migration (Lindell et al., 2019) Phasor-fields (Liu et al., 2019)

𝑓 –𝑘 migration Phasor-fields

Offline reference Offline reference

Original offline implementation Original offline implementation

3. Background: non-line-of-sight imaging This section formalizes the active time-of-flight NLOS imaging problem, and introduces the two wave-based NLOS imaging methods used in this work: 𝑓 –𝑘 migration (Lindell et al., 2019) and phasor-fields-based reconstruction (Liu et al., 2019), focusing on the parts that are relevant for our GPU implementation. Figure 2 illustrates a typical NLOS imaging setup. A laser emits ultrashort pulses toward points 𝐱𝑙 ∈  on a visible relay wall. A small fraction of the light scatters from the wall toward the hidden scene and then returns to the relay wall after interacting with the hidden object. An ultrafast sensor records the indirect illumination at points 𝐱𝑠 ∈ , yielding a time-resolved impulse response 𝐻(𝐱𝑙 , 𝐱𝑠 , 𝑡). Here, 𝑡 denotes the time of flight from 𝐱𝑙 to 𝐱𝑠 through the hidden scene. NLOS reconstruction then operates on the impulse response 𝐻(𝐱𝑙 , 𝐱𝑠 , 𝑡) to estimate a hidden-scene response 𝑓 (𝐱𝑣 ) at points 𝐱𝑣 in the reconstruction volume. There are two capture modalities: confocal capture, where the laser and sensor are directed to the same relaywall point (thus 𝐱𝑙 = 𝐱𝑠 ), and non-confocal capture, where illumination and sensing points can differ. Confocal capture is slower, since laser and sensor must move together, but it leads to a more structured reconstruction problem, since every light path starts and ends at the same point. The standard 𝑓 –𝑘 migration formulation assumes confocal data, whereas phasor-fields-based reconstruction is formulated for both modalities. Still, non-confocal measurements can also be approximately converted to a confocal form through temporal shifts of 𝐻(𝐱𝑙 , 𝐱𝑠 , 𝑡) (Lindell et al., 2019). Thus, one can apply both reconstruction methods under any capture modality. Next, we briefly review the two methods.

3.1. NLOS imaging methods 𝑓 –𝑘 migration. Lindell et al. (2019) interpret the hidden scene as a wave field Ψ(𝑥, 𝑦, 𝑧, 𝑡) in which each scene point emits a spherical wave at 𝑡 = 0. The confocal impulse response 𝐻(𝐱𝑠 , 𝐱𝑠 , 𝑡) samples this wave field at relay-wall positions 𝐱𝑠 = (𝑥𝑠 , 𝑦𝑠 , 𝑧𝑠 ), with the wall plane placed at 𝑧𝑠 = 0, and at later times 𝑡 > 0. The reconstruction can then be written as a boundaryvalue problem for this wave field. Stolt interpolation is used in the frequency domain to migrate the measured field from the relay wall, where 𝑧 = 0, back to the time slice 𝑡 = 0, A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

which corresponds to the hidden scene: Stolt

Ψ(𝑥, 𝑦, 𝑧 = 0, 𝑡)

⟹

Ψ(𝑥, 𝑦, 𝑧, 𝑡 = 0);

(1)

|2

(2)

𝑓 (𝐱𝑣 ) = ||Ψ(𝑥𝑣 , 𝑦𝑣 , 𝑧𝑣 , 𝑡 = 0)|

Phasor-fields method. Liu et al. (2019) treat the relay wall as the virtual aperture of a camera with line-of-sight toward the hidden scene. They modulate an illumination signal (𝐱𝑙 , 𝑡) on response, yielding { top of the impulse } ̂ 𝜔 (𝐱𝑙 , 𝐱𝑠 ) =  (𝐱𝑙 , 𝑡) ∗𝑡 𝐻(𝐱𝑙 , 𝐱𝑠 , 𝑡) (𝜔) which represents out-of-focus light waves at the relay wall of angular frequency 𝜔 via a convolution ∗𝑡 and Fourier transform  . Then, NLOS reconstruction is formulated as a planeto-plane propagation from the relay wall to the hidden scene, analogous to how a camera focuses light, to generate an image: ̂ 𝜔 (𝐱𝑣 ) =

∫

𝑒𝑖𝜔𝑡𝑠

∫

𝑒𝑖𝜔𝑡𝑙 ̂ 𝜔 (𝐱𝑙 , 𝐱𝑠 ) d𝐱𝑙 d𝐱𝑠 ,

(3)

where 𝑡𝑙 = |𝐱𝑙 − 𝐱𝑣 |∕𝑐 and 𝑡𝑠 = |𝐱𝑣 − 𝐱𝑠 |∕𝑐 represent times of flight, with 𝑐 the speed of light. The resulting volumetric reconstruction 𝑓 (𝐱𝑣 ) is the sum over all frequencies 𝑓 (𝐱𝑣 ) = +∞ ∫−∞ ̂ 𝜔 (𝐱𝑣 ) d𝜔. When the relay wall and reconstruction planes are parallel, we can apply two optimizations. Liu et al. (2020) formulate Equation 3 using a spatial convolution. When both planes are also aligned, the times of flight 𝑡𝑙 and 𝑡𝑠 have radial symmetry. This structure enables more efficient phasorfields formulations based on ring- or radius-wise kernels (Jiang et al., 2022), reducing computation and memory use while preserving the same propagation model. Since the two-dimensional Fourier transform of a radial function is itself radial, a single one-dimensional profile carries the same information as a full two-dimensional map.

4. Proposed GPU implementation In this section, we describe the design of our NLOS imaging implementations. We first discuss general processing steps (Section 4.1), and later present our optimizations for 𝑓 –𝑘 migration (Section 4.2) and phasor-fields (Section 4.3) methods. Next, we present several techniques to combine consecutive NLOS reconstructions to reduce background noise (Section 4.4). Most of the section is written from the Page 4 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction Laser illuminates the relay wall

SPAD measures indirect photons

Photon detection events Relay wall point

Arrival time

Point 2 Point 8 ...

3.25 ns 4.23 ns ...

Binning

Histograms

Temporal Fourier transform

Frequency-domain phasors

Figure 2: Illustration of the capture setup and the collected data. An ultrafast laser emits a pulse onto the visible relay wall, from which light scatters toward the hidden scene. The scattered light then returns to the relay wall and is captured by the SPAD array as third-bounce photons. The recorded measurements consist of individual photon arrival events, which can be aggregated into transient histograms or directly accumulated into temporal-frequency phasors for wave-based reconstruction methods such as phasor-fields.

streaming point of view, where frames arrive repeatedly and launch overhead matters. Finally, we describe how the offline path differs (Section 4.5). Unless stated otherwise, the results presented throughout this section are generated using the dataset of Nam et al. (2021), with dimensions 190 × 190 × 3002 (𝑥 × 𝑦 × 𝑡). The 3002 temporal bins are transformed into 208 frequency components and propagated over 52 depth planes.

4.1. Producer-consumer data processing SPAD devices output a list of timestamps corresponding to photon arrival times, which are processed to form a transient histogram 𝐻 and later used for reconstruction (see Figure 2). Our work accomplishes this following the producer-consumer (PC) structure from Figure 3. The PC pattern is well suited for simultaneously solving the multiple tasks required by a SPAD asynchronously in a multi-threaded environment. Briefly, PC systems are composed of threads that (i) Produce data and push them into a queue, (ii) pop data for Consuming, or (iii) perform both operations (PC). In a PC system, multiple instances of each worker may exist; however, for simplicity in this explanation, we assume a single instance of each component. We adapted the work of Nam et al. (2021) to implement the PC data processing. Our implementation supports reading data from three different sources. These include photon timestamps read from a SPAD device (live data), or read from a binary file (photon disk data). We can also read existing NLOS datasets (Galindo et al., 2019; Lindell et al., 2019), where photons already come pre-binned into histograms (preprocessed disk data). The Data reader & parser module reads raw photon data and streams the complete record of photons captured by the SPAD array after iterating over every relay wall target. Next, the Photon binner module accumulates the photons into histograms, one per relay wall target, yielding the impulse A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

response 𝐻. For our phasor-fields pipeline, this stage can also perform a time-domain Fourier transform to compute the phasor field required by the method. The Reconstructor works either over the impulse response 𝐻 (𝑓 –𝑘 migration) or its corresponding frequency version (phasor-fields), and writes the result into a texture from the Image pool. For preprocessed disk data, there is no Data reader & parser. The Photon binner instead copies the already available transient data into the layout expected by the selected reconstruction path: index-major temporal bins for 𝑓 –𝑘 migration, and index-major compact phasors for our phasor-fields pipeline.

4.1.1. Photon binner Processing data from the SPAD requires binning raw photon records before reconstruction. We perform this operation on the GPU, whereas Nam et al. (2021) accumulate the photon events directly into a Fourier-domain histogram on the CPU using multi-threading. As discussed in the supplementary material, this stage becomes increasingly heavy when the photon count grows. In our default path with large photon counts, photon batches are copied into pinned host memory allocated with CUDA. This memory is locked in physical RAM and can be transferred efficiently through Direct Memory Access (DMA). Figure 4 compares the performance recorded when using DMA and GPU buffers. At low photon counts, the extra indirection during computations through host-pinned memory can outweigh the transfer benefit. At higher counts, the DMA path is more robust because moving the photon batch becomes a larger part of the frame. We therefore use it as the default for high-throughput streaming. Another important speed-up comes from organizing phasor accumulation at the warp level. Histogramming photons into phasors still requires atomic additions to the output histogram, because many photons may contribute to the same relay-wall point and frequency. Instead of assigning Page 5 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction Section 4.1. Producer-consumer data processing

8

{ float2 sinCos = &cosSinLUT[lutIndex];

9

Offline data

SPAD data

Real-time data

Simulated data

SPAD array

Photon data Point 42 Point 8

...

Data reader & parser

...

H(xl, xs, t)

atomicAdd(&(hist[histIndex].x), sinCos.x); atomicAdd(&(hist[histIndex].y), sinCos.y);

11 12 13

}

Listing 1: Warp-per-photon phasor binning. One warp owns one photon, and its lanes stride over temporal frequencies before atomically updating the compact phasor histogram.

PC

3.25 ns 4.23 ns

10

PC

Section 4.1.1. Photon binner

4.1.2. Display The purpose of the Image pool and Renderer modules in Figure 3 is to enable interoperability between OpenGL (rendering) and CUDA (compute). The Image pool maintains a set of buffers that are written by the reconstructor and read by the Renderer. The pipeline operates as follows: (i) the reconstructor waits for an available image buffer; once written, (ii) the image is pushed into the presentation queue and removed from the writing queue. Asynchronously, (iii) the renderer waits for an image in the presentation queue; once displayed, (iv) the image is returned to the writing queue for reuse. The Renderer writes the contents of an available image into a CUDA surface, which can be both written from CUDA kernels and sampled from OpenGL shaders.

Real-time reconstruction

Section 4.2. f-k migration

Section 4.3. Phasor fields PC

Reconstructor Section 4.4. Combining frames

Section 4.1.2. Display C

Image pool

Renderer

Figure 3: Overview of our producer-consumer pipeline. Labels P and C indicate whether a module is a producer, a consumer, or both.

4.1.3. Optimization of CUDA kernels In streaming mode, the same reconstruction pipeline is executed once per frame. Kernel launch overhead therefore becomes visible, especially after the main kernels have been optimized. CUDA graphs let us record a fixed sequence of operations and replay it with lower scheduling overhead. The captured work includes cuFFT calls, custom kernels, reductions, normalization, and optional band-pass filtering. Presentation through the CUDA/OpenGL display path remains outside this reconstruction graph. Section 5.2 measures how graph replay affects performance.

4.2. 𝑓 –𝑘 migration optimizations The 𝑓 –𝑘 migration method consists of four main stages: (i) input scaling and padding to measure the wave field Figure 4: Histogram construction time for DMA and GPUΨ(𝑥, 𝑦, 𝑧, 𝑡) at 𝑧 = 0, 𝑡 > 0, (ii) a 3D Fast Fourier Transform resident buffers using different photon counts. By default, (FFT), (iii) Stolt interpolation, and (iv) an inverse 3D FFT histograms are built with 1.2M photons, matching the binary to finally obtain Ψ at 𝑧 > 0, 𝑡 = 0 which corresponds datasets of Nam et al. (2021); the 20M and 100M cases are generated by uniformly oversampling this baseline. to the reconstructed volume 𝑓 (𝐱𝑣 ). The Fourier transforms dominate runtime, and the buffers they pass between stages dominate memory. That footprint grows cubically with the side length of the padded reconstruction volume. Existing a thread to every photon-frequency pair, our phasor path implementations therefore become too expensive when they assigns a warp to one photon and lets the lanes cover different allocate separate arrays for padding, remapping coordinates, temporal frequencies, as illustrated in Listing 1. This keeps shifted spectra, and FFT (see the supplementary material). In the writes atomic, but reduces scheduling overhead and our implementation, the input transient volume is stored with produces a more regular access pattern for the compact index-major layout, thus keeping temporal samples adjacent phasor layout.  and giving cuFFT a regular 3D volume to process. 1 2 3

uint thread = blockIdx.x * blockDim.x + threadIdx.x; uint warpId = thread >> 5; // / 32 uint lane = threadIdx.x & 31; // % 32

4 5

// Omitted indexing and data recovery lines

6 7

for (uint freq = lane; freq < numFreqs; freq += 32)

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

4.2.1. Video memory footprint We reduce the overall memory footprint by solving the Stolt interpolation on the fly, avoiding auxiliary remapping buffers. Each thread handles one output voxel and evaluates Page 6 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction Padding (2×)

Without padding

Padding (2×)

Performance (FPS)

Without padding

Video memory (MB)

↑ higher is better

2500

485

500

1500

295

100

Figure 5: Two scenes reconstructed with and without padding using 𝑓 –𝑘. The reconstruction on the left side does not show artifacts without padding, whereas the right one does; however, the hidden scene remains recognizable and the versions without padding were computed at higher speed.

its source position directly from its discrete frequency coordinates. In our implementation, the spatial coordinates map back to shifted Cartesian grid locations, so only the temporal/depth coordinate requires interpolation. The remapping kernel therefore reads two neighboring samples along the contiguous third dimension and performs a linear interpolation. Listing 2 shows the interpolation and write step of the CUDA kernel. We also evaluated a shared-memory variant that pre-caches the block data into shared memory, but we observed marginal gains, as shown later in Section 5.2. 1 2

uint z0u = static_cast<uint>(z0); uint z1u = static_cast<uint>(z1);

3 4 5

uint idx0 = getVolumeIdx(x0, y0, z0u, fftRes); uint idx1 = getVolumeIdx(x0, y0, z1u, fftRes);

6 7 8

cufftComplex c0 = &H[idx0], c1 = &H[idx1]; cufftComplex res = complexLerp(c0, c1, dz);

9 10 11

uint outZ = wrapShiftedIndex(z, shift.z, fftRes.z); uint outIdx = getVolumeIdx(x0, y0, outZ, fftRes);

12 13 14

float scale = fabs(fz) * safeRCP(sqrt_term); result[outIdx] = complexMulScalar(res, scale);

Listing 2: Partial CUDA code of Stolt interpolation. Moreover, our 𝑓 –𝑘 algorithm uses two working FFT buffers, Ψ(0) and Ψ(1) , in addition to the input transient volume Ψ ∈ ℝ𝑁𝑥 ×𝑁𝑦 ×𝑁𝑡 and the reconstructed image 𝑓 (𝐱𝑣 ). Each processing stage alternately reads from one working buffer and writes to the other. In the padded configuration used in our experiments, the working buffers have size 2𝑁𝑥 × 2𝑁𝑦 × 2𝑁𝑧 . This zero-padding reduces FFT wrap-around artifacts, as illustrated in Figure 5, but also increases memory use and runtime. We therefore include an unpadded variant as a faster alternative. This trade-off is summarized quantitatively in Figure 6, where PAD ×1 denotes the unpadded configuration and PAD ×2 denotes ×2 zero-padding. Prior real-time approaches (Nam et al., 2021) have likewise omitted padding, presumably for similar efficiency reasons.

4.2.2. Data quantization Although 𝑓 –𝑘 consumes much less video memory than standard phasor-fields, its complex-valued working data can still be stored in reduced precision. In the FP32 path, Ψ(0) and Ψ(1) are stored as 64-bit complex values. In the A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

1216

1000

200

0

2176

2000

400 300

↓ lower is better

38 Pad ×1 FP32

Pad ×1 FP16

Pad ×2 FP32

500

384

320

Pad ×1 FP32

Pad ×1 FP16

23 Pad ×2 FP16

0

Pad ×2 FP32

Pad ×2 FP16

Figure 6: Performance and video memory footprint of 𝑓 –𝑘 with different padding and precision configurations, measured on a confocal dataset shaped as 256 × 256 × 256.

FP16 path, they are stored as 32-bit complex values and transformed with cuFFT Xt; however, cuFFT Xt operates over

data with power-of-two dimensions, so non-power-of-two padded volumes are rounded up before planning the FFT, which can cancel the memory saved by quantization. FP16 also has a limited numeric range, so the input is scaled before the transform and the spectrum is scaled before the inverse FFT to keep the reconstruction stable. As shown in Figure 6, FP16 reduces video memory consumption, but its performance depends on the padding configuration and the conversion overhead. For instance, the FP16-PAD (×2) setting  exhibited worse performance due to FP32-FP16 conversions, and could have been even worse if the dataset was not shaped with power-of-two dimensions. Therefore, we do not treat reduced precision as a general optimization. In practice, 𝑓 –𝑘 benefits most from reducing temporary allocations and avoiding unnecessary remapping buffers.

4.2.3. Algorithm details A key optimization arises from the observation that the original 𝑓 –𝑘 formulation keeps only one half of the interpolated spectrum (see the supplementary material). As illustrated in Figure 7, this reduces the active Stolt domain to half of the padded frequency volume. The destination buffer is first cleared, and the remapping kernel writes only this active half. Additional acceleration comes from fusing operations that were originally separate. Input scaling, which multiplies each temporal bin by its normalized time index, is fused with zero padding. The cyclic shift is applied directly during the remapping write (see Listing 2). After the inverse FFT, a warp-level reduction computes the maximum response along depth for each spatial location while discarding the padded region. The resulting image is then normalized and optionally band-pass filtered.

4.3. Phasor-fields optimizations The phasor-field framework of Liu et al. (2019) derives several reconstruction methods, each modelled after a lineof-sight imaging system. We focus on the one modelled after a confocal camera. The name refers to the virtual imaging system the operator emulates, not to the capture configuration: as Equation 4 makes explicit through (𝜎𝑠 , 𝜎𝑡 ), the same operator serves confocal and non-confocal captures. As explained before, preprocessing performs a time-domain Fourier transform over photon timestamps to obtain a phasor Page 7 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction Input transient volume Ψ 𝑥×𝑦×𝑡

Distance falloff compensation + zero-padding 2𝑥 × 2𝑦 × 2𝑡

2𝑥 × 2𝑦 × 2𝑡

3D FFT 2𝑥 × 2𝑦 × 2𝑓

2𝑥 × 2𝑦 × 2𝑓

Stolt remapping + cyclic shift 2𝑥 × 2𝑦 × 2𝑓

2𝑥 × 2𝑦 × 𝑓 written

3D inverse FFT 2𝑥 × 2𝑦 × 2𝑧

Figure 8: Mean radial magnitude of the dense phasor-fields Fourier-domain kernel and of our radial approximation, averaged over temporal frequencies and depths. The comparison is one of shape rather than of exact values.

2𝑥 × 2𝑦 × 2𝑧

Max over depth + crop padded region 𝑥×𝑦

Final image 𝑓 (𝐱𝑣 ) Relative runtime cost lower medium

higher

Figure 7: Overview of our optimized 𝑓 –𝑘 migration pipeline with 2× padding. Each stage shows the corresponding input and output dimensions. Stolt remapping only writes the active half of the padded spectrum, while the final stage selects the maximum response along depth for each spatial location. The third dimension is denoted as 𝑡, 𝑓 , or 𝑧 depending on the domain, but all three have the same resolution.

field ̂ 𝜔 (𝐱𝑙 , 𝐱𝑠 ) ∈ ℂ𝑁𝑥 ×𝑁𝑦 ×𝑁𝑓 . The algorithm consists of three main stages: (i) batched 2D spatial FFT, (ii) convolution with plane-to-plane propagation operators, and (iii) batched 2D inverse FFT. In this case, runtime and memory footprint are both dominated by propagation operators. For this reason, real-time phasor-fields implementations precompute a complex propagation kernel 𝐾 ∈ ℂ𝑁𝑥 ×𝑁𝑦 ×𝑁𝑓 ×𝑁𝑑 and keep it for the whole session. Kernel size grows with padding, temporal frequencies and depths from plane-toplane propagations. Additional optimizations exploit the radial symmetry of the phasor-fields kernel in the Fourier domain. Jiang et al. (2022) described this idea through ring and radiusbased kernel representations. Our implementation follows the same intuition, but stores a representation more suitable for GPU runtime. We precompute a unit-ring Fourier basis analytically, then assemble a compact radial kernel for each frequency-depth pair. The supplementary material covers the radial sampling, cutoff, and interpolation used by the implementation. The stored layout is [𝜌, 𝑓 , 𝑑], where 𝜌 is the radial coordinate in the Fourier domain, 𝑓 is temporal frequency, and 𝑑 is propagated depth. During propagation, each Cartesian Fourier pixel accesses this radial kernel through a precomputed pixel-to-radius map. Figure 9 compares our pipeline with the real-time phasor-fields implementations of Nam et al. (2021) and Mu et al. (2025). Dense 2D kernel slices are therefore never constructed at runtime. Constructing those slices would either require

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

allocating a dense buffer again or would reduce concurrency by constructing only one slice at a time. Jiang et al. (2022)’s algorithm rebuilds the dense Cartesian kernel from the compact basis before every use, through an index map that assigns each Cartesian pixel its nearest stored sample. They report this mapping as the heaviest function in their implementation. Instead, we never rebuild it by storing the final radial kernels directly and only perform radius lookup, interpolation, low-pass weighting as well as multiply and accumulate operations. This makes our phasor-fields implementation especially useful in memory-constrained real-time settings. The approximation does not exactly reproduce the dense phasor-fields kernel (see Figure 8). The continuous operator between parallel planes is exactly radial, as Equation 4 later makes explicit, but the dense kernel discretizes it on the finite rectangular grid of the relay wall, whose window and Cartesian sampling break that symmetry; the discretized kernel is therefore only approximately radial, and a radial representation cannot reproduce it pointwise by construction. What the radial form discards is mostly angular structure in the low-energy tail, which is why the reconstructions of Figure 16 remain comparable.

4.3.1. Video memory footprint The dense kernel 𝐾 described at the start of this section dominates the memory budget of the original phasor-fields. Our compact representation consists of: (i) the Fourierdomain phasor data, (ii) the compact radial kernel, (iii) the pixel-to-radius map, and (iv) the working buffers needed for accumulation and inverse FFT. This reduces memory usage substantially and removes dense kernel reconstruction from runtime. 4.3.2. Data quantization The most direct way to further reduce video memory is to quantize the stored propagation data, which here means the compact radial kernel rather than a dense Cartesian buffer. Radial kernels are accumulated in FP32 during precomputation and only downcast when they are stored, which keeps precomputation stable. The FFT planning and the main reconstruction stages are left unchanged.

Page 8 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction a) Nam et al.

b) Physics to the Rescue (PyTorch)

c) Our phasor fields

many explicit GPU stages

many tensor operations

fewer GPU stages, batched, radial kernel

𝑥×𝑦×𝑡

𝑥×𝑦×𝑓

Input phasors ̂ 𝜔

𝑥×𝑦×𝑓

Input transient histograms

𝑥×𝑦×𝑡

𝑥×𝑦×𝑓

𝑥×𝑦×𝑓

Temporal RFFT

2D FFTs

Input phasors ̂ 𝜔

𝑥 × 𝑦 × (𝑡∕2 + 1)

𝑥×𝑦×𝑓

Pad and batched 2D forward FFTs𝑥 × 𝑦 × 𝑓

for each depth 𝑥×𝑦

Frequency filter

Initialize depth buffer

𝑥×𝑦×𝑓

𝑥×𝑦

Frequency weighting

Apply harmonic phase

𝑥×𝑦×𝑓

𝑥×𝑦×𝑑

Sample radial kernel and accumulate frequencies 𝑥×𝑦×𝑑

for each frequency 𝑥×𝑦×𝑓

𝑥×𝑦×𝑓

Multiply data and kernel

𝑥×𝑦×𝑑

Spatial FFT

𝑥×𝑦×𝑓

Add weighted frequency contribution

𝑥×𝑦×𝑓

𝑥×𝑦×𝑑

𝑥×𝑦×𝑑×𝑛

Apply phase term

2D inverse FFT

Multi-frame enhancement

𝑥×𝑦×𝑓 ×𝑑

𝑥×𝑦×𝑑

𝑥×𝑦×𝑓 ×𝑑

Extract magnitude

Extract magnitude

𝑥×𝑦×𝑓 ×𝑑

𝑥×𝑦×𝑓 ×𝑑

𝑥×𝑦×𝑑

𝑥×𝑦×𝑑

𝑥×𝑦×𝑑

Multiply data and kernel

𝑥×𝑦

𝑥×𝑦

Batched 2D inverse FFTs 𝑥×𝑦×𝑓

𝑥×𝑦×𝑓

𝑥×𝑦×𝑑

𝑥×𝑦×𝑑

Sum over frequencies

Maximum Z

𝑥×𝑦×𝑑

𝑥×𝑦×𝑑

𝑥×𝑦

𝑥×𝑦×𝑑

𝑥×𝑦×𝑑×𝑛

Inverse FFT

Multi-frame enhancement

𝑥×𝑦×𝑑

𝑥×𝑦×𝑑

𝑟2 +𝑑 2

,

𝑥×𝑦×𝑑 𝑥×𝑦×𝑓

𝑥×𝑦×𝑓

𝑥×𝑦×𝑑

𝐴(𝑟, 𝑑, 𝑓 ) = 𝑒 𝑖𝜎𝑠 𝜅𝑓

𝑥×𝑦×𝑑

𝑥 × 𝑦 × (𝑡∕2 + 1)

Initialize depth buffer

√ radius 𝑟 = 𝑥2 + 𝑦2 , so it can be written as a weighted sum of unit rings 𝛿𝑟 of radius 𝑟, ∑ 𝑘(⋅, 𝑑, 𝑓 ) = 𝑒 𝑖𝜅𝑓 (𝜎𝑡 𝑑+Δ) 𝐴(𝑟, 𝑑, 𝑓 ) 𝛿𝑟 , 𝑟 (4) √ where 𝜅𝑓 is the wavenumber of temporal frequency 𝑓 and Δ a calibration offset of the optical path. The pair (𝜎𝑠 , 𝜎𝑡 ) selects the capture modality: (2, 0) for confocal data, where illumination and detection share a wall point and the light covers the same distance twice, and (1, 1) for nonconfocal data, where only one leg is propagated and the other is approximated by the plane √ depth. The kernel carries phase only: we drop the 1∕ 𝑟2 + 𝑑 2 amplitude falloff, which varies slowly across the aperture; discarding it mildly reweights contributions across the aperture rather than adding structure, and the per-frame normalization before display absorbs the remaining overall scale. The two-dimensional Fourier transform of a unit ring is analytic,

𝑥×𝑦×𝑑

𝑥×𝑦×𝑑

Shift output volume

Maximum Z 𝑥×𝑦

𝑥×𝑦

𝑥×𝑦×𝑑

𝐵(𝑟, 𝜌) = 2𝜋𝑟 𝐽0 (2𝜋𝑟𝜌),

(5)

𝑥×𝑦×𝑑

Extract magnitude

Normalize and shift 𝑥×𝑦

𝑥×𝑦

𝑥×𝑦×𝑑

𝑥×𝑦×𝑑

Band-pass filter

Maximum Z 𝑥×𝑦

𝑥×𝑦

𝑥×𝑦

Normalize 𝑥×𝑦

Relative runtime cost low medium high

Figure 9: Comparison of three phasor-fields reconstructions. Nam et al. (2021)’s phasor-fields requires many explicit GPU stages, Mu et al. (2025)’s phasor-fields utilizes many individual tensor operations, and our phasor-fields implementation compresses the main work into fewer stages. Darker boxes indicate heavier stages. Dimensions are annotated at the top-left and bottom-right of each stage, with 𝑓 temporal frequencies, 𝑑 propagated depths and 𝑛 the frames used by multi-frame enhancement.

4.3.3. Our improved ring-and-radius construction Jiang et al. (2022) keep a single radial profile and rebuild the dense kernel at reconstruction time. We remove that rebuild: since the Fourier transform of a ring is analytic (Equation 5), the radial kernel can be assembled once offline and sampled directly during propagation. After padding and forward FFTs, the input phasors are represented as ̂ 𝜔 (𝑢, 𝑣, 𝑓 ) ∈ ℂ𝑁𝑥 ×𝑁𝑦 ×𝑁𝑓 , where (𝑢, 𝑣) index the relay-wall spatial frequencies produced by the batched 2D FFT. This is the phasor field of Section 3.1, sampled after the spatial transform rather than at the wall. Let 𝑘(⋅, 𝑑, 𝑓 ) denote the phasor-fields propagation kernel that takes the relay wall to a plane at depth 𝑑 for temporal frequency 𝑓 , written in the spatial domain; the dense kernel 𝐾 introduced above is its two-dimensional Fourier transform, 𝐾 = 2𝐷 {𝑘}. For parallel relay-wall and reconstruction planes it depends only on the wall coordinates through the A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

where 𝜌 is the sampled radial coordinate in the Fourier domain and 𝐽0 is the zeroth-order Bessel function of the first kind. Only the zeroth order appears because a unit ring is constant in the angular coordinate, so every harmonic 𝑛 ≠ 0 2𝜋 carries a factor ∫0 𝑒𝑖𝑛𝜃 d𝜃 = 0 and drops out. Equation 5 is therefore exact rather than a truncated expansion, which leaves no accuracy parameter to tune at this step. Integrating the ring coefficients of Equation 4 against this basis yields the stored radial kernel 𝐾𝑟 (𝜌, 𝑓 , 𝑑) = 𝑒 𝑖𝜅𝑓 (𝜎𝑡 𝑑+Δ) ≈ 𝑒 𝑖𝜅𝑓 (𝜎𝑡 𝑑+Δ)

∫ ∑

𝐴(𝑟, 𝑑, 𝑓 ) 𝐵(𝑟, 𝜌) d𝑟 (6) 𝑞𝑗 𝐴(𝑟𝑗 , 𝑑, 𝑓 ) 𝐵(𝑟𝑗 , 𝜌),

𝑗

where 𝑞𝑗 are trapezoidal weights over the sampled radii 𝑟𝑗 . The modality phase does not depend on 𝑟, so it leaves the quadrature untouched and is applied once per (𝑓 , 𝑑) pair after the sum. We compute this integral offline just once for each (𝑓 , 𝑑) pair. This precomputation introduces two numerical approximations: using a discrete sum (quadrature) and stopping at a maximum radius. Our precomputed kernel differs from the exact discrete kernel in two other ways: we ignore the amplitude falloff (as noted above), and we enforce perfect radial symmetry on the square relay-wall grid. At runtime, each Cartesian Fourier pixel (𝑢, 𝑣) is associated with a continuous radial coordinate 𝜌(𝑢, 𝑣) through a precomputed lookup map. The propagation stage can then be written as 𝑁

̂ 𝑣, 𝑑) = (𝑢,

𝑓 ∑

̃𝑟 (𝜌(𝑢, 𝑣), 𝑓 , 𝑑), (7) ̂ 𝜔 (𝑢, 𝑣, 𝑓 ) 𝑤(𝑓 ) 𝐾

𝑓 =1

̃𝑟 denotes the sampled radial kernel evaluated at the where 𝐾 radius associated with pixel (𝑢, 𝑣). Page 9 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

In practice, 𝜌(𝑢, 𝑣) usually does not fall exactly on an integer radial sample, so the kernel is evaluated by linearly interpolating between neighboring radial bins to avoid the extra artifacts that would appear with a nearest-neighbor lookup. The access pattern is shown in Listing 3: one lane computes the radius-dependent indices and weights, broadcasts them to the lanes covering neighboring depths, and each lane then loads the corresponding FP16 radial-kernel values. 1 2 3 4 5

uint lane = threadIdx.x; uint fullMask = 0xffffffffu; int i0 = 0, i1 = 1; float a = 0.0f, cutoffWeight = 1.0f; int outsideCutoff = 0;

6 7 8 9 10 11

if (lane == 0) { float rf = radiusMap[xy]; cutoffWeight = lowPassWeight(rf, rhoCutoffCoord, ringSize); outsideCutoff = (cutoffWeight <= 0.0f) ? 1 : 0;



Without DDA

Nw = 3

Nw = 5

Nw = 15

294 MB

+30 MB

+45 MB

+120 MB

Figure 10: Denoising results from DDA on a dynamic capture, for temporal windows of 𝑁𝑤 frames. Each row is one frame of the sequence, with the corresponding pose rendered at the left; the leftmost column reconstructs that frame without denoising. All panels of a row show the same frame, and the inset magnifies the same region in every panel.

12

int lastBase = max(0, min(int(rhoCutoffIndex) - 1, int(ringSize) - 2)); i0 = max(0, min(int(floorf(rf)), lastBase)); i1 = min(i0 + 1, int(rhoCutoffIndex)); a = rf - static_cast<float>(i0);

13 14 15 16 17

Single reconstruction

Naïve mean

Coherenceweighted mean

Ground truth

}

18 19 20 21 22 23

i0 = __shfl_sync(fullMask, i0, 0); i1 = __shfl_sync(fullMask, i1, 0); a = __shfl_sync(fullMask, a, 0); cutoffWeight = __shfl_sync(fullMask, cutoffWeight, 0); outsideCutoff = __shfl_sync(fullMask, outsideCutoff, 0);

24 25 26

if (outsideCutoff || !validDepth) return;

27 28 29 30 31 32 33 34

__half2 h0 = rsdKernel[(i0 * numFreqs + wave) * numDepths + depth]; __half2 h1 = rsdKernel[(i1 * numFreqs + wave) * numDepths + depth]; float2 k0 = __half22float2(h0); float2 k1 = __half22float2(h1); float2 k = make_float2( cutoffWeight * fmaf(a, k1.x - k0.x, k0.x), cutoffWeight * fmaf(a, k1.y - k0.y, k0.y));

Listing 3: Symmetry-aware phasor-fields lookup. Lane zero computes the radial interpolation and cutoff terms, broadcasts them across the warp, and each lane loads the two FP16 radialkernel samples for its depth. A frequency-parallel implementation would require atomics on the output buffer. Instead, our propagation kernel assigns one thread to each output element (𝑢, 𝑣, 𝑑) and accumulates all temporal frequencies locally. Threads sharing the same Fourier pixel broadcast the radius lookup, interpolation weights, input phasors, and frequency weights across the warp while neighboring lanes cover adjacent depths. Once all temporal frequencies have been accumulated, the depth slices are transformed back through a batched inverse FFT, performing 𝑁𝑑 two-dimensional inverse FFTs in parallel.

4.4. Post-processing for background noise removal The raw output of NLOS reconstruction often shows small artifacts, especially on the background region. We present and evaluate three noise removal strategies. Two exploit temporal consistency across nearby frames: an extension of the depth-dependent average (DDA) of Nam et al. (2021), and a new averaging strategy driven by our coherence metric. The third, our Mixture of Experts (MoE), A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

Figure 11: Merging reconstructed frames in the complex domain on a simulated scene. Weighting the coherent mean by 𝐶𝑖 (𝐱𝑣 ) attenuates the incoherent background speckle while preserving the hidden surface, whose phase is consistent across the window.

instead combines several reconstructions of the same frame, which is only practical because our pipeline can reconstruct each one in a fraction of the frame budget.

4.4.1. Depth-Dependent Average (DDA) This algorithm follows the idea of Nam et al. (2021), but we adapt it to our data and enable a variable number of frames, instead of only three. The Depth-Dependent Average (DDA) stores the last 𝑁𝑤 reconstructed depth volumes in a circular buffer and computes a weighted average according to the depth, where the window length 𝑁𝑤 is a runtime parameter. This allows temporal denoising to operate before the 3D information is collapsed into a 2D image. The weights vary with depth. After temporal averaging, the rest of the pipeline remains as in previous sections; the output is reduced over depth, normalized for display, and optionally band-pass filtered. Figure 10 illustrates the trade-off introduced by temporal averaging. The effect is mostly temporal rather than spatial. On a moving sequence, a short window leaves the per-frame appearance nearly unchanged while reducing frame-to-frame flicker, without visibly smearing the moving object. The benefit saturates quickly with window length, so we keep short windows by default; longer ones mainly add smoothing and runtime overhead.

Page 10 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

4.4.2. Coherence-weighted merging Both of our reconstruction operators are wave-based, so the volume they return is complex valued and each voxel carries a phase as well as an amplitude. The depthdependent average works on |𝑓𝑖 (𝐱𝑣 )| and discards that phase. It therefore cannot tell whether the successive estimates of a voxel actually agree, or merely have similar magnitudes. Signal from a static surface keeps a consistent phase from frame to frame, whereas reconstruction noise does not, so the phase is exactly the quantity that separates them. We exploit this by measuring, over a window of 𝑁𝑤 frames, how coherently the complex estimates add up, |∑𝑖 | | 𝑗=𝑖−𝑁 +1 𝑓𝑗 (𝐱𝑣 )| 𝑤 | , (8) 𝐶𝑖 (𝐱𝑣 ) = | ∑𝑖 𝜖 + 𝑗=𝑖−𝑁 +1 |𝑓𝑗 (𝐱𝑣 )| 𝑤

with 𝜖 a small constant for numerical stability. The ratio compares the magnitude of the coherent sum with the sum of magnitudes, so 0 ≤ 𝐶𝑖 (𝐱𝑣 ) < 1; the upper bound is approached but never reached, since 𝜖 > 0. It approaches one where the contributions share a phase, and collapses toward zero where they cancel. The merged frame is the coherent mean scaled by this factor, 1 | 𝑓̄𝑖 (𝐱𝑣 ) = 𝐶𝑖 (𝐱𝑣 ) | 𝑁𝑤 |

𝑖 ∑

𝑗=𝑖−𝑁𝑤 +1

| 𝑓𝑗 (𝐱𝑣 )|. |

(9)

Figure 11 shows the effect on a simulated scene: the background speckle, which is incoherent across the window, is strongly attenuated, while the surface of the hidden object survives the weighting. The assumption behind Equation 8 is that the underlying signal is static over the window.

4.4.3. Mixture of experts (MoE) Given the low latency of our reconstruction methods, a direct approach for reducing background noise is to mix complementary reconstructions and keep the structures on which they agree. Intuitively, phasor-fields recovers finer detail at the price of more background noise, while 𝑓 –𝑘 migration trades resolution for a cleaner background. Our goal is to keep high-frequency detail while retaining the noise suppression of the low-pass response. Our Mixture of experts (MoE) shares part of the spectral work between the two branches. First, the input is padded into a shared frequency-minor layout and a batched spatial FFT is applied once over the relay-wall dimensions. The phasor-fields branch reads this shared spectrum, applies its temporal-frequency transform and propagates it with the radial kernel representation, and the resulting depth volume is reduced to an image. The 𝑓 –𝑘 migration branch expands the same shared spectrum into the padded 𝑓 –𝑘 migration volume, performs the temporal FFT, Stolt remapping, inverse FFT and depth maximum. Since both branches only read the shared spectrum and write to disjoint buffers, the captured CUDA graph executes them concurrently on separate streams with dedicated FFT work areas, joining right before mixing. Figure 12 compares each branch with the mixture. The two are not redundant, even though both are filtered A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

Our phasor fields

Our f-k

Our phasor fields (DDA, Nw = 5)

Our MoE

Our MoE (DDA, Nw = 5)

294 MB

361 MB

339 MB

474 MB

519 MB

Figure 12: Merging phasor-fields and 𝑓 –𝑘 migration with MoE. The inset magnifies the same region through a display tone curve (𝛾 = 2.2).

backprojections of the same measurement. What separates them is the filter each applies along the temporal-frequency axis (Marco et al., 2026). 𝑓 –𝑘 migration implements the exact inverse, so it keeps the full band and implicitly favors high frequencies. The phasor-fields branch instead applies a narrow band-pass centered on the virtual wavelength, with almost no response at zero frequency. Neither filter is preferable everywhere in a scene, so we combine them rather than choosing one: ( )𝛽 𝑓 (𝐱𝑣 ) = 𝑓f-k (𝐱𝑣 ) 𝑓phasor-fields (𝐱𝑣 ) + 𝜖 , (10) with 𝜖 a small constant for numerical stability. Formally it is a product of experts, in which 𝑓 –𝑘 migration supplies the amplitude and phasor-fields acts as a confidence gate of weight 𝛽. Setting 𝛽 = 1 recovers the plain product of the two reconstructions. We use 𝛽 = 0.3; larger values keep suppressing background but start to erode dim geometry with it. The mixed result is normalized and optionally band-pass filtered.

4.5. Offline reconstructions Offline reconstruction uses separate implementations from the streaming path. Since each offline reconstruction is typically executed only once for a given input volume, it does not benefit from CUDA graphs. For 𝑓 –𝑘 migration, the offline implementation follows the same main ideas as the streaming version: fused input scaling and padding, on-the-fly Stolt remapping, two FFT work buffers, and optional FP16 execution. For phasor-fields, however, we do not use the streaming phasor-fields implementation. Instead, for fair comparison with existing offline baselines, we accelerate the phasor-fields implementation released by Lindell et al. (2019), which follows work from Liu et al. (2019). This offline path precomputes the plane-toplane propagation operators and the transform matrices for the given volume, uses cuBLAS/cuFFT for the main operations, and can use many more temporal samples and propagated depths than the real-time configuration.

5. Results and evaluation We evaluate the proposed implementations from two viewpoints. First, we measure the complete streaming path Page 11 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

used for real-time NLOS imaging, where photon records are binned while the previous frame is reconstructed. This is the setting in which frame rate and GPU memory pressure matter at the same time. Second, we evaluate offline reconstruction, where all the captured data is already available and the experiments mainly reflect reconstruction time. We then isolate the main implementation choices: padding, reduced precision, CUDA graphs, photon binning, and Stolt interpolation.

Execution hardware. All experiments were run on a desktop system with an Intel(R) Core(TM) i7-14700KF CPU at 3.40 GHz, 64 GB RAM, an NVIDIA RTX 4080 SUPER GPU with 16 GB VRAM, and Windows 11. The implementation is written in C++23 with CUDA 13.0 and OpenMP for CPU multi-threading. Real-time visualization uses OpenGL 4.6 together with CUDA/OpenGL interoperability.

Table 2 Average frame time and peak RAM and VRAM usage on dynamic data, for our implementations and Nam et al. (2021). The capture is a raw SPAD recording binned frame by frame, so it measures whether a method keeps interactive rates under frame-by-frame streaming processing. PttR (Mu et al., 2025) is absent here because it does not handle dynamic captures. Dataset dimensions are shown above the results. Dynamic data 190 × 190 × 3002 ⇒ 208 frequencies, 52 depths, nlosbox1 dataset Approach

↓ Frame time

↑ FPS

↓ Peak VRAM

Our 𝑓 –𝑘 Our 𝑓 –𝑘 (w/o padding)

12.282 ± 0.475 ms 1.087 ± 0.061 ms

81.42 919.85

1, 104.1 MB 302.1 MB

Phasor-fields (Nam et al., 2021) Our phasor-fields

93.769 ± 4.951 ms 2.230 ± 0.126 ms

10.66 448.48

3, 425.3 MB 280.4 MB

↓ Peak RAM

1, 154.5 MB

Adaptation to real SPAD arrays. Our experiments pro-

(a) Photon binning Padded SM 1.01× slower (1.1% fewer FPS)

77 76

900 800

73

700

71

5.1. Real-time performance

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

SM 1.05× slower (4.6% fewer FPS)

74

72

f-k

f-k + SM

f-k

f-k + SM

(b) Stolt shared-memory variant of 𝑓 –𝑘 migration Our f-k 1.11

Frame time (ms)

Our real-time comparison uses the dynamic SPAD data released by Nam et al. (2021). The data consists of raw photon timestamps, so both implementations start from the same input and include photon processing before reconstruction. The reported frame times are not isolated kernel timings, but the times needed by the streaming pipeline to produce one reconstructed frame. Table 2 reports average frame time, frame rate, peak GPU memory, and peak CPU memory. We also include an unpadded 𝑓 –𝑘 migration variant because prior real-time systems commonly avoid FFT padding to reduce cost. The fastest configuration is unpadded 𝑓 –𝑘 migration, which reaches 919.9 FPS on this dataset. This number should be interpreted together with the padding discussion in Section 4.2: removing padding is a speed-memory tradeoff and may introduce wrap-around artifacts. With padding enabled, 𝑓 –𝑘 migration still reaches 81.4 FPS. Our phasorfields variant reaches 448.5 FPS in 0.27 GB of video memory, the least of any configuration here. The reference phasorfields runs at 10.7 FPS and allocates 3.34 GB, an order of magnitude more, mainly because it uses more intermediate buffers. The supplementary material breaks this capture down across every precision and padding configuration, and shows the same general trend.

No padding

1100 1000

75

FPS

cess raw photon records captured by a SPAD array, but they do not stream directly from live acquisition hardware. Instead, we read the recorded photons from disk as fast as possible. This lets us stress the processing pipeline without being limited by the data interface of a particular camera. The supplementary material examines how our frame rate varies with the photon budget, and finds that the pipeline stays above interactive rates at photon counts far larger than those of the dynamic capture. We will therefore release our implementation with the hope that future SPAD devices and data interfaces will accommodate our faster pipeline.

Our phasor fields

Graphs 1.013× faster (1.31% less time)

2.24

1.10

2.23

1.09

2.22

Graphs 1.003× faster (0.28% less time)

2.21

1.08

2.20

1.07

2.19 No graph

Graph

No graph

Graph

(c) CUDA graph replay Figure 13: Ablation of implementation decisions: (a) the photon-binning path, (b) the shared-memory Stolt-remapping variant of 𝑓 –𝑘 migration, and (c) replaying a captured CUDA graph instead of launching each operation separately.

5.2. Implementation choices and ablations Photon binning. Photon binning becomes visible in the frame time once the number of photons grows. Subfigure (a) of Figure 13 compares the relevant variants. The warp-perphoton phasor path improves throughput by distributing the frequency work for one photon across a warp, while still using atomic additions for the final histogram updates. DMA through pinned host memory is not always faster for small batches, but Figure 4 shows why it is useful in the largerphoton-count regime: transfer cost eventually dominates, Page 12 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

and avoiding an additional copy through pageable memory becomes worthwhile. Stolt interpolation. Subfigure (b) shows a case where a more complex memory scheme is not automatically better. Preloading the interpolation neighborhood into shared memory reduces some global reads, but the optimized Stolt kernel already performs only a simple two-sample interpolation along the contiguous third dimension. The extra synchronization therefore cancels most of the expected gain. CUDA graphs. Graphs address a different source of overhead. Once the main kernels are short, the launch sequence itself becomes a measurable part of the frame time. Figure 13c shows the effect of replaying a captured reconstruction graph instead of launching each operation separately: 1.3% for 𝑓 –𝑘 and 0.3% for our phasor-fields. The gain is small because our reconstruction submits only about fifteen operations per frame, so there is little launch overhead once the kernels have been fused; reported speed-ups of 1.4× and above come from iterative solvers that launch hundreds of small kernels per step (Ekelund et al., 2025; Li et al., 2026). The benefit is most relevant for repeated streaming reconstruction, where the same sequence of FFTs and kernels is executed every frame.

5.3. Comparison with other methods We next compare throughput against other fast NLOS reconstruction approaches. Table 3 reports millions of voxels processed per second, average reconstruction time, and peak video memory. The CPU and GPU baselines include Fast Backprojection (Arellano et al., 2017), CUDA-based backprojection (Sun et al., 2026), PttR (Mu et al., 2025), and the phasor-fields pipeline of Nam et al. (2021). PttR is the relevant reference here: it is the fastest published GPU pipeline on preprocessed data, whereas Nam et al. (2021) is our reference for streaming reconstruction and is included for completeness. We also include the FPGA result of Liao et al. (2022); since no source code is available, we report their timing for their stated configuration instead of re-running it. The comparison mainly shows where each family of methods scales. Backprojection methods avoid the dense Fourier-domain kernel, but their throughput remains low for the tested volumes. Our phasor-fields variant gives the best throughput in the reported real-time comparisons because it keeps the propagation kernel small. 𝑓 –𝑘 remains the most stable option across larger temporal dimensions, since its cost is controlled only by FFT volume size. The memory column is what separates the two phasorfields implementations. Dense propagation kernels grow with the number of retained frequencies, so on the office scene PttR needs 12.7 GB of video memory at the 115 frequencies of its block, and Nam et al. (2021) 4.6 GB, against 0.22 GB for our symmetry-aware variant. The gap is a property of the representation rather than of the operating point: a dense kernel is stored per spatial pixel, ours per radius.

Table 3 Performance of fast reconstruction algorithms. Throughput is the processed data volume divided by reconstruction time: for the wave-based methods this is the input transient cube of each block header, of which the phasor-fields family propagates the listed frequency and depth counts, while backprojection uses its reconstructed volume, 256 × 256 × 491 voxels on Z and 360 × 260 × 512 on office. The frequency counts are chosen from the virtual wavelength and capped at the Nyquist limit of the temporal sampling. Green and red cells mark the best and worst value per metric within each dataset. Approach

Mvox/s ↑

Time (s) ↓ VRAM (GB) ↓

RTX 4080 SUPER - Z: 2562 × 512, 255 frequencies, 30 depths Our 𝑓 –𝑘 4,846.80 0.0069 ± 0.0000 9,577.03 0.0035 ± 0.0000 Our phasor-fields Phasor fields (Nam et al., 2021) 521.67 0.0643 ± 0.0042 Phasor fields (Mu et al., 2025) 1,242.38 0.0270 ± 0.0023 CUDA backprojection 1.21 26.593 ± 0.377 0.37 91.238 ± 0.781 Fast backprojection

0.750 0.293 4.389 8.851 0.366 0.689

RTX 4080 SUPER - office: 360 × 260 × 512, 115 frequencies, 51 depths Our 𝑓 –𝑘 3,406.01 0.0141 ± 0.0001 1.496 Our phasor-fields 17,745.27 0.0027 ± 0.0000 0.219 Phasor fields (Nam et al., 2021) 855.57 0.0560 ± 0.0029 4.565 Phasor fields (Mu et al., 2025) 1,269.18 0.0378 ± 0.0023 12.707 CUDA backprojection 1.74 27.586 ± 0.229 0.538 0.10 467.194 ± 3.895 4.164 Fast backprojection Stratix 10 FPGA: 1282 × 𝑡 (N.A.), 69 frequencies, 51 depths FPGA-acc. (Liao et al., 2022) N.A. 0.04

N.A.

not include the CPU memory needed to load the full dataset, which is common to all methods. To keep the tables readable, the input volume is listed separately from the phasor-fields frequency count ΩPF and the phasor-fields depth count 𝑑PF . There is no corresponding column for 𝑓 –𝑘 migration: it migrates the entire wave field rather than selecting a band of frequencies, so its footprint is governed directly by the volume in the first column. The tables use the reducedprecision configuration that is most relevant for memory pressure: our phasor-fields and MoE store their propagation data in FP16, while 𝑓 –𝑘 remains in FP32 in these runs because the FP16 cuFFT path imposes power-of-two constraints. Padding increases the memory footprint of all methods. MoE remains the largest consumer because it runs both branches, allocating the padded 𝑓 –𝑘 migration volume alongside the phasor-fields working buffers; however, since its phasor-fields expert shares the radial kernel representation of Section 4.3, its stored propagation data no longer grows with the spatial dimensions of the relay wall. Our phasorfields implementation needs 4.5 GB in the largest padded configuration and is the only method that completes it: the padded 𝑓 –𝑘 transform exceeds what cuFFT accepts at that size, and the MoE runs out of memory. That headroom is the practical reason for adopting the radial kernel representation in the real-time pipeline despite the slight loss in visual quality it introduces.

5.3.1. Memory scalability Table 4 and Table 5 report peak GPU memory for increasingly large reconstruction volumes. These measurements do A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

Page 13 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

Volume

ΩPF

𝑑PF

𝑓 –𝑘 Phasor fields

MoE

90 × 65 × 128 180 × 130 × 256 360 × 260 × 512 720 × 520 × 1024

29 58 115 229

13 0.031 26 0.188 51 1.496 103 11.960

0.005 0.026 0.208 1.658

0.035 0.234 1.872 14.976

Table 5 Peak GPU memory usage for our reconstruction methods with padding. ΩPF and 𝑑PF are the number of temporal frequencies and propagated depths used by the phasor-fields runs. Missing entries are reported as OOM. Precision: 𝑓 –𝑘 in FP32; phasor fields and MoE in FP16. Peak GPU memory usage (GB) Volume

ΩPF

𝑑PF

90 × 65 × 128 180 × 130 × 256 360 × 260 × 512 720 × 520 × 1024

29 58 115 229

13 0.133 26 1.064 51 8.477 103 OOM

𝑓 –𝑘 Phasor fields

MoE

0.009 0.072 0.567 4.530

0.143 1.147 9.170 OOM

Figure 14: Stage-time decomposition for the main reconstruction pipelines.

5.4. Decomposition of stage time Figure 14 breaks the frame time into the processing stages. This view is useful because the fastest algorithm is not always limited by the same part of the pipeline. In 𝑓 –𝑘 migration, the FFT and inverse FFT stages dominate. For our phasor-fields variant, the compact kernel reduces memory traffic in propagation, so launch overhead and common preprocessing steps become more visible. This explains why there is not a single optimization workflow for all the reconstruction methods.

5.5. Data quantization Reduced precision is useful only when it reduces the actual bottleneck. Figure 15 therefore reports both time and GPU memory for FP16 and FP32 variants. For our phasorfields, FP16 mainly affects the stored propagation kernel, so the expected benefit is lower VRAM usage with limited impact on the FFT stages. For 𝑓 –𝑘, the trade-off is less direct: A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

FP16

FP32

7.5 (b) Our phasor fields 6.0 4.5 3.0 1.5 0.0 56 12 256 8×2 8×5 64× ×12 ×12 64× 128 128

3000 2400 1800 1200 600 0

024 8×1 ×12 128

(x , y , t)

024 6×1 ×25 256

048 6×2 ×25 256

600 450 300 150 0

VRAM (MB)

GPU VRAM

30 (a) Our f-k migration 24 18 12 6 0

VRAM (MB)

Time (ms)

Peak GPU memory usage (GB)

Time (ms)

Time

Table 4 Peak GPU memory usage for our reconstruction methods without padding. ΩPF and 𝑑PF are the number of temporal frequencies and propagated depths used by the phasor-fields runs. Missing entries are reported as OOM. Precision: 𝑓 –𝑘 in FP32; phasor fields and MoE in FP16.

Figure 15: Frame time and peak GPU memory for FP16 and FP32 configurations.

FP16 reduces element size, but cuFFT Xt requires power-oftwo dimensions. As a result, FP16 is not automatically faster

for 𝑓 –𝑘 migration even when it reduces memory. This is why the memory-scalability tables use FP16 storage for our phasor-fields propagation data but keep 𝑓 –𝑘 in FP32 for large data volumes.

5.6. Equivalence with reference implementations We compare our reconstructions using the implementation from Nam et al. (2021) as a reference to show that our results are practically equivalent across quantization levels, propagation kernel formulations and consecutive reconstructed frames. Figure 16 illustrates this on a dynamic dataset. The top panel tracks the mean absolute difference for each frame, and the bottom panel shows the corresponding reference, reconstruction, and error map for one representative frame. We show FLIP error (Andersson et al., 2020), and the absolute difference between reference and candidate images. For the tested methods, the mean error remains low throughout the sequence, and our reconstruction qualitatively looks very similar to the reference frame. Over the whole sequence the mean absolute difference stays at 0.016 for our phasorfields and 0.008 for our 𝑓 –𝑘. The FP16 and FP32 phasorfields curves overlap in the plot. The residual difference with respect to Nam et al. (2021) is therefore attributable to the radial approximation of Section 4.3 rather than to reducedprecision storage, which is consistent with the kernel-domain comparison of Figure 8.

5.7. Background noise suppression The mixture described in Section 4.4.3 is meant to reduce background response without removing scene content, so we measure both. Every row of Table 6 applies the same single normalization and band-pass step. The object level is the mean of the brightest 5% of pixels, and the background level is the mean over an outer border where the hidden object never appears. Their ratio is the contrast, and the mean gradient magnitude over the object serves as a sharpness proxy. Figure 17 draws both regions on a measured frame, so each column of the table can be read against the part of

Page 14 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

Mean abs. difference

Our phasor fields (FP32)

1.00

Our phasor fields (FP16)

zoomed view

0.75

Our f − k (FP16)

frame 106

0.02 0.01

0.50

0.00

0

50

100

150

Table 6 Background suppression on nlosbox1, in FP32. Each method is reported with and without padding. Green and bold mark the best value in each column. The sweep over the pooling weight 𝛽 is reported in the supplementary material.

0.25

Method

Padding

↑ Object

↓ Background

↑ Contrast

↑ Sharpness

0.00

Our 𝑓 –𝑘

no yes

0.615 0.656

0.0379 0.0475

16.2 13.8

0.043 0.037

Our phasor-fields

no yes

0.405 0.410

0.0115 0.0079

35.1 51.6

0.056 0.056

Our MoE (𝛽 = 0.3)

no yes

0.449 0.479

7.4 ⋅ 10−6 0.0009

60,981 512

0.053 0.052

0

25

50

75

100 Frame

125

150

175

(a) Sequence-wide visual error. Reference

Ours

FLIP mean 0.082

Abs. difference mean 0.019

Phasor fields (FP32)

mean 0.082

mean 0.019

mean 0.038

mean 0.009

Phasor fields (FP16)

f − k migration (FP16)

(b) Representative frame. Figure 16: Equivalence with reference implementation on the dynamic capture of Nam et al. (2021). The sequence plot reports the mean absolute difference over time, while the frame panel shows the reference image, our reconstruction, FLIP metric and absolute difference. The reference for phasorfields-based methods is Nam et al. (2021)’s reconstruction; for 𝑓 –𝑘, the reference is 𝑓 –𝑘 using FP32. The two phasor-fields curves overlap.

grow and the mixture’s contrast falls from 60,981 to 512. The ordering between methods holds: phasor-fields is the exception, its background falls and its contrast in fact improves with padding, from 35.1 to 51.6, whereas 𝑓 –𝑘 loses both contrast and sharpness. We also compared the two temporal strategies of Section 4.4.1 and Section 4.4.2 on real dynamic captures, using the same contrast measure. The subject of nlosbox1 moves slowly relative to the frame rate. There, merging in the complex domain over 𝑁𝑤 = 3 frames raises contrast from 36 to 48 and halves the frame-to-frame flicker, well ahead of the depth-dependent average. The coherence weight of Equation 8 is not what produces this gain. Applying it lowers contrast to 42, and it degrades further as the window grows. On a faster capture the ordering reverses altogether, with complex merging dropping contrast from 27 to 22 while the depth-dependent average raises it to 32. Coherence weighting assumes a static signal over the window, and a moving surface violates it, so we use the depth-dependent average by default and use complex merging for slow or static content.

5.8. Offline reconstruction (a) Object and background

(b) Sharpness: |∇| in object

brightest 5% = object

border = background

Figure 17: The regions used in Table 6. (a) The brightest 5% of pixels form the object (red contour), and the outer border forms the background (blue band). Their ratio is the contrast. (b) Sharpness is the mean gradient magnitude inside the same object contour.

the image it summarizes. Values are averaged over the full 100-frame sequence of the nlosbox1 capture. Unpadded, the mixture leaves a background of 7.4 ⋅ 10−6 , three orders of magnitude below either branch alone, while keeping 73% of the 𝑓 –𝑘 object level and more than the phasor-fields branch retains on its own. Padding changes the picture in a way worth stating. It raises every object level slightly, but it also admits real signal into the outer border, so the 𝑓 –𝑘 and mixture backgrounds A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

Table 7 reports one representative configuration per dataset, and Figure 18 summarizes offline throughput in millions of voxels processed per second. Both are extracts: the complete sweep over every resolution is in the supplementary material. We evaluate datasets from Lindell et al. (2019) and Galindo et al. (2019), including confocal and exhaustive captures. Exhaustive captures are first converted to confocal measurements using Normal Moveout Correction. In this setting, the input transient volume is already available, so the reported time focuses on reconstruction rather than on streaming photon records. The reported time is not the reconstruction kernels alone. It covers the whole perrun sequence: allocating and zeroing the output volume, reconstructing, and releasing every intermediate buffer. We measure it this way because those stages are a real cost of running a reconstruction once. Across most experiments, the optimized 𝑓 –𝑘 implementation gives the highest throughput. The optimized offline phasor-fields path is faster by a few milliseconds in only a few tests, and whether it wins depends on the selected dimensions. Both optimized implementations are notably faster than their original MATLAB counterparts, which run entirely on the CPU. Across the full offline benchmark Page 15 of 27

Millions of voxels / second

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction Bike dataset (confocal) 300 250 200

Our f-k f-k

Our phasor fields Phasor fields 230.2

181.3

152.2

128.5

150 100

63.5

50 0

5.1

3.2

128×128×128

4.9

3.2

256×256×256

0.9

35.4

4.3

512×512×512

Millions of voxels / second

Dataset size Office dataset (non-confocal) 250

229.1

227.8

200

153.9

150

153.3

100

47.7

50 0

2.6

2.5

360×260×375

2.4

1.1

360×260×750

0.3

25.9

1.0

360×260×1500

Millions of voxels / second

Dataset size Bunny dataset (non-confocal) 100.3

100

76.5 60.9

56.8

3.0 0

48.3

45.5

50

31×31×512

3.0

2.7

1.9

31×31×1024

2.8

0.9

31×31×2048

Dataset size

Figure 18: Offline throughput, in millions of voxels processed per second, for representative confocal and non-confocal datasets. Hatched bars are our implementations and solid bars the original MATLAB ones. Full timing and memory results are reported in the supplementary material.

table, reconstruction time is reduced by 98.3% for 𝑓 –𝑘 and 96.1% for phasor-fields on average, and peak memory by 52.5%. The curves in Figure 18 are not strictly monotonic: small datasets are affected by allocation and setup overhead, medium-sized datasets make better use of the GPU, and the largest datasets become limited by FFT size and memory traffic. In our measurements, increasing spatial resolution is more expensive than increasing the number of temporal bins, because the dominant FFT and inverse FFT stages grow over the spatial dimensions as well.

6. Conclusions and future work We have presented CUDA implementations of wavebased NLOS reconstruction pipelines for both streaming and offline processing. The proposed system reorganizes their execution around GPU dataflow: fused kernels, compact memory layouts, batched transforms, CUDA graphs, and reduced intermediate allocations. On the phasor-fields side, building the ring-and-radius kernels of Jiang et al. (2022) offline, from the analytic Fourier transform of a ring, keeps the propagation data compact and removes dense kernel reconstruction from the runtime. On the dynamic SPAD dataset of Nam et al. (2021), the padded 𝑓 –𝑘 path reaches 81.4 FPS and our phasor-fields variant reaches 448.5 FPS A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

in 0.27 GB of video memory, compared with 10.7 FPS in 3.34 GB for the reference streaming phasor-fields pipeline of Nam et al. (2021): 42× the frame rate in 8.2% of the memory. Against the fastest published GPU baseline, Physics to the Rescue (Mu et al., 2025), we reconstruct full multi-plane volumes roughly an order of magnitude faster (7.7× on the Z scene and 14.0× on the office scene) while using on average 2.5% of its peak video memory; on the single-plane setting it is built for, the margin is 1.4×. In the offline setting, where the reference implementations are CPU-bound MATLAB code, reconstruction time is reduced by 98.3% for 𝑓 –𝑘 migration and 96.1% for phasor-fields on average, and by 99.4% and 96.1% on the largest office dataset. From the reported results we can conclude that no single optimization dominates all setups. Unpadded 𝑓 –𝑘 migration gives the highest frame rate on small dynamic inputs, but padding remains an important quality and robustness tradeoff. Our phasor-fields variant replaces the dense propagation kernel with a radial representation, reducing memory traffic and video-memory usage while introducing a small visual error. Similarly, FP16 is useful when it targets our stored phasor-fields propagation data, but it is not a universal replacement for FP32; for 𝑓 –𝑘, cuFFT planning constraints and conversion overhead can even outweigh the reduced element size. Our streaming experiments process raw photon records from disk rather than from a live SPAD interface, so future systems should evaluate the same reconstruction path together with camera-side transfer constraints. Photon binning and data movement also remain important costs at high photon counts, motivating more asynchronous uploads and better integration between acquisition and reconstruction. Finally, both 𝑓 –𝑘 migration and phasor-fields still depend heavily on Fourier-domain operators, so future improvements will likely come from better transform scheduling, more memory-efficient propagation representations, and denoising strategies that preserve temporal detail while exploiting the high frame rates enabled by the GPU pipeline.

Declaration of generative AI and AI-assisted technologies in the manuscript preparation process During the preparation of this work the authors used Claude to revise and improve parts of the manuscript text and to generate the scripts that produce the figures and tables from human-measured data. The authors reviewed and edited the output as needed and take full responsibility for the content of the published article.

References Ahn, B., Dave, A., Veeraraghavan, A., Gkioulekas, I., Sankaranarayanan, A., 2019. Convolutional Approximations to the General Non-Lineof-Sight Imaging Operator, in: 2019 IEEE/CVF International Conference on Computer Vision (ICCV), pp. 7888–7898. URL: https:// ieeexplore.ieee.org/abstract/document/9009533, doi:10.1109/ICCV.2019. 00798. iSSN: 2380-7504.

Page 16 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction Table 7 Offline comparison against the original MATLAB implementations, showing the largest configuration measured for each dataset. This is a representative extract: every dataset was also evaluated at lower spatial and temporal resolutions, and those results, together with the per-resolution trends, are reported in full in the supplementary material. Throughput is the reconstructed volume divided by reconstruction time. For exhaustive captures, the dimensions column lists the original measurement grid before normal moveout correction; normal moveout correction produces a 31 × 31 × 𝑇 confocal transient, which is the volume used for the throughput calculation. Max. tracked memory reports the peak footprint of each implementation’s main buffers: host RAM for the original MATLAB code, and the peak of all tracked CUDA allocations for ours, which counts device buffers and pinned host staging together and is therefore not a pure video-memory figure. ↓ Time (s) Dataset

Dimensions

Algorithm

Ours

Original

↑ Throughput (Mvox/s)

↓ Max. tracked memory (GB)

Ours

Original

Ours

Original

Confocal – 𝑓 –𝑘 migration dataset (Lindell et al., 2019) bike

512×512×512

𝑓 –𝑘 migration Phasor fields

2.114 ± 0.03 3.791 ± 0.03

146.624 ± 8.47 30.872 ± 0.67

63.50 35.40

0.92 4.35

17.000 22.001

52.430 33.579

teaser

512×512×512

𝑓 –𝑘 migration Phasor fields

2.016 ± 0.01 3.916 ± 0.08

154.730 ± 14.77 31.458 ± 2.65

66.59 34.28

0.87 4.27

17.000 22.001

52.433 33.579

statue

512×512×512

𝑓 –𝑘 migration Phasor fields

2.018 ± 0.01 3.758 ± 0.02

152.567 ± 8.22 32.195 ± 1.86

66.51 35.72

0.88 4.17

17.000 22.001

52.430 33.579

Confocal – Zaragoza dataset (Galindo et al., 2019) usaf

256×256×1964

𝑓 –𝑘 migration Phasor fields

1.747 ± 0.00 4.172 ± 0.01

139.592 ± 22.22 31.655 ± 1.26

73.70 30.85

0.92 4.07

16.303 21.112

52.429 33.621

bunny

256×256×1964

𝑓 –𝑘 migration Phasor fields

1.734 ± 0.00 4.169 ± 0.02

138.419 ± 8.19 33.418 ± 0.67

74.24 30.88

0.93 3.85

16.303 21.112

52.429 33.621

z

256×256×1964

𝑓 –𝑘 migration Phasor fields

1.780 ± 0.08 4.208 ± 0.02

142.930 ± 6.27 32.743 ± 1.42

72.29 30.58

0.90 3.93

16.303 21.112

52.429 33.636

Tal-generated dataset (Royo and Luesia-Lahoz, 2024), inspired by a Phasor Fields scene (Liu et al., 2020) 2.945 ± 0.05 5.414 ± 0.02

0.28 1.02

17.784 23.023

75.938 66.424

107.02 68.74

3.21 0.90

0.352 0.436

0.770 0.558

0.744 ± 0.05 4.263 ± 0.15

122.69 78.85

2.53 0.44

0.351 0.436

0.773 0.562

0.685 ± 0.10 2.187 ± 0.03

100.31 48.33

2.76 0.86

0.352 0.436

0.769 0.560

360×260×1500

𝑓 –𝑘 migration Phasor fields

concavities

16×16×16×16×2048

𝑓 –𝑘 migration Phasor fields

0.018 ± 0.00 0.027 ± 0.00

0.588 ± 0.04 2.107 ± 0.01

t (in a box)

16×16×16×16×2048

𝑓 –𝑘 migration Phasor fields

0.015 ± 0.00 0.024 ± 0.00

bunny

16×16×16×16×2048

𝑓 –𝑘 migration Phasor fields

0.019 ± 0.00 0.039 ± 0.02

office

507.808 ± 47.68 138.281 ± 13.87

47.67 25.93

Exhaustive – Zaragoza dataset (Galindo et al., 2019)

Andersson, P., Nilsson, J., Akenine-Möller, T., Oskarsson, M., Åström, K., Fairchild, M.D., 2020. FLIP: A Difference Evaluator for Alternating Images. Proceedings of the ACM on Computer Graphics and Interactive Techniques 3, 15:1–15:23. doi:10.1145/3406183. Arellano, V., Gutierrez, D., Jarabo, A., 2017. Fast back-projection for non-line of sight reconstruction. Optics Express 25, 11574–11583. URL: https://opg.optica.org/oe/abstract.cfm?uri=oe-25-10-11574, doi:10.1364/OE.25.011574. Badawood, D.Y., Aldosari, F.M., 2026. An Enhanced LatencyBounded GPU-Resident Pipeline for Real-Time Market Stream Visualization. Computation 14, 140. URL: https://doi.org/10.3390/ computation14060140, doi:10.3390/computation14060140. Buttafava, M., Zeman, J., Tosi, A., Eliceiri, K., Velten, A., 2015. Nonline-of-sight imaging using a time-gated single photon avalanche diode. Optics Express 23, 20997–21011. Conca, E., Riccardo, S., Sesta, V., Portaluppi, D., Zappa, F., Tosi, A., 2019. Design of a 16 x 16 fast-gated SPAD imager with 16 integrated shared picosecond TDCs for non-line-of-sight imaging, in: Emerging Imaging and Sensing Technologies for Security and Defence IV, SPIE. pp. 25–32. doi:10.1117/12.2536352. Efe, B.C., Padmanabhan, P., Zhang, C., Ximenes, A.R., Lee, M.J., Charbon, E., 2026. Characterization of a 3-D stacked CMOS SPAD depth sensor in 45/22-nm technology. Optics Express 34, 11630–11651. URL: https://opg.optica.org/oe/fulltext.cfm?uri= oe-34-7-11630, doi:10.1364/OE.583545. Ekelund, J., Markidis, S., Peng, I., 2025. Boosting Performance of Iterative Applications on GPUs: Kernel Batching with CUDA

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

Graphs, in: 2025 33rd Euromicro International Conference on Parallel, Distributed, and Network-Based Processing (PDP), pp. 70– 77. URL: http://arxiv.org/abs/2501.09398, doi:10.1109/PDP66500.2025. 00019. arXiv:2501.09398 [cs]. Faccio, D., Velten, A., Wetzstein, G., 2020. Non-line-of-sight imaging. Nature Reviews Physics 2, 318–327. Galindo, M., Marco, J., O’Toole, M., Wetzstein, G., Gutierrez, D., Jarabo, A., 2019. A dataset for benchmarking time-resolved non-line-of-sight imaging, in: ACM SIGGRAPH 2019 Posters, ACM. pp. 1–2. URL: https://graphics.unizar.es/nlos, doi:10.1145/3306214.3338583. Jaros, M., Riha, L., Strakos, P., Kozubek, T., 2026. Out-of-core aware multiGPU rendering for large-scale scene visualization. Future Generation Computer Systems 175, 108109. URL: https://www.sciencedirect. com/science/article/pii/S0167739X25004030, doi:10.1016/j.future.2025. 108109. Jiang, D., Liu, X., Luo, J., Liao, Z., Velten, A., Lou, X., 2022. Ring and Radius Sampling Based Phasor Field Diffraction Algorithm for NonLine-of-Sight Reconstruction. IEEE Transactions on Pattern Analysis and Machine Intelligence 44, 7841–7853. URL: https://ieeexplore. ieee.org/document/9560039, doi:10.1109/TPAMI.2021.3117962. Jurado, J.M., Padrón, E.J., Jiménez, J.R., Ortega, L., 2022. An out-ofcore method for GPU image mapping on large 3D scenarios of the real world. Future Generation Computer Systems 134, 66–77. URL: https:// www.sciencedirect.com/science/article/pii/S0167739X22000978, doi:10. 1016/j.future.2022.03.022. Li, Z., Huang, T.W., Ogras, U., 2026. SET: Stream-Event-Triggered Scheduling for Efficient CUDA Graph Pipelines. URL: http://arxiv.

Page 17 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction org/abs/2606.05495, arXiv:2606.05495. Liao, Z., Jiang, D., Liu, X., Velten, A., Ha, Y., Lou, X., 2022. FPGA Accelerator for Real-Time Non-Line-of-Sight Imaging. IEEE Transactions on Circuits and Systems I: Regular Papers 69, 721–734. URL: https:// ieeexplore.ieee.org/document/9599476, doi:10.1109/TCSI.2021.3122309. Lindell, D.B., Wetzstein, G., O’Toole, M., 2019. Wave-based non-lineof-sight imaging using fast f-k migration. ACM Trans. Graph. 38, 116:1–116:13. URL: https://dl.acm.org/doi/10.1145/3306346.3322937, doi:10.1145/3306346.3322937. Liu, S., Dong, X., Cao, J., Chu, R., Wang, Z., Wang, Q., Bai, X., 2026. Mapo: Performance model driven GPU memory access code optimization. Future Generation Computer Systems 184, 108606. URL: https:// www.sciencedirect.com/science/article/pii/S0167739X26002402, doi:10. 1016/j.future.2026.108606. Liu, X., Bauer, S., Velten, A., 2020. Phasor field diffraction based reconstruction for fast non-line-of-sight imaging systems. Nature Communications 11, 1645. URL: https://www.nature.com/articles/ s41467-020-15157-4, doi:10.1038/s41467-020-15157-4. Liu, X., Guillén, I., La Manna, M., Nam, J.H., Reza, S.A., Huu Le, T., Jarabo, A., Gutierrez, D., Velten, A., 2019. Non-line-of-sight imaging using phasor-field virtual wave optics. Nature 572, 620– 623. URL: https://www.nature.com/articles/s41586-019-1461-3, doi:10. 1038/s41586-019-1461-3. Luesia-Lahoz, P., Gutierrez, D., Muñoz, A., 2023. Zone Plate Virtual Lenses for Memory-Constrained NLOS Imaging, in: ICASSP 2023 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 1–5. URL: https://ieeexplore.ieee.org/ document/10094929, doi:10.1109/ICASSP49357.2023.10094929. Luesia-Lahoz, P., Sultan, T., Peterson, F.B., Velten, A., Gutierrez, D., Muñoz, A., 2025. Zero-Phase Phasor Fields for Non-Line-of-Sight Imaging, in: 2025 IEEE International Conference on Computational Photography (ICCP), IEEE Computer Society. pp. 1–9. URL: https://www.computer. org/csdl/proceedings-article/iccp/2025/11143850/29JEeeYXDIk, doi:10. 1109/ICCP64821.2025.11143850. Luszczek, P., Abdelfattah, A., Anzt, H., Suzuki, A., Tomov, S., 2024. Batched sparse and mixed-precision linear algebra interface for efficient use of GPU hardware accelerators in scientific applications. Future Generation Computer Systems 160, 359–374. URL: https:// www.sciencedirect.com/science/article/pii/S0167739X24003017, doi:10. 1016/j.future.2024.06.004. Magro, A., 2014. A Real-Time, GPU-Based, Non-Imaging Back-End for Radio Telescopes. Ph.D. thesis. University of Malta. URL: http://arxiv.org/abs/1401.8258, arXiv:1401.8258. Marco, J., Jarabo, A., Nam, J.H., Tosi, A., Gutierrez, D., Velten, A., 2026. A comprehensive study of time-of-flight non-line-of-sight imaging. URL: http://arxiv.org/abs/2603.09548, doi:10.48550/arXiv.2603.09548. arXiv:2603.09548. Mu, F., Mo, S., Peng, J., Liu, X., Nam, J.H., Raghavan, S., Velten, A., Li, Y., 2025. Physics to the Rescue: Deep Non-Line-of-Sight Reconstruction for High-Speed Imaging. IEEE Transactions on Pattern Analysis & Machine Intelligence 47, 6146–6158. URL: https://doi.ieeecomputersociety. org/10.1109/TPAMI.2022.3203383, doi:10.1109/TPAMI.2022.3203383. place: Los Alamitos, CA, USA. Nam, J.H., Brandt, E., Bauer, S., Liu, X., Renna, M., Tosi, A., Sifakis, E., Velten, A., 2021. Low-latency time-of-flight non-lineof-sight imaging at 5 frames per second. Nature Communications 12, 6526. URL: https://www.nature.com/articles/s41467-021-26721-x, doi:10.1038/s41467-021-26721-x. Royo, D., Luesia-Lahoz, P., 2024. diegoroyo/tal: y-tal v0.16.1. URL: https://zenodo.org/records/11197745, doi:10.5281/zenodo.11197745. Schwartz, D.E., Charbon, E., Shepard, K.L., 2008. A Single-Photon Avalanche Diode Array for Fluorescence Lifetime Imaging Microscopy. IEEE Journal of Solid-State Circuits 43, 2546–2557. URL: https:// ieeexplore.ieee.org/document/4685420, doi:10.1109/JSSC.2008.2005818. Sun, Y., Hong, Y., Qiu, Z., Li, W., Li, W., Sun, Q., Xu, F., 2026. CUDAaccelerated Non-line-of-sight imaging with irregular relay surfaces. Optics and Lasers in Engineering 200, 109591. URL: https:// www.sciencedirect.com/science/article/pii/S0143816625007754, doi:10.

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

1016/j.optlaseng.2025.109591. Tsai, Y.H.M., Beams, N., Anzt, H., 2023. Three-precision algebraic multigrid on GPUs. Future Generation Computer Systems 149, 280–293. URL: https://www.sciencedirect.com/science/article/pii/ S0167739X23002741, doi:10.1016/j.future.2023.07.024. Velten, A., Willwacher, T., Gupta, O., Veeraraghavan, A., Bawendi, M.G., Raskar, R., 2012. Recovering three-dimensional shape around a corner using ultrafast time-of-flight imaging. Nature Communications 3, 745. Velten, A., Wu, D., Jarabo, A., Masia, B., Barsi, C., Joshi, C., Lawson, E., Bawendi, M., Gutierrez, D., Raskar, R., 2013. Femto-Photography: Capturing and Visualizing the Propagation of Light. ACM Transactions on Graphics 32, 1–8. West, R., Golchin, A., Njavro, A., 2023. Real-Time USB Networking and Device I/O. ACM Trans. Embed. Comput. Syst. 22, 67:1–67:38. URL: https://dl.acm.org/doi/10.1145/3604429, doi:10.1145/3604429. Xue, W., Wang, H., Roy, C.J., 2024. CPU–GPU heterogeneous code acceleration of a finite volume Computational Fluid Dynamics solver. Future Generation Computer Systems 158, 367–377. URL: https:// www.sciencedirect.com/science/article/pii/S0167739X24001900, doi:10. 1016/j.future.2024.04.049. Ye, J., Hong, Y., Su, X., Yuan, X., Xu, F., 2024. Plug-and-Play Algorithms for Dynamic Non-line-of-sight Imaging. ACM Trans. Graph. 43, 155:1– 155:12. URL: https://dl.acm.org/doi/10.1145/3665139, doi:10.1145/ 3665139. Yuan, L., Yan, D., Adhikari, S., Ahmad, A., Wang, F., 2026. 3DPipe: A Pipelined GPU Framework for Scalable Generalized Spatial Join over Polyhedral Objects. URL: http://arxiv.org/abs/2604.19982, arXiv:2604.19982.

Page 18 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

A. Integration in a SPAD array In this section, we evaluate the potential results that could be achieved by running our code on current and next-generation live SPAD arrays. The experiments still use recorded photons rather than a physical camera stream, but they answer two practical questions: the number of photons our work can process every frame, and the image quality to expect from these photons. These experiments use the dynamic dataset released by Nam et al. (2021) for streamed reconstruction. Each frame contains approximately 1.2 ⋅ 106 photon records. Reading those records from disk lets us stress the processing pipeline without being limited by the bandwidth or buffering behavior of a specific acquisition device. With this, we take the opposite perspective: instead of fixing the photon count and reporting the resulting frame rate, we fix the frame rate and determine how many photons our system can process per second. Figure 19 shows that our pipeline remains above 60 FPS for photon counts up to approximately 47M per frame, and above 24 FPS up to approximately 96M photons per frame. These counts are much larger than those of the dynamic dataset used in the main streamed comparison. While our work handles several million photons until surpassing the 60 FPS threshold, we also aim to demonstrate that, for a lower number of photons, the hidden scene remains recognizable using both 𝑓 –𝑘 and phasor-fields algorithms. Figure 20 shows reconstructions of a scene with progressively fewer photons. Both our 𝑓 –𝑘 and phasor-fields preserve the main hidden shape at reduced counts, although phasor-fields is more robust to noise in this example. The third column uses fewer photon records than a typical frame of the dynamic dataset from Nam et al. (2021), and the fourth column uses less than half of that amount. For context, the dynamic dataset of Nam et al. (2021) is processed at approximately 400 FPS in this benchmark. The teaser dataset of Lindell et al. (2019), which contains 178M photons captured over 180 min, is processed at slightly below 16 FPS. The reconstructions with fewer photons in Figure 20 suggest that much shorter effective exposures or

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

9,931,460 photons

1,000,000 photons

500,000 photons

33.02

25.04

22.06

29.69

26.78

24.60

Ground truth

Figure 19: Photon binning time for a dynamic NLOS capture of shape 190 × 190 × 3002 (𝑥 × 𝑦 × 𝑡), whose photons are accumulated directly into 208 temporal frequencies. The photon count varies slightly from frame to frame and typically stays near 1.2 ⋅ 106 .

Our phasor fields

Our f-k migration

178,630,366 photons

Figure 20: Teaser dataset from (Lindell et al., 2019). The first and second columns correspond to measurements of 180 and 10 minutes, respectively. We downsample the former to 106 and 5 ⋅ 105 photons to show that the hidden scene remains recognizable with fewer captured photons. The first and second rows show our 𝑓 –𝑘 and phasor-fields reconstructions, respectively. The numbers within each reconstruction indicate the Peak Signal-to-Noise Ratio with respect to the most informed reconstruction, shown in the first column.

higher acquisition rates can still be useful for interactive operation. Finally, we estimate the acquisition data rate. Nam et al. (2021) report photon events streamed from the hardware queue over USB 3.0, whose speed rate is 5 Gbit s−1 , i.e. at most about 500 MB s−1 of payload (West et al., 2023). Since each photon record occupies 4 bytes (Nam et al., 2021), the link is bounded at about 125M photon events per second. With the photon counts of the dynamic sequence, that bandwidth gives an ideal upper bound of about 104 frames per second; protocol and device overheads probably place the practical rate below that. If fewer photons are sufficient for the target scene, as in Figure 20, the acquisition-side frame rate can be higher.

B. Further implementation details B.1. f–k migration implementation The following algorithms show the pseudocode of the original 𝑓 –𝑘 migration implementation (Algorithm 1) and our version (Algorithm 2). We mainly use them to highlight the dataflow difference. The reference pipeline allocates separate arrays for scaling, padding, Fourier shifts, and Stolt remapping. Our CUDA path keeps the input volume and uses two complex working buffers of size (2𝑁𝑥 , 2𝑁𝑦 , 2𝑁𝑧 ) for the FFT-domain stages. The initial time-bin scaling is fused with padding, cyclic shifts are handled inside the kernels, and Stolt remapping is evaluated on the fly instead of being stored as an auxiliary coordinate volume. Finally, after the inverse FFT we take the maximum magnitude over depth, rather than the squared amplitude of Equation 2 of the main

Page 19 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

paper, which in our experiments improves the tone mapping of the final NLOS image. Algorithm 1: 𝑓 –𝑘 migration (Lindell et al., 2019). Reconstruction pipeline following prior work. Pseudocode shown in a Python-like style. 𝑁𝑥 , 𝑁𝑦 , and 𝑁𝑧 denote the volume dimensions. 𝑑max refers to the maximum distance, whereas aptwidth is the physical size in meters of the scanned relay wall. 1 def 𝑓 –𝑘_original (Ψ, 𝑁𝑥 , 𝑁𝑦 , 𝑁𝑧 , 𝑑max , aptwidth ): 2 3 4 5 6 7 8 9 10 11 12 13 14

(𝑋𝑤 , 𝑌𝑤 , 𝑍𝑤 ) ← mgrid[−𝑁𝑥 ∶𝑁𝑥 , −𝑁𝑦 ∶𝑁𝑦 , −𝑁𝑧 ∶𝑁𝑧 ] (𝑋𝑤 , 𝑌𝑤 , 𝑍𝑤 ) ← (𝑋𝑤 ∕𝑁𝑥 , 𝑌𝑤 ∕𝑁𝑦 , 𝑍𝑤 ∕𝑁𝑧 )  ← tile(linspace(0, 1, 𝑁𝑧 ), (𝑁𝑥 , 𝑁𝑦 , 𝑁𝑧 )) Ψ←Ψ⊙ ⊳ Initial scaling and zero-padding Ψ′ ← zeros(2𝑁𝑥 , 2𝑁𝑦 , 2𝑁𝑧 ) Ψ′ [∶𝑁𝑥 , ∶𝑁𝑦 , ∶𝑁𝑧 ] ← Ψ ̂ ←  {f f tshif t(Ψ′ )} Ψ /( ) 4𝑁𝑧 ⋅ (aptwidth ∕2) 𝑠 ← 𝑁𝑥 ⋅ 𝑑max √ ′ ← 2 + 𝑌 2) + 𝑍2 𝑠2 (𝑋𝑤 𝑍𝑤 𝑤 𝑤

⊳ Forward FFT ⊳ Stolt remapping

̂ (𝑋𝑤 , 𝑌𝑤 , 𝑍 ′ )) Ψ′ ← interpn(𝑋𝑤 , 𝑌𝑤 , 𝑍𝑤 , Ψ, 𝑤 ′ ̂ ← Ψ ⊙ 𝟏𝑍 >0 Ψ ⊳ Spectral filtering/compensation 𝑤 ̂ ←Ψ ̂ ⊙ |𝑍𝑤 |′ Ψ max(𝑍𝑤 ) ̂ Ψ ← if f tshif t( −1 {Ψ}) ⊳ Inverse FFT

15

𝑓 (𝐱𝑣 ) ← max_z(|Ψ2 |)

16

return 𝑓 (𝐱𝑣 )

⊳ Final reconstruction

Algorithm 2: Our CUDA-based 𝑓 –𝑘 migration pipeline. 1 def our_𝑓 –𝑘 (Ψ, 𝑁𝑥 , 𝑁𝑦 , 𝑁𝑧 , 𝑑max , aptwidth ):

̂ ⋅, ⋅) ← 0 Ψ(⋅, ̂ ← scale_pad_f f tshif t(Ψ) Ψ ⊳ Time-bin scaling fused with padding ̂ ←  {Ψ} ̂ 4 Ψ ( ) 5 𝑠 ← 𝑁𝑥 ⋅ 𝑑max ∕ 4𝑁𝑧 ⋅ aptwidth ∕2 ̂ 𝑠) 6 Ψ′ ← stolt(Ψ, ⊳ On-the-fly remapping −1 ̂ 7 Ψ ←  {Ψ′ } ̂ 8 𝑓 (𝐱𝑣 ) ← if f tshif t_max_magnitude(Ψ) 9 return 𝑓 (𝐱𝑣 ) 2 3

Algorithm 3: Phasor-fields as implemented by Nam et al. (2021). Reconstruction pipeline following prior work. 𝑃 holds the sampled phasor field ̂ 𝜔 and 𝐾 the dense propagation kernel, both as defined in Section 4.3 of the main paper; 𝑃̂ denotes the stack of transformed slices 𝑃̂𝑐 produced by the first loop. 𝑁𝑥 , 𝑁𝑦 , and 𝑁𝑓 denote the relay-wall spatial dimensions and the number of temporal-frequency slices; 𝑑min , 𝑑max , and Δ𝑑 define the propagated-depth range, and 𝑤 weights each frequency slice. 1 def phasor_fields_nam (𝑃 , 𝐾, 𝑁𝑥 , 𝑁𝑦 , 𝑁𝑓 , 𝑑min , 𝑑max , Δ𝑑, 𝑤): 2 3

for 𝑐 ← 1 to 𝑁𝑓 : 𝑃̂𝑐 ← 2𝐷 {𝑃 (⋅, ⋅, 𝑐)}

𝑖←0 𝑃 (⋅, ⋅, ⋅) ← 0 6 for 𝑑 ← 𝑑min to 𝑑max , step Δ𝑑 : 7  ← 𝑃̂ ⊙ 𝐾(⋅, ⋅, ⋅, 𝑖) ⊳ Convolution over the whole frequency stack 8 for 𝑐 ← 1 to 𝑁𝑓 : 9 𝑃 (⋅, ⋅, 𝑖) ← 𝑃 (⋅, ⋅, 𝑖) + 𝑤𝑐 (⋅, ⋅, 𝑐) −1 {𝑃 (⋅, ⋅, 𝑖)} 10 𝑃 (⋅, ⋅, 𝑖) ← 2𝐷 ⊳ Inverse FFT 11 𝑃̂(⋅, ⋅, 𝑖) ← |𝑃 (⋅, ⋅, 𝑖)| ⊳ Magnitude at depth 𝑖 12 𝑖←𝑖+1 13 𝑓 (𝐱𝑣 ) ← max_z(𝑃̂) 14 return 𝑓 (𝐱𝑣 ) 4 5

B.2.1. Our improved ring-and-radius construction Our symmetry-aware phasor-fields variant replaces the dense Cartesian propagation kernel used in streamed reconstructions with a compact radial representation. We store the radial kernel as a 3D array indexed by radial Fourier bin 𝜌, temporal frequency 𝑓 , and depth 𝑑. A precomputed lookup table maps each Cartesian Fourier pixel to its continuous radius. At runtime, the GPU linearly interpolates this radial representation and accumulates all temporal frequencies locally for each (𝑢, 𝑣, 𝑑) output voxel. We evaluate the analytic Bessel function representation only during precomputation. For a circular ring of spatial radius 𝑟, the angular component of its 2D Fourier transform depends solely on the Fourier radius 𝜌, reducing to the zerothorder Hankel term: 𝐵(𝑟, 𝜌) = 2𝜋𝑟 𝐽0 (2𝜋𝑟𝜌).

B.2. Phasor-fields implementation Similarly to the previous section, Algorithm 3 presents the reference phasor-fields implementation used by Nam et al. (2021). In the reference implementation, the initial Fourier transforms are not batched: although they could be executed as a group, they are instead handled sequentially. The convolved volume  is also processed sequentially, weighting each slice (⋅, ⋅, 𝑐) by 𝑤𝑐 . Moreover, the Fouriertransformed matrices are not padded, which reduces memory consumption at the expense of introducing artifacts.

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

⊳ 2D forward FFTs

(11)

This expression is exact, as noted in Section 4.3 of the main paper. Because a unit ring is uniform across the angular coordinate, all higher-order harmonics (𝑛 ≠ 0) integrate to zero. This formulation is exact and involves no series truncation. The only approximations are in the finite trapezoidal quadrature integration over sampled radii, the high-frequency radial cutoff described below, and enforcing radial symmetry on the discrete, rectangular relay-wall grid. During precomputation, we evaluate this basis once over a grid of spatial and Fourier radii. The pipeline then multiplies these real-valued basis samples by the complex propagation phase for each (𝑓 , 𝑑) pair, integrates across the sampled spatial rings, and stores the resulting radial kernel 𝐾𝑟 (𝜌, 𝑓 , 𝑑). Page 20 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

The sampling is chosen in two stages, both derived from the capture geometry. Let Δ𝑥 and Δ𝑦 be the relay-wall sample spacings along 𝑥 and 𝑦, and let 𝑁𝑢 and 𝑁𝑣 be the padded FFT grid extents along 𝑦 and 𝑥 respectively, so that 𝑁𝑣 Δ𝑥 and 𝑁𝑢 Δ𝑦 are the physical side lengths of the padded aperture. The largest spatial radius covered by that grid and the radial Nyquist limit of the Cartesian spectrum are then: √ )2 ( )2 ( 1 1 𝑟max = 𝑁 Δ𝑥 + 𝑁 Δ𝑦 , 𝑣 𝑢 2 2 (12) √ )2 ( )2 ( 1 1 + 2Δ𝑦 . 𝜌max = 2Δ𝑥 First, spatial radii are sampled uniformly over [0, 𝑟max ]. To resolve the phase oscillations of the phasor-fields kernel 𝐴(𝑟, 𝑑, 𝑓 ) (Equation 4 of the main paper), we differentiate its phase √ to obtain the instantaneous radial frequency 𝜎𝑠 𝑟∕(𝜆 𝑟2 + 𝑑 2 ). This quantity is upper-bounded by 𝜎𝑠 ∕𝜆min cycles per meter, where 𝜆min is the shortest reconstructed wavelength and 𝜎𝑠 is the modality factor (2 for confocal, 1 for non-confocal). Rounding this Nyquist sampling limit up to the nearest power of two determines the number of rings and their uniform spacing: 𝑁𝑟 = 2⌈log2 (⌈2𝜎𝑠 𝑟max ∕𝜆min ⌉+1)⌉ ,

Δ𝑟 =

𝑟max , (13) 𝑁𝑟 − 1

floored at 𝑁𝑟 = 2. The radial integral of Equation 6 of the main paper uses trapezoidal weights over these samples, which for uniform spacing are: 𝑞0 = 𝑞𝑁𝑟 −1 = 12 Δ𝑟,

𝑞𝑗 = Δ𝑟

(0 < 𝑗 < 𝑁𝑟 −1). (14)

Second, the radial spectral domain spans [0, 𝜌max ], using the same Nyquist criterion to similarly determine the sample count: 𝑁𝜌 = 2⌈log2 (⌈2𝑟max 𝜌max ⌉+1)⌉ .

(15)

Because discrete spatial rings cannot resolve arbitrarily high spatial frequencies, the ring sampling rate 𝑓𝑠 = (𝑁𝑟 −1)∕𝑟max caps the usable bandwidth at: ( ) 𝜌cutof f = min 𝜌max , 21 𝑓𝑠 . (16) Samples beyond 𝜌cutof f are set to zero during precomputation. The runtime lookup works in the discrete bin coordinate of the stored axis: 𝜌̃ = 𝜌

𝑁𝜌 − 1 𝜌max

(17)

.

To prevent boundary artifacts, we apply a smoothstep taper over the final Δtaper = 8 bins below this cutoff: ⎧1, ⎪ 𝑎cutof f (𝜌) ̃ = ⎨𝑡2 (3 − 2𝑡), ⎪ ⎩0,

𝜌̃ ≤ 𝜌̃cutof f − Δtaper , 𝜌̃cutof f − Δtaper < 𝜌̃ < 𝜌̃cutof f , 𝜌̃ ≥ 𝜌̃cutof f ,

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

Sampled spatial radii 𝑟 uniform spacing, trapezoidal weights, sized to the shortest wavelength used

Unit ring basis 𝐵(𝑟, 𝜌) = 2𝜋𝑟 𝐽0 (2𝜋𝑟𝜌) evaluated over spatial radius 𝑟 and Fourier radius 𝜌

Multiply by propagation phase (𝑓 , 𝑑) and integrate over rings

Stored radial kernel 𝐾𝑟 (𝜌, 𝑓 , 𝑑) size 𝜌 × 𝑓 × 𝑑; samples beyond the radial cutoff are zero

Figure 21: Offline precomputation of our symmetry-aware phasor-fields radial kernel. Spatial radii are sampled and turned into a Bessel-function ring basis 𝐵(𝑟, 𝜌), which is multiplied by the propagation phase and integrated over rings to form the stored radial kernel 𝐾𝑟 (𝜌, 𝑓 , 𝑑); both are defined in Section 4.3 of the main paper. This stage runs once, ahead of runtime. Section B.2.1 describes how the runtime lookup then samples this kernel per pixel.

(18) with 𝑡 = (𝜌̃cutof f − 𝜌)∕Δ ̃ taper , where 𝜌̃cutof f is 𝜌cutof f expressed in the same bin coordinate. For the 190 × 190 nonconfocal capture of Nam et al. (2021) used throughout this appendix, where 𝜎𝑠 = 1, these expressions give 𝑁𝑟 = 128 rings (from a Nyquist minimum of 92) and 𝑁𝜌 = 256 stored radii (from a minimum of 191), with 𝜌max = 70.34 m−1 and 𝜌cutof f = 47.02 m−1 . A confocal capture of the same geometry doubles the ring count through 𝜎𝑠 = 2. At runtime, each Cartesian FFT pixel reads a precomputed continuous radius from the pixel-to-radius map, applies the taper weight of Equation 18, and linearly interpolates between the two neighboring 𝜌 bins before accumulating over temporal frequencies. Figure 21 summarizes this offline precomputation phase. Algorithm 4 details the runtime execution loop for this streamed radial phasor-fields kernel. For each temporal frequency, 𝑎cutof f evaluates the low-pass taper (Equation 18) near the radial sampling boundary, while linear_sample interpolates between adjacent radial bins. The taper factor 𝑎 operates independently of the frequency weights 𝑤(𝑓 ), which are applied directly during local accumulation. By bypassing the generation of dense 2D kernel slices at runtime, this approach substantially reduces VRAM consumption and memory bandwidth.

B.2.2. Offline phasor-fields simplification This subsection discusses our separate offline phasorfields implementation used for the preprocessed experiments, different from our streamed pipeline from Section 4.3. We use as a baseline the implementation by Lindell et al. (2019), which is slightly different from that of Nam et al. (2021) Page 21 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

Algorithm 4: Our symmetry-aware phasor-fields dataflow. The dense kernel 𝐾 from Algorithm 3 is replaced by a radial kernel 𝐾𝑟 and a pixel-to-radius map 𝑅. 1 def ours_phasor_fields (𝑃 , 𝐾𝑟 , 𝑅, 𝑁𝑥 , 𝑁𝑦 , 𝑁𝑑 , 𝑁𝑓 , 𝑤): 2 3 4 5 6 7 8 9 10 11 12 13

𝑃̂ ←  {padtl (𝑃 )} ⊳ Top-left padding and batched 2D FFTs for (𝑢, 𝑣, 𝑑) in parallel : (𝑢, 𝑣, 𝑑) ← 0 𝜌 ← 𝑅(𝑢, 𝑣) for 𝑓 ← 1 to 𝑁𝑓 : 𝑎 ← 𝑎cutof f (𝜌) ̃ ⊳ Taper of Equation 18 ℎ ← 𝑎 ⋅ linear_sample(𝐾𝑟 (⋅, 𝑓 , 𝑑), 𝜌) (𝑢, 𝑣, 𝑑) ← (𝑢, 𝑣, 𝑑) + 𝑤(𝑓 ) 𝑃̂(𝑢, 𝑣, 𝑓 ) ℎ 𝑃̂ ←  −1 {} ⊳ Batched inverse FFTs over depth 𝑃 ′ ← |𝑃̂| 𝑓 (𝐱𝑣 ) ← max_z(𝑃 ′ ) return 𝑓 (𝐱𝑣 )

and follows the algorithm of Liu et al. (2019). We use this baseline as reference for all our offline reconstruction experiments. The first memory issue we addressed was the convolution of the propagation kernel 𝐾 of Section 4.3 of the main paper, which acts here as the point spread function of the plane-to-plane propagation, with the signal’s phasor representation. The original implementation allocates two real-valued buffers for the phasor representation, 𝑃cos and 𝑃sin , which hold the cosine and sine components of the phasor field ̂ 𝜔 of the main paper. In this formulation the two buffers are treated as if they were independent signals: each one is Fourier-transformed, multiplied by 𝐾 in the frequency domain, and then inversely transformed, and the final complex result is obtained by combining both outputs. However, these two components are in fact the real and imaginary parts of a single complex phasor. Since the Fourier transform, the pointwise multiplication by 𝐾, and the inverse transform are all linear operations, applying them separately to 𝑃cos and 𝑃sin is equivalent to applying them once to their complex combination. Therefore, we found that the following two expressions are equivalent: 𝑃 ′ =  −1 {𝑃̂cos ⋅ 𝐾} + 𝑖  −1 {𝑃̂sin ⋅ 𝐾} =  −1 {(𝑃̂cos + 𝑖 𝑃̂sin ) ⋅ 𝐾}

(19) (20)

with 𝑃̂cos =  {𝑃cos } and 𝑃̂sin =  {𝑃sin }. Additionally, another time-consuming step is the construction of the transform operator that maps phasor data from a Cartesian layout to a ring-based layout. In the original implementation, this operator is built as a sparse matrix of size 𝑀 2 × 𝑀 with 𝑀 ← max(𝑥, 𝑦) and 𝑀 2 nonzero entries. We index samples from zero throughout, 𝑖 = 0, … ,√ 𝑀 2 − 1, so that the one-based radial index of sample 𝑖 is ⌈ 𝑖 + 1⌉. The reference implementation is written in one-based √ MATLAB indexing, where the same quantity reads ⌈ √ 𝑖⌉. Each sample contributes a value of √1 at position (𝑖, ⌈ 𝑖 + 1⌉), and its row is then scaled by 1∕ 𝑖 + 1, which accounts for the number of samples that fall on the same ring. A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

The matrix is then iteratively collapsed over log2 (max(𝑥, 𝑦)) iterations; in each iteration, pairs of non-contiguous rows are averaged. After log2 𝑀 iterations, the matrix is reduced to size 𝑀 × 𝑀, and the accumulated weights have been scaled 1 by a further factor of ( 21 )log2 𝑀 = 𝑀 . Composing the two scalings, this hierarchical construction is equivalent to directly creating an 𝑀 × 𝑀 matrix and mapping√each sample index 𝑖 to its corresponding ring index ⌈ 𝑖 + 1⌉, with a final weight of √1 : the 𝑀 𝑖+1 √ 1∕ 𝑖 + 1 comes from the per-row scaling and the 1∕𝑀 from the repeated averaging. In other words, the hierarchical averaging performed in the original implementation can be collapsed into a single kernel. Constructing the transform operators as in the original formulation required building large sparse matrices (e.g., using the Eigen library) and repeatedly collapsing them on the CPU. The simplified construction removes this sparse-matrix setup from the critical path and makes the offline implementation easier to run on the GPU. Finally, note that it is not necessary to construct the inverse operator explicitly; the kernel can access it by simply using transposed indices of the forward transform operator. In our offline implementation, the forward and backward matrix multiplications are solved with the cuBLAS module on top of CUDA.

C. Additional results and benchmarks We provide four additional experiments to support our main findings. First, a breakdown of the dynamic capture from Nam et al. (2021) across different configurations (specifically varying the algorithm, precision, and padding). Second, more per-frame error maps (extending Section 5.6, which only shows the mean absolute difference). Third, a sweep of our Mixture of Experts’ (MoE) pooling weight 𝛽. Finally, the complete offline benchmark tables expanding on Section 5.8.

C.1. Dynamic reconstruction across configurations Figure 22 breaks the nlosbox1 capture of Section 5.1 of the main paper down across every configuration. The comparison includes 𝑓 –𝑘 migration, our phasor-fields and MoE, each with and without padding and at both precisions, together with the phasor-fields implementation of Nam et al. (2021), which does not offers a half-precision path nor a padded reconstruction volume and therefore appears as a single configuration. Frame rates here are derived from mean frame times, the same statistic as Table 2 of the main paper, so the configurations the two share report identical values. As in the main results, unpadded 𝑓 –𝑘 migration gives the highest frame rate, but it is also the configuration most exposed to padding artifacts. Our phasor-fields produces sharper letters in this sequence than 𝑓 –𝑘 migration while using far less GPU memory than the reference implementation. The useful reading of this figure is the speed-memory frontier. Unpadded 𝑓 –𝑘 migration sits at the high-FPS end, Page 22 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

Our f-k

Phasor fields (Nam et al.'21) unpadded

unpadded

Our phasor fields

padded

unpadded

Our MoE

padded

unpadded

padded

(a) Qualitative comparison of reconstruction results. For visualization only, the displayed intensities are gamma-corrected (𝛾 = 2) to make low-intensity background noise more visible. Our f-k

Our phasor fields

Method

Prec.

Peak video memory (GB) ↓

3

Padding

↑ FPS

↓ GB

no yes no yes

464.34 68.36 448.48 67.89

0.26 0.48 0.27 0.50

no yes no yes

524.16 22.00 919.85 81.42

0.34 1.21 0.30 1.08

no yes no yes

233.31 34.98 230.90 34.21

0.43 1.60 0.44 1.62

no

10.69

3.34

FP16

Our phasor fields

2

FP32

1 FP16

0

Our 𝑓 –𝑘 FP32

Our MoE

Phasor fields (Nam et al.'21) FP16

3

Our MoE FP32

2 Phasor fields (Nam et al.’21)

1 0

0

200

400

600

800

1000

0

200

400

600

800

1000

(c) Numerical runtime and memory summary.

FPS ↑ Precision FP32

FP16

FP32

FPS: from mean frame time. GB: peak VRAM. Pad.: padded reconstruction volume.

Variant no padding

padding

(b) Response time and memory trade-off. Markers show FPS and peak VRAM. Figure 22: Per-configuration real-time comparison on the nlosbox1 dynamic capture. The top row shows a representative reconstruction, while the bottom row summarizes the runtime-memory trade-off across methods. Frame rates are derived from mean frame times, as in Table 2 of the main paper.

our phasor-fields keeps most of the visual benefit in the low-memory corner, and the reference implementation sits an order of magnitude away on both axes. MoE is therefore better interpreted as a visual-combination option than as the best point on the performance frontier. One configuration runs against the trend: 𝑓 –𝑘 migration in FP 16 is both slower and slightly larger than its FP 32 counterpart, so half precision is not a useful operating point for that branch.

C.2. Additional image-space quality Section 5.6 of the main paper reports the sequence-wide mean absolute difference. Here, we show the per-frame error maps, together with the corresponding FLIP errors. Figure 23 shows the office scene, which contains several structures at different depths and is therefore useful for stress testing, although that also makes local error patterns harder to read. The FP16 and FP32 rows are indistinguishable, so the residual difference is the radial kernel representation rather than the lack of precision from FP16. Our phasor-fields reaches a peak signal-to-noise ratio of 37.0 dB against the dense-kernel reference, and 𝑓 –𝑘 reaches 57.0 dB against its FP32 reference. A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

Reference

Ours

FLIP

Abs. difference

mean 0.048

mean 0.009

mean 0.048

mean 0.009

mean 0.005

mean 0.001

Phasor fields (FP16)

Phasor fields (FP32)

f-k migration (FP16)

Figure 23: Image-space error on the office scene. Each row compares our reconstruction with its reference and reports FLIP and absolute error maps. The reference for the phasor-fields rows is Nam et al. (2021)’s phasor-fields implementation, so the comparison isolates the radial kernel representation. For 𝑓 –𝑘, the reference is 𝑓 –𝑘 using FP32.

C.3. Pooling weight of the mixture Table 8 sweeps the pooling weight 𝛽 (Equation 10) on the nlosbox1 capture using the post-processing setup from Page 23 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction Table 8 Impact of the pooling weight 𝛽 on the 100-frame nlosbox1 sequence (FP32, unpadded). Because every metric changes monotonically with 𝛽, there is no single optimal setting: 𝛽 = 0.5 provides the strongest background rejection but also the weakest object signal. Bold text highlights our chosen operating point (𝛽 = 0.3), which matches the MoE row in Table 6 of the main paper. 𝛽 0.15 0.20 0.25 0.30 0.35 0.40 0.50

↑ Object

↓ Background

↑ Contrast

↑ Sharpness

0.520 0.495 0.471 0.449 0.428 0.408 0.372

8.7 ⋅ 10−5 3.8 ⋅ 10−5 1.7 ⋅ 10−5 7.4 ⋅ 10−6 3.2 ⋅ 10−6 1.4 ⋅ 10−6 2.5 ⋅ 10−7

6,008 12,993 28,186 60,981 132,666 291,566 1,478,556

0.050 0.051 0.052 0.053 0.054 0.055 0.056

video memory, the largest configurations report totals above the 16 GB of our GPU. Across these tables, the largest speed-ups occur when the original implementations allocate and transform large intermediate volumes. The optimized versions are still affected by setup overhead at very small sizes and by memory traffic at the largest sizes, but the throughput remains high across confocal, simulated, and exhaustive inputs. This supports our claim from the main paper that the speed-up comes from dataflow and layout changes rather than from a dataset-specific shortcut.

Section 5.7 of the main paper. As 𝛽 increases, background noise decreases exponentially by roughly 2.3× per step, whereas object intensity drops linearly. Image sharpness increases slightly over the same range. Because the response is monotonic without a distinct inflection point, 𝛽 controls a direct tradeoff between background suppression and object signal retention. We select 𝛽 = 0.3 across all experiments: this setting reduces background noise by three orders of magnitude relative to either standalone branch while preserving 73% of the 𝑓 –𝑘 object intensity, reproducing the MoE row in Table 6 of the main paper. Higher weights yield additional background suppression at the expense of object visibility: at 𝛽 = 1, where pooling simplifies to a plain product, object intensity drops to 36% of the 𝑓 –𝑘 baseline.

C.4. Offline reconstruction The remaining tables give the full offline reconstruction timings. They cover several spatial and temporal resolutions, starting from the original dataset sizes and then downsampling the spatial or temporal dimensions. The 𝑓 –𝑘 migration datasets from Lindell et al. (2019) are 5123 , while the largest Zaragoza datasets (Galindo et al., 2019) are 256×256×4096. We limit the temporal dimension to 2048 bins in these experiments because 4096 bins put substantial memory pressure on both the original and optimized implementations. Besides reconstruction time, we report throughput in millions of voxels per second. For the exhaustive captures, the dimensions column lists the original measurement grid. Normal moveout correction turns it into 31 × 31 × 𝑇 confocal data, which is the format used by the throughput figures. The last two columns report the maximum observed memory footprint of each implementation’s main buffers. We label the metric tracked memory rather than VRAM because the two sides measure different resources: the original MATLAB implementations allocate their main buffers in CPU memory, whereas for our implementations the reported value is the peak of a process-scoped counter over live CUDA allocations, which counts pinned host staging buffers alongside device memory. Since that counter is not restricted to dedicated A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

Page 24 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction Table 9 Offline performance and memory comparison between our accelerated 𝑓 –𝑘 and phasor-fields implementations and their original MATLAB counterparts, across every measured resolution. Throughput is the reconstructed volume divided by reconstruction time. For exhaustive captures, the dimensions column lists the original measurement grid before normal moveout correction, which produces a 31×31×𝑇 confocal transient volume. Max. tracked memory reports the peak footprint of each implementation’s main buffers: host RAM for the original MATLAB code, and the peak of all tracked CUDA allocations for ours, counting device buffers and pinned host staging together (see the discussion preceding these tables). Bold marks our value when it beats the baseline. ↓ Time (s) Dataset

↑ Throughput (Mvox/s)

↓ Max. tracked memory (GB)

Ours

Original

Ours

Original

63.50 35.40

0.92 4.35

17.000 22.001

52.430 33.579

20.815 ± 1.51 14.793 ± 0.35

237.72 157.83

3.22 4.54

8.500 11.000

26.214 16.778

0.143 ± 0.00 0.217 ± 0.00

10.695 ± 0.66 6.894 ± 0.20

234.41 154.81

3.14 4.87

4.250 5.501

13.107 8.393

𝑓 –𝑘 migration Phasor fields

0.141 ± 0.00 0.213 ± 0.00

10.364 ± 0.37 7.537 ± 0.27

237.99 157.66

3.24 4.45

4.250 5.500

13.107 8.389

256×256×256

𝑓 –𝑘 migration Phasor fields

0.073 ± 0.00 0.110 ± 0.00

5.243 ± 0.14 3.428 ± 0.15

230.17 152.22

3.20 4.89

2.125 2.750

6.554 4.194

128×128×512

𝑓 –𝑘 migration Phasor fields

0.038 ± 0.00 0.057 ± 0.00

2.552 ± 0.11 1.607 ± 0.02

219.35 147.91

3.29 5.22

1.062 1.376

3.277 2.101

256×256×128

𝑓 –𝑘 migration Phasor fields

0.038 ± 0.00 0.057 ± 0.00

2.650 ± 0.18 1.751 ± 0.06

222.97 146.93

3.17 4.79

1.062 1.375

3.277 2.097

128×128×256

𝑓 –𝑘 migration Phasor fields

0.020 ± 0.00 0.030 ± 0.00

1.253 ± 0.03 0.801 ± 0.03

207.55 138.92

3.35 5.24

0.531 0.688

1.638 1.049

128×128×128

𝑓 –𝑘 migration Phasor fields

0.012 ± 0.00 0.016 ± 0.00

0.652 ± 0.02 0.408 ± 0.02

181.34 128.48

3.22 5.14

0.266 0.344

0.819 0.524

512×512×512

𝑓 –𝑘 migration Phasor fields

2.016 ± 0.01 3.916 ± 0.08

154.730 ± 14.77 31.458 ± 2.65

66.59 34.28

0.87 4.27

17.000 22.001

52.433 33.579

512×512×256

𝑓 –𝑘 migration Phasor fields

0.280 ± 0.00 0.427 ± 0.00

21.723 ± 0.88 15.453 ± 0.38

239.68 157.32

3.09 4.34

8.500 11.000

26.214 16.777

256×256×512

𝑓 –𝑘 migration Phasor fields

0.144 ± 0.00 0.218 ± 0.00

11.462 ± 0.65 7.155 ± 0.27

232.97 153.91

2.93 4.69

4.250 5.501

13.109 8.393

512×512×128

𝑓 –𝑘 migration Phasor fields

0.142 ± 0.00 0.212 ± 0.00

10.528 ± 0.08 7.910 ± 0.17

236.05 158.02

3.19 4.24

4.250 5.500

13.107 8.389

256×256×256

𝑓 –𝑘 migration Phasor fields

0.073 ± 0.00 0.111 ± 0.00

5.597 ± 0.20 3.568 ± 0.11

228.38 151.14

3.00 4.70

2.125 2.750

6.554 4.194

128×128×512

𝑓 –𝑘 migration Phasor fields

0.039 ± 0.00 0.057 ± 0.00

2.713 ± 0.21 1.621 ± 0.06

217.19 146.98

3.09 5.17

1.062 1.376

3.277 2.101

256×256×128

𝑓 –𝑘 migration Phasor fields

0.038 ± 0.00 0.057 ± 0.00

2.857 ± 0.07 1.917 ± 0.07

218.52 146.45

2.94 4.38

1.062 1.375

3.277 2.097

128×128×256

𝑓 –𝑘 migration Phasor fields

0.021 ± 0.00 0.031 ± 0.00

1.435 ± 0.18 0.822 ± 0.06

200.98 136.51

2.92 5.10

0.531 0.688

1.638 1.049

128×128×128

𝑓 –𝑘 migration Phasor fields

0.012 ± 0.00 0.017 ± 0.00

0.672 ± 0.04 0.409 ± 0.01

180.10 125.27

3.12 5.13

0.266 0.344

0.819 0.524

512×512×512

𝑓 –𝑘 migration Phasor fields

2.018 ± 0.01 3.758 ± 0.02

152.567 ± 8.22 32.195 ± 1.86

66.51 35.72

0.88 4.17

17.000 22.001

52.430 33.579

512×512×256

𝑓 –𝑘 migration Phasor fields

0.280 ± 0.00 0.427 ± 0.00

21.298 ± 0.59 14.895 ± 0.40

239.83 157.20

3.15 4.51

8.500 11.000

26.214 16.778

256×256×512

𝑓 –𝑘 migration Phasor fields

0.144 ± 0.00 0.217 ± 0.00

11.118 ± 0.87 7.063 ± 0.17

233.54 154.43

3.02 4.75

4.250 5.501

13.107 8.393

512×512×128

𝑓 –𝑘 migration Phasor fields

0.141 ± 0.00 0.213 ± 0.00

10.614 ± 0.21 7.774 ± 0.14

237.14 157.64

3.16 4.32

4.250 5.500

13.107 8.389

256×256×256

𝑓 –𝑘 migration Phasor fields

0.073 ± 0.00 0.111 ± 0.00

5.340 ± 0.21 3.476 ± 0.10

229.77 151.73

3.14 4.83

2.125 2.750

6.554 4.195

128×128×512

𝑓 –𝑘 migration Phasor fields

0.038 ± 0.00 0.057 ± 0.00

2.747 ± 0.16 1.668 ± 0.04

221.81 147.39

3.05 5.03

1.062 1.376

3.277 2.101

256×256×128

𝑓 –𝑘 migration Phasor fields

0.038 ± 0.00 0.057 ± 0.00

2.856 ± 0.25 1.820 ± 0.06

220.54 146.61

2.94 4.61

1.062 1.375

3.277 2.097

128×128×256

𝑓 –𝑘 migration Phasor fields

0.020 ± 0.00 0.030 ± 0.00

1.316 ± 0.06 0.830 ± 0.04

204.85 139.19

3.19 5.05

0.531 0.688

1.638 1.049

128×128×128

𝑓 –𝑘 migration Phasor fields

0.011 ± 0.00 0.016 ± 0.00

0.686 ± 0.14 0.392 ± 0.00

182.86 128.49

3.06 5.35

0.266 0.344

0.819 0.524

Dimensions

Algorithm

Ours

Original

512×512×512

𝑓 –𝑘 migration Phasor fields

2.114 ± 0.03 3.791 ± 0.03

146.624 ± 8.47 30.872 ± 0.67

512×512×256

𝑓 –𝑘 migration Phasor fields

0.282 ± 0.00 0.425 ± 0.00

256×256×512

𝑓 –𝑘 migration Phasor fields

512×512×128

Confocal – 𝑓 –𝑘 migration dataset (Lindell et al., 2019)

bike

teaser

statue

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

Page 25 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

Table 10 Continuation of Table 9. ↓ Time (s) Dataset

Dimensions

Algorithm

Ours

Original

↑ Throughput (Mvox/s)

↓ Max. tracked memory (GB)

Ours

Original

Ours

Original

Confocal – Zaragoza dataset (Galindo et al., 2019) 256×256×1964

𝑓 –𝑘 migration Phasor fields

1.747 ± 0.00 4.172 ± 0.01

139.592 ± 22.22 31.655 ± 1.26

73.70 30.85

0.92 4.07

16.303 21.112

52.429 33.621

256×256×982

𝑓 –𝑘 migration Phasor fields

0.280 ± 0.00 0.428 ± 0.00

20.740 ± 1.52 14.870 ± 0.24

230.05 150.21

3.10 4.33

8.151 10.552

26.214 16.794

128×128×1964

𝑓 –𝑘 migration Phasor fields

0.142 ± 0.00 0.222 ± 0.00

10.361 ± 0.48 8.209 ± 0.05

226.13 144.67

3.11 3.92

4.076 5.289

13.107 8.454

256×256×491

𝑓 –𝑘 migration Phasor fields

0.141 ± 0.00 0.220 ± 0.01

10.714 ± 0.29 6.665 ± 0.28

228.16 146.26

3.00 4.83

4.076 5.275

13.107 8.393

128×128×982

𝑓 –𝑘 migration Phasor fields

0.073 ± 0.00 0.109 ± 0.00

4.976 ± 0.32 3.424 ± 0.17

221.74 147.40

3.23 4.70

2.038 2.641

6.554 4.211

128×128×491

𝑓 –𝑘 migration Phasor fields

0.039 ± 0.00 0.057 ± 0.00

2.468 ± 0.18 1.611 ± 0.04

206.95 140.95

3.26 4.99

1.019 1.320

3.277 2.101

64×64×1964

𝑓 –𝑘 migration Phasor fields

0.039 ± 0.00 0.059 ± 0.00

2.590 ± 0.13 3.285 ± 0.09

206.38 136.22

3.11 2.45

1.019 1.333

3.277 2.163

64×64×982

𝑓 –𝑘 migration Phasor fields

0.022 ± 0.00 0.031 ± 0.00

1.509 ± 0.23 0.987 ± 0.03

181.21 131.32

2.67 4.08

0.509 0.663

1.638 1.065

64×64×491

𝑓 –𝑘 migration Phasor fields

0.014 ± 0.00 0.018 ± 0.00

0.605 ± 0.05 0.449 ± 0.04

148.94 110.21

3.32 4.48

0.255 0.331

0.819 0.528

256×256×1964

𝑓 –𝑘 migration Phasor fields

1.734 ± 0.00 4.169 ± 0.02

138.419 ± 8.19 33.418 ± 0.67

74.24 30.88

0.93 3.85

16.303 21.112

52.429 33.621

256×256×982

𝑓 –𝑘 migration Phasor fields

0.280 ± 0.00 0.430 ± 0.00

22.058 ± 0.97 14.709 ± 0.33

229.92 149.75

2.92 4.38

8.151 10.552

26.214 16.793

128×128×1964

𝑓 –𝑘 migration Phasor fields

0.143 ± 0.00 0.223 ± 0.00

10.740 ± 0.57 8.290 ± 0.14

225.81 144.02

3.00 3.88

4.076 5.289

13.107 8.454

256×256×491

𝑓 –𝑘 migration Phasor fields

0.140 ± 0.00 0.213 ± 0.00

10.384 ± 0.18 7.006 ± 0.26

229.36 151.12

3.10 4.59

4.076 5.275

13.107 8.393

128×128×982

𝑓 –𝑘 migration Phasor fields

0.072 ± 0.00 0.109 ± 0.00

5.020 ± 0.15 3.410 ± 0.07

222.54 147.34

3.20 4.72

2.038 2.641

6.554 4.211

128×128×491

𝑓 –𝑘 migration Phasor fields

0.038 ± 0.00 0.057 ± 0.00

2.552 ± 0.03 1.556 ± 0.06

209.20 141.99

3.15 5.17

1.019 1.320

3.277 2.101

64×64×1964

𝑓 –𝑘 migration Phasor fields

0.039 ± 0.00 0.059 ± 0.00

2.633 ± 0.23 3.293 ± 0.10

207.03 136.40

3.06 2.44

1.019 1.333

3.277 2.163

64×64×982

𝑓 –𝑘 migration Phasor fields

0.021 ± 0.00 0.031 ± 0.00

1.232 ± 0.08 1.026 ± 0.04

187.23 129.02

3.26 3.92

0.509 0.663

1.638 1.065

64×64×491

𝑓 –𝑘 migration Phasor fields

0.014 ± 0.00 0.018 ± 0.00

0.679 ± 0.08 0.420 ± 0.01

148.54 110.09

2.96 4.79

0.255 0.331

0.819 0.528

256×256×1964

𝑓 –𝑘 migration Phasor fields

1.780 ± 0.08 4.208 ± 0.02

142.930 ± 6.27 32.743 ± 1.42

72.29 30.58

0.90 3.93

16.303 21.112

52.429 33.636

256×256×982

𝑓 –𝑘 migration Phasor fields

0.290 ± 0.01 0.430 ± 0.00

23.453 ± 0.80 14.576 ± 0.38

222.01 149.71

2.74 4.42

8.151 10.552

26.214 16.794

128×128×1964

𝑓 –𝑘 migration Phasor fields

0.147 ± 0.00 0.223 ± 0.00

10.414 ± 0.36 8.356 ± 0.20

218.17 144.47

3.09 3.85

4.076 5.289

13.107 8.454

256×256×491

𝑓 –𝑘 migration Phasor fields

0.144 ± 0.00 0.212 ± 0.00

10.379 ± 0.20 6.895 ± 0.36

222.93 151.78

3.10 4.67

4.076 5.275

13.107 8.400

128×128×982

𝑓 –𝑘 migration Phasor fields

0.077 ± 0.00 0.109 ± 0.00

5.169 ± 0.11 3.464 ± 0.04

208.14 147.30

3.11 4.64

2.038 2.641

6.554 4.211

128×128×491

𝑓 –𝑘 migration Phasor fields

0.041 ± 0.00 0.057 ± 0.00

2.933 ± 0.21 1.567 ± 0.06

196.97 140.21

2.74 5.13

1.019 1.320

3.277 2.101

64×64×1964

𝑓 –𝑘 migration Phasor fields

0.039 ± 0.00 0.059 ± 0.00

2.525 ± 0.16 3.268 ± 0.04

206.03 135.43

3.19 2.46

1.019 1.333

3.277 2.163

64×64×982

𝑓 –𝑘 migration Phasor fields

0.022 ± 0.00 0.031 ± 0.00

1.270 ± 0.08 1.028 ± 0.01

185.57 128.08

3.17 3.91

0.509 0.663

1.638 1.065

64×64×491

𝑓 –𝑘 migration Phasor fields

0.014 ± 0.00 0.018 ± 0.00

0.642 ± 0.07 0.416 ± 0.01

144.58 111.45

3.13 4.83

0.255 0.331

0.819 0.528

usaf

bunny

z

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

Page 26 of 27

Memory-efficient GPU pipelines for real-time non-line-of-sight reconstruction

Table 11 Continuation of Table 10. ↓ Time (s) Dataset

Dimensions

Algorithm

Ours

Original

↑ Throughput (Mvox/s)

↓ Max. tracked memory (GB)

Ours

Ours

Original

Original

Tal-generated dataset (Royo and Luesia-Lahoz, 2024), inspired by a Phasor Fields scene (Liu et al., 2020) 360×260×1500

𝑓 –𝑘 migration Phasor fields

2.945 ± 0.05 5.414 ± 0.02

507.808 ± 47.68 138.281 ± 13.87

47.67 25.93

0.28 1.02

17.784 23.023

75.938 66.424

360×260×750

𝑓 –𝑘 migration Phasor fields

0.308 ± 0.01 0.458 ± 0.00

64.483 ± 7.22 28.830 ± 1.34

227.77 153.27

1.09 2.43

8.893 11.510

37.969 33.194

180×130×1500

𝑓 –𝑘 migration Phasor fields

0.150 ± 0.00 0.232 ± 0.00

15.883 ± 1.75 14.810 ± 0.46

234.54 151.47

2.21 2.37

4.446 5.762

18.984 16.654

360×260×375

𝑓 –𝑘 migration Phasor fields

0.153 ± 0.00 0.228 ± 0.00

14.011 ± 0.08 13.549 ± 0.21

229.12 153.85

2.51 2.59

4.447 5.755

18.984 16.593

180×130×750

𝑓 –𝑘 migration Phasor fields

0.076 ± 0.00 0.115 ± 0.00

7.340 ± 0.39 6.564 ± 0.12

229.82 152.06

2.39 2.67

2.223 2.879

9.492 8.311

180×130×375

𝑓 –𝑘 migration Phasor fields

0.040 ± 0.00 0.059 ± 0.00

3.350 ± 0.07 3.101 ± 0.02

221.61 148.09

2.62 2.83

1.112 1.439

4.746 4.151

90×65×1500

𝑓 –𝑘 migration Phasor fields

0.041 ± 0.00 0.061 ± 0.00

3.818 ± 0.14 4.814 ± 0.02

214.11 142.68

2.30 1.82

1.112 1.447

4.746 4.213

90×65×750

𝑓 –𝑘 migration Phasor fields

0.022 ± 0.00 0.032 ± 0.00

1.809 ± 0.04 1.782 ± 0.02

203.46 137.10

2.43 2.46

0.556 0.721

2.373 2.090

90×65×375

𝑓 –𝑘 migration Phasor fields

0.012 ± 0.00 0.018 ± 0.00

0.912 ± 0.06 0.806 ± 0.01

176.53 122.81

2.41 2.72

0.278 0.360

1.187 1.041

16×16×16×16×2048

𝑓 –𝑘 migration Phasor fields

0.018 ± 0.00 0.027 ± 0.00

0.588 ± 0.04 2.107 ± 0.01

107.02 68.74

3.21 0.90

0.352 0.436

0.770 0.558

16×16×16×16×1024

𝑓 –𝑘 migration Phasor fields

0.011 ± 0.00 0.016 ± 0.00

0.305 ± 0.02 0.479 ± 0.01

88.41 60.09

3.09 1.97

0.176 0.215

0.384 0.263

16×16×16×16×512

𝑓 –𝑘 migration Phasor fields

0.007 ± 0.00 0.012 ± 0.00

0.157 ± 0.02 0.153 ± 0.00

67.26 38.93

3.01 3.08

0.088 0.106

0.192 0.127

16×16×16×16×2048

𝑓 –𝑘 migration Phasor fields

0.015 ± 0.00 0.024 ± 0.00

0.744 ± 0.05 4.263 ± 0.15

122.69 78.85

2.53 0.44

0.351 0.436

0.773 0.562

16×16×16×16×1024

𝑓 –𝑘 migration Phasor fields

0.008 ± 0.00 0.015 ± 0.00

0.372 ± 0.09 0.491 ± 0.00

112.12 63.70

2.53 1.92

0.176 0.214

0.384 0.269

16×16×16×16×512

𝑓 –𝑘 migration Phasor fields

0.008 ± 0.00 0.010 ± 0.00

0.159 ± 0.01 0.158 ± 0.01

59.02 48.51

2.97 2.99

0.088 0.106

0.192 0.127

16×16×16×16×2048

𝑓 –𝑘 migration Phasor fields

0.019 ± 0.00 0.039 ± 0.02

0.685 ± 0.10 2.187 ± 0.03

100.31 48.33

2.76 0.86

0.352 0.436

0.769 0.560

16×16×16×16×1024

𝑓 –𝑘 migration Phasor fields

0.012 ± 0.00 0.016 ± 0.00

0.348 ± 0.03 0.490 ± 0.01

76.51 60.86

2.71 1.93

0.176 0.215

0.384 0.262

16×16×16×16×512

𝑓 –𝑘 migration Phasor fields

0.008 ± 0.00 0.010 ± 0.00

0.159 ± 0.02 0.158 ± 0.00

56.80 45.54

2.97 2.99

0.088 0.106

0.192 0.127

office

Exhaustive – Zaragoza dataset (Galindo et al., 2019) concavities

t (in a box)

bunny

A. López-Ruiz and D. Royo: Preprint submitted to Elsevier

Page 27 of 27

Related documents

• Non-common path aberration compensation and a dark hole loop with a pyramid adaptive optics system: Application to SAXO+#405867
• 4D beam matrix reconstruction in particle accelerators#495726
• Accelerating Presto with GPUs#303302
• AsyncSparse: Accelerating Sparse Matrix-Matrix Multiplication on Asynchronous GPU Architectures#120471
• GPU-RMQ: Accelerating Range Minimum Queries on Modern GPUs#2775