Efficient Monte Carlo simulation of spatio-temporal speckles and their correlations

,

Despite much research on temporal or spatial statistics, there 30 is little usage of joint spatio-temporal statistics, in part due to 31 limited understanding of such statistics, and due to the lack of 32 Third, the motions of different scatterers are statistically indepen-106 dent, and their average displacement is an order of magnitude 107 smaller than the distance between scatterers. Finally, we ignore 108 refraction and reflection events at the interface of volume V.  Scatterers are moving independently 124 with displacements following a distribution T (∆ t ). We model it 125 as a Brownian motion with a drift, so that where w is sampled from a 3-dimensional unit normal distribu-127 tion, and the scalar t is the time interval of the displacement. The  tuations with a semi-random noise structure known as speckles 162 (see flatland speckles in Fig. 1(b,c)). A camera usually only 163 measures the intensity of these fields To capture speckle statistics, we can begin with the speckle 165 mean, Assuming scatterer density is stationary, the speckle mean is 167 time-invariant and we shorten notation and denote it as m i v . 168 We can similarly define higher-order statistics of speckles. Of 169 particular importance will be the speckle covariance, where (·) * denotes complex conjugation.
In this case, is present). Therefore, when characterizing speckle statistics, the 180 most challenging part is computing the covariance.

181
Note that our goal in this paper is to compute covariances between complex speckle fields. Such covariances can be easily translated to covariances of positive intensity images. Assuming the speckles follow a zero-mean multi-variate Gaussian distribution, the Siegert relation can be used: Computing speckle statistics.
A straightforward approach for computing the speckle covariance is to sample N different scatterer configurations O 1 (t), . . . O N (t), solve the wave equation for each configuration at each time instance , and then compute the empirical covariance: , note the strong correlation along a diagonal shifted from the center. This is the memory effect correlation, as speckles generated under nearby illuminations are correlated but shifted. (e) Temporal-only covariance for a fixed illumination i o and a fixed sensor point . At the same time instance the speckle fields are most similar, which is manifested by a strong diagonal in the correlation matrix. As the time difference increases, the speckles start to change. (f) Spatio-temporal covariance, where one wave is fixed and the other varies at both space and time: . The correlation is strongest at the central spot t = t o , v = v o where the two waves are identical, but it has a strong elongated diagonal corresponding to a non-zero constant velocity component U in the scatterer motion.

189
The basis for a Monte Carlo evaluation of speckle covariance  For ease of notation, w.l.o.g we consider covariances between 201 two time instances of the form t 1 = −t/2, t 2 = t/2. We denote 202 the mean particle position by We start by defining a mean path through the volume as an Our main result, proven in supplementary, is in showing that 211 the covariance can be expressed as the integral over the space of 212 such joint path pairs where each path has a contribution that can be expressed as a

214
Markovian product of its segments The contribution of central segments is identical to the contribution considered by classical Monte Carlo algorithms for evaluating temporal-only correlations [42]: and the formula for the start end/segments (the b = 1, b = B cases) is provided in supplement. Here ρ(·) is the phase function of the material denoting the probability of turning from one direction to the next one upon interacting with a scatterer, andα is the exponential attenuation along paths in the volume times its redial decay:α(o 1 , is the integral of phase variations over all possible motions that

224
The main difference between our derivation and the classi-225 cal temporal-only correlation is that we evaluate covariances 226 between two different source-sensor pairs, and therefore the con- The throughput is a complex number with a phase ξ propor-243 tional to the path length ℓ(⃗ x): where k = 2π/λ is the wavenumber at wavelength λ. 245 As a result, the covariance between two speckle fields can be 246 expressed as an integral over many pairs of paths: ⃗ x 1 from i 1 to 247 v 1 at time −t/2; and path⃗ x 2 from i 2 to v 2 at time t/2. Each path 248 has a contribution whose phase is proportional to the difference 249 between the two path lengths: considering joint pairs of paths that share their central segments where we denote the displacement sequence of the path as Note that in the above equation we use 260 the same central sequence, but at two different time instances.

261
Such joint path sequences are illustrated in Fig. 2

266
By considering the space sub-paths ⃗ x s and displacement 267 paths ⃗ ∆ t , we can evaluate the covariance as the integral of con-268 tributions from joint path pairs: Unfortunately, the path integral in Eq.

283
Our goal is to compute the spatio-temporal covariance of Eq. (11).

284
Given the absence of analytical solution, we want to approximate 285 the integral using importance sampling. For that we sample N 286 paths from a distribution q(⃗ x s ) of choice, and approximate the 287 covariance as While any distribution q(⃗ x s ) would provide an unbiased esti-289 mate in Eq. (21), the variance of this estimator can be largely 290 reduced using importance sampling rather than naive uniform 291 sampling. The estimation largely improves when we can define 292 a sampling strategy for which q is a good approximation of c ⃗ x s,n . 293 We review below two sampling strategies for the temporal-only 294 case, which are defined by different choices of q. We then adapt 295 them to the case of spatio-temporal covariance.
Previous approaches com-297 pute temporal correlation alone without spatial variation of 298 the source and/or sensor position. In this case i 1 = i 2 and 299 v 1 = v 2 , and the paths ⃗ x 1 ,⃗ x 2 in the derivation of the pre-300 vious section share all their segments. In particular, the for- 13)) is equivalent for all b values, including 302 the first and end segments. We also note that apart from the 3. If the last node on pathō B is on sensor v, update using 313 Eq. (20) Overall a path in this approach is sampled from a distribution The contribution of the path c ⃗ x s,n is non zero, only if it reaches 316 the area of the sensor v, for such cases we note that the ratio is only the γ(·) terms. Thus the estimate 318 can be expressed as  In a next-event estimation approach all paths contribute to the estimate, and moreover, every node on the path is explicitly connected to the sensor. We weight each connection by the probability that a path-segment leaving that node will actually hit the sensor. (c) Next-event estimation of path pairs: we only sample the central segment of the path, and connect it to two sources and two sensors. The covariance estimate sums the throughput contributed from such connections.
As a result, the covariance estimate can be expressed as Since we do not sample the last segment, the ratio leaves the term f I,n B in Eq. (28). Note that effectively, the for- node to the source as in Fig. 3(c), and add to the summation its 367 throughput f I 0 : The complete path sampling algorithm is provided in supple-

378
In this section our goal is to sample a speckle field with the cor-379 rect spatio-temporal statistics. For that we assume we are given section. That is, for every j, k, A straightforward way to do that would be to use the algo-  We follow the strategy of the previous section and sample 392 N subpaths ⃗ x s,n . For each subpath we sample a sequence of 393 temporal displacements ⃗ ∆ n t j from T . We define N × J paths by 394 concatenating the same sub-paths to all sources and sensors, as 395 illustrated in Fig. 4.
We define the sampled fields as the sum of contributions from 397 these paths. Each path has a phase proportional to its length and 398 we also need to take into account the attenuation and scattering Fig. 4. Sampling speckle fields: We sample subpaths ⃗ x s,n . To achieve a field with consistent entries, each subpath is connected to all source and sensors.
amplitude function at the first and last segments, as those are 400 not sampled by q. This leads to the fields A detailed pseudo-code for the field-sampling algorithm is pro-  We emphasize that the fields generated by the algorithm de-
In Fig. 5 we start by validating our ap-432 proach on 3D scenes using the publicly available MCX pack- where l t is the transport MFP, r 1 = δ 2 + l 2 t , and r 2 = 445 δ 2 + (2.33 · l t ) 2 (ignoring refraction at surface interface).

448
That is, we plot the correlation between a field uˆiv illuminated by  Fig. 6. Spatial-only memory effect correlation between the speckle fields generated from a sample under two directional illuminations, plotted as a function of the tilt ∆ θ between the illumination directions. We plot such correlations for media of two different thicknesses. We keep the mean free path fixed at 0.3 cm and vary the anisotropy parameter g of the phase function (so that while mean free path is fixed, the transport mean free path varies). For the thinner volume our correlation decay agrees with the theoretical C1 prediction at g = 0, but differs from it at larger g values as the diffusion limit is not yet achieved. For the thicker sample the diffusion approximation is also valid with g = 0.9. The simulation uses λ = 500 nm.
path varies). For the thinner volume our correlation decay agrees 459 with theoretical prediction at g = 0, but differs from it at larger 460 g values as the diffusion limit is not yet achieved. For the thicker 461 sample the diffusion approximation is also valid with g = 0.9.
In Fig. 7 we test the accuracy of 463 the spatio-temporal correlations computed by our algorithm.
The first row of each correlation image corresponds to the ∆ x = 0 476 case, which is the usual temporal-only correlation. As can be 477 seen, there are also other correlations in the data which depend 478 on the spatial positions.

479
In particular, note that if we have a purely linear motion so 480 that for all nodes on a path ∆ tb = tU, and we also select i 2 − i 1 =  Fig. 7(a) the illumination is fixed. In Fig. 7(b) we repeat the 489 experiment, this time when the illumination is also shifting in a 490 way that matches the scatterer velocity, With this tracking configuration the correlation along the diago-492 nal is even stronger. Finally, in Fig. 7(c) we repeat the simulation we also shift the illumination source to track particle motion, demonstrating an even higher correlation. This evaluates c(∆ . Finally in (c) we simulate correlations under directional illumination, which can detect strong correlations without moving the light source. We simulate three motion types illustrated in (d-f): Brownian only motion, mixture of Brownian and linear motion as well as a linear only motion. For fully linear motion, strong spatio-temporal correlations can be detected over long time instances. The simulation uses an homogeneous MFP= 40λ, D = 0.015λ 2 /s, U = 1λ/s, ∆ x ranging from 0 to 2λ, and t ranging from 0 to 2s. correlation of a tracking system, when the illumination and sen-552 sor point are shifting with the same velocity as the particles. To 553 see this, note that in a standard temporal-only correlation we 554 would measure: and in a tracking system we measure Views v The setup of our simulation is illustrated in Fig. 10(a)

t[µs] t[µs]
Brownian-Linear motion with tracking Brownian-Linear motion without tracking Isotropic Brownian motion Fig. 9. Separating motion components. We demonstrate the temporal decay in correlation for a static source and sensor, compared with the case where they are both shifting to track the linear part of the motion velocity. Without tracking the correlation decay mixes both a Brownian and a linear motion components. When tracking is used, the remaining correlation is only a function of the Brownian component. It precisely matches the correlation observed in a volume with the same D parameter and no linear term (U = 0). The simulation uses isotropic scattering, MFP = 0.1 cm, λ = 500 nm, a linear motion velocity of U x = 5 cm/s, and two different Brownian motions simulated in the two sub-figures.
where g 1 (t) is the normalized temporal covariance of the com-633 plex fields, which is expressed earlier in this paper as .

634
Intuitively, w is the dynamic portion of the volume, and when 635 this portion is larger the temporal correlation decays faster. In