Reflection-mode diffraction tomography of multiple-scattering samples on a reflective substrate from intensity images

Our reflection-mode DT technique aims to bridge the gap between transmission-mode DT and OCT, by simultaneously and quantitatively characterizing the 3D FS and BS information of a sample placed on a highly reflective substrate, as illustrated in Fig. 1(a). With the assistance of the reflective substrate, such as a silver mirror, the measured scattered field contains both forward- and back-scattered signals, each conveying complementary information about the sample, as shown in the zoomed-in inset of Fig. 1(a). To perform the reconstruction, we first take brightfield intensity images under varying illumination angles $I_{\text{m},\text{i}}$.

Once the intensities are collected, the 3D FS and BS information is retrieved by solving the inverse scattering problem to reconstruct the RI distribution, as shown in the iterative loop in Fig. 1(b). This process involves two main stages: the forward and backward processes.

In the forward process, MBS serves as an accurate and fast-converging forward model, which has demonstrated its state-of-the-art performance in transmission-mode DT for complex biological samples [26]. Starting with an initial guess of the RI (typically the background RI), MBS iteratively refines the electric field $\boldsymbol{\psi}$, incorporating multiple scattering and handling reflections from the substrate. After convergence, the field on the simulation domain’s surface ($z=z_{0}$) is projected onto the camera plane through a detection operator $\hat{D}$. The simulated intensity $I_{\text{s,i}}$ is compared with the measured data, and the loss function $\mathcal{L}_{\text{i}}$ is computed.

In the backward process, AM is used to compute the RI gradient for updating its distribution. An adjoint simulation is performed, generating an adjoint field $\boldsymbol{\psi}_{a}$ with a mirrored incident direction. The RI gradient is obtained by element-wise multiplication of $\boldsymbol{\psi}_{a}$ and $\boldsymbol{\psi}$. An example of the calculated gradient in a 2D simulation is shown in Fig. 1(c). Due to the capability of MBS to accurately model the complex multiple scattering processes, the Fourier spectrum of the computed gradient demonstrates that both FS and BS information are efficiently extracted during the reconstruction process, inherently concentrating in the regions around FS and BS supports predicted by the single-scattering model.

After each gradient descent step, the updated RI is used for the next iteration with new measurements under different illumination conditions, continuing until the loss function converges. Following reconstruction, an additional processing step is applied to the Fourier spectrum of the reconstructed RI to isolate the corresponding regions of FS and BS information for further analysis. This outlines the main workflow of our reflection-mode DT reconstruction strategy, with further details elaborated in the following sections.

MBS is an efficient and rigorous forward model for solving the inhomogeneous vectorial Helmholtz equations, which can be expressed as,

$$\nabla\times\nabla\times\boldsymbol{\psi}(\boldsymbol{r})-k^{2}(\boldsymbol{r})\boldsymbol{\psi}(\boldsymbol{r})=\boldsymbol{S}(\boldsymbol{r}),$$

(1)

where $\boldsymbol{\psi}(\boldsymbol{r})$ denotes the electric field at spatial position $\boldsymbol{r}\equiv(\boldsymbol{r}_{\parallel},z)$, $\boldsymbol{S}(\boldsymbol{r})$ is the electric current source at wavelength $\lambda$ and $k(\boldsymbol{r})=k_{0}n(\boldsymbol{r})$ is the wavenumber in an isotropic medium with $n(\boldsymbol{r})$ representing the RI distribution, and $k_{0}=2\pi/\lambda$ is the vacuum wavenumber. By using Green’s function, its solution can be formulated as,

$$\boldsymbol{\psi}(\boldsymbol{r})=\overset{\leftrightarrow}{\mathscr{G}}V(\boldsymbol{r})\boldsymbol{\psi}(\boldsymbol{r})+\overset{\leftrightarrow}{\mathscr{G}}\boldsymbol{S}(\boldsymbol{r}),$$

(2)

where $\overset{\leftrightarrow}{\mathscr{G}}\equiv\mathscr{F}^{-1}\big{[}\tilde{\overset{\leftrightarrow}{\boldsymbol{g}}}(\boldsymbol{q})\cdot\mathscr{F}\left(\cdot\right)\big{]}$ is the Green’s function operator, implemented via the angular spectrum method. $V(\boldsymbol{r})\equiv k(\boldsymbol{r})-k_{b}^{2}-i\eta$ is the scattering potential, where $k_{b}\equiv n_{b}k_{0}$ is the wavenumber of the background medium with RI $n_{b}$, and $\eta$ denotes its loss. Notably, the choice of $n_{b}$ and $\eta$ does not affect the solution of Eq. (1). $\mathscr{F}$ and $\mathscr{F}^{-1}$ denote the Fourier and inverse Fourier transforms. The dyadic Green’s function in the background medium has the Fourier transform $\tilde{\overset{\leftrightarrow}{\boldsymbol{g}}}=[\overset{\leftrightarrow}{\boldsymbol{I}}-\boldsymbol{q}\boldsymbol{q}^{T}/(k_{b}^{2}+i\eta)]/(|\boldsymbol{q}|^{2}-k_{b}^{2}-i\eta)$, where $\overset{\leftrightarrow}{\boldsymbol{I}}$ is the identity tensor, and $\boldsymbol{q}$ is the spatial frequency coordinate.

Expanding $\boldsymbol{\psi}(\boldsymbol{r})$ based on Eq. (2) recursively leads to the conventional Born series,

$$\boldsymbol{\psi}_{\mathrm{B}}(\boldsymbol{r})=[\overset{\leftrightarrow}{\boldsymbol{I}}+\overset{\leftrightarrow}{\mathscr{G}}V(\boldsymbol{r})+\overset{\leftrightarrow}{\mathscr{G}}V(\boldsymbol{r})\overset{\leftrightarrow}{\mathscr{G}}V(\boldsymbol{r})+\cdots]\overset{\leftrightarrow}{\mathscr{G}}\boldsymbol{S}(\boldsymbol{r}).$$

(3)

However, the arbitrary choice of $\eta$ in conventional Born series causes serious divergence issues during recursions [25]. To address this, G. Osnabrugge et al. proposed a selection strategy for $\eta$, where $\eta$ is chosen as $\eta\geq\max|k^{2}(\boldsymbol{r})-k_{b}^{2}|$ to reduce the influence radius of the Green’s function from the original entire simulation domain to a spherically decaying region, ensuring the stability across all scattering regions. The effect of $\eta$ is illustrated by the red band matching as the series order increases in Fig. 1(b). Accordingly, the energy loss is compensated by the $\eta$ term in the scattering potential $V(\boldsymbol{r})$. A preconditioner $\gamma(\boldsymbol{r})=iV(\boldsymbol{r})/\eta$ is then applied to modify Eq. (3), yielding the recursive form of MBS:

$$\boldsymbol{\psi}_{\mathrm{M}}(\boldsymbol{r})=[\overset{\leftrightarrow}{\boldsymbol{I}}+\overset{\leftrightarrow}{\mathscr{M}}+\overset{\leftrightarrow}{\mathscr{M}}\overset{\leftrightarrow}{\mathscr{M}}+\cdots]\gamma(\boldsymbol{r})\overset{\leftrightarrow}{\mathscr{G}}\boldsymbol{S}(\boldsymbol{r}),$$

(4)

where $\overset{\leftrightarrow}{\mathscr{M}}\equiv\gamma(\boldsymbol{r})\overset{\leftrightarrow}{\mathscr{G}}V(\boldsymbol{r})-\gamma(\boldsymbol{r})+\overset{\leftrightarrow}{\boldsymbol{I}}$ (see mathematical details in Section 1 in Supplement 1).

To replicate real experimental conditions within a finite simulation domain, BCs are commonly applied to prevent artificial reflections and unwanted scattering. A widely used method is the absorbing BC, which absorbs the light waves as they exit the simulation domain in an open system. In MBS, absorbing BCs have been implemented using ultra-thin boundary layers [2], offering low computational and memory costs. To extend MBS for reflection-mode conditions, we further introduce two additional BCs, including Bloch BC (Fig. 2(a)) and PEC BC (Fig. 2(c)).

Bloch BC [36] is commonly used in numerical methods like FDTD to handle the oblique incidence wave $e^{i\boldsymbol{k}_{\mathrm{in}}\cdot\boldsymbol{r}}$, where $\boldsymbol{k}_{\mathrm{in}}=[\boldsymbol{k}_{\parallel},k_{z}]$. In FDTD, implementing Bloch BC is straightforward: a constant Bloch phase factor, $e^{\pm i\boldsymbol{k}_{\parallel}\cdot\boldsymbol{L}}$, is first applied to the boundary field before copying it to the opposite side. Here, $\boldsymbol{L}$ represents the dimensions of the simulation domain, and the $+$ and $-$ signs correspond to the right and left boundaries, respectively. This ensures the correct phase shift for waves exiting and re-entering the simulation domain.

The most straightforward way of applying the Bloch BC based on the angular spectrum method, with a similar way as in FDTD, requires first padding regions at least the same size as the simulation domain before applying the Bloch phase factor to the padded region, in order to both ensure the correct phase shift for waves exiting and re-entering the simulation domain (see details in Section 2 in Supplement 1). However, this approach results in too much memory overhead.

To overcome this, we propose adapting the acyclic convolution (ACC) method [11, 2] to achieve a padding-free Bloch BC in MBS, as illustrated in Fig. 2(a). First, a phase ramp $e^{-i\boldsymbol{k}_{\parallel}\cdot\boldsymbol{r}_{\parallel}}$ is applied to the entire simulation domain in the real space to effectively demodulate the field before performing discrete Fourier transform (DFT). Correspondingly, we use a frequency-shifted Green’s function $\tilde{\overset{\leftrightarrow}{\boldsymbol{g}}}(\boldsymbol{q}+\boldsymbol{k}_{\mathrm{\parallel}})$ to compute the scattered field. Finally, the scattered field is re-modulated by $e^{i\boldsymbol{k}_{\parallel}\cdot\boldsymbol{r}_{\parallel}}$ in the real space. Accordingly, the modified Green’s function operator is

$$\overset{\leftrightarrow}{\mathscr{G}}_{\mathrm{Bl}}\equiv\mathscr{F}_{\mathrm{Bl}}^{-1}\big{[}\tilde{\overset{\leftrightarrow}{\boldsymbol{g}}}(\boldsymbol{q}+\boldsymbol{k}_{\parallel})\cdot\mathscr{F}_{\mathrm{Bl}}\left(\cdot\right)\big{]},$$

(5)

where $\mathscr{F}_{\mathrm{Bl}}\left(\phi\right)\equiv\mathscr{F}\left(\phi e^{-i\boldsymbol{k}_{\parallel}\cdot\boldsymbol{r}_{\parallel}}\right)$ and $\mathscr{F}_{\mathrm{Bl}}^{-1}(\tilde{\phi})\equiv\mathscr{F}^{-1}(\tilde{\phi})e^{i\boldsymbol{k}_{\parallel}\cdot\boldsymbol{r}_{\parallel}}$. This method ensures that the fields leaking from the boundaries automatically satisfy the Bloch BC, without the need for additional padding or calculations. In Fig. 2(b), 2D MBS simulations on a circular scatterer with Bloch and periodic BCs are compared. The periodic BC introduces strong artifacts due to the mismatch between the incident field at the two ends. In contrast, our proposed padding-free Bloch BC ensures that the re-entering field seamlessly matches the internal field. Note that ACC is not limited to implementing Bloch BC in MBS; it can also be easily extended to other forward models based on the angular spectrum method, such as angular spectrum diffraction, single-scattering, and multi-slice methods.

We implemented a self-adaptive PEC boundary in MBS to simulate strong reflections from a substrate. While MBS can accurately model reflection at an interface, incorporating a high RI substrate typically requires finer grids and smaller propagation steps, significantly increasing computational complexity. To simplify this, we use a PEC boundary as a lossless reflection model for the substrate interface, providing insights into reflection-mode DT without excessive computational demands.

Based on the equivalence principle [35], a PEC boundary is equivalent to a surface electric current source $\boldsymbol{S}_{\text{surf}}$ on the boundary surface $\partial\Omega$. Therefore, applying a PEC BC is converted to solving for $\boldsymbol{S}_{\text{surf}}$ such that $\hat{\boldsymbol{n}}\times\boldsymbol{\psi}(\boldsymbol{r})\big{|}_{\partial\Omega}=0$.

At each MBS iteration, we incrementally add a surface current source $\boldsymbol{S}^{\prime}_{\text{surf,j}}\propto\hat{\boldsymbol{n}}\times\boldsymbol{\psi}(\boldsymbol{r})\big{|}_{\partial\Omega}$, as shown in Fig. 2(c). The recursive formula for the $j$-th MBS scattering field $\boldsymbol{\phi}_{\text{j}}$ can be written as

$$\boldsymbol{\phi}_{\text{j}+1}=\overset{\leftrightarrow}{\mathscr{M}}\boldsymbol{\phi}_{\text{j}}+i\boldsymbol{\gamma}\overset{\leftrightarrow}{\mathscr{G}}\left[\boldsymbol{\psi}_{\text{j}}\big{|}_{\partial\Omega}-\hat{\boldsymbol{n}}\cdot\boldsymbol{\psi}_{\text{j}}\big{|}_{\partial\Omega}\right],$$

(6)

where $\hat{\boldsymbol{n}}$ is the normal vector of $\partial\Omega$ and $\boldsymbol{\psi}_{\text{j}}$ is the total field at the $j$-th iteration. By defining the added surface current source as $\boldsymbol{S}_{\text{surf,j}}(\boldsymbol{r})=i\left[\boldsymbol{\psi}_{\text{j}}\big{|}_{\partial\Omega}-\hat{\boldsymbol{n}}\cdot\boldsymbol{\psi}_{\text{j}}\big{|}_{\partial\Omega}\right]$, the induced electric field cancels the tangential components of the total field estimated on the PEC surface. Moreover, when the PEC BC is fully satisfied, the added source will vanish, making this a self-adaptive mechanism for integrating the PEC BC in MBS (see mathematical details in Section 2 in Supplement 1).

A typical simulation domain in our framework is shown in Fig. 2(d). Four side walls are configured with Bloch BCs (purple lines), the bottom plane is assigned a PEC BC (gray line), and the top plane is set with an absorbing BC (double dashed line).

In reflection mode, unlike transmission mode, incident and scattered wave components intermingle within the simulation domain, creating an axial standing wave, as illustrated in Fig. 2(d). This overlap complicates field analysis, making it essential to isolate the scattered field from the total field for accurate interpretation. To achieve this, the simulation is divided into distinct regions for the total and scattered fields.

A one-way source is introduced to isolate the scattered field by generating a plane wave propagating in a single direction. This is achieved using a dual-layer source (red line), consisting of two staggered planar sources placed one grid point apart in the $z$-direction with an interval $\Delta z$. In the lateral direction, each planar source has a phase distribution $e^{i\boldsymbol{k}_{\parallel}\cdot\boldsymbol{r}_{\parallel}}$ to produce the desired oblique incidence. The bottom source also includes an additional phase factor $-\exp(-ik_{z}\Delta z)$, canceling any backward-propagating incident waves. The planar source is positioned two wavelengths below the absorbing boundary, with its $k_{z}$ aligned along the positive $z$-axis.

The region between the source and the PEC boundary serves as the total-field region, where structures are placed to simulate scattering. The area between the source and the absorbing boundary is the scattered-field region, containing only the pure scattered field from the structure and substrate. The $z_{0}$ plane is positioned in the middle of this gap to capture the scattered field, which is then projected onto the camera plane. The two-wavelength gap effectively eliminates potential evanescent waves near the source or boundary, allowing for accurate analysis of the scattered field.

To validate our method, we simulate the scattering field of a 3D sphere positioned near a PEC boundary, and benchmark the result against FDTD. A spherical scatterer, with a radius of 1.5 $\mu$m and a RI of 1.461, is placed 6 $\mu$m above the PEC boundary within a background medium of RI 1.336. The simulation domain, measuring 12 $\mu$m, is discretized into 20 nm voxels for both MBS and FDTD simulations. The scatterer is illuminated by a circularly polarized light with a wavelength of 532 nm at an incident angle of $30^{\circ}$. Fig. 2(e) compares the $x$-$z$ profiles of the simulated electric fields obtained from MBS and FDTD (additional comparisons are provided in Section 2 in Supplement 1). The simulations show strong agreement, capturing the full scattering lobes, reflections from the substrate, half-wave loss, and standing wave patterns in both cases.

After the forward process of MBS, the calculated scattered field at the top slice, $\boldsymbol{\psi}_{M}(z_{0})$, is projected onto the camera plane using a detection operator $\hat{D}$, as illustrated in Fig. 1(c). The detection operator $\hat{D}$ models the imaging process in the microscopy system to simulate the captured intensity $I_{\text{s,i}}$. First, $\boldsymbol{\psi}_{M}(z_{0})$ is numerically propagated to the focal plane of the objective lens and modulated by its pupil function, which ideally acts as a hard aperture with a radius defined by the NA. Finally, the intensities of each polarized component are summed to produce the simulated intensity $I_{\text{s,i}}$ (see details in Section 3 in Supplement 1).

To solve the inverse scattering problem, we define a loss function $\mathcal{L}_{\text{i}}$ that quantifies the difference between the simulated intensity $I_{\text{s,i}}$ and the measured intensity $I_{\text{m}}$ for each LED illumination:

$$\mathcal{L}_{\text{i}}=\left|\left|I_{\text{s,i}}-I_{\text{m,i}}\right|\right|_{2}^{2},$$

(7)

where $\parallel\cdot\parallel_{2}$ denotes the $\ell_{2}$-norm. Since MBS is a nonlinear forward model, we minimize $\mathcal{L}_{\text{i}}$ using a gradient descent approach.

In Wirtinger calculus [51], the gradient $\nabla_{n^{*}}\mathcal{L}_{\text{i}}$ is computed to update the RI distribution in order to minimize the real-valued loss function $\mathcal{L}_{\text{i}}$. To reduce memory costs associated with storing intermediate fields, we employ AM rather than traditional chain-rule-based methods, which has been applied in areas like nanophotonic inverse design [22] and DT [23]. The gradient with AM in the discrete form can be expressed as

$$\nabla_{n^{*}}\mathcal{L}_{\text{i}}=\left[2k^{2}_{0}\text{diag}(\boldsymbol{n})\cdot\boldsymbol{\psi}_{\mathrm{a}}^{*}\odot\boldsymbol{\psi}_{\mathrm{M}}\right]^{\dagger},$$

(8)

where $\boldsymbol{\psi}_{\mathrm{a}}$ is the solution to the adjoint simulation, given by

$$\nabla\times\nabla\times\boldsymbol{\psi}-k_{0}^{2}\text{diag}(\boldsymbol{\epsilon}^{*})\boldsymbol{\psi}=\left(\frac{\partial\mathcal{L}_{\text{i}}}{\partial\boldsymbol{\psi}_{\mathrm{M}}}\right)^{\dagger},$$

(9)

This adjoint simulation of MBS (see details in Section 4 in Supplement 1) only requires the field at the $z_{0}$-plane to define $\mathcal{L}_{\text{i}}$, making $(\partial\mathcal{L}_{\text{i}}/\partial\boldsymbol{\psi}_{\mathrm{M}})^{\dagger}\in\mathbb{C}^{N\times 3}$ a plane-shaped source located at the $z_{0}$-plane. If Bloch BCs are used in the forward simulation, the Bloch wave vector must be reversed in the adjoint simulation.

Since AM is derived from the Helmholtz equation (Eq. (1)), Eq. (8) serves as a universal formula for calculating gradients across all forward models. Moreover, Eq. (8) defines the minimal information required to calculate the gradient in a scattering problem. While gradients in MBS can be computed using chain-rule-based backpropagation, this approach requires storing all intermediate fields generated during each MBS iteration, resulting in significant memory overhead, especially given MBS’s volumetric iterative nature. In contrast, using AM for gradient calculation reduces memory demands, as only the final scattering field $\boldsymbol{\psi}_{\mathrm{M}}$ needs to be retained. This leads to substantial memory savings.

For validation, we use AM to compute the gradient in a 2D simulation. In the simulation, a 2D circular scatterer, with a radius of 1.5 $\mu$m and a RI of 1.361, is placed 5 $\mu$m above the PEC boundary within a background RI $n_{b}=1.336$. The scatterer is illuminated by a circularly polarized light with a wavelength of 632 nm and an incident angle of $40^{\circ}$. A pupil with a NA of 0.95 in air is used to capture the scattered field. The simulation domain has a size of $L=10$ $\mu$m, and is discretized into 20 nm pixels to simulate both the forward and backward process. For comparison, the gradient is also calculated using the chain rule, implemented via PyTorch’s automatic differentiation module (see details in Section 4 in Supplement 1). Both simulations were run on an Nvidia RTX 4090 24 GB graphics card to compare memory and time consumption. The gradient results indicate that the gradients calculated using AM closely align with those from the chain rule. However, the difference in computational cost is significant: AM requires only 504 MB of memory and completes in 157 ms, while the chain-rule-based method, which must store all intermediate fields, consumes 20.29 GB and takes 527 ms.

Before performing reconstruction tasks, we first analyze the information extracted by reflection-mode DT through interpreting the gradient derived from the loss function. In real space, the gradient reveals lateral stripes with strong oscillations along the axial direction. These axial oscillations correspond to two bright regions in the high axial-frequency domain around $k_{z}=\pm 2k_{b}$ in the Fourier spectrum, as shown in Fig. 1(c).

To guide our interpretation, we plot the Fourier support based on the single-scattering model [37], with orange and blue masks representing FS and BS contributions, respectively. Although our model includes multiple scattering, which includes higher-order effects, the alignment of the high axial frequency components with the BS region predicted by the single-scattering model indicates that the observed axial oscillations originate from captured backscattered signals. Notably, our multiple scattering model also predicts frequency components extending beyond the single-scattering region, as seen in Fig. 1(c). In contrast, FS corresponds to the slowly varying components in the gradient, aligning with the low-frequency region predicted by the single-scattering model, while multiple scattering also introduces additional frequency components extending beyond this region.

We first validate our reflection-mode DT technique in simulation by using gradient descent to reconstruct the 3D RI distribution from simulated intensity measurements. In our simulations, we used a light with a wavelength of 632 nm and a pupil with a finite NA of 0.95. Circular polarization was used in both simulation and experiments to mimic the isotropic and unpolarized illumination from the LED.

The reflection-mode DT was first tested on beads, with ground truth and reconstruction shown in the top and bottom panels of Fig. 3, respectively. Since both FS and BS features are present in the gradient (Fig. 1(c)), the reconstructed 3D RI contains both slowly varying features and lateral dense stripes in its $y$-$z$ profile. We perform the Fourier transform on the recovered RI distribution to show the 3D Fourier spectrum of the reconstructed FS and BS components. In the $k_{y}$-$k_{z}$ cross-section of the 3D Fourier spectrum, the FS and BS components are naturally separated along the $k_{z}$ direction, concentrating in the low- and high-frequency axial regions. These regions correspond to the FS and BS supports predicted by the single-scattering model, outlined by white curves, with identical lateral and axial bandwidths. To analyze the information embedded in these two components, they are separately extracted from the Fourier spectrum, followed by an inverse Fourier transform to visualize their real space distributions in the right panel of Fig. 3.

We first analyze the reconstructed low-frequency axial component. The ground truth of the low-frequency RI distribution ($\text{RI}_{\text{LF}}$) is obtained by isolating the region corresponding to the FS from the 3D Fourier spectrum. In extracting $\text{RI}_{\text{LF}}$ from the reconstructed RI, a low-pass filter along the axial dimension is applied in its 3D Fourier spectrum to isolate only the low-frequency components around the FS support. This filter slightly extends beyond the single-scattering support to account for extra contributions from multiple scattering (see details in Section 6 in Supplement 1). The extracted $\text{RI}_{\text{LF}}$ from the reconstruction shows a smooth distribution, with the dense stripe features removed. The $y$-$z$ profile of $\text{RI}_{\text{LF}}$ illustrates similar reconstruction results as the conventional transmission-mode DT and shows good agreement with the ground truth. Due to the missing cone problem in FS, the boundaries of the reconstructed $\text{RI}_{\text{LF}}$ are sharp and distinctly resolved in the lateral direction but elongated and blurred in the axial direction. To evaluate the reconstruction accuracy after incorporating multiple scattering through MBS, the reconstruction is also performed using the single-scattering model (1st Born approximation) for comparison, as illustrated in Section 6 of Supplement 1. The single-scattering model performs poorly in reconstructing the RI, with particularly degraded performance in samples with high RI contrast and complex 3D geometries.

Similar to the $\text{RI}_{\text{LF}}$, the high-frequency RI distribution ($\text{RI}_{\text{HF}}$) in the ground truth is obtained by isolating the BS region in its 3D Fourier spectrum. The $\text{RI}_{\text{HF}}$ in the reconstruction is extracted by applying a high-pass filter at the BS support to isolate the recovered high-frequency components. High-frequency oscillations with an approximate frequency of $2n_{b}/\lambda$ along the $z$-direction appear in both the real and imaginary parts of $\text{RI}_{\text{HF}}$ due to its shift from the zero frequency. Since our method provides the complex-valued reconstruction of $\text{RI}_{\text{HF}}$, the envelope of $\text{RI}_{\text{HF}}$ along $z$ can be obtained by calculating its absolute value, as shown in Fig. 3. In both the ground truth and reconstruction of $\text{RI}_{\text{HF}}$, strong features are present near the top and bottom surfaces of the beads, while they diminish inside the beads (see details in Section 6 in Supplement 1). This result suggests that the recovered information from BS is more sensitive to the lateral interfaces (with a normal vector along the $z$-direction) of the scatterers.

Our experimental setup is based on our previously developed reflection-mode LED microscope [33], as shown in Fig. 4(a). A 10$\times$/0.28 NA (Mitutoyo Plan Apo Infinity Corrected Long WD) objective lens collects the scattered signals. For oblique illumination, a 25-LED array (Kingbright) is positioned at the focal plane of a 4-$f$ system and relayed to the back focal plane of the objective. The LED array consists of a central LED and two concentric rings, with the outermost ring matching the objective NA, providing illumination NAs of 0.14 and 0.28 for the rings. Each LED illuminates the sample as a plane wave with a central wavelength of 632 nm and a 20 nm bandwidth. The intensity is captured by a camera (Imaging Source, DMK38UX541, 2.74 $\mu$m pixel size). The illumination angle is calibrated using the method from [34, 33].