{'fraction_non_alphanumeric': 0.03614561027837259, 'fraction_numerical': 0.016274089935760173, 'mean_word_length': 4.595734483584951, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Automatic Differentiation for All Photons Imaging to See Inside Volumetric Scattering Media

Imaging through dense scattering media - such as biological tissue, fog, and smoke - has applications in the medical and robotics fields. We propose a new framework using automatic differentiation for All Photons Imaging through homogeneous scattering media with unknown optical properties for non-invasive sensing and diagnostics. We overcome the need for the imaging target to be visible to the illumination source in All Photons Imaging, enabling practical and non-invasive imaging through turbid media with a simple optical setup. Our method does not require calibration to acquire the sensor position or optical properties of the media.


Introduction
Imaging through dense, volumetric, scattering media has been a significant challenge for various applications, such as imaging through skin tissue and industrial remote sensing through the smoke. In conventional imaging, photons directly travel from the target to the sensor. However, in scattering media, photons randomly move, blurring the imaging target both in time and space (Figure 1, right). This blur makes it challenging to resolve small features of the target from the measurement.
Several approaches have been developed to tackle the challenge of imaging through a scattering medium. When the medium is not dense or thick, ballistic photons, that do not undergo scattering, can be used, effectively ignoring the effect of scattering. Optical coherence tomography [1], confocal microscopy [2], and two-photon imaging [3] optically reject scattered photons to capture only ballistic photons. Time-of-flight measurements enable the time-gating of ballistic photons [4], so that the measurement only relies upon the first photon to hit the Reconstruction Complexity Simplicity of Imaging System (Less Need of Calibrations) All Photons Imaging [26] Traditional DOT Computational DOT [24] This Work (Simple Imaging Setup, Calibration-free) Figure 2: Towards imaging more complex scene with simple imaging system. Our technique exploits auto-differential diffusion approximation to simultaneously recover the target and the optical properties of the scattering media.
remote sensing. In this research, we capture the time-resolved measurement of the spatio-temporal profile of the scattered photons and recover the optical properties and 2D target embedded in the dense scattering media. While the diffusion approximation we use is well studied in the past [19,23], our adaptation of a modern automatic differentiation framework [24] enables the robust estimation of the optical properties in complex scenes with a 2D target embedded within the scattering media.
Relation to closely related prior art: Our method uses the same imaging system as All Photons Imaging (API) [22], but extends API to the a realistic scene (a 2D target embedded inside a phantom, instead of the target on top of a phantom) with the use of automatic differentiation to the forward model. We use the same optical setup as Lyons et al. [17]. While Lyons et al. [17] required a prior knowledge on the optical properties of the scattering media, our technique jointly recovers the unknown optical properties and embedded 2D target. In contrast to traditional DOT, our technique requires fewer priors, such the position of illumination sources and detectors, and is capable of single-shot measurements. Fig 2 summarizes the relation of this work to previous works.

Forward Model
The radiative transfer function (RTE) describes the random walk of photons through a homogeneous medium: The left two terms represent the change of the energy (per unit volume per unit solid angle) and the divergence of the beam, while the right three terms represent the extinction, scattered energy from the surrounding, and light source. µ a and µ s respectively represent the expected number of of absorption and scattering per unit light path, and µ t is a sum of µ a and µ s . The terms r, ω, t denote the location, direction and time respectively. c is a speed of light in the media. L(r, ω, t), P (ω, ω ) denote the radiance and phase function respectively. The phase function is often represented by the anisotropy parameter g =< cos θ >, where θ is the angle between the old and new trajectories.
The RTE is hard to solve, so two approaches are often used to approximate. When the number of scatterings is sufficiently large and µ a << µ s , the diffusion equation approximates the RTE, where the photon fluence rate φ(r, t) = L(r, ω, t)dΩ satisfies the diffusion equation. The diffusion approximation becomes less accurate as the number of scatterings becomes small, as the propagation of photons cannot be considered diffusion. In situations with limited scattering, a Monte Carlo approximation is used to approximate the RTE. However, Monte Carlo approximations are computationally expensive as they involve the simulation of the random walks of millions of photons.
Here, we consider dense and thick scattering media, where the diffusion approximation can provide a reasonably accurate solution to the RTE. The diffusion equation is where D = [3(µ a + µ s (1 − g))] −1 is the diffusion coefficient. For an isotropic light source with short time pulse, the photon fluence rate in the infinite scattering media is written as .
As illustrated in Fig. 3, our imaging setup is in transmission mode, where a target is between scattering media. Unlike traditional DOT, the goal of our technique is to reconstruction a 2D target within a scattering media with high contrast. We express the measurement as the result of a 3D convolution with K 1 (x, y, t) (scattering between illumination and target), followed by a 2D multiplication with T (x, y) (2D mask) and another 3D convolution with x t x t Figure 3: A sketch of the forward model and reconstruction. (a) The illumination (I(x, y)) is blurred as photons scatter in the first scattering media (K 1 (x, y, t; D, d 1 )) before passing through the 2D target (T (x, y)). Photons scatter in the second scattering media (K 2 (x, y, t; D, d 2 )) before reaching the detector. We illustrate the forward model in 2D for simplicity. Our reconstruction method takes the 3D measurement over x, y, and t for input. (b) Our technique estimates the unknown optical properties of the media by minimizing the observed measurement m(x, y, t) and estimated measurement (m(x, y, t)). I(x, y, t) denotes the illumination on the surface of the scattering media. Our forward model is similar to the model proposed by Lyons et al. [17]. In this paper, we consider a 2D target with sparse opening pixels. Therefore, the geometry of the phantoms before and after the phantom can be approximated as slobs. We follow diffusion approximation model proposed by Patterson et al. [23] to write the 3D kernel K 1 (x, y, t; D, d 1 ), K 2 (x, y, t; D, d 2 ) that represent the propagation of photons through the scattering media with thickness d as follows.
where z 0 = [µ s (1 − g)] −1 and |r| 2 = x 2 + y 2 + d 2 . This is different from Eq. 2 because we consider the photon transmittance of a homogeneous slab media that is bounded along optical axis, and infinite in the other two axes. Applying two boundary conditions for surfaces of the slab geometry results in Eq. 5.

Reconstruction Algorithm
Our method estimates the 2D target as well as the optical properties of the homogeneous scattering media from the measurement. We write the objective of the algorithm as follows.
where Γ(T (x, y)) is a regularization on the target. We solve this problem in the model by minimizing with respect to the optical property and the target in an alternate manner. Scattering coefficient µ s , anisotropy g, and absorption coefficient µ a are commonly used to represent the scattering properties of media. In the diffusion paradigm, the scattering coefficient and anisotropy are interchangeable, so reduced scattering coefficient µ s = µ s (1 − g) is used instead. While we assume µ a << µ s , µ a still has a significant impact on the 3D kernel K(x, y, t). Hence, our goal is to estimate µ s and µ a from the measurement.
Estimation of µ a : As Chance et al. [25] suggest, the absorption coefficient µ a can be estimated with the asymptotic slope of the log of the measurement along the time dimension. Let us consider the blur kernel K(x, y, t; D, d) for large t (See Supplemental Material for the derivation): Hence. t is sufficiently large µ a can be estimated aŝ Our forward model is more complex than the expression of the 3D kernel in Eq. 5. However, the insight from Eq. 7 shows that the slope of the log of the measurement is dominantly determined by µ a . Hence, we estimate the absorption coefficient as written in Eq. 8.
Estimation of µ s : While µ a can be estimated because absorption dominates the intensity fall-off for sufficiently large t, an analytical expression of µ s in terms of the measurement is challenging. We numerically estimate µ s with gradient descent. As shown in Eq. 4 and Eq. 5, the forward model is fairly complex. Our work uses recent advancements in auto-differentiation to calculate the gradient of the objective function (Eq. 6) with respect to µ s . Automatic differentiation efficiently and accurately computes the derivative of functions from the forward model. It has been exploited for robust design optimization [26] computational fluid dynamics [27], and has recently gained significant interest from machine learning community [28]. In imaging, automatic differentiation has been used in DOT with more complex forward model [29,30], fluorescent lifetime imaging [31]. We use JAX [24] to implement the forward model in Eq. 4 and compute the gradient with respect to µ s , and RMSProp as the optimizer [32]. While Patterson et al. [23] suggest a closed-form estimate of the reduced scattering coefficient from the measurement on slab geometry, our imaging setup contains a 2D target in the media, which introduces model mismatch. Moreover, such a closed-form solution depends heavily on the estimate of the absorption coefficient, which can vary in low SNR regime. Therefore, a reduced scattering coefficient estimated by this method is not accurate, resulting in an inferior reconstruction of the target (See Supplemental Material Section 3). The use of automatic differentiation on the forward model with diffusion approximation enables robust recovery of the target, as our technique aims to minimize the loss between the expected measurement and the actual measurement.
Estimation of the target: Given µ s and µ a , the measurement can be computed by Eq. 4. The target can be recovered by gradient descent with the gradient computed from automatic differentiation. However, the model can be simplified as a linear system to solve the minimization problem efficiently: where b and v denote the vectorized 3D measurement m(x, y, t) and 2D target T (x, y) respectively. A represents the mapping between the measurement target, and can be obtained by computing the expected measurement for each pixel of the target. For example, ith column of A is the expected measurement for the case v i = 1, and given as where T i (x, y) is 1 for the ith pixel and 0 for all the other pixels. vec(·) denotes an operator that vectorizes 3D matrix into a normalized column vector. Eq. 9 is ill-posed because A is a blurring operator. Therefore, we solve this problem with l1 and total variation regularizer: The prior term is written as follows: with Algorithm 1: ADMM algorithm to solve Eq. 11 Input: Input: v 0 , z 0 = vec (T (x, y)), u 0 = 0, ρ, λ 1 , λ 2 while not converged do where D x , D y are matrices that compute the gradient of the target image in x and y directions, respectively. We use the alternating direction method of multipliers (ADMM) [33] to solve this minimization problem. The iterative steps of ADMM are summarized in Alg. 3.

Number of iterations
While we use l1 sparsity prior to our target, our formulation could extend to types of common priors such as total-variation or sparsity over the discrete some bases (e.g., discrete cosine transform). When such priors are not available, Tikhonov (l2) regularization can be used (See Supplemental Material).

Experimental Demonstration
We demonstrate our technique with both simulated and experimental data.  the target. The simulated measurement does not contain any sensor noise, such as dark noise. The scattering coefficient, anisotropy, and absorption coefficient of the simulated phantom were 20.0 mm -1 , 0.9, and 0.01 mm -1 . The reduced scattering coefficient is 2.0 mm -1 . Fig. 4 shows the recovered target and reduced scattering coefficient for the first 50 iterations of µ s estimation. The estimated reduced scattering coefficient and the absorption coefficient was 2.04 mm -1 and 0.011 mm -1 . Fig. 4 shows that different initialization converges close to the ground truth. Fig. 5 illustrates the experimental setup. The thickness of the phantom was 26 mm in total, where the target is placed in the middle. Fig. 6 shows the reconstruction results in 32 by 32 resolution with a pixel size of 2 mm, compared with time-averaged and ballistic-photon imaging. We integrated 2D frames with respect time for the time-averaged frame, which is equivalent to an image taken without the temporal resolution. To generate a frame of early-arriving photons, we extract the first photons above the noise level for each pixel. (See supplemental video for the raw 3D measurement.) The time-averaged frame shows blur due to the scattering, and the ballistic frame is noisy because the number of ballistic photons is small. Our method successfully recovers the shape of the 2D target. The estimated reduced scattering and absorption coefficient are 2.45 mm -1 and 0.0082 mm -1 for "A" target, and 2.39 mm -1 and 0.0088 mm -1 for the "wedge" target. Fig. 7 shows reconstruction results for different acquisition time. Our method is robust even when the SNR of the measurement is so low that the time-averaged and ballistic frames are noisy. Our technique reconstructed the "A" target with 150 ms exposure time. This shows the benefit of our simple optical setup. Fast acquisition is challenging with traditional DOT, which requires raster scanning of the illumination source. The accuracy of the estimation of the optical properties decreases as the acquisition time decreases. For example, for the acquisition time of 50ms, the estimated reduced scattering coefficient and absorption coefficient are 2.04 mm -1 and 0.0031 mm -1 . This is due to the estimation of the absorption coefficient from the noisy temporal histogram from a small number of photons. However, as shown in Fig. 8 (a)-(e), the estimated blur kernels K(x, y, t; D, d) for different acquisition times are still similar to the ground truth. Fig. 8 (f) plots the normalized intensity over time. The estimated kernel roughly follows the ground truth blur kernel, which we generated from the estimate of the optical properties with phantoms with no target.

Discussion
Practical limitations: In this paper, we demonstrate reconstruction of a 2D mask in an unknown homogeneous volumetric scattering media. Our technique extends the limitations of previous works in complexity [22] or ease of use without calibration [17]. The use of automatic differentiation provides a solution towards practical use of previously proposed works. However, our technique still requires prior knowledge of the target depth, and the transmission-mode imaging system limits the reconstruction to 2D shapes. We also only consider flat surfaces for the illumination and detection.

Resolution:
We used 2D targets that consist of two 2 mm by 10 mm bars that are 5 mm, 10 mm, 15 mm away from each other to evaluate the resolution of our system. Figure 9 shows the 1D reconstruction of the targets. This illustrates  Figure 9: Reconstruction of bars separated by 15 mm, 10 mm, 5 mm space to evaluate the resolution of our methods compared to timeaveraged and ballistic reconstruction. Time-averaged reconstruction suffers from blur due to scattering, and ballistic reconstruction suffers from low SNR. Our method resolves two points separated by 15mm and 10mm.
the resolution of our method for the same scattering media described in the experimental demonstration as compared to time-averaged and ballistic photon imaging. The resolution of imaging is determined by the blur from the scattering media. The illustrated resolution should be worse as the scattering coefficient or the thickness of the phantom increase. Traditional DOT places illumination sources and detectors at different locations to capture more information for better reconstruction quality. In contrast, the goal of our technique is to improve the practicality of the imaging system.
Need for Time-of-Flight Measurement: Time-of-flight information plays an important role in the estimation of the diffusion coefficient. Our algorithm initializes an estimation of the target with the time-averaged frame. Without the time of flight, the estimation of the blur kernels K 1 (x, y, t; D, d 1 ) and K 2 (x, y, t; D, d 2 ) becomes an impulse response as the estimate becomes equivalent to the measurement.
The Use of Diffusion Approximation and Linear Model: While the linear model provides an approximation to RTE, the accuracy of the reconstruction is limited. Furthermore, the linear model may introduce an error because it does not model photon's random walk from the second phantom back to the first phantom across the target plane. Therefore, Monte Carlo-based reconstruction should produce more accurate reconstruction result [18]. However, such a method requires many iterations of Monte-Carlo rendering, which is computationally expensive. In contrast, diffusion-approximation-based reconstruction is more computationally efficient. In practice, our method can be used for the initialization of the Monte Carlo-based approach for faster convergence.

Conclusion:
We propose a novel algorithm to reconstruct a mask in dense, scattering media with unknown optical properties and demonstrate that it can reconstruct targets embedded in thick, scattering media with limited calibration.
The automatic differentiation to the forward model overcomes the challenge of All Photons Imaging, and enables recovery of object embedded in the dense volumetric scattering media. We hope that our approach inspires more practical imaging through scattering media without extensive calibrations for biomedical and remote sensing applications.

Methods
Optical setup: Fig. 5 describes the optical setup. The SPAD array is focused on the surface of the tissue phantom, and one pixel corresponds to a 2 mm square patch on the surface. The illumination source is a SuperK super-continuum laser, where white light is filtered at 628 nm with acousto-optic tunable filters followed by a 50-degree circle diffuser (ThorLab, ED1-C50-MD), which was placed 6.4cm away from the phantom. The instrument response has roughly 150 ps full-width-half-max (FWHM). Each measurement is 32 x 32 x 70 tensor with a 50ps time bin. Exposure time is 0.5 ms for each frame of the SPAD array, and we measure 200000 frames for the complete measurement. The total acquisition time to generate reconstruction in Fig. 6 is 100 seconds. The biomimetic phantom consists of epoxy resin (East Coast Resin) as a base and TiO 2 powder (Pure Organic Ingredients) as a scattering agent [34]. The refractive index of the media is 1.4 [34] and the total thickness of the phantom is 26mm (13mm on each side of the target). We measured the absorption coefficient of the phantom as 0.080 mm -1 using Eq. 8 and the reduced scattering coefficient as 2.32 mm -1 by fitting the measurement with Eq. 5. For simulated measurement, we use the same parameters for the detector, and the illumination source was placed 10 cm away from the phantom. Reconstruction was performed on a 32 x 32 x 60 3D x-y-t measurement with a 50 ps time bin.
Algorithm parameters: All the reconstruction on experimental data was performed with ρ = 1, λ 1 = 8 × 10 −4 , and λ 2 = 10 −3 λ 1 . The update of µ s was performed with a step size of 5 × 10 −2 , and terminated when a certain number of iterations are performed or when the change of value is less than 0.01% in a single iteration. Each µ s and T (x, y) has 300 and 50 iterations, alternated three times. The initial estimation of the reduced scattering coefficient was 1.5mm −1 for all of our reconstruction. We used the time-averaged frame as the initialization of the target estimate. The reconstruction on simulated data was performed with ρ = 1, λ 1 = 1 × 10 −4 , and λ 2 = 10 −3 λ 1 . Each µ s and T (x, y) has 300 and 50 iterations, alternated three times. We used larger λ 1 for the experimental data because of higher noise level.

Supplementary Videos
Video S1 shows the 3D measurement of the "A" target for different exposure time (100 s, 500 ms, 150 ms, 50 ms) that are used for the reconstruction. Video S2 shows the 3D measurement for the "wedge" target. Video S3 shows the3D measurement for the bars that we used to evaluate the resolution of our imaging system.

Estimation of Absorption Coefficient
The left side of Eq.7 in the main manuscript can be expanded as follows: We define f (t) as following: Since log(f (t)) = f (t)/f (t) and Since f (t) has 1 t 2 term while f (x) doesn't, lim t=∞ log(f (t)) = 0. Therefore, the term on the second and third line of Eq. 14 becomes 0. Reconstruction with closed-form estimation of reduced scattering coefficient Patterson et al. [23] proposed a closed-form estimate of the reduced scattering coefficient from the time-resolved measurement. The forward model is complex, and the previously proposed method of the optical estimation suffers from model mismatch. Our framework is more robust than closed-form estimation, as our method estimates reduced scattering coefficient by minimizing the L2 loss between the forward model and actual measurement. Fig. 10 illustrates the reconstruction results with closed-form estimation of the optical properties. Because the recovery of the 3D blur kernel is not robust, the quality of the reconstruction is worse than our method (Fig. 6 in the main manuscript).

Reconstruction with Tikhonov Regularizartion
Our reconstruction result exploited l1 sparsity and small total variation for the prior. This prior is valid for various practical scenarios, such as fluorescence imaging. Here, we also show reconstruction results with Tikhonov regularization, which is more general:v With Tikhonov regularization, the target can be estimated aŝ Fig. 11 shows the reconstruction results of "A" and Wedge targets with λ = 0.1. Though the reconstruction is more blurry than the results with l1 + total variation regularization, it still reveals details in the target that are not obvious from time-averaged and ballistic frames. For example, the horizontal line of the "A" target is visible in the reconstruction with the Tikhonov regularization.  Figure 11: Reconstruction results with Tikhonov regularization. While Tikhonov regularization does not make the reconstruction as sharp as those shown in Fig.6 in the main text, it stills reveals details (e.g., the horizontal line in "A") that are not apparent from the time-averaged image or ballistic frame.