Measuring light transport properties using speckle patterns as structured illumination

The measurement of light absorption and scattering properties of biological materials has several diagnostic and therapeutic applications. We can measure these properties for skin without contact using structured illumination and imaging. However, building simple handheld devices remains challenging due to motion artefacts and moving targets. To overcome this limitation, we project random speckle patterns instead of discrete spatial frequencies on the target. Since random patterns are spatially broadband, they capture more information per image, enabling frame-by-frame analysis. In this paper, we describe the statistics of objective speckles and demonstrate how the optical system is designed for spatially bandlimited illumination. Next, we use a diverse set of liquid tissue phantom to validate the method. We successfully demonstrate that a calibrated instrument can independently predict the two primary light transport properties of a homogeneous turbid system. This work is a starting point for analysing skin and other heterogeneous biological media in the future.


Supplementary Material
Fully Developed Speckles as White Noise Figure S1 illustrates a simple optical configuration to generate objective speckle patterns. A collimated beam of coherent light is incident on a rough surface with a randomly varying surface texture height f (x , y ). The scattered light forms the characteristic granular speckles on the parallel screen kept at a distance D. The intensity u of the light pattern at any point (x, y) or ρ on the screen is proportional the square of the complex electric field amplitude E: Here, η is the characteristic impedance of the medium, and usually suppressed for brevity. Since singly-scattered light retains its coherence, wavefronts that originate from the many scattering points on the rough surface interfere on the screen. The electric field E ( ρ) is therefore expressed as the sum of multiple phasors: Here, E s ( ρ , ρ) is the electric field amplitude of the scattered wavefront at a point ρ on the rough surface, towards the point ρ on the screen. λ is the wavelength of incident light, r is the scalar distance between points ρ and ρ. The phasor E s ( ρ , ρ) has an additional independent random phase angle proportional to f (x, y) because of the corresponding path length difference. If f (x, y) is of the order of the λ , which is usually the case, electric field E ( ρ) is effectively a sum of random phasors.
It is possible to describe the first-order statistics of speckle patterns by analysing sums of random phasors. However, for an intuitive explanation of the second-order statistics, we proceed to simplify the expression in eq. (S2) for far-field speckles using the paraxial approximation and the plane wave assumption. The paraxial approximation (D | ρ − ρ|) is used in Fresnel Diffraction Formula for free space propagation. Using it we can approximate distance r and therefore the Electric field E( ρ) as follows: We rewrite ( ρ /λD) as spatial frequency k. This reduces the integral in eq. (S4) to a Fourier Transform: The plane wave assumption ( D /| ρ − ρ| π | ρ − ρ| /λ) is used in far-field Fraunhofer Diffraction to further extend the paraxial approximation. It allows us to ignore the quadratic phase terms in eq. (S5) and hence simplify the expressions for Electric Field E( ρ) and Intensity u( ρ): Figure S1. Objective speckles generated due to interference of coherent light reflected of a rough surface The Intensity distribution u on the screen is therefore approximately proportional to the square of the Fourier transform of the scattered electric field E s .
We assume that the surface profile f (x , y ) is a white noise random process in 2-dimensions. If f is uniformly distributed over a range larger than wavelength λ , phasor angle ∠E s varies uniformly between −π and π. Assuming the scattering is isotropic, or the phasor magnitude |E s | is constant, we rewrite the Fourier transform in eq. (S7) as: Here A is the area of the illuminated region on scattering surface. The cosine of a random variable is another random variable, varying between −1 and 1. In eq. (S9), ∠E s (x , y ) is a uniform white noise random process, and it's cosine is another such process with a non-uniform distribution. Regardless of the distribution, the Fourier transform of a random cosine function is random itself. This is illustrated in Figures S2a to S2e with the example of a discrete function is a random variable uniformly distributed between −π and π, independent of n. The square of the Fourier transform |X[Ω]| 2 is a white noise random number with a χ 2 distribution with 2 degrees of freedom. Similarly, the square of the 2D Fourier transform in eq. (S9) also results in a white noise random variable with a χ 2 distribution with 2 degrees of freedom. The speckle pattern, given the above assumptions is therefore a white noise pattern with a vanishing specular field. The size of individual grains is diffraction limited to wavelength λ . Such speckles are known as Fully Developed Speckle Patterns.

Stationarity and Ergodicity of Speckle Patterns
In the present work, we are interpreting random speckle patterns as broadband spatial input for system identification. This relies on the assumption that input speckle patterns are stationary and ergodic. To experimentally validate these assumptions, we collected sets of N = 10 independent speckle pattern images for three different aperture sizes. Each optical configuration represents a random process U(x, y) and each image is a finite area realization u i (x, y); i = 1 . . . N of this random process. A cropped portion of one observed speckle pattern for each configuration is shown in Figures S3a to S3c (1) for illustration. To experimentally characterize a random process is challenging with finite data. Therefore, we limit the scope to validating that the process is Wide-Sense Stationary (WSS) and ergodic in the mean. We first discuss the stationarity because ergodicity is not defined for non-stationary processes.
For a time-series WSS process, the first and second order statistics remain invariant with time. A necessary and sufficient condition for WSS stationarity is that the Auto-Correlation Function (ACF) is dependent only on time delay and independent of position in time. We extend this to 2D processes in XY space. Figures S3a to S3c (2) give the estimated ACFR UU (∆x, ∆y) of U(x, y) for each configuration. We arrive at the estimate by averaging the biased ACF estimates for all u i (x, y) realizations: Here P = 1280 and Q = 1024 are the number of pixels captured per realization along X and Y directions. The ACF estimate has elliptical contours because of the elliptical shape of the aperture. For clarity, Figures S3a to S3c (3) give the individual and averaged ACF along the X direction (∆y = 0). Figures S3a to S3c (4) give the same along the Y direction (∆x = 0). Since the ACF estimate decays much faster than the maximum observed spatial delay, we can confidently claim that the process is WSS. Note that the images were mean-centred and normalized before estimating the ACFs, hence auto-correlation is the same as auto-covariance.
A time-series process is ergodic in the mean if the average along time almost always equals the average across multiple realizations. In essence, it ascertains that each independent realization of a stationary random process is statistically identical. Figures S3a to S3c (5) give the cross-correlationsR u i u j (∆x, ∆y); i = j, averaged for all N C 2 = 45 possible pairs of realizations. Since the cross-correlation is always below 5 × 10 −3 , it is statistically insignificant. Therefore, we can consider each realization independent of all others. Figures S3a to S3c (6) gives a histogram of the means of pixel readings across realizations and the mean estimated across XY space for each realization. The average mean estimate is close to the estimate from each realization. Figures S3a to S3c (6) gives a histogram of the standard deviations of pixel readings across realizations and that estimated across XY space for each realization. We arrive at the same conclusion as earlier. We can therefore confidently claim that the process is ergodic in the mean.

4/9
(1) Light flux density in the semi-infinite medium and the radially-symmetric backscatter impulse response h(r), predicted using 100,000 photon scattering simulations

Backscatter Impulse and Frequency Response
When a narrow beam of light with unit power is incident at any point (x 0 , y 0 ) on a turbid medium, the diffuse backscatter is the linear 2D impulse response h(x, y; x 0 , y 0 ). For a uniform semi-infinite medium, the impulse response is independent of position (x 0 , y 0 ) and reduces to a radially symmetric function h (r), where r is the scalar distance from point (x 0 , y 0 ). We estimate the impulse response h(r) for a given turbid medium with known light transport properties using White Monte Carlo (WMC) Ray Tracing simulations. The results for one such simulation are shown in Figure S4 for illustration. Next, we define the spatial frequency response H(k x , k y ) as the Fourier transform of 2D impulse response h(x, y). Since h(x, y) reduces to a radially symmetric function h (r), the 2D Fourier transform also reduces to a zeroth-order 1D Hankel transform. Figure S5 presents the impulse and frequency responses estimated for a range of light transport properties (µ a , µ s and g).