Probing diffusive media through speckle differencing

Temporally varying speckle patterns, produced by light-matter interaction encode valuable information about inhomogeneities embedded within a scattering medium. These speckle fluctuations arise either from the tuning of the emission frequency of a laser illuminating a static scattering medium or from the microscopic motion of scatterers within a dynamically scattering medium. In this work, we detect embedded inhomogeneities by probing static and dynamic scattering media with coherent light and leveraging the statistical distribution of temporal speckle differences. In addition, we utilize the insights from the speckle differencing paradigm, to present the first experimental results of detecting inhomogeneities embedded within a scattering medium using bio-inspired neuromorphic sensors. The proposed neuromorphic approach simplifies the optical and electronic design, and significantly reduces data throughput by capturing only the differential information in the form of 1-bit spikes.


Introduction
Speckle patterns arising from the interaction of coherent light with a scattering medium encode a wealth of information about the objects embedded within them.For instance, speckle patterns can be used to recover images of objects embedded within scattering layers by leveraging angular [1,2] or spectral correlations [3] in scattered light.Temporal fluctuations in speckle produced by rough surfaces enable noninvasive detection of vital signals such as heartbeat or EEG [4].Methods like laser speckle contrast imaging (LSCI) utilize information within speckle patterns to visualize blood flow in biological tissue [5,6].In diffusely scattering media, speckle methods are employed to recover the scattering or absorption properties [7][8][9][10][11][12][13][14] or the statistical properties of the motion of scatterers [5,6,[15][16][17][18][19][20].The embedded inhomogeneities often have distinct scattering or absorption signatures, which serve as endogenous contrast mechanisms that facilitate their detection [21,22].Depending on the type of scattering medium, various acquisition and processing strategies are employed to extract the information encoded within speckle.
For instance, Diffuse Correlation Spectroscopy (DCS) methods use a continuous wave laser operating at a fixed wavelength to recover motion statistics of the scatterers within a dynamic scattering medium [15][16][17][18].In static scattering scenarios, speckle observed at distinct optical wavelengths is used to recover the transport properties of the medium [7][8][9][10][11][12][13][14].Furthermore, several techniques like LSCI [5,6] measure the local contrast in time-integrated speckle patterns to recover information about variations in dynamics or transport properties over a wide field of regard.Thus, various illumination, detection, and processing schemes are used to retrieve information about inhomogeneities from a scattering medium.
This work aims to develop a holistic approach that can detect inhomogeneities-embedded within both static and dynamic media over a wide field of view-by exclusively relying on endogenous contrast mechanisms.To this end, we rely on the statistics of the temporal gradients of speckle patterns observed at the distal end of scattering medium.Broadly speaking, the information encoded in the speckle patterns can be retrieved using two distinct constructs: one derived from signal processing concepts and the other random processes.A vast majority of existing methods employ signal processing constructs, wherein the temporal auto-correlation or cross-correlation of the speckle fluctuations in each pixel is computed.The shape and extent of these correlation curves are used to extract the properties of the medium.The primary disadvantage of the correlation methods arises from the need to capture substantial amounts of data that is uniformly sampled in time or optical frequency.Speckle contrast methods, which utilize random process constructs to detect inhomogeneities, operate by calculating the ratio of the standard deviation of intensity to the mean intensity across a spatial neighborhood of pixels (typically 7 × 7 pixels, and in some cases 50 × 50 pixels to enhance signal-to-noise ratio [44).Although these methods are effective in detecting inhomogeneities, they do so at the expense of spatial detail.
This work employs random processes constructs (detailed in Section 3) to detect inhomogeneities through speckle differences.In the proposed speckle differencing scheme, illustrated in Fig. 1, the static or dynamic scattering media is probed using coherent light to produce a time varying speckle data cube.The temporal or spectral gradients at each pixel of the acquired speckle data cube are computed.Since the rate of speckle fluctuations is modulated by the embedded inhomogeneities, the distribution of speckle differences varies across different spatial regions.Specifically, the population of gradients with a particular amplitude will differ by spatial region.A thresholding procedure is applied to enhance these differences and generate an inhomogeneity map.
Furthermore, this study introduces a bio-inspired neuromorphic sensor-based detection scheme to identify inhomogeneities within a scattering medium.The necessity for a bio-inspired detection scheme arises from the need to detect minor changes in speckle intensity at high rates.This demands an imaging system with high dynamic range, bit-depth, and frame rate.Traditional imaging systems, with their limited space-time-bandwidth product, are not optimized for such high data transfer rates.Biological vision systems address this elegantly by responding exclusively to change.Neuromorphic sensors mimic this behavior by responding exclusively to changes in the light distribution incident at each pixel.Additionally, the log-companding circuitry in these cameras ensures they respond to changes in temporal contrast (∆I/I) of the incident light, similar to human vision.These sensors generate a 1-bit spike whenever the temporal contrast (∆I/I) exceeds a preset threshold (see Section 4.3).The contrast stretching characteristics of the logarithm enable these sensors to operate in large dynamic range scenes (>120 dB compared to 60 dB with traditional cameras).
The high speed (≈ 1 µs resolution), large dynamic range, low data throughputs, and high pixel density (780 × 1920) of the neuromorphic sensors is furthering their adoption in several fields, including computer vision [23] and microscopy [24].Typically, these sensors are used in scenes where the temporal changes are confined to sparse spatial locations.Initially, it may appear challenging to use these sensors for monitoring spatially dense speckle patterns.However, the thresholding functionality intrinsic to neuromorphic sensors effectively reduces the number of pixels that generate an event at any given time.In the case of the 100ℓ * medium shown in Fig. 1, only 30% of pixels exceed a spectral gradient amplitude (used as a threshold here) value that is half the maximum value.Therefore, the application of these sensors to speckle patterns aligns with their conventional operation.
In this work, we utilize neuromorphic sensors to monitor the speckle pattern emerging from a diffusive medium.The number of events generated in a spatial neighborhood provides an estimate of the fluctuation rate of speckle patterns, which is related to the local transport properties of the diffusive medium.Thus, the number of events generated at each pixel serves as a proxy for the transport properties, enabling the detection of inhomogeneities.This neuromorphic strategy simplifies the acquisition and processing complexity by eliminating the need to capture hundreds Fig. 1.Probing diffusive media through speckle.Time varying speckle patterns are produced from a scattering medium either by tuning a laser (static media) or scatterer motion (dynamic media).The presence of inhomogeneities produces a change in statistics of speckle fluctuations.In static scattering media, speckle fluctuations are produced by tuning the illumination wavelength, whereas in dynamic media, the motion of scatterers produces these fluctuations.A) Illustration of proposed speckle differencing approach.The gradients of speckle exhibit variations based on transport properties.This can be leveraged to detect inhomogeneities.Simulation details are provided in Section 3 B) A neuromorphic sensor directly records the gradient information in the form of spike, enabling the detection of inhomogeneities from event activity. of 8 − bit or 12 − bit images.Instead, it directly detects inhomogeneous regions using a 1 − bit spike or event information produced at each pixel.In Section 4, we present, to our knowledge, the first experimental results of using a neuromorphic sensor to probe diffusely scattering media.
The following are the key contributions of this work: • We introduce a speckle differencing scheme, that utilizes constructs from random processes to detect static or dynamic inhomogeneities within a diffusive medium.We demonstrate the ability of the proposed approach to operate in both static and dynamic media, allowing it to detect changes in the transport properties and dynamics.
• Using the insights gleaned from the speckle differencing paradigm, we present, to our knowledge, the first experimental results of probing inhomogeneities embedded within a scattering medium using a bio-inspired neuromorphic sensor.The proposed approach does not require complex optical or electronic design and significantly reduces data throughput by exclusively recording the gradient information in the form of 1 − bit spikes.
The remainder of this manuscript is organized as follows: Section 2 provides a brief overview of the statistics of speckle patterns emerging from a scattering medium.Section 3 provides analytic expressions describing the proposed speckle differencing approach.Section 4 details the experiments conducted to validate the proposed methods, and Section 5 provides a discussion of these approaches.

Background
This Section provides background information on the statistical properties of speckle patterns emerging from a scattering medium.The expressions presented here and in the following sections utilize the notation and symbols specified in Table 1.
In a volumetric scattering medium, the incident light is scattered multiple times such that it travels in complex, convoluted paths that resemble a random walk.Under coherent illumination, such as the output from a narrowband laser, these randomly phased light rays interfere to produce a stochastic pattern, called a speckle pattern at the output of the medium [25].The intensity, I of a fully developed speckle pattern follows a negative exponential probability distribution [25].In several applications, the limitations imposed by the spatiotemporal sampling density of the detectors result in the recording of a blurred and integrated speckle pattern.The integrated speckle can be modeled as a sum of N independent speckles, each following a negative exponential probability distribution.Consequently, the probability density function of the detected speckle pattern becomes a Gamma distribution [26], Here, ⟨I⟩ is the mean intensity of the integrated speckle pattern, and Γ(N) is the Gamma function evaluated at N. The variance (σ 2 ) of the distribution is given as ⟨I⟩ 2 N −1 .The speckle contrast (σ/µ) becomes 1/ √ N. Since the superposition of time-delayed wavefronts produces speckle patterns, they also encode information about the time-of-flight distribution of light emerging from a medium [7][8][9][10][11][12][13][14].This information is typically extracted by varying the illumination frequency of the laser.As the illumination frequency varies, the phase relationship between the interfering wavefronts changes, causing the speckle to fluctuate.Specifically, the field-field correlation between the speckle fields E( ν, r s , r d ) and E( ν + ∆ν, r s , r d ) observed at optical frequencies ν Hz and ν + ∆ν Hz, respectively, is proportional to the Fourier Transform of the temporal response, p(t; r s , r d ) of the scattering medium [7,8], Here p(t; r s , r d ) represents the temporal response or the time-of-flight distribution of light emerging from the scattering medium for an ultra-short pulse source positioned at r s : (x s , y s , z s ) and a detector positioned at r d : (x d , y d , z d ).P(∆ν; r s , r d ) represents the complex-valued frequency response or the Fourier transform of p(t; r s , r d ), evaluated at ∆ν.A ⟨I⟩ denotes the average intensity over an ensemble of measurements.The temporal response p(t; r s , r d ) is also called the Temporal Point Spread Function (TPSF) of the scattering medium.
Since speckle intensity measurements are more straightforward to acquire experimentally, the normalized intensity correlation function, also called the Spectral Speckle Intensity Cumulant Correlation function (SICCF), C I (∆ν; r s , r d ) is typically computed to characterize the scattering media.Typically, C I (∆ν; r s , r d ) is defined as [8], Here, Ĩ = (I− ⟨I ⟩) where ⟨I⟩ and I σ represent the mean and the standard deviation of the observed speckle intensity.Consequently, speckle patterns acquired at a series of closely spaced optical frequencies provide indirect access to the TPSF of the scattering medium.This approach bypasses the need for ultra-short pulse sources and time-resolved detectors to measure p(t; r s , r d ).Moreover, the time-of-flight distribution, p(t; r s , r d ) can be parameterized as a function of the transport properties: the reduced scattering coefficient, µ ′ s and the absorption coefficient, µ a of the scattering medium [27].These insights were used to develop methods that estimate the transport properties of scattering media [9][10][11][12][13][14].In these methods, a series of speckle intensity images are captured at closely spaced wavelengths.These speckle images were used to compute the SICCF disclosed in Eq. ( 3).The transport properties of the medium are then estimated by computing the TPSF from the power spectrum using bi-spectrum techniques [9][10][11][12][13][14].Here, the speckle intensity images are acquired either by stepping through discrete optical frequencies or by continuously tuning the emission frequency of the laser.The latter approach, adopted in this work (Section 1 in Supplement 1), results in a temporally varying speckle pattern that can be monitored using a framing or neuromorphic camera.
In dynamic scattering media, such as colloidal suspensions, the motion of scatterers changes the relative phases of the scattered wavefronts, causing the speckle intensity to fluctuate [15].This is typically encountered in applications such as cerebral blood flow imaging, where the blood cells flowing through vessels deep inside tissue cause the speckle to fluctuate.Here, the normalized temporal field correlation of speckle, provides access to the mean squared displacement of the scattering particles, ⟨∆r 2 (τ)⟩ in time τ [15,16].In the diffusive regime, the normalized field correlation inside the medium is related to the field correlation at the boundary through the correlation diffusion equation.The expression for normalized field correlation of the speckle fields for an infinite extent, homogeneous medium is given as [15,16], where k 0 represents the wave number, |r s − r d | represents the source-detector distance, µ ′ s and µ a represent the reduced scattering and absorption coefficient, respectively.The normalized field correlation function, g 1 (τ), can be related to the normalized intensity correlation function, , using the relation [15,16], From Eq. ( 4) and Eq. ( 5), it is evident that the changes in the mean squared displacement of the scatterers, ⟨∆r 2 (τ)⟩ or the transport properties, or a combination thereof, lead to local variations in the temporal evolution of speckle patterns.In techniques such as Laser Speckle Contrast Imaging (LSCI), the local speckle contrast of a time-integrated speckle pattern is used to identify regions with moving scatterers, such as blood cells [5,6].Techniques like DCS directly obtain the intensity autocorrelation profile by sampling the temporally evolving speckle patterns using a high-speed detector [15][16][17][18].This profile is then fitted to analytical expressions to estimate parameters such as ⟨∆r 2 (τ)⟩, Brownian diffusion coefficient (D B ), or physiologically relevant properties like the blood flow index (BFi) [18].

Statistics of speckle differences
From the discussion in the previous Section, it is evident that the fluctuations in speckle patterns encode valuable information about the properties of the scattering medium.The fluctuations induced by changing the emission frequency of the laser (ref.Equation ( 3)) provide insights into the transport properties of the medium.In contrast, the speckle fluctuations observed in a dynamic scattering medium reveal information about the motion of scatterers deep within the scattering medium (ref.Equation ( 5)).Therefore, any spatial non-uniformity in the speckle fluctuations observed at the boundary of the medium indicates the presence of inhomogeneities within the medium.
To mathematically describe the proposed approach, consider the two discrete random variables, X = I(ν 0 ) and Y = I(ν 0 + ∆ν), where I represents the intensity at a given pixel.Since X and Y are identically distributed (representing the speckle differences observed at the exact location), the mean value of the difference between these two random variables, Z = X − Y is given as The variance of the random variable Z can further be obtained as, where ρ XY denotes the cross-correlation between the two random variables, which is given by g 2 (τ; r s , r d ) for dynamic media and C I (∆ν; r s , r d ) for static media.Since X and Y are identically distributed, σ 2 X = σ 2 Y = σ 2 .Thus, the variance of Z can be simplified to be, It is evident from Eq. ( 8) that the variance and standard deviation of the distribution of the differences are related to the cross-correlation of the two random variables, which is given by g 2 (τ; r s , r d ) for dynamic media and C I (∆ν; r s , r d ) for static media.Additionally, an upper bound on the probabilities can be determined using Chebyshev's Inequality, which is stated as follows [28]: i.e., the probability that |Z − µ z | exceeds kσ z is upper bounded by k −2 .The above expression can be simplified by substituting the values of µ z and σ z , such that, , the above expression can be rewritten as, From Eq. ( 11), the probability of |Z| exceeding nσ is upper bounded by 2n −2 (1 − ρ XY ).An increase in the correlation ρ XY reduces this bound, whereas a reduction in the correlation increases it.Any spatial variations in ρ XY i.e., g 2 (τ; r s , r d ) for dynamic media and C I (∆ν; r s , r d ) for static media, are thus reflected in the probability distribution of Z = X − Y. Therefore, nσ can be used as a threshold to emphasize the differences between the inhomogeneous regions and the background medium, thereby enhancing the visibility of the inhomogeneities.Section 3 in the Supplement 1 delves into the impact of noise terms on the above formalism.
It is worth noting that the proposed approach is inspired by the use of similar methods in fields such as food inspection and metrology [29][30][31][32].In metrology applications [29,30], zero crossings in speckle patterns have been used as a proxy for measuring velocity and related parameters.Similarly, in food inspection applications [31,32], generalized differences or weighted generalized differences of speckle patterns have been used to monitor activity over specific regions in fruit.The threshold can be viewed as a weighting function, similar to that used in the Weighted General Difference in bio-speckle applications, to emphasize regions based on necessity [31].
A series of simulations were performed to evaluate the effectiveness of the distribution of speckle differences in revealing variations in static (µ a or µ ′ S ) or dynamic properties (D B ) of the medium.These simulations are described in detail in the subsequent sections.

Simulating the spectrally evolving speckle pattern
The approach described here simulates the speckle field observed at a series of optical wavelengths by leveraging the relationship between the TPSF of the medium and the spectral evolution of speckle, as outlined in Eq. ( 3).For this simulation, a three-dimensional, complex-valued speckle field z was generated, with its real and imaginary parts independently and identically distributed as delta-correlated Gaussian random variables.This ensures that the simulated speckle field exhibits complex circular Gaussian statistics.
The spatio-temporal speckle field thus generated was multiplied by the TPSF of the scattering medium, p(t; r D ).Here r D : [x, y, L] T denotes the detector coordinates at the exit fact of the medium, with L being the thickness of the medium.The source, r D : [x, y, 0] T is assumed to be confocal with the detector.In the simulations shown in Fig. 2, the TPSF p(t; r D ) is obtained from the analytic expression for the TPSF of an infinite extent, homogeneous medium with transport properties specified subsequently [27].The spectral speckle field, E ′ (x, y, ν; L) is then estimated by computing the Fourier transform of the scaled three-dimensional stack, i.e., In addition to a finite spectral correlation, the speckle fields observed in experiments also exhibit a finite correlation along the spatial dimension due to the blurring introduced by the aperture of the imaging optics.To simulate this effect, the spatial Fourier transform of E ′ (x, y, ν; L), i.e., F (E ′ (x, y, ν; L)) is computed.This is then multiplied by a window function w uv (u, v, ν), where u and v represent the spatial frequency in the Fourier domain.The window function w uv (u, v, ν) multiplies each (u, v) slice in F (E ′ (x, y, ν; L)) by a circular aperture function.The spatially and spectrally blurred speckle field speckle field E(x, y, ν) is then obtained by computing an inverse Fourier transform along the spatial dimensions, The above procedure was used to simulate the spectral speckle in steps of 25 MHz for 250 GHz, producing a three-dimensional stack of size 51 × 51 × 10000.The speckle patterns were generated for an infinite extent, homogeneous scattering medium with a µ ′ s = 6.7 cm −1 and µ a = 0.02 cm −1 , mimicking the optical properties of gray matter in the human brain [33].The speckle patterns were simulated for two different transport mean free paths (Lµ ′ s = 12 and 100) by varying the thickness of the medium L from 1.8 cm to 15 cm.A speed of light (c) of 2 × 10 8 ms −1 in used in these simulations.
Since the field-field correlation, G E (∆ν; r s , r d ) and p(t; r s , r d ) form a Fourier transform pair (refer to Eq. ( 2)), a wider time-of-flight distribution (more scattering) must result in a speckle pattern that decorrelates faster.More specifically, the speckle pattern observed within a spectral window must fluctuate more rapidly with the increasing severity of scattering.This effect can be seen in Fig. 2(A), which shows a simulated speckle pattern over a frequency range of 10 GHz for the scattering media with increasing Lµ ′ s .Furthermore, the histograms of the spectral derivatives of the simulated speckle accumulated over a 5 × 5 neighborhood of pixels (representing the coherence area of the speckle pattern) for each of these scattering scenarios are displayed in Fig. 2(B).These plots indicate the ability of the spectral gradients of speckle patterns to discern variations in the transport properties of the medium.

Simulating temporal speckle
The procedure adopted for simulating spectrally evolving speckle is inapplicable for simulating dynamic speckle seen in applications like Laser Speckle Contrast Imaging (LSCI) or Diffuse Correlation Spectroscopy (DCS), as only the autocorrelation profile of the speckle is known in these applications.Various methods have been proposed in the literature to address this.The approach proposed in [33] is used here for generating spatio-temporal speckle patterns with the  8)) desired temporal correlation.This method employs two independent and identically distributed speckle fields W ′ and W ′′ to generate the desired speckle field W n using the expression below: where ρ n denotes the correlation coefficient between the W n and W ′ .For ρ n = 1 the speckle fields are fully correlated, whereas for ρ n = 0, they are uncorrelated.The three-dimensional speckle stack with a pre-defined temporal correlation can thus be generated by varying the value of ρ n in the expression disclosed in Eq. ( 4).The speckle field W ′ is generated using the procedure disclosed in Section 3.1.The resulting speckle field W ′ δ (x, y) in the spatial domain exhibits circular Gaussian statistics and is δ-correlated.To introduce a finite correlation width along the spatial dimensions caused by the finite aperture size of the camera, the complex-valued field, F (W ′ δ (x, y)) is multiplied by a circular aperture function w uv (u, v).The speckle field W ′ (x, y) is then obtained by computing an inverse Fourier transform.
The same procedure is repeated with a different realization of random numbers to generate W ′′ (x, y), which is then used to generate the speckle patterns with the desired temporal correlation using Eq. ( 4).This method was employed to simulate the spectral speckle in steps of 1 µs increments over a span of 0.3 ms, producing a three-dimensional stack of size 51 × 51 × 3000.The speckle patterns were generated for an infinite extent, homogeneous scattering medium with a µ ′ s = 8 cm −1 , µ a = 0.05 cm −1 , and a D B = 5 × 10 −8 cm 2 /s.For the second scenario, the µ ′ s = 4 cm −1 , µ a = 0.05 cm −1 and D B = 2.5 × 10 −8 cm 2 /s were used.A thickness of 1 cm for the medium was assumed in these simulations.From the simulations shown in Fig. 3(A), it can be seen that the speckle pattern fluctuates more rapidly with increasing Brownian diffusion coefficient (D B ).The histogram plots shown in Fig. 3(B) indicate the ability of the gradients of speckle patterns to discern variations in the dynamic properties of the medium.

Experiments
The schematic setup shown in Fig. 4(A) is used to validate the proposed speckle differencing approach experimentally.In this setup, a narrowband tunable Distributed Bragg Reflector (DBR) laser (Part number: PH780DBR040BF-ISO, linewidth <1MHz) operating at 780 nm is used to illuminate the scattering medium.As shown in Fig. 1, the emission wavelength of the laser is tuned to produce a time-varying speckle pattern at the distal end of a static scattering medium.For dynamic media, the wavelength of the source is held constant, and the motion of scatterers within the medium is used to produce a time-varying speckle pattern.This time-varying speckle pattern is relayed to a Focal Plane Array (FPA) using an imaging lens.The FPA is used to record the temporal variations in these speckle patterns.
The experiments presented here utilize two types of FPAs: a traditional CMOS FPA and a neuromorphic FPA.The CMOS FPA records the speckle patterns at different time instances to produce a spatio-temporal speckle stack.This spatio-temporal is then processed using the procedure shown in Fig. 4(B) to recover the spatial maps of the embedded inhomogeneities.The neuromorphic sensor, however, directly records the gradient information by generating a sequence of events in response to the incident time varying intensity.The subsequent sections provide a detailed description of these experiments.

Static medium
The 3d printed scattering medium, shown in Fig. 4, was used to experimentally validate the ability of the proposed approach to detect inhomogeneities within a static medium.The scattering medium was fabricated using a mixture of white pigment (RS-F2-CRWH-01) and Color Base resin (RS-F2-GPCB-01), which is commercially available from Formlabs.The standard recipe and mixing protocol prescribed by Formlabs were used to synthesize the white resin and print the targets.The reduced scattering coefficient (µ ′ s ) of the 3d printed medium was estimated to be 15 cm −1 using oblique incidence reflectometry [34].The absorption coefficient: µ a = 0.01 cm −1 and the refractive index: n = 1.5 were obtained from the literature [35].The scattering medium has the dimensions of 165 mm × 150 mm × 40 mm, giving an overall scattering thickness of 60ℓ * .Various shapes were cut out of 60 lb white matte paper and placed midway within a thick scattering medium to generate scattering inhomogeneities.The light output from the tunable DBR laser was amplified using a tapered amplifier (Toptica BoosTA Pro, 780 nm) to produce an output power of 180 mW.The output from the tapered amplifier was expanded to a diameter of 150 mm to illuminate the scattering medium.The speckle distribution emerging from the medium in a transillumination geometry was relayed to a focal plane array (IDS UI-3000SE) using an imaging optic (Sigma 35 mm f/1.4 DG HSM Art lens).
The spectral correlation techniques for characterizing scattering media typically acquire speckle measurements by stepping through discrete optical frequencies [12].For each speckle image, the emission frequency of the laser is either logged or actively controlled using a reference Fig. 5. A) Results from imaging through a 60ℓ * thick static scattering medium using the proposed approach.B) Images acquired with a conventional wide field imager using an 850 nm LED source.cavity, such as a scanning Fabry-Perot interferometer.However, very small frequency increments, on the order of 10 MHz used in this work, necessitate the use of a highly sensitive reference cavity or wavemeter, with high stability throughout the experiment.Since this is challenging to realize, a computational procedure described in Section 1 in Supplement 1 was used to thermally tune the emission frequency of the laser at a rate of 0.4 GHz/s.The sensor captured 100 frames of size 1992 × 2484 at an exposure time of 50 ms.Thus, these images were captured at 20 MHz increments.The corresponding speckle evolution over a small temporal window of size 51 × 51 × 100 is shown in Fig. 4(A).
The processing workflow depicted in Fig. 4(B) is used to assemble the spatial map of the inhomogeneities embedded within the medium.The acquired images are independently rescaled to be in range [0, 1], following which the temporal gradients of the speckle intensity stack are computed.From Sections 2 and 3, it is evident that the distribution of the temporal gradients exhibits differences that can be exploited to assemble the spatial map of the inhomogeneity.To exaggerate these differences, we adopt a thresholding procedure, where only the gradients above a certain threshold value θ are retained.The thresholded gradients are accumulated at each pixel to produce an inhomogeneity map.The above procedure exploits the fact that the population of gradients within the threshold value will differ for different regions (see the histogram plot in Fig. 4(B)).In the extreme case, the gradients corresponding exclusively to either the background or the inhomogeneity can be retained to produce a map with high contrast.
Here, the standard deviation, σ b of the gradients computed over the entire data stack is used as a proxy for the spread of gradients of the background medium.In this work, the gradients exceeding a threshold θ>4σ b were retained after the thresholding procedure.Here, it is assumed that the pixels containing the inhomogeneity are significantly smaller than the background medium.If this condition is violated, an alternate strategy must be adopted to identify the threshold.In the datasets acquired from static scattering medium, the gradient magnitude is in range [0, 0.16], and θ ≈ 0.035.The thresholded gradients were accumulated and denoised using a 3 × 3 median filter to generate the wide-field images of inhomogeneities shown in Fig. 5(A).In Fig. 5(B), the images of these scattering media in a transillumination geometry using an 850 nm LED source is provided.From these results, it is evident that the proposed approach significantly enhances the detectability of embedded objects.

Dynamic medium
A compressed breast liquid phantom, shown in Fig. 6(A) was used to experimentally validate the ability of the proposed approach to detect inhomogeneities within a dynamic scattering medium.The hollow compressed breast phantom mold was modeled as a scaled replica of the DigiBreast phantom [36].This model was 3D printed with clear resin (RS-F2-GPCL-04) on a 3D printer (Formlabs Form 3).The mold had the dimensions of 150 mm × 90 mm × 30 mm was filled with a dilute intralipid solution (Intralipid 20, Sigma Aldrich).The solution was created by adding 2.5 % intralipid (v/v) to 200 ml of water to create a phantom with a µ ′ s = 5 cm −1 and µ a = 0.03 cm −1 at 780 nm [37].The inhomogeneity was created by creating another intralipid solution with µ ′ s = 2.5 cm −1 into a cuvette of dimensions 10 mm × 10 mm × 35 mm.The cuvette was suspended in the middle of the liquid phantom i.e., half way between the source and detector planes.The Brownian diffusion coefficient (D B ) varies linearly with the concentration of the intralipid [38].Moreover, the speckle decorrelation time depends on the transport properties and the Brownian diffusion coefficient, as disclosed in Eq. ( 4), causing a variability in the temporal evolution statistics of the speckle.The experimental setup was similar to the setup used for the static scattering scenario.The liquid phantom was illuminated using a laser beam of 1 cm diameter in a transillumination geometry.The speckle patterns at the distal end of the medium were relayed to a focal plane array (IDS UI-3370) using an imaging optic (Thorlabs MVL50TM23).The sensor was used to capture 1500 frames of size 528 × 528 with an exposure time of 3 ms, ≈8× Gain (analog gain of 2 and digital gain of 4) and frame rate of 140 fps.Due to a large D B associated with the intralipid, the speckle fluctuations in this experiment far exceeds the frame rate of the sensor.
However, the integrated speckle patterns exhibited residual correlations sufficient to identify the inhomogeneous regions within the sample using the proposed approach.The recorded image stack was processed using the workflow outlined in Fig. 4(A).The standard deviation of the temporal gradients computed over the entire stack was used as an aggregate measure representing the fluctuations observed in the background medium.This assumption is valid so long as the inhomogeneous regions are sparse.The gradient values exceeding this threshold were accumulated at each pixel to assemble the spatial map of the inhomogeneity.In the datasets acquired from static scattering medium, the gradient magnitude is in range [0, 0.18], and θ ≈ 0.06. Figure 6(B) depicts the results for a single inhomogeneity that was suspended at a lateral offset of 5 mm from the illumination beam center.Whereas Fig. 6(C) shows the results for two scattering inhomogeneities separated by ≈1 cm.The overlaid regions displayed in these images were normalized to be in range [0, 1].
In the dynamic media used here, the change in intralipid concentration produces a change in both the dynamic properties (D B ) and transport properties (i.e., µ ′ s and µ a ) of the medium.This is also generally the case with most biological media.However, an experiment to examine the ability of the proposed approach to detect changes exclusively arising due to the dynamic properties of the medium is also performed.This is described in Section 2 of the Supplement 1.Additionally, the dynamic scattering experiments showcased in this work were limited to operating at ≈15ℓ * thick media and covering a ≈1cm 2 area.This is due to the use of low frame rate of the camera used in this work and not a limitation of the approach.

Neuromorphic sensing
The experimental evidence presented in the previous sections demonstrates the capability of the speckle-differencing approach introduced in this work for wide-field probing of scattering media.However, it is important to note that this approach necessitated capturing over 100 multi-megapixel, 12 − bit quantized image frames, of which a significant portion was discarded during the thresholding procedure.This raises a critical question: "Is there a more efficient method to retrieve the same information?"Indeed, there is a more efficient method.By monitoring temporal differences or gradients in speckle intensity, the speckle-differencing approach shares similarities with the functioning of neuromorphic sensors.
These sensors are designed to mimic the mechanisms employed by the human visual system to process visual information.Unlike pixels in a traditional FPA that synchronously record intensity over an array of pixels at regular intervals, neuromorphic pixels respond exclusively to changes in intensity.They asynchronously generate 1 − bit spikes or events only when the intensity at a pixel exceeds a threshold.More specifically, each neuromorphic pixel independently generates a "spike" when the natural log intensity (∆ logI) changes by a value that exceeds a preset threshold (Θ on or Θ off ), as shown in Fig. 7(A).
Since these sensors exclusively record change information, they generate significantly less data.Moreover, the logarithm operation allows the neuromorphic sensor to exaggerate differences, in a manner similar to contrast stretching in image processing.However, here, the digitization happens after the logarithm operation, thereby allowing these sensors to operate at >120 dB dynamic range compared to the ≈60 dB dynamic range of traditional image sensors [23].These unique features have been leveraged to enhance the performance of imaging techniques across varied field, including computer vision [23], microscopy [24], velocimetry [39], and x-ray imaging [40].Additionally, these cameras are increasingly used for applications involving speckle patterns [41][42][43] due to their high dynamic range, necessitating highly sensitive detectors.
The simulation, shown in Fig. 7(A), was performed to examine the use of neuromorphic sensors for detecting inhomogeneous regions within scattering media.In this simulation, the speckle patterns at a pixel in the transillumination geometry were generated for two different transport mean free paths (Lµ ′ s = 12 and 100) using the procedure outlined in Section 3.1.The intensity was scaled by factor of 10, followed by an addition of a bias value of 10, such that the speckle intensity is in range [10,100].The response of a neuromorphic pixel was simulated by generating an event whenever the change in the logarithm of speckle intensity exceeded θ = 10.It can be seen from Fig. 7(A) that the number of events generated by the neuromorphic pixel depends on the scattering thickness of the medium.Thus, the inhomogeneity can be detected by pointing an event sensor at the scattering medium and accumulating the spikes generated in response to time-varying intensity, significantly simplifying the acquisition and processing complexity.This insight forms the basis of the approach presented here.

A. Response of a neuromorphic pixel to incident speckle intensity: Simulation
The approach is experimentally validated by replacing a traditional camera with a neuromorphic camera in Fig. 4(A).Here, the neuromorphic camera was used to image a scattering inhomogeneity, shown in the inset, embedded inside the scattering medium used in Section 4.1.The laser beam was expanded to a diameter of 150 mm to illuminate the scattering medium.The emission frequency of the laser was tuned at a rate of 4 GHz/s over a range of 400 GHz.The speckle patterns were monitored in a transillumination geometry using a neuromorphic FPA (Prophesee EVK4).The embedded objects of various shapes were cut out in white matte paper and placed halfway inside a 60ℓ * thick scattering medium.These cutouts represent a scattering or absorbing inhomogeneity in an otherwise homogeneous medium.The positive and negative events triggered by the neuromorphic camera (θ on = 100, θ off = 100) during the sweep were accumulated to produce an event map.The accumulated event map was then denoised using a 3 × 3 median filter to produce the wide field activity maps shown in Fig. 3-8(B).These results demonstrate the fact that neuromorphic sensors can indeed be used to detect inhomogeneous regions within scattering media.
The accumulated event map can be physically interpreted as being related to the spectral correlation of speckle pattern at each pixel, as described by the Chebyshev's inequality in Eq. ( 12), Thus, insights into the factors influencing the detectability of the inhomogeneity can be gleaned by examining Eq. ( 8), which upper bounds the probability of |X − Y |>nσ to be 2n −2 (1 − ρ XY ).In the current case, X = Ĩ(ν 0 ), Y = Ĩ(ν 0 + ∆ν) and ρ XY = C I (∆ν).Here, the spatial dependence of these variables is omitted.Thus, the contrast accumulated events at any pixel is proportional to (1 − C I (∆ν) )/n 2 .Here, the sweep rate of the factor C I (∆ν) is determined by the sweep rate of the laser, whereas n is determined by the threshold parameter of the event sensor.Thus, these two factors influence the detectability or the spatial contrast within the topographic maps obtained from the abovementioned procedure.

Discussion
The increasing trend towards the use of coherent sensing methods for probing scattering media has resulted in the development of a suite of methods to retrieve information encoded within speckle patterns.For instance, methods like SCOS [19], that shares similarities with LSCI, allows the retrieval of dynamics of scatterers by measuring the speckle contrast of intensity images acquired at various integration times.Alternatively, methods like Fourier domain DCS [20] utilizes a heterodyne interferometric arrangement to measure the dynamic properties from frequency domain measurements.Another set of methods leverage the advancements in detection technology, such as the availability of megapixel SPAD arrays, to develop parallelized detection schemes with a significantly large number of detection channels [45,46].For example, in [46], a 0.25 MPix SPAD pixel array operating at 100 kfps was used to achieve an order of magnitude improvement in SNR associated with DCS measurements.However, further improvements can be increasingly hard to realize due to the challenges emerging from high data transfer rates.
This work introduces the speckle differencing paradigm and combines the insights from this paradigm with a bio-inspired neuromorphic sensor to develop an approach that detects inhomogeneities.The designed approach does not require complex optical or electronic design, preserves spatial detail.Moreover, the economic data representation in the form of 1 − bit spikes, significantly reduces data throughput while retaining pixel density.It is hoped that the methods and analyses presented in this work can further the adoption of these neuromorphic devices in biomedical applications.Additionally, the experiments presented in Section 4 demonstrate the ability of the speckle differencing paradigm to operate in both static and dynamic media, enabling the detection of changes in transport properties and dynamics.This method provides a holistic way to detect inhomogeneities within static and dynamic scattering media by utilizing temporal gradients of speckle.While state-of-the-art speckle-based methods can also be adapted to detect variations in transport properties, to our knowledge, this is the first demonstration of such a capability.
However, the proposed methods are not without limitations.The proposed approach in its current form cannot provide quantitative information about variations in transport properties or dynamics.However, it is important to note that the accumulated gradients are related to g 2 (τ) or C I (∆ν) (see Eq. ( 11)).Therefore, it may be possible to extend the method to infer quantitative details from the thresholded gradients or event activity.One potential solution is to train a spiking neural network to extract this information from event activity [47].It is also worth noting that the neuromorphic sensor-based approach is not single photon sensitive as a SPAD or interferometric methods.Nonetheless, combining computational imaging techniques [48] such as patterned illumination schemes [49][50][51][52][53], or foveated sensing with point illumination can help enhance the sensitivity to inhomogeneities, or allow probing through larger scattering volumes.Finally, the analog bandwidth of the existing neuromorphic sensors restricts their operating speed to <10 kHz.This can limit the utility of the approach in application like CBF monitoring.This challenge can be mitigated by using emerging sensors, such as asynchronous SPAD arrays that possess functionality similar to a neuromorphic sensor [54].

Fig. 2 .
Fig. 2. Simulation of the spectral evolution of speckle as a function of the number of transport mean free paths (A) As the number of transport mean free path increase the speckle intensity fluctuates more rapidly, indicating a faster decorrelation.B) The histogram of spectral gradients of the simulated speckle intensity stack.As the scattering severity increases, speckle decorrelates faster resulting in an increase in the spread in the distribution of gradients (see Eq.(8))

Fig. 3 .
Fig. 3. Simulation of the temporal evolution of speckle as a function of the changing dynamic properties (A) As the Brownian diffusion coefficient (D B ) increases the speckle intensity fluctuates more rapidly, indicating a faster decorrelation.B) The histogram of spectral gradients of the simulated speckle intensity stack.As D B increases, speckle decorrelates faster resulting in an increase in the spread in the distribution of gradients (see Eq. (8))

Fig. 4 .
Fig. 4. The speckle differencing system used to perform wide field detection of inhomogeneities embedded within a scattering medium.A) Schematic of the setup used in the experiments.The scattering medium shown here is used for static scattering demonstrations B) Processing workflow of the speckle differencing approach.

Fig. 6 .
Fig. 6.A) Liquid phantom used for demonstrating the detection of inhomogeneities within a dynamic scattering medium B) With a single inhomogeneous region, C) Two inhomogeneous regions, with lower rate of speckle fluctuation compared to background

Fig. 7 .
Fig. 7. Probing a scattering medium using neuromorphic sensors.A) Simulation of the response of a neuromorphic pixel to spectrally evolving speckle from scattering media with different transport properties.The spiking output is proportional to the rate of fluctuation of the speckle pattern.B) Experimental results of imaging inhomogeneities embedded within a 60ℓ * thick static scattering medium using neuromorphic sensor.

Table 1 . List of symbols and notation Symbol/Notation Description a :
[a 1 , b 1 , c 1 ]Position vector associated with a point within the medium I (∆ν; r s , r d ) normalized intensity correlation between the speckle intensities observed at optical frequencies ν Hz and ν + ∆ν Hz g 1 (τ; r s , r d ) normalized field correlation of the speckle fields observed at times t seconds and (t + τ) seconds g 2 (τ; r s , r d ) normalized intensity correlation of the speckle intensities observed at times t seconds and (t + τ) seconds θ Intensity threshold N Number of speckle cells at a pixel within the integration time of the camera