ECCV 2026 3D Vision Autonomous Driving single-photon LiDAR SPAD dead time Markov-renewal process photon count statistics histogram cube simulation

Markov-Renewal Single-Photon LiDAR Simulator¶

Conference: ECCV 2026
Paper: ECCV 2026 / Project page
Area: 3D Vision / Autonomous Driving
Keywords: single-photon LiDAR, SPAD dead time, Markov-renewal process, photon count statistics, histogram cube simulation

TL;DR¶

MaRS formulates photon registration under SPAD dead time as a Markov-renewal process (MRP), analytically predicting each pixel's photon-count distribution (Gaussian mean and variance) and its temporal PDF; delay invariance plus a precomputed lookup table then synthesize a whole H×W×n_b histogram cube in one parallel sampling pass — statistically indistinguishable from per-photon rejection-sampling gold standard, yet generating a 128×128×1024 cube in 1.01 s versus 40.8 s (>40×).

Background & Motivation¶

Single-photon LiDAR (SP-LiDAR) pairs SPAD arrays with TCSPC timing to accumulate per-laser-period photon timestamps into a histogram, achieving single-photon-sensitivity, high-resolution 3D imaging. That sensitivity comes with awkward measurement statistics: photon arrivals are themselves a random point process, and the detector state feeds back into that process — every registered photon puts the SPAD into a dead time of length \(t_d\) during which subsequent arrivals are simply lost (in asynchronous operation the detector only reopens once \(t_d\) has elapsed). The effect is nonlinear: it destroys statistical independence between detections and creates temporal correlations, and it simultaneously suppresses the registration rate and distorts the histogram shape (the classic nonlinear pile-up). Any Poisson / multinomial simulator that ignores dead time is therefore physically wrong.

Existing simulators are trapped in a dilemma. On one side, Poisson-style models are fast but physically incorrect: the high-flux histograms they produce show no pile-up distortion and their background photons spread uniformly across the period. On the other side, sequential "gold-standard" simulators first draw all photon arrivals via multinomial Poisson sampling and then reject, photon by photon, any detection falling inside a dead-time window. They are physically correct but need per-photon, per-pixel loops, so they collapse as soon as flux, spatial resolution, or the number of realizations grows — a 128×128×1024 cube takes 40.8 s, while deep-learning training needs thousands of such cubes. That is precisely the bottleneck this paper removes: because scalable high-fidelity simulation is unavailable, SP-LiDAR reconstruction networks are trained on simplified low-flux simulators or on physically wrong models. Recently Zhang et al. tried to sidestep the per-photon loop by learning the dead-time-distorted statistics with a neural network, but their temporal-PDF estimator must be retrained for every dead-time setting, their count prediction is visibly inaccurate in variance, and the pipeline remains pixel-wise rather than fully parallel.

The key tension is that simulating one histogram needs two things — a temporal PDF and a photon count — and while the former has been solved, the latter has never been done right. The Markov chain model of Rapp et al. (later accelerated by Zhang et al. through a "dead time as row shifting" reformulation) accurately yields the stationary temporal PDF, but it only describes which bin the next registration lands in; the absolute times \(T_k\) are discarded, so it cannot answer how many photons were registered over an exposure, and therefore cannot fix the histogram's scale. Classical renewal theory does not help either, because it requires i.i.d. inter-registration times whereas dead time makes successive intervals strongly dependent, and self-exciting point processes are computationally intractable. The core idea of this paper is to recast registration as a Markov-renewal process (MRP): state transitions supply the relative timestamps while holding times restore the absolute time, so that "count × shape" can be solved separately — the count is the inverse of the cumulative holding time, and mixing guarantees a central limit theorem, giving a Gaussian limiting law with closed-form mean and variance; delay invariance then compresses a three-dimensional lookup table into two dimensions, and the entire cube is drawn in one pass of batched memory reads plus cheap shifts.

Method¶

Overall Architecture¶

MaRS takes as input a system configuration \(\theta=(t_r, N_r, \sigma_t, t_d, n_b, \tau, S, B, N_{\text{iter}})\), where \(t_r\) is the laser repetition period, \(N_r\) the number of pulses (total exposure \(t=N_r t_r\)), \(\sigma_t\) the pulse width, \(t_d\) the dead time and \(n_b\) the number of histogram bins, together with scene parameters — a depth image \(Z\), a reflectivity image \(R\), and a background level \(B\). It outputs \(N_{\text{iter}}\) histogram cubes of size \(H\times W\times n_b\). Unlike a sequential simulator that generates every photon and then rejects it, MaRS factors a histogram into "count × shape": the shape comes from the stationary temporal PDF of a Markov chain, the count from the newly proposed MRP count model. Both depend only on per-pixel \((S,B,t_d)\) and the delay \(\tau\), so they can be precomputed, looked up, and sampled in parallel. The pipeline is: a statistical engine computes each pixel's count distribution and temporal PDF, spectral truncation makes the variance in that distribution tractable, delay invariance collapses the lookup table from 3D to 2D, and batched lookup plus shifting plus two-stage sampling produces the whole cube at once.

%%{init: {'flowchart': {'rankSpacing': 24, 'nodeSpacing': 28, 'padding': 6, 'wrappingWidth': 400}}}%%
flowchart TD
    A["Input: depth Z / reflectivity R / background B<br/>+ system parameters θ"] --> B["Statistical engine: shape + count<br/>Markov-chain stationary PDF ＋ MRP count"]
    B --> C["Spectral truncation<br/>long-lag covariance closed form"]
    C --> D["Lookup table + delay invariance<br/>statistics decoupled from delay"]
    D --> E["Parallel cube sampling<br/>count × shape in one pass"]
    E --> F["Histogram cube<br/>→ depth / reflectivity reconstruction"]

Key Designs¶

1. MRP photon-count model: recovering the absolute time the Markov chain threw away

The Markov chain model above solves the temporal PDF accurately, but its state is the relative timestamp \(X_k = T_k \bmod t_r\); the absolute time \(T_k\) is dropped entirely, so "how many photons arrived in total" is unanswerable — the histogram shape is right while its radiometric scale is still wrong, and that scale directly biases reflectivity estimation. The fix is to add an absolute-time dimension to that chain: registration becomes an MRP \(\{(X_k, W_k)\}\), where \(X_k\) remains the Markov state (relative timestamp) and \(W_k = T_k - T_{k-1}\) is the holding time on the absolute scale, the two being modeled jointly. The transition kernel is the old transition matrix plus one conditional holding-time distribution:

\[Q_{ij}(t) = \Pr\!\left(X_{k+1}=j,\; W_{k+1}\le t \mid X_k = i\right) = \tilde{P}_{ij}\,F_{ij}(t)\]

where \(\tilde{P}_{ij}\) is the existing Markov transition probability and \(F_{ij}(t)\) is the conditional CDF of the holding time for transition \(i\to j\). The clever part is that the count \(N(t)\) is the inverse of the cumulative holding time \(T_n=\sum_{k\le n}W_k\): denser registration means smaller increments \(W_k\). Although the \(W_k\) are dependent (because of dead time) and therefore violate the i.i.d. requirement of classical renewal theory, the underlying Markov states mix fast enough that \(T_n\) still obeys a central limit theorem; inverting that relation yields an asymptotically Gaussian law for the photon count:

\[N(t) \xrightarrow{\;d\;} \mathcal{N}\!\left(\frac{t}{\mu},\; \frac{t\,\sigma^2}{\mu^3}\right)\quad (t\to\infty),\qquad \mu \triangleq \lim_{n\to\infty}\frac{\mathbb{E}[T_n]}{n},\quad \sigma^2 \triangleq \lim_{n\to\infty}\frac{\mathrm{Var}[T_n]}{n}\]

(The \(\mu^3\) denominator in the variance comes from inverting the counting process; ⚠️ the cached text is corrupted at this equation, refer to the original paper.) What remains is to write \(\mu,\sigma^2\) in closed form. Conditioned on a transition \(i\to j\), the holding time has two parts — the minimum awake interval \(t_{i'j}\) from the reactivation instant \(x_i'\) to the next registration \(x_j\), plus a number of complete laser periods of waiting; since the probability that at least one photon arrives within a period is \(1-e^{-\Lambda}\) (with \(\Lambda=S+B\) the total photon flux per period), that number of full periods is geometrically distributed, giving

\[\mu_{ij} = t_{i'j} + t_r\frac{e^{-\Lambda}}{1-e^{-\Lambda}},\qquad \sigma^2_{ij} = t_r^2\frac{e^{-\Lambda}}{\left(1-e^{-\Lambda}\right)^2}\]

(⚠️ both expressions of this lemma are corrupted in the cached text; they are reconstructed here from the geometric-waiting physics — refer to the original paper.) From these one builds the first-moment kernel \(Q_\mu=\tilde P\odot M\), the second-moment kernel \(Q_{\mu^2}=\tilde P\odot M\odot M\) and the variance kernel \(Q_{\sigma^2}=\tilde P\odot\Sigma\) (\(\odot\) is element-wise multiplication and \(M,\Sigma\) collect all \(\mu_{ij},\sigma^2_{ij}\)); the effective mean is then \(\mu=\pi Q_\mu \mathbf{1}\), where \(\pi\) is the stationary distribution of the Markov chain from \(\pi\tilde P=\pi\), and the effective variance splits into a steady-state term plus covariance terms, \(\sigma^2=\sigma^2_{ss}+\gamma_{ss}\) with \(\sigma^2_{ss}=\pi(Q_{\mu^2}+Q_{\sigma^2})\mathbf{1}\). The count thus has an analytical answer for the first time, and it shares the same \(\tilde P\) and the same stationary distribution as the shape — the two halves of a histogram are consistent by construction rather than glued together.

2. Spectral truncation: a convergent closed-form tail for infinitely long correlations

\(\gamma_{ss}=\sum_{l\ge1}\gamma_l\) is an infinite series that must be truncated to be computable, yet the decay profile of \(\gamma_l\) varies enormously across \((S,B,t_d)\) (Fig. 4(a)): a fixed cutoff \(c\) is either too large (slow) or too small (biased), and no universal \(c\) exists. The paper borrows the observation of Zhang et al. that the long-term dynamics of \(\tilde P\) are governed by a few dominant eigenvalues: short-range covariance terms depend on all eigenmodes but there are few of them and they are cheap to sum exactly, whereas long-range correlations are dictated by the slowest-decaying modes and can be captured analytically by the low-rank approximation \(\tilde P^{\,l-1}\approx V_p\Lambda_p^{\,l-1}U_p^\top\) (with \(\{(\lambda_i,v_i,u_i)\}_{i=1}^p\) the leading eigenpairs satisfying the biorthogonality \(U_p^\top V_p=I_p\)):

\[\hat\gamma_{ss}(L,p) = \sum_{l=1}^{L}\gamma_l \;+\; \sum_{i=2}^{p}\frac{\alpha_i\beta_i\,\lambda_i^{\,L}}{1-\lambda_i},\qquad \alpha_i = (\pi Q_\mu V_p)_i,\quad \beta_i=(U_p^\top Q_{\mu^2})_i\]

The first term covers exactly computed short-range correlations, the second is a closed-form estimate of the tail from the dominant modes. The payoff is not a better approximation but a change of complexity class: an uncomputable infinite sum becomes constant-cost. Empirically, eigenvalues beyond the fifth have magnitude \(\lesssim 0.7\), so their contribution is negligible after \(L\approx6\) steps, and the paper fixes \(L=6\), \(p=5\) across all tested flux / dead-time regimes. The control experiment is convincing: in several difficult settings, even summing 20 covariance terms deterministically still fails to capture the long-range tail (the variance does not converge), whereas the spectral rule recovers it at negligible cost.

3. Lookup table + delay invariance: collapsing a 3D table into 2D with one physical symmetry

Even with both designs above making the statistics computable, evaluating them on the fly does not scale: accurate count statistics at native temporal resolution (e.g. \(n_b=2^{14}\)) are expensive, and since every pixel has its own \((S,B,\tau)\), recomputing per pixel simply relocates the disease of the sequential simulator. The natural idea is a precomputed lookup table (LUT), but tabulating all three variables explodes: with \(n_S=256\), \(n_B=100\), \(n_\tau=1024\) and \(n_b=1024\), each entry stores \((\mu,\sigma^2)\) plus a length-\(n_b\) PDF \(\pi\) (\(n_b+2=1026\) floats), i.e. \(256\cdot100\cdot1024\cdot1026\cdot4\approx 108\) GB — impractical.

The turning point is a physical symmetry: delay only shifts time, it does not change the renewal structure. Writing the per-pixel arrival rate as \(\lambda_\tau(t)=S\,s(t-\tau)+B/t_r\), changing \(\tau\) merely moves the time origin of the flux; the stochastic arrival mechanism and the way dead time acts are untouched, so the registered count's mean and variance are independent of \(\tau\) while the temporal PDF is merely a cyclic shift:

\[\mu_\tau=\mu_{\tau'},\qquad \sigma^2_\tau=\sigma^2_{\tau'},\qquad \boldsymbol{\pi}_{\tau'}=\mathrm{shift}_{(\tau'-\tau)}\left(\boldsymbol{\pi}_\tau\right)\]

(Distance-dependent attenuation is absorbed into \(S\) (or \(B\)) as a preprocessing step before applying the theorem.) The LUT can therefore be built over \((S,B)\) alone at a canonical delay (pulse centered at \(t_r/2\)), leaving \(\tau\) to a single shift at sampling time: storage drops from 108 GB to \(256\cdot100\cdot1026\cdot4\text{ B}\approx100\) MB, small enough to stay resident and be accessed in batches. Cube generation becomes three steps — batched lookup of \((\mu,\sigma^2,\pi_0)\), a shift of the PDF by \(Z-t_r/2\), then drawing count maps \(n_i\sim\mathcal{N}(\mu,\sigma^2)\) followed by histograms \(H_i\sim\text{Multinomial}(n_i,\Pi)\) — all vectorized memory reads, shifts and samples, with no per-photon or per-pixel loop. The table is not universal: it is bound to a given system configuration \((t_r,n_b,\sigma_t,t_d)\) and pulse shape \(s(\cdot)\), and must be regenerated for other operating regimes (or other \((S,B)\) ranges and discretizations).

A Worked Example¶

Take a 128×128 depth and reflectivity image, a repetition period \(t_r\) corresponding to a 7.5 MHz pulse train, \(n_b=1024\), and a lookup table built over \(S\in[0.01,10]\) (256 points) and \(B\in[0.01,10]\) (100 points), occupying 100.4 MB. For one pixel, its reflectivity \(R\) and background \(B\) determine \((S,B)\), which either hits a table grid point or falls between neighbors; either way the pixel retrieves its \((\mu,\sigma^2,\pi_0)\) — and \(\pi_0\) is still the "pulse-centered" shape. The pixel's depth \(Z\) converts to a delay \(\tau=2Z/c\), which cyclically shifts \(\pi_0\) by \(\tau-t_r/2\) into that pixel's true echo-time distribution. To generate one realization, an integer photon count \(n\) is drawn from the pixel's Gaussian count distribution, then a single multinomial draw against the shifted \(\pi\) produces its 1024-bin histogram. All 16384 pixels do this simultaneously, repeated \(N_{\text{iter}}\) times. The three measured stages take \(0.020\pm0.005\) s (batched lookup), \(0.144\pm0.008\) s (temporal alignment) and \(0.849\pm0.070\) s (histogram sampling), \(1.013\pm0.071\) s in total, against \(40.831\pm8.045\) s for the sequential gold standard on the same cube.

Loss & Training¶

MaRS is a simulator with no trainable parameters and no loss; its relevance to training lies in supplying data for downstream learned reconstruction. The experiments use a DDPM U-Net modified to ingest the full spatio-temporal histogram cube, and the same architecture is trained separately for depth and for reflectivity recovery until convergence. Four versions (Poisson, Renewal, Zhang et al., MaRS) share that architecture and the same LUT specification, differing only in which simulator generated their training data. Evaluation runs twice: first on each simulator's own data to confirm stable learning, then all models on gold-standard data, which best reflects real sensor behavior.

Key Experimental Results¶

Main Results¶

Accuracy of the photon-count distribution itself (averaged over a grid of \((S,B,t_d)\); ground truth is the empirical distribution from the sequential simulator). This is the paper's central table, since it directly tests the statistic prior work left unsolved.

Method	Wasserstein dist. ↓	KL div. ↓	Mean diff. ↓	Var. diff. ↓
Poisson	7966.130	24.208	7966.130	9257.464
Renewal (equal-energy constant-rate approx.)	207.565	18.821	207.564	45.927
Zhang et al. [39]	2.337	0.409	0.282	64.558
MaRS	0.309	0.013	0.269	1.588

Downstream depth and reflectivity recovery (same U-Net architecture, trained on data from each of the four simulators; "own" = tested on its own simulator's data, "GS" = tested on gold-standard data, which best reflects real sensor behavior).

Training simulator	Test domain	Reflectivity PSNR ↑	Reflectivity RMSE ↓	Depth PSNR ↑	Depth RMSE ↓
Poisson	own	36.4946	0.0722	33.6827	0.0606
Poisson	GS	7.0237	1.6537	7.2490	0.7853
Renewal	own	29.7230	0.0543	30.9704	0.0514
Renewal	GS	12.1210	0.0959	11.6988	0.0967
Zhang et al. [39]	own	22.7444	0.0526	21.9506	0.0506
Zhang et al. [39]	GS	22.5421	0.0820	22.2137	0.1193
MaRS	own	24.0321	0.0532	24.2821	0.0553
MaRS	GS	24.2488	0.0739	23.5776	0.0949

(The two column groups correspond to the original's Reflectivity and Depth blocks; the original table also reports an SBR=0.6 and an SBR=1.0 setting, while the cached text preserves only 10 numbers per row — ⚠️ for the full grouping and the SBR=1.0 values refer to Table 5 of the original.)

Ablation Study¶

Which statistic each simulator actually gets right (\(\mu\) count mean, \(\sigma^2\) count variance, \(\pi\) temporal PDF, use of a LUT). This table explains the origin of every performance difference above.

Simulator	Dead-time-correct \(\mu\)	Dead-time-correct \(\sigma^2\)	Dead-time-correct \(\pi\)	Uses LUT
Poisson	✗	✗	✗	✗
Renewal	✗	✗	✓	✓
Zhang et al. [39]	✓	✗	✓	✓
MaRS	✓	✓	✓	✓

Runtime breakdown for generating a 128×128×1024 cube, together with the LUT-versus-sequential comparison.

Stage / method	Time (s)
Batched lookup	0.020 ± 0.005
Temporal alignment (PDF shifted by \(\tau\))	0.144 ± 0.008
Histogram sampling	0.849 ± 0.070
MaRS total	1.013 ± 0.071
Sequential gold standard	40.831 ± 8.045

Key Findings¶

The missing count variance, not the mean, is Zhang et al.'s weak point: their mean difference (0.282) is essentially on par with MaRS (0.269), but the variance difference is 64.558 versus 1.588 — roughly 40× — showing that getting the mean right is far from enough to fix a histogram's radiometric scale.
Networks trained on a physically wrong simulator do not transfer at all: Poisson reaches 36.49 reflectivity PSNR on its own data but collapses to 7.02 on GS with a negative SSIM (−0.075); Renewal drops from 29.72 to 12.12. MaRS shows almost no cross-domain gap (reflectivity 24.03 → 24.25, depth 24.28 → 23.58) and achieves the best GS-domain depth and reflectivity recovery, turning "dead-time-correct photon counts" from a distributional improvement into a downstream-task gain.
In-domain scores alone are misleading: Poisson has the highest in-domain PSNR, but only because its own simulated data is artificially clean; the cross-simulator evaluation is what exposes the real gap.
The LUT is what makes it efficient: MaRS with lookup runs in near-constant time along all three axes (flux, realizations, pixels), the no-lookup variant degrades because of repeated large matrix operations, and the sequential simulator worsens sharply in every setting. The single-pixel picture is equally telling: Poisson 0.2 s, MaRS 0.3 s, gold standard 16.1 s.
The spectral truncation parameters are empirically grounded: across \((S,B,t_d)\) settings, eigenvalues beyond the fifth have magnitude \(\lesssim 0.7\), hence \(L=6\), \(p=5\); the control experiment shows that even 20 covariance terms fail to capture the long-range tail in difficult regimes, while the spectral rule recovers it at negligible cost. A scene's \((S,B,\tau)\) need not land on the LUT grid, but the resulting quantization error is minor and is absorbed by the reconstruction network.
Real-hardware validation goes beyond "matching another simulator": using a 7.5 MHz pulsed laser and a SPAD with 85 ns dead time, with detections over 9975 laser cycles forming one realization and 10000 realizations total, the authors estimate signal and background levels following Rapp et al.'s calibration procedure and predict the count distribution; MaRS matches the measured distributions for both a bright and a dark target at fixed distance and ambient illumination. The MRP count model thus captures real detector dead-time statistics rather than merely reproducing a sequential simulator.

Highlights & Insights¶

The "histogram = count × shape" decoupling is the paper's lever. Prior work (Rapp et al., Zhang et al.) devoted all its effort to the shape (temporal PDF) and did it well; this paper points out that the unresolved part is the scale (count), and that a wrong scale corrupts reflectivity estimation as a radiometric bias. Re-cutting the problem beats optimizing the old cut.
The MRP formulation is elegant: the Markov chain handles relative state, renewal theory handles absolute time, and one conditional holding-time distribution stitches them together; the observation that the count is the inverse of the cumulative holding time converts the hard \(N(t)\) into "CLT plus inversion." Classical renewal theory fails here (intervals are not i.i.d.) yet mixing is strong enough for the limit to stay Gaussian — that "one step back still works" judgment is the most beautiful part of the theory.
One physical symmetry buys three orders of magnitude in storage: delay invariance turns a 108 GB 3D table into a 100 MB 2D table plus a shift. Any simulator parameterized by a global shift (different ranges, different trigger instants) can reuse this recipe: precompute the shift-invariant part, defer the shift to sampling time.
A transferable template: for any state-dependent point process where an absolute count is required — dead-time-limited detector arrays, Poisson-regime optical communication, counting-type sensors beyond single-photon imaging — the chain "state chain + holding times + CLT inversion + low-rank spectral tail" applies directly.

Limitations & Future Work¶

The authors explicitly acknowledge that the LUT is not universal: it is bound to a given system configuration \((t_r,n_b,\sigma_t,t_d)\) and pulse shape, and must be regenerated for other regimes or other \((S,B)\) ranges and discretizations.
The tabulated flux range is limited to \(S,B\in[0.01,10]\) (roughly indoor levels); strongly reflective or long-range low-flux scenes would need a regenerated table or extrapolation, and extrapolation behavior is not reported.
The Gaussian approximation rests on the asymptotic CLT and thus needs \(t=N_r t_r\) to be large enough (high enough \(\Lambda\)). The experimental grid concentrates on the high-flux dead-time regime; whether the "Gaussian × multinomial" combination still holds for very few photons (small \(N(t)\)), and where its failure boundary lies, is not addressed.
All efficiency experiments run in NumPy on CPU (AMD Ryzen 7 9700X, 32 GB DDR5) with no GPU implementation, and the baseline is a same-machine sequential implementation; how the advantage holds against a heavily optimized GPU photon-level simulator remains open.
Real-hardware validation covers only two pixel distributions (bright and dark) at fixed distance and fixed ambient illumination, leaving distance sweeps, multiple targets, and extreme dynamic range untested.
Plausible improvements: make the LUT differentiable or adaptively interpolated to shrink the quantization error; model the breakdown regime of the Gaussian approximation explicitly and fall back to the exact distribution at low flux; replace the empirical \((L,p)\) constants with a choice adapted to the spectral gap of \(\tilde P\).

vs Poisson / multinomial simulators: they ignore dead time — fast (0.2 s) but physically wrong, with no pile-up distortion and uniformly spread background photons. MaRS delivers physically correct statistics at a comparable cost (0.3 s).
vs the sequential "gold-standard" simulator: per-photon rejection sampling is correct but does not scale (40.8 s for a 128×128×1024 cube, degrading with flux, pixels and realizations). MaRS is statistically indistinguishable from it while running >40× faster.
vs Zhang et al. [39] (neural-mapping simulator): an autoencoder learns the distorted temporal PDF, requiring retraining for each dead-time setting, count variance is visibly inaccurate (variance difference 64.6 versus 1.6), and the pipeline stays pixel-wise. MaRS is a unified closed-form model requiring no training and fully parallel.
vs the Markov chain models of Rapp et al. / Zhang et al. [41]: they provide only the temporal PDF, never the count. MaRS reuses that line as the "shape" component and supplies the missing "count," sharing the same transition matrix so the combination is consistent rather than bolted on.
vs classical renewal theory and self-exciting point processes: the former's i.i.d. interval assumption fails here and the latter is computationally intractable; the MRP carves an analytically tractable path between them.

Rating¶

Novelty: ⭐⭐⭐⭐ First analytical characterization of the photon-count distribution under dead time (MRP + CLT + spectral tail), genuinely breaking the simulator's speed–accuracy dilemma.
Experimental Thoroughness: ⭐⭐⭐⭐ Covers simulation accuracy, real SPAD hardware validation, downstream depth/reflectivity, and 3-axis efficiency scaling; hardware validation is small-scale and the efficiency comparison lacks a GPU baseline.
Writing Quality: ⭐⭐⭐⭐ The "count × shape" reframing is explained clearly and the physics-to-formula flow is smooth; the cost is high mathematical density, with most derivations deferred to the supplement.
Value: ⭐⭐⭐⭐ Provides scalable, physically faithful data generation for learned SP-LiDAR and plug-in completion for existing Markov-chain tooling.