Paraxial Theory of Phasor-Field Imaging

TL;DR

This paper provides a formal paraxial wave optics analysis of phasor-field imaging, deriving propagation primitives that extend its capabilities for non-line-of-sight imaging scenarios.

Contribution

It introduces a set of propagation primitives based on the two-frequency, spatial Wigner distribution for enhanced phasor-field imaging analysis.

Findings

Primitive analysis of occluded and unoccluded geometries

Application of primitives to practical imaging scenarios

Insights into light modulation effects on imaging

Abstract

The phasor field has been shown to be a valuable tool for non-line-of-sight imaging. We present a formal analysis of phasor-field imaging using paraxial wave optics. Then, we derive a set of propagation primitives---using the two-frequency, spatial Wigner distribution---that extend the purview of phasor-field imaging. We use these primitives to analyze a set of simple imaging scenarios involving occluded and unoccluded geometries with modulated and unmodulated light. These scenarios demonstrate how to apply the primitives in practice and reveal what kind of insights can be expected from them.

Equations

$${\mathcal{E}}_0({\mathbf{\rho }}_0,\omega )\equiv \int \mathrm{d}t{E}_0({\mathbf{\rho }}_0,t)e^{i\omega t},$$

$${\mathcal{E}}_0^{\prime }({\mathbf{\rho }}_0,\omega )={\mathcal{E}}_0({\mathbf{\rho }}_0,\omega )\exp [i({\omega }_0+\omega )h_0({\mathbf{\rho }}_0)/c],$$

$${\mathcal{E}}_1({\mathbf{\rho }}_1,\omega )=\int {\mathrm{d}}^2{\mathbf{\rho }}_0{\mathcal{E}}_0^{\prime }({\mathbf{\rho }}_0,\omega )\frac{\exp [i({\omega }_0+\omega )L_1/c+i({\omega }_0+\omega )|{\mathbf{\rho }}_1-{\mathbf{\rho }}_0{|}^2/2c{L}_1]({\omega }_0+\omega )}{i2\pi c{L}_1},$$

$${\mathcal{E}}_1^{\prime }({\mathbf{\rho }}_1,\omega )={\mathcal{E}}_1({\mathbf{\rho }}_1,\omega )\sqrt{T({\mathbf{\rho }}_1)}\exp [i({\omega }_0+\omega )h_1({\mathbf{\rho }}_1)/c],$$

$${\mathcal{E}}_2({\mathbf{\rho }}_2,\omega )=\int {\mathrm{d}}^2{\mathbf{\rho }}_1{\mathcal{E}}_1^{\prime }({\mathbf{\rho }}_1,\omega )\frac{\exp [i({\omega }_0+\omega )L_2/c+i({\omega }_0+\omega )|{\mathbf{\rho }}_2-{\mathbf{\rho }}_1{|}^2/2c{L}_2]({\omega }_0+\omega )}{i2\pi c{L}_2},$$

$${\mathcal{E}}_2^{\prime }({\mathbf{\rho }}_2,\omega )={\mathcal{E}}_2({\mathbf{\rho }}_2,\omega )\exp [i({\omega }_0+\omega )h_2({\mathbf{\rho }}_2)/c],$$

$$⟨I_1({\mathbf{\rho }}_1,t)⟩$$

$$=\int \frac{\mathrm{d}{\omega }_{-}}{2\pi }\left[\int \frac{\mathrm{d}{\omega }_{+}}{2\pi }⟨{\mathcal{E}}_1({\mathbf{\rho }}_1,{\omega }_{+}+{\omega }_{-}/2){\mathcal{E}}_1^{\ast }({\mathbf{\rho }}_1,{\omega }_{+}-{\omega }_{-}/2)⟩\right]e^{-i{\omega }_{-}t}$$

$$=\int \frac{\mathrm{d}{\omega }_{-}}{2\pi }{\mathcal{P}}_1({\mathbf{\rho }}_1,{\omega }_{-})e^{-i{\omega }_{-}t},$$

$${\mathcal{P}}_1({\mathbf{\rho }}_1,{\omega }_{-})=\int {\mathrm{d}}^2{\mathbf{\rho }}_0\int {\mathrm{d}}^2{\mathbf{\rho }}_0^{\prime }\int \frac{\mathrm{d}{\omega }_{+}}{2\pi }{\mathcal{E}}_0({\mathbf{\rho }}_0,\omega ){\mathcal{E}}_0^{\ast }({\mathbf{\rho }}_0^{\prime },{\omega }^{\prime })⟨e^{i[({\omega }_0+\omega )h_0({\mathbf{\rho }}_0)-({\omega }_0+{\omega }^{\prime })h_0({\mathbf{\rho }}_0^{\prime })]/c}⟩$$

$$\times ({\omega }_0+\omega )({\omega }_0+{\omega }^{\prime })e^{i(\omega -{\omega }^{\prime })L_1/c+i[({\omega }_0+\omega )|{\mathbf{\rho }}_1-{\mathbf{\rho }}_0{|}^2-({\omega }_0+{\omega }^{\prime })|{\mathbf{\rho }}_1-{\mathbf{\rho }}_0^{\prime }{|}^2]/2c{L}_1}/(2\pi c{L}_1{)}^2,$$

$${\mathcal{P}}_1({\mathbf{\rho }}_1,{\omega }_{-})$$

$$\times e^{i(\omega -{\omega }^{\prime })L_1/c+i[({\omega }_0+\omega )|{\mathbf{\rho }}_1-{\mathbf{\rho }}_0{|}^2-({\omega }_0+{\omega }^{\prime })|{\mathbf{\rho }}_1-{\mathbf{\rho }}_0^{\prime }{|}^2]/2c{L}_1}.$$

$$⟨e^{i{\omega }_0[h_0({\mathbf{\rho }}_0)-h_0({\mathbf{\rho }}_0^{\prime })]/c}⟩=\exp \{-{\omega }_0^2[{\sigma }_h^2-K_h(|{\mathbf{\rho }}_0-{\mathbf{\rho }}_0^{\prime }|)]/c^2\}.$$

$$⟨e^{i{\omega }_0[h_0({\mathbf{\rho }}_0)-h_0({\mathbf{\rho }}_0^{\prime })]/c}⟩\approx {\lambda }_0^2\delta ({\mathbf{\rho }}_0-{\mathbf{\rho }}_0^{\prime }),$$

$${\mathcal{P}}_1({\mathbf{\rho }}_1,{\omega }_{-})$$

$$=\int {\mathrm{d}}^2{\mathbf{\rho }}_0{\mathcal{P}}_0({\mathbf{\rho }}_0,{\omega }_{-})e^{i{\omega }_{-}{L}_1/c+i{\omega }_{-}|{\mathbf{\rho }}_1-{\mathbf{\rho }}_0{|}^2/2c{L}_1}/L_1^2.$$

$${\mathcal{P}}_0({\mathbf{\rho }}_0,{\omega }_{-})=\int \frac{\mathrm{d}{\omega }_{+}}{2\pi }{\mathcal{E}}_0({\mathbf{\rho }}_0,{\omega }_{+}+{\omega }_{-}/2){\mathcal{E}}_0^{\ast }({\mathbf{\rho }}_0,{\omega }_{+}-{\omega }_{-}/2),$$

$$⟨I_1({\mathbf{\rho }}_1,t)⟩=\int {\mathrm{d}}^2{\mathbf{\rho }}_0{I}_0({\mathbf{\rho }}_0,t-L_1/c-|{\mathbf{\rho }}_1-{\mathbf{\rho }}_0{|}^2/2c{L}_1)/L_1^2,$$

$$\frac{\displaystyle\exp\!\left[i\omega\sqrt{L_{1}^{2}+|\bm{\rho}_{1}-\bm{\rho}_{0}|^{2}}/c\right]}{\displaystyle\sqrt{L_{1}^{2}+|\bm{\rho}_{1}-\bm{\rho}_{0}|^{2}}}\approx\frac{\displaystyle\exp(i\omega L_{1}/c+i\omega|\bm{\rho}_{1}-\bm{\rho}_{0}|^{2}/2cL_{1})}{\displaystyle L_{1}},\mbox{ for $|\omega|\leq\Delta\omega$}$$

$${\mathcal{P}}_2({\mathbf{\rho }}_2,{\omega }_{-})$$

$$=\int {\mathrm{d}}^2{\mathbf{\rho }}_1{\mathcal{P}}_1({\mathbf{\rho }}_1,{\omega }_{-})T({\mathbf{\rho }}_1)\exp (i{\omega }_{-}{L}_2/c+i{\omega }_{-}|{\mathbf{\rho }}_2-{\mathbf{\rho }}_1{|}^2/2c{L}_2)/L_2^2,$$

$${\mathcal{P}}_2({\mathbf{\rho }}_2,{\omega }_{-})$$

$$\times T({\mathbf{\rho }}_1)\exp (i{\omega }_{-}{L}_2/c+i{\omega }_{-}|{\mathbf{\rho }}_2-{\mathbf{\rho }}_1{|}^2/2c{L}_2)/L_2^2.$$

$${\mathcal{P}}_2({\mathbf{\rho }}_2,0)=\int {\mathrm{d}}^2{\mathbf{\rho }}_{1}T({\mathbf{\rho }}_1)\int {\mathrm{d}}^2{\mathbf{\rho }}_0{\mathcal{P}}_0({\mathbf{\rho }}_0,0)/(L_1{L}_2{)}^2,$$

$${\mathcal{E}}_{\mathrm{im}}({\mathbf{\rho }}_{\mathrm{im}},\omega )$$

$$\times \int {\mathrm{d}}^2{\mathbf{\rho }}_2{\mathcal{E}}_2^{\prime }({\mathbf{\rho }}_2,\omega )\frac{e^{i({\omega }_0+\omega )L_3/c+i({\omega }_0+\omega )|{\mathbf{\rho }}_3-{\mathbf{\rho }}_2{|}^2/2c{L}_3}}{i{\lambda }_0{L}_3}$$

$$=e^{i({\omega }_0+\omega )|{\mathbf{\rho }}_{\mathrm{im}}{|}^2/2c{L}_{\mathrm{im}}}\int {\mathrm{d}}^2{\mathbf{\rho }}_2{\mathcal{E}}_2^{\prime }({\mathbf{\rho }}_2,\omega )\frac{e^{i({\omega }_0+\omega )(L_3+L_{\mathrm{im}})/c+i({\omega }_0+\omega )|{\mathbf{\rho }}_2{|}^2/2c{L}_3}}{i{\lambda }_0{L}_3}$$

$$\times \int {\text{displaylimits}}_{|{\mathbf{\rho }}_3|\le D/2}{\mathrm{d}}^2{\mathbf{\rho }}_3\frac{e^{-i(\omega +{\omega }_0){\mathbf{\rho }}_3\cdot ({\mathbf{\rho }}_2/L_3+{\mathbf{\rho }}_{\mathrm{im}}/L_{\mathrm{im}})/c}}{i{\lambda }_0{L}_{\mathrm{im}}}.$$

$${\mathcal{E}}_{\mathrm{im}}({\mathbf{\rho }}_{\mathrm{im}},\omega )=e^{i({\omega }_0+\omega )|{\mathbf{\rho }}_{\mathrm{im}}{|}^2/2c{L}_{\mathrm{im}}}$$

$$⟨{\mathcal{E}}_2^{\prime }({\mathbf{\rho }}_2,\omega ){\mathcal{E}}_2^{\prime \ast }({\mathbf{\rho }}_2^{\prime },{\omega }^{\prime })⟩\approx {\lambda }_0^2⟨{\mathcal{E}}_2({\mathbf{\rho }}_2,\omega ){\mathcal{E}}_2^{\ast }({\mathbf{\rho }}_2,{\omega }^{\prime })⟩\delta ({\mathbf{\rho }}_2-{\mathbf{\rho }}_2^{\prime }),$$

$${\mathcal{P}}_{\mathrm{im}}$$

$$\times e^{i{\omega }_{-}(L_3+L_{\mathrm{im}})/c+i{\omega }_{-}|{\mathbf{\rho }}_2{|}^2/2c{L}_3+i{\omega }_{-}|{\mathbf{\rho }}_{\mathrm{im}}{|}^2/2c{L}_{\mathrm{im}}}{\left[\frac{\pi D^2}{4{\lambda }_0{L}_3{L}_{\mathrm{im}}}\frac{J_1\left(\frac{\pi D}{{\lambda }_0}\left|\frac{{\mathbf{\rho }}_2}{L_3}+\frac{{\mathbf{\rho }}_{\mathrm{im}}}{L_{\mathrm{im}}}\right|\right)}{\frac{\pi D}{2{\lambda }_0}\left|\frac{{\mathbf{\rho }}_2}{L_3}+\frac{{\mathbf{\rho }}_{\mathrm{im}}}{L_{\mathrm{im}}}\right|}\right]}^2.$$

$$⟨I_{\mathrm{im}}({\mathbf{\rho }}_{\mathrm{im}},t)⟩$$

$$\times {\left[\frac{\pi D^2}{4{\lambda }_0{L}_3{L}_{\mathrm{im}}}\frac{J_1\left(\frac{\pi D}{{\lambda }_0}\left|\frac{{\mathbf{\rho }}_2}{L_3}+\frac{{\mathbf{\rho }}_{\mathrm{im}}}{L_{\mathrm{im}}}\right|\right)}{\frac{\pi D}{2{\lambda }_0}\left|\frac{{\mathbf{\rho }}_2}{L_3}+\frac{{\mathbf{\rho }}_{\mathrm{im}}}{L_{\mathrm{im}}}\right|}\right]}^2.$$

$$E_{0}(\bm{\rho}_{0},t)=\left\{\begin{array}[]{ll}\sqrt{\frac{\displaystyle 8P_{0}}{\displaystyle\pi d^{2}}}\,e^{-4|\bm{\rho}_{0}|^{2}/d^{2}}\cos(\Delta\omega t/2),&\mbox{for $|t|\leq t_{0}/2$,}\\[7.22743pt] 0,&\mbox{otherwise,}\end{array}\right.$$

Full text

I Introduction

Non-line-of-sight (NLoS) imaging, colloquially known as imaging around corners, is an important and growing area of research in the imaging community. Kirmani et al. Kirmani2011 introduced the concept of transient NLoS imaging by using short pulses and time-resolved detection together with multipath analysis to recover the geometry of simple, occluded scenes. Their approach was independent of bidirectional reflectance distribution function (BRDF) and albedo, and they demonstrated its experimental feasibility. Velten et al. Velten2012 revisited the problem, focusing on the case of diffuse reflection, using ultrafast streak cameras and computational backprojection. With these more powerful and developed tools, they were able to demonstrate human-identifiable reconstructions of relatively detailed geometry from around a corner. A major obstacle to applying Velten et al.’s approach in practice is the relative expense of their advanced equipment. This barrier was addressed by Heide et al. 6909808 who applied similar techniques with success to data collected by relatively inexpensive photonic-mixer-device (PMD) time-of-flight sensors. Buttafava et al. Buttafava:15 also improved upon the practical feasibility—bearing in mind cost, power, size, etc.— of implementing these approaches by demonstrating NLoS imaging with single-photon avalanche diode (SPAD) detectors. Whereas all of this work had focused on static geometry reconstruction, Gariepy et al. Gariepy2015 extended these techniques using SPAD detectors to detect motion and track moving objects around corners. With an awareness of the depth of the preceding work, Kadambi et al. kadambi provided a unified theoretical framework for the problem of occluded geometry reconstruction and motion tracking, including an analysis of expected performance and a consideration of commercially available equipment. They also generalized their theory to deal with imaging through diffusers, in addition to the around-the-corner scenario, and offered experimental demonstration of the effectiveness of their framework. Pointing out that the experimentally collected data in the previous literature had quality and resource issues owing to experimental practicalities, Klein et al. doi:10.1117/12.2241179 developed a simulation engine fit for thinking more broadly about NLoS imaging tasks without the limitations of real data. Additionally, leveraging their newfound ability to quickly simulate NLoS scenarios, they developed and demonstrated a new simulation-based inversion technique as an alternative to the computational backprojection methods that had been used in most of the prior work. Making further improvements in the area of reconstruction techniques and coping with practical resource limitations, O’Toole et al. O'Toole2018 demonstrated a confocal NLoS imaging system which facilitated the development and use of a closed-form inversion formula.

With the goal of further advancing the field of NLoS imaging, Reza et al. Reza2018 recently introduced the phasor-field ($\mathcal{P}$-field) representation for light transport that involves diffuse reflection (such as occurs in NLoS imaging) or diffuse transmission. Attempting to apply their light transport model to NLoS geometries that include intermediate occluding objects or non-Lambertian reflections will reveal that the $\mathcal{P}$ field is an insufficient representation of the underlying field at the site of such features. Nevertheless, Liu et al. Liu2018 used the $\mathcal{P}$-field approach to propose and demonstrate that line-of-sight imaging techniques can be fruitfully applied, in a computational manner, to NLoS operation, even in the presence of intermediate occluders and non-Lambertian reflections. In doing so, they presented what may be the most robust and detailed reconstructions of NLoS scenes to date. Their success in this endeavor is due to their development of reconstruction techniques that obviate the need for a full light transport model by relying on there being initial and final Lambertian reflections. These techniques are fortunately, and somewhat surprisingly, not burdened by the limitations inherent in applying $\mathcal{P}$-field propagation to scenarios more general than purely Lambertian reflections. Very recently, Reza et al. Reza2019 reported an elegant series of experiments that verify the $\mathcal{P}$ field’s legitimacy. These experiments clearly demonstrate the $\mathcal{P}$ field’s wave-like properties, which offer the possibility of NLoS imaging without the need for computational reconstructions by using a $\mathcal{P}$-field lens instead.

The success of Liu et al.’s experiments is impressive, and Reza et al.’s $\mathcal{P}$-field lens is quite promising. However, we believe that even greater performance might be possible if afforded a complete transport model that can account for all features that might be encountered in NLoS imaging. At the very least, such a transport model would facilitate anticipatory preparation and analysis for particular scenarios of interest. The argument could be made that the propagation rules for the optical-frequency field—not those for the $\mathcal{P}$ field—already provide such a transport model, but the aforementioned works have demonstrated the intuitive utility of the $\mathcal{P}$-field approach. Consequently, we believe it is worthwhile to pursue propagation primitives that can readily establish the $\mathcal{P}$-field input-output relation for the initial and final Lambertian reflections when occluders and non-Lambertian reflectors are present in the intervening space.

In this paper, we develop a set of propagation primitives that extend the $\mathcal{P}$-field formalism to scenarios that go beyond what was considered in Reza2018 by Reza, et al. For convenience, we assume a transmissive geometry (without reflections) that is an unfolded proxy for occlusion-aided, three-bounce NLoS imaging Xu2018 ; Thrampoulidis2018 and use scalar-wave, paraxial optics although these restrictions are not essential. In Sec. II we present our own development and analysis of the $\mathcal{P}$-field notion. We begin by tracing light propagation through an example transmissive geometry wherein a natural definition for the $\mathcal{P}$ field presents itself. Continuing this analysis, we arrive at a paraxial $\mathcal{P}$-field propagator analogous to that reported by Reza et al. Reza2018 . Using this result, we analyze the performance of $\mathcal{P}$-field imaging for unoccluded transmissive geometries. Next, moving beyond the $\mathcal{P}$ field, in Sec. III we introduce the two-frequency spatial Wigner distribution and present primitives for its propagation through a diffuser, through a deterministic occluder, through a specular-plus-diffuser mask, and through Fresnel diffraction. With these primitives, we then derive the $\mathcal{P}$-field input-output relation for occlusion-aided, diffuse-object, transmissive imaging. With that analysis in hand, we compare the $\mathcal{P}$-field point-spread function for diffuse-object imaging using modulated light in the absence of an occluder with those for diffuse-object imaging using unmodulated light that is aided by the presence of either a Gaussian-pinhole occluder or a Gaussian-pinspeck occluder. Finally, in Sec. IV we summarize our results and consider directions for further research.

II $\bm{\mathcal{P}}$-Field Propagation and Imaging

In this section we consider electromagnetic field propagation through a paraxial, transmissive geometry that serves as a surrogate for an around-the-corner imaging configuration. As was done by Reza et al. Reza2018 , we define the $\mathcal{P}$ field as the Fourier transform of the short-time average irradiance. Using this definition, we derive a formula for paraxial propagation of the $\mathcal{P}$ field, which we find to be similar to the traditional Fresnel-diffraction formula for the propagation of the electromagnetic field, as reported by Reza et al. in Reza2018 . We then apply this understanding of the $\mathcal{P}$ field to the task of imaging through diffusers and analyze the associated performance.

II.1 Setup for Paraxial Propagation through Multiple Diffusers

Figure 1 shows the transmissive geometry we shall address in this paper for $\mathcal{P}$-field propagation within the paraxial regime, i.e., wherein Fresnel diffraction applies. Here, $E_{0}(\bm{\rho}_{0},t)$ is the baseband, complex-field envelope for a quasimonochromatic, scalar-wave, modulated laser field entering the $z=0$ plane, expressed as a function of the transverse spatial coordinates, $\bm{\rho}_{0}=(x_{0},y_{0})$, and time, $t$. This field has center frequency $\omega_{0}$ and bandwidth $\Delta\omega\ll\omega_{0}$, so that the optical-frequency field is ${\rm Re}[E_{0}(\bm{\rho}_{0},t)e^{-i\omega_{o}t}]$. Its units are $\sqrt{{\rm W/m}^{2}}$, making $I_{0}(\bm{\rho}_{0},t)=|E_{0}(\bm{\rho}_{0},t)|^{2}$ the short-time average irradiance footnote1 illuminating the $z=0$ plane. It will be assumed, in all that follows, that $\Delta\omega$ is such that available photodetectors can fully resolve the time dependence of $I_{0}(\bm{\rho}_{0},t)$. As soon will be seen, it will be valuable to employ the time-domain Fourier transform of $E_{0}(\bm{\rho}_{0},t)$, viz. limits ,

[TABLE]

for use analyzing the Fig. 1 configuration.

After propagating through the thin diffuser $h_{0}(\bm{\rho}_{0})$, the Fourier-domain field at $z=0_{+}$ is

[TABLE]

where $c$ is light speed and we have normalized away the diffuser’s refractive index. Physically, we are modeling this diffuser as a space-dependent $h_{0}(\bm{\rho}_{0})/c$ time delay. Because it is unreasonable to presume we can accurately account for this delay as a deterministic quantity, we shall suppress its average value—across an ensemble of statistically identical diffusers—and consider $h_{0}(\bm{\rho}_{0})$ to be a zero-mean, homogeneous, isotropic, Gaussian random function of $\bm{\rho}_{0}$, with covariance function $K_{h}(|\Delta\bm{\rho}|)=\langle h_{0}(\bm{\rho}_{0}+\Delta\bm{\rho})h_{0}(\bm{\rho}_{0})\rangle$, where angle brackets denote ensemble average. Moreover, in keeping with $h_{0}(\bm{\rho}_{0})$’s being a diffuser, we shall take its standard deviation, $\sigma_{h}=\sqrt{K_{h}(0)}$ to be much greater than the center wavelength, $\lambda_{0}=2\pi c/\omega_{0}$, and its coherence length $\rho_{c}$—the transverse distance beyond which $K_{h}(|\Delta\bm{\rho}|)$ vanishes—to be at most a few $\lambda_{0}$. Furthermore—and this condition is essential to there being a useful $\mathcal{P}$-field propagator—we shall assume that $\sigma_{h}$ is much smaller than the wavelength of the modulation bandwidth, $\Delta\lambda=2\pi c/\Delta\omega$.

Within the paraxial (Fresnel-diffraction) propagation regime we have that

[TABLE]

is the time-domain Fourier transform of $E_{1}(\bm{\rho}_{1},t)$, the field illuminating the $z=L_{1}$ plane. This illumination results in

[TABLE]

being the time-domain Fourier transform of $E^{\prime}_{1}(\bm{\rho}_{1},t)$, the field that emerges at $z=L_{1_{+}}$, after propagation through a deterministic thin transmission screen with intensity transmission pattern $T(\bm{\rho}_{1})$, and a thin diffuser, $h_{1}(\bm{\rho}_{1})$, that we will take to be statistically independent of, but identically distributed as, $h_{0}(\bm{\rho}_{0})$.

Paraxial propagation to $z=L_{1}+L_{2}$, now gives us

[TABLE]

and propagation through the thin diffuser at $z=L_{1}+L_{2}$ results in

[TABLE]

being the time-domain Fourier transform of $E^{\prime}_{2}(\bm{\rho}_{2},t)$, the field that emerges at $z=(L_{1}+L_{2})_{+}$. We will assume that $h_{2}(\bm{\rho}_{2})$ is statistically independent of, but identically distributed as, $h_{0}(\bm{\rho}_{0})$ and $h_{1}(\bm{\rho}_{1})$.

Before proceeding further, let us briefly comment on how the Fig. 1 geometry relates to three-bounce NLoS active imaging. The $z=0$ diffuser, which is illuminated by modulated laser light, represents a Lambertian-reflecting visible wall with a uniform albedo. The combination of the intensity transmission pattern $T(\bm{\rho}_{1})$ and the $z=L_{1}$ diffuser represent a Lambertian-reflecting hidden wall with spatially-varying albedo $T(\bm{\rho}_{1})$. The $z=L_{1}+L_{2}$ diffuser represents a second Lambertian reflection at the visible wall, where statistical independence from the first visible-wall reflection can be ensured by the NLoS imaging sensor’s viewing a different section of that wall than what the laser illuminates. The goal of three-bounce NLoS active imaging in this setting is to use the third-bounce light returned from the visible wall to reconstruct the hidden wall’s albedo $T(\bm{\rho}_{1})$. In the next section, we will derive the $\mathcal{P}$-field propagator for the preceding transmission geometry.

II.2 $\mathcal{P}$-Field Propagator in the Paraxial Regime

To start our derivation, consider $\langle I_{1}(\bm{\rho}_{1},t)\rangle$, where $I_{1}(\bm{\rho}_{1},t)\equiv|E_{1}(\bm{\rho}_{1},t)|^{2}$ is the short-time average irradiance illuminating the $z=L_{1}$ plane and angle brackets denote averaging over the statistics of $h_{0}(\bm{\rho}_{0})$. Going to the temporal-frequency domain, we have that

[TABLE]

where ∗ denotes complex conjugate, $\omega_{+}\equiv(\omega+\omega^{\prime})/2$, $\omega_{-}\equiv\omega-\omega^{\prime}$, and we have introduced the $\mathcal{P}$ field at the $z=L_{1}$ plane as the Fourier transform of $\langle I_{1}(\bm{\rho}_{1},t)\rangle$. Next, employing Eqs. (2) and (3), we get

[TABLE]

where, as before, $\omega_{+}\equiv(\omega+\omega^{\prime})/2$ and $\omega_{-}\equiv\omega-\omega^{\prime}$. Because $\Delta\omega\ll\omega_{0}$ and $\sigma_{h}\ll\Delta\lambda$, the preceding result can be reduced to

[TABLE]

A standard result for Gaussian random functions gives us Gaussian

[TABLE]

Then, because $\sigma_{h}\gg\lambda_{0}$ and $\rho_{c}\sim\lambda_{0}$ we can use an impulse approximation, viz.,

[TABLE]

in Eq. (11) to obtain

[TABLE]

Here, the $\mathcal{P}$ field at $z=0$ is

[TABLE]

with no averaging brackets required, because the laser illumination of the $z=0$ plane is deterministic.

Equation (15)—which coincides with the result of applying the Fresnel approximation to Reza et al.’s Rayleigh-Sommerfeld $\mathcal{P}$-field propagator Reza2018 —is our essential result for paraxial $\mathcal{P}$-field propagation over a distance $L_{1}$. It shows that the field emerging from a diffuser that imposes complete spatial incoherence at the optical frequency, but is smooth at the modulation frequency, leads to paraxial $\mathcal{P}$-field propagation at frequency $\omega_{-}$ over a distance $L_{1}$ that is governed by a modified version of the $\mathcal{E}$-field’s Fresnel-diffraction formula, viz., one in which the exponent’s optical frequency in the $\mathcal{E}$-field Fresnel formula is replaced by the $\mathcal{P}$ field’s modulation frequency and the $\mathcal{E}$-field formula’s $\omega_{0}/i2\pi cL_{1}$ factor is replaced by the $\mathcal{P}$ field’s $1/L_{1}^{2}$ factor. By inverse Fourier transformation of Eq. (15), we see that irradiance propagation from the diffuser at $z=0$ to the $z=L_{1}$ plane is governed by

[TABLE]

which has the following pleasing physical interpretation: Paraxial propagation of the short-time average irradiance from the diffuser’s output to the $z=L_{1}$ presumes that

[TABLE]

can be employed, and results in $\langle I_{1}(\bm{\rho}_{1},t)\rangle$ being governed by the paraxial form of geometric optics, viz., the differential contribution of $I_{0}(\bm{\rho}_{0},t)$ to $\langle I_{1}(\bm{\rho}_{1},t)\rangle$ is time delayed by $L_{1}/c+|\bm{\rho}_{1}-\bm{\rho}_{0}|^{2}/2cL_{1}$ and attenuated by the inverse-square-law factor $1/L_{1}^{2}$.

Paralleling the previous development, it is now easy to show that

[TABLE]

where the averaging brackets in Eq. (19) represent averaging over the $h_{0}(\bm{\rho}_{0})$ and the $h_{1}(\bm{\rho}_{1})$ ensembles. Combining this result with what we have already obtained for relating $\mathcal{P}_{1}(\bm{\rho}_{1},\omega_{-})$ to $\mathcal{P}_{0}(\bm{\rho}_{0},\omega_{-})$ we get

[TABLE]

Before continuing, it is crucial to note the behavior of $\mathcal{P}_{2}(\bm{\rho}_{2},0)$. From Eq. (21) we immediately find that

[TABLE]

indicating that there is no spatial information about $T(\bm{\rho}_{1})$ available in $\mathcal{P}_{2}(\bm{\rho}_{2},0)$. This behavior is a consequence of using the paraxial approximation. Going beyond the paraxial-propagation regime—to Rayleigh-Sommerfeld diffraction—will yield a $\mathcal{P}_{2}(\bm{\rho}_{2},0)$ containing some spatial information about $T(\bm{\rho}_{1})$, but the inverse problem for recovering $T(\bm{\rho}_{1})$ from $\mathcal{P}_{2}(\bm{\rho}_{2},0)$ will still be poorly conditioned in the Fig. 1 configuration. This behavior has been seen by Xu et al. Xu2018 and Thrampoulidis et al. Thrampoulidis2018 in their work on NLoS active imaging with pulsed illumination, in which occlusion-aided operation was needed to obtain useful albedo reconstructions when transient behavior was ignored.

II.3 $T(\bm{\rho}_{1})$ Reconstruction in the Paraxial Regime $\mathcal{P}$-Field Formalism

Equation (21) shows that the intensity transmission pattern, $T(\bm{\rho}_{1})$, we wish to reconstruct is illuminated by $\mathcal{P}_{1}(\bm{\rho}_{1},\omega_{-})$, the $\mathcal{P}$ field that results from propagation of the laser illumination’s $\mathcal{P}_{0}(\bm{\rho}_{0},\omega_{-})$ from $z=0$ to $z=L_{1}$. After transmission through $T(\bm{\rho}_{1})$ and the diffuser $h_{1}(\bm{\rho}_{1})$, $\mathcal{P}$-field propagation from to $z=L_{1}+L_{2}$ results in $\mathcal{P}_{2}(\bm{\rho}_{2},\omega_{-})$, which encounters another diffuser. Because that last diffuser will render the field emerging from it spatially incoherent, we will use the conventional thin-lens imaging system, shown in Fig. 2, to gather the data needed to reconstruct $T(\bm{\rho}_{1})$.

Let $E^{\prime}_{2}(\bm{\rho}_{2},t)$ be the baseband, complex-field envelope emerging from the diffuser in the $z=L_{1}+L_{2}$ plane, and let $\mathcal{E}^{\prime}_{2}(\bm{\rho}_{2},\omega)$ be its time-domain Fourier transform. After Fresnel propagation from $z=L_{1}+L_{2}$ to $z=L_{1}+L_{2}+L_{3}$, propagation through the diameter-$D$ circular-pupil, focal-length-$f$, thin lens, and Fresnel propagation over an additional $L_{\rm im}$ distance where $1/f=1/L_{3}+1/L_{\rm im}$, the resulting image-plane field $E_{\rm im}(\bm{\rho},t)$ has time-domain Fourier transform given by

[TABLE]

Performing the integration over $\bm{\rho}_{3}$ results in

[TABLE]

where $J_{1}(\cdot)$ is the first-order Bessel function of the first kind, and we have used $\pi D/\lambda_{0}$ in lieu of $(\omega_{0}+\omega)D/2c$ in the Airy pattern because $\Delta\omega\ll\omega_{0}$.

The presence of the diffuser $h_{2}(\bm{\rho}_{2})$ makes

[TABLE]

which together with Eq. (25) yields

[TABLE]

and hence

[TABLE]

So, by measuring $\langle I_{\rm im}(\bm{\rho}_{\rm im},t)\rangle$, i.e., the diffuser-averaged, short-time average, image-plane irradiance, we obtain a $1.22\lambda_{0}/D$-angular-resolution, image of $\langle I_{2}(\bm{\rho}_{2},t-(L_{3}+L_{\rm im})/c-|\bm{\rho}_{2}|^{2}/2cL_{3})\rangle$. From that irradiance image we can then compute a $1.22\lambda_{0}/D$-angular-resolution image of $\mathcal{P}_{2}(\bm{\rho}_{2},\omega_{-})$ at any modulation frequency of interest.

For reconstructing $T(\bm{\rho}_{1})$, let us suppose that the $z=0$ illumination is a duration $t_{0}$, cosinusoidally-modulated, collimated Gaussian-beam laser field where $\Delta\omega t_{0}\gg 1$, i.e.,

[TABLE]

with $P_{0}t_{0}/2$ being the energy illuminating the $z=0$ plane. This field’s short-time average irradiance is then

[TABLE]

which leads to

[TABLE]

and hence

[TABLE]

because $\Delta\omega t_{0}\gg 1$. Although this expression can be evaluated analytically, we shall not bother. We just note that with $\Delta\omega/2\pi\sim 1\,$GHz, $d\sim 1$ mm, and $L_{1}\sim 1$ m, we have $cL_{1}/\Delta\omega d^{2}\gg 1$ from which it follows that the spatial extent of $\mathcal{P}_{1}(\bm{\rho}_{1},\Delta\omega)$ will be \sim$$cL_{1}/\Delta\omega d\gg d. In other words, the effect of the diffuser $h_{0}(\bm{\rho}_{0})$ is to ensure that a finite, but much larger than diameter-$d$, region of the $z=L_{1}$ plane is illuminated by the frequency-$\Delta\omega$ $\mathcal{P}$ field.

To proceed further, assume we have generated the computed image,

[TABLE]

of $\mathcal{P}_{2}(\bm{\rho}_{2},\Delta\omega)$ from the $\langle I_{\rm im}(\bm{\rho}_{\rm im},t)\rangle$ measurement. We can computationally invert Eq. (20) to obtain a reconstruction of $T(\bm{\rho}_{1})\mathcal{P}_{1}(\bm{\rho}_{1},\Delta\omega)$ and use our knowledge of $\mathcal{P}_{1}(\bm{\rho}_{1},\Delta\omega)$ to obtain a $T(\bm{\rho}_{1})$ image. In particular, suppose we measure $\langle I_{\rm im}(\bm{\rho}_{\rm im},t)\rangle$ for $|\bm{\rho}_{\rm im}|\leq d_{\rm im}/2$, and then define $\tilde{T}(\tilde{\bm{\rho}}_{1})$ by

[TABLE]

where $D^{\prime}\equiv d_{\rm im}L_{3}/L_{\rm im}$. Neglecting noise, and assuming that the $1.22\lambda_{0}/D$ angular resolution is sufficient to make

[TABLE]

for $|\bm{\rho}_{2}|\leq D^{\prime}/2$, Eq. (34) leads to

[TABLE]

Thus, over the region in the $z=L_{1}$ plane wherein $|\mathcal{P}_{1}(\bm{\rho}_{1},\Delta\omega)|$ has an appreciable value, the $\mathcal{P}$-field imager using cosinusoidal $E$-field modulation at frequency $\Delta\omega/2$ achieves a spatial resolution of $1.22\Delta\lambda L_{2}/D^{\prime}$, where: $\Delta\lambda=2\pi c/\Delta\omega$; $L_{2}$ is the distance from the transparency-containing plane to the plane visible to the sensor; and $D^{\prime}=d_{\rm im}L_{3}/L_{\rm im}$, with $L_{3}$ being the distance from the plane visible to the sensor to the sensor’s entrance pupil, $L_{\rm im}$ being the distance from that entrance pupil to the image plane where irradiance measurements are made, and $d_{\rm im}$ being the diameter of the image-plane region over which those measurements are made.

III Two-Frequency Spatial Wigner Distribution and Occlusion-Aided Imaging

In this section, we consider a generalized version of our paraxial, transmissive geometry which allows for the presence of deterministic occluders in the light’s path and a more general target transmissivity mask. The $\mathcal{P}$ field alone does not suffice to track the evolution of the light through all intermediate planes of this geometry, so we go beyond this quantity to define a more comprehensive one: the two-frequency spatial Wigner distribution. We demonstrate how the two-frequency spatial Wigner distribution relates to other better-known quantities for characterizing propagation through random media and present a set of propagation primitives for it, relevant to our transmissive geometry. Finally, we use these propagation primitives to analyze occlusion-aided imaging scenarios and demonstrate that the presence of intermediate occluders has the potential to improve performance, as seen previously in Xu et al. Xu2018 and Thrampoulidis et al. Thrampoulidis2018 .

III.1 Setup for Paraxial Propagation through Multiple Diffusers with Occlusion

Figure 3 shows a generalized setup for transmissive $\mathcal{P}$-field imaging. Here, two occluders, having field-transmission functions $P(\bm{\rho}_{d})$ and $P^{\prime}(\bm{\rho}^{\prime}_{d})$, have been introduced in the $z=L_{1}-L_{d}$ and $z=L_{1}+L^{\prime}_{d}$ planes, and the $z=L_{1}$ plane contains a field-transmission mask $F(\bm{\rho}_{1})$ that has both specular and diffuse components. In the NLoS analogy, the two occluders represent objects in the hidden space—encountered by the light as it propagates towards and returns from the hidden wall, respectively—and the generalized field-transmission mask accounts for more general, non-Lambertian hidden walls. This configuration—if $F(\bm{\rho}_{1})$ is purely diffuse with a space-varying albedo that is to be imaged, i.e., equivalent to the stacked intensity-transmission mask and thin diffuser from Fig. 1—is our unfolded proxy for Xu et al.’s experiments Xu2018 .

The ultimate goal of a phasor-field transport model is to provide the short-time average irradiance at the output of some system—or equivalently, its Fourier transform: the $\mathcal{P}$ field—given the short-time average irradiance, or its associated $\mathcal{P}$ field, at the input of the system. This is possible in NLoS or diffuse transmissive-imaging scenarios—provided that the system can be summarized by a linear transformation of the underlying electromagnetic field—when the input and output facets of the systems in question are Lambertian walls (NLoS case) or diffusers (transmissive case). Such facets destroy all directionality information, viz., all spatial coherence, so that $\mathcal{P}$-fields fully characterize the light they reflect (NLoS case) or transmit (transmissive case). Free-space propagation increases spatial coherence, but provided we only care about the short-time average irradiances at input and output planes containing pure diffusers, a $\mathcal{P}$-field input-output model propagation is possible as those diffusers will, respectively, destroy the initial and propagation-created coherence. If, however, as at $z=L_{1}-L_{d}$, $z=L_{1}$, or $z=L_{1}+L^{\prime}_{d}$ in Fig. 3, we are interested in planes that do not contain pure diffusers, the $\mathcal{P}$ field is insufficient to fully characterize the electromagnetic field emerging from them. Thus, owing to what can be viewed as a lack of directionality information, the $\mathcal{P}$ field at those output planes fails to provide enough information to determine the increased spatial coherence that will accrue from subsequent free-space diffraction. Accordingly, we find the $\mathcal{P}$ field insufficient for the task of building a complete light-transport model for scenarios including occluders and specular-plus-diffuser masks. Indeed, although omitted for brevity, carrying out a Fig. 3 propagation analysis—like that done for Fig. 1—confirms that a $\mathcal{P}$-field input-output relation built up from propagating the $\mathcal{P}$ field from each plane containing an optical element to the next such plane is impossible.

To tackle these scenarios, we start from the beginning, and instead of considering the short-time average irradiance we consider a variant with directionality information—the time-dependent specific irradiance from small-angle-approximation linear transport theory Ishimaru1978a :

[TABLE]

In computer vision, this quantity is known as the 5D light field lf1 ; lf2 ; lf3 . By replacing $2\pi{\bf s}/\lambda_{0}$ with ${\bf k}$, the time-dependent specific irradiance can be seen to be a time-indexed spatial Wigner distribution, cf. the spatial Wigner distribution of a monochromatic scalar wave, viz.,

[TABLE]

which has long been recognized as a useful tool in optics, see, e.g., Walther1968 ; Bastiaans1980 ; Alonso2011 . The short-time average irradiance is obtained from $I_{z}(\bm{\rho}_{+},{\bf s},t)$ by integrating out its directionality information,

[TABLE]

and the $\mathcal{P}$ field is then obtained by time-domain Fourier transformation.

As before, we find it to convenient to carry out our analysis in the temporal-frequency domain. Paralleling the development in Eqs. (7)–(9) we have:

[TABLE]

where $\omega_{+}\equiv(\omega+\omega^{\prime})/2$ and $\omega_{-}\equiv\omega-\omega^{\prime}$ as we employed in Sec. II. The bracketed quantity in Eq. (41) is the Fourier transform of the time-dependent specific irradiance, so it contains equivalent information. Comparing to our Sec. II analysis, this quantity is the directionality-augmented analog of the $\mathcal{P}$ field, and as it turns out would be sufficient to build a transport model for the Fig. 3 scenario. Out of prudence though, having learned from the insufficient generality of the $\mathcal{P}$ field, we feel it is wise to build our Fig. 3 analysis on the quantity in parentheses within Eq. (41), the two-frequency spatial Wigner distribution (TFSWD):

[TABLE]

from which the time-dependent specific irradiance can be obtained via

[TABLE]

The merit of the TFSWD’s added generality can be seen by considering the space-time autocorrelation function,

[TABLE]

that is used in parabolic-approximation propagation theory through random media Ishimaru1978b . The time-dependent specific irradiance can be found from the space-time autocorrelation function, viz., we have that

[TABLE]

but the converse is not true, i.e., the space-time autocorrelation function cannot in general be found from knowledge of the time-dependent specific irradiance alone. However, the space-time autocorrelation function is equivalent to the TFSWD because we have that

[TABLE]

where $t_{+}\equiv(t_{1}+t_{2})/2$, $t_{-}\equiv t_{1}-t_{2}$, and

[TABLE]

For $E$-field propagation through an arbitrary linear transformation of the form

[TABLE]

the input’s space-time autocorrelation function suffices to determine the output’s space-time autocorrelation function, and hence the output-plane $\mathcal{P}$ field. Morevoer, the same must be true for the TFSWD. Because knowledge of the time-dependent specific irradiance alone does not in general determine the space-time autocorrelation function, it does not suffice to characterize second-moment propagation through an arbitrary linear transformation of the form given in Eq. (48), i.e., it cannot determine the output $\mathcal{P}$ field. For example, the time-dependent specific irradiance cannot account for propagation that involves a linear time-invariant filtering in time, e.g., through a transparency that has a frequency-dependent transmissivity. So, although this capability is not fully exploited in this paper, by building our theory around the TFSWD we are prepared to handle arbitrary linear transformations of the $E$ field, rather than just those that can be characterized by the time-dependent specific irradiance. Note that the 6D light field,

[TABLE]

would also suffice in this regard, as it is the time-domain inverse Fourier transform of the TFSWD.

The $z$-plane $\mathcal{P}$ field can be found from that plane’s TFSWD as follows:

[TABLE]

From this result we see that the TFSWD allows us to realize the goal of analyzing occluded phasor-field imaging if we can: (1) propagate $W_{\mathcal{E}_{z}}(\bm{\rho}_{+},{\bf k},\omega_{+},\omega_{-})$ through a $z$-plane field-transmission mask, whether that be a diffuser, deterministic occluder, or specular-plus-diffuser mask; and (2) propagate $W_{\mathcal{E}_{z}}(\bm{\rho}_{+},{\bf k},\omega_{+},\omega_{-})$ through a distance $L$ of Fresnel diffraction. All of these propagation calculations are done Appendix A. For convenience, we summarize these results below:

**Propagation through a diffuser:

**For propagation through a diffuser characterized by the impulse approximation in Eq. (13), we have

[TABLE]

.

**Propagation through a deterministic occluder:

**With $W_{P}(\bm{\rho}_{+},{\bf k})\equiv$ $\int\!{\rm d}^{2}\bm{\rho}_{-}\,P(\bm{\rho}_{+}+\bm{\rho}_{-}/2)P^{*}(\bm{\rho}_{+}-\bm{\rho}_{-}/2)e^{-i{\bf k}\cdot\bm{\rho}_{-}}$, we have

[TABLE]

**Propagation through a specular-plus-diffuser mask: **

With $F(\bm{\rho}_{1})$ having nonzero mean $\langle F(\bm{\rho}_{1})\rangle\neq 0$, and covariance, $\langle\Delta F(\bm{\rho}_{+}+\bm{\rho}_{-}/2)\Delta F^{*}(\bm{\rho}_{+}-\bm{\rho}_{-}/2)\rangle\approx\lambda_{0}^{2}\mathcal{F}(\bm{\rho}_{+})\delta(\bm{\rho}_{-})$ where $0\leq\mathcal{F}(\bm{\rho}_{+})\leq 1$ , we get

[TABLE]

**Fresnel diffraction:

**For Fresnel diffraction from the $z=0_{+}$ plane to the $z=L_{1}-L_{d}$ plane, we get

[TABLE]

III.2 Occlusion-Aided Imaging

In Sec. II we noted that, in the paraxial limit, unoccluded imaging configurations without modulated light are unconditioned with respect to reconstructing the target mask’s albedo. Moreover, we showed that the addition of modulation enabled reconstruction of the target mask’s albedo at a resolution limited by the bandwidth of that modulation. What remains then is to examine the unmodulated and modulated cases for occluded geometries. For clarity and convenience, we will consider a simplified version of Fig. 3 in which the first occluder is absent and the screen at $z=L_{1}$ is purely diffuse. In the NLoS analogy, this corresponds to a geometry in which a single occluding object is encountered in the hidden space only on the light’s return trip from a Lambertian hidden wall. Further convenience, without appreciable loss of generality, is afforded by our assuming that the laser light incident on the $z=0$ plane is a $+z$-going plane wave of short-time average irradiance $I_{0}(t)$, and that the distances in Fig. 3 satisfy $L_{1}=L_{2}=L$, and $L_{d}=L/2$

The TFSWD of the plane-wave laser light is easily shown to be

[TABLE]

where

[TABLE]

After the diffuser in the $z=0$ plane we get

[TABLE]

and after propagation to the $z=L$ plane, we find

[TABLE]

At $z=L_{1}$ this Wigner distribution encounters a diffuse target mask, i.e., one whose field-transmission function $F(\bm{\rho}_{1})$ has zero mean and covariance $\langle\Delta F(\bm{\rho}_{1})\Delta F^{*}(\bm{\rho}_{2})\rangle=\lambda_{0}^{2}\mathcal{F}[(\bm{\rho}_{1}+\bm{\rho}_{2})/2]\delta(\bm{\rho}_{1}-\bm{\rho}_{2})$, which results in

[TABLE]

Fresnel propagation to $z=3L/2$ now gives us

[TABLE]

and passage through the occluder in that plane leads to

[TABLE]

Fresnel propagation over another $L/2$ distance then gives

[TABLE]

from which we get

[TABLE]

Now, using

[TABLE]

and changing variables to ${\bf k}_{-}={\bf k}-{\bf k}^{\prime}$ and ${\bf k}_{+}=({\bf k}+{\bf k}^{\prime})/2$ we have

[TABLE]

where we have suppressed the $\bm{\rho}_{+}$ argument of $\mathcal{P}_{0}(\bm{\rho}_{+},\omega_{-})$ because that field has no such dependence for the plane-wave source we have assumed. We define a new function

[TABLE]

With this definition we have

[TABLE]

Changing variables again, $\tilde{\bm{\rho}}=\bm{\rho}_{+}-cL{\bf k}_{+}/\omega_{0}$, we get our final result

[TABLE]

Owing to the Fresnel-propagation kernel in Eq. (68), this result is a superposition integral with image inversion, rather than a convolution integral with image inversion.

To get to a simpler result that will afford us insight into the advantage of occlusion-aided imaging, we shall assume that the initial laser illumination is monochromatic, i.e., the optical-frequency field that illuminates the $z=0$ plane is ${\rm Re}[E_{0}(\bm{\rho}_{0})e^{-i\omega_{0}t}]$. In this unmodulated case we can use the usual spatial Wigner distribution, i.e.,

[TABLE]

of the $z=0$-plane field, in lieu of the TFSWD. The propagation primitives given earlier for the TFSWD all apply to the spatial Wigner distribution function for the unmodulated case with the only difference being that we set $\omega_{-}=0$ in the Fresnel-diffraction primitive. Paralleling the development that led to Eq. (68) assuming that $E_{0}(\bm{\rho}_{0})=\sqrt{I_{0}}$ is a constant, we get

[TABLE]

where

[TABLE]

and we have used the evanescence cutoff, $|{\bf k}|\leq 2\pi/\lambda_{0}$, to justify replacing $\int\!{\rm d}^{2}{\bf k}\,I_{0}/(2\pi)^{2}$ with $\pi I_{0}/\lambda_{0}^{2}$.

Equations (70) and (71) show that this unmodulated case offers no spatial information about $\mathcal{F}(\bm{\rho})$ in the absence of an occluder, i.e., we get $G(\bm{\rho})=\pi/L^{2}$ when $P(\bm{\rho})=1$, as seen previously in Eq. (22). To quantify the spatial information afforded by the presence of an occluder in the unmodulated scenario, we consider two simple cases: the Gaussian pinhole

[TABLE]

and the Gaussian pinspeck,

[TABLE]

where $\rho_{0}$ is the $e^{-1/2}$-attenuation radius of the Gaussian functions. The Gaussian-pinhole camera can be analyzed with far less complication than our approach to obtaining Eqs. (70) and (71), but (after accounting for image inversion) its point-spread function (psf) $G_{\rm ph}(\bm{\rho})$ is revealing. The Gaussian-pinspeck camera, on the other hand, is more relevant to the experiments of Xu et al. Xu2018 , but its psf $G_{\rm ph}(\bm{\rho})$ is more complicated. In both cases, however, the Gaussian functions involved enable us to get closed-form psf results.

For the Gaussian pinhole, we find that

[TABLE]

where $k_{0}\equiv\omega_{0}/c=2\pi/\lambda_{0}$ is the wave number at the optical frequency and $\Omega\equiv 4k_{0}\rho_{0}^{2}/L$ is the Fresnel number for the pinhole’s propagation geometry. The spatial resolution of $G_{\rm ph}(\bm{\rho})$ improves with decreasing $\rho_{0}$ when $\Omega>1$, and degrades with decreasing $\rho_{0}$ when $\Omega<1$. Thus the Gaussian pinhole’s resolution-optimized psf,

[TABLE]

is obtained when $\rho_{0}=\sqrt{L/4k_{0}}=\sqrt{\lambda_{0}L/8\pi}$. The optimized psf’s spatial resolution—taken to be its $e^{-\pi}$-attenuation radius—is then $\sqrt{\lambda_{0}L}$, which is far superior to the $1.22\Delta\lambda L/D^{\prime}$ for the unoccluded, modulated case governed by Eq. (36). For example, with $\lambda=1\,\mu$m and $L=1\,$m the optimum spatial resolution of occlusion-aided unmodulated imaging is 1 mm, while that of unoccluded modulated imaging, with $\Delta\lambda=3\,$cm ($\Delta\omega/2\pi=10\,$GHz) and $D^{\prime}=10$ cm, is 37 cm at $L=1$ m. For comparison with the Gaussian pinspeck’s psf, it is worth noting that the Gaussian pinhole’s psf maintains its Gaussian shape for all values of its Fresnel number $\Omega$, with only its overall amplitude $G_{\rm ph}({\bf 0})$ and its spatial resolution $\rho_{\rm res}(\Omega)\equiv\sqrt{4\pi(1+\Omega^{2})}\,\rho_{0}/\Omega$ changing, i.e., we have that

[TABLE]

for the Gaussian pinhole.

For the Gaussian pinspeck, we get

[TABLE]

This psf is a bit more complicated than what we found for the Gaussian pinhole. Nevertheless, it shows the expected result for a pinspeck camera, viz., that the image-bearing part of the psf is embedded in a uniform background term whose presence creates photodetection shot noise that degrades signal-to-noise ratio. As was the case for the Gaussian pinhole, we see that optimum spatial resolution occurs when $\Omega=1$, in which case we get

[TABLE]

On the other hand, unlike the Gaussian pinhole’s psf, the Gaussian pinspeck’s psf does not preserve its shape as the Fresnel number is varied. This is illustrated in Fig. 4, where we have plotted $G_{\rm ps}(\bm{\rho})/G_{\rm ps}({\bm{\infty}})$ versus $\bm{\rho}/\rho_{\rm res}(\Omega)$ for $\bm{\rho}=(x,0)$ and $\Omega=0.1,1,$ and 10, where $\rho_{\rm res}(\Omega)$ is the Gaussian pinhole’s spatial resolution.

Note that in the near-field region, wherein $\Omega\gg 1$, Eq. (77) reduces to the geometric optics result,

[TABLE]

which is analogous to the geometric optics treatment used by Xu et al. Xu2018 and Thrampoulidis et al. Thrampoulidis2018 for the hard-aperture, circular occluder

[TABLE]

IV Discussion

In summary, we have presented a complete light transport model, in phasor-field terms, capable of describing propagation through a transmissive, paraxial geometry—including intermediate occluders and a specular-plus-diffuser mask—that serves as an unfolded proxy for occlusion-aided, three-bounce NLoS imaging. For imaging purely diffuse objects without intermediate occluders, we phrased our analysis in terms of the $\mathcal{P}$ field and provided a straightforward derivation of its behavior, analogous to that reported by Reza et al. Reza2018 . To handle more general scenarios, we introduced and presented propagation primitives for the two-frequency spatial Wigner distribution (TFSWD). With these in hand, we turned our attention to the task of diffuse-object, occlusion-aided imaging and arrived at closed-form results for occlusion-aided imaging with unmodulated light using either a Gaussian-pinhole occluder or a Gaussian-pinspeck occluder. Our results show that imaging unoccluded diffuse objects with unmodulated light is not possible in the paraxial regime, but phasor-field imaging provides techniques for image construction if modulated light is used or object occlusion can be exploited. For imaging non-occluded diffuse objects with modulated light, spatial resolution is the diffraction limit at the modulation frequency. For occlusion-aided imaging of the same object with unmodulated light, spatial resolution is set by the optical-frequency diffraction limit of the occluder. Although the latter can be far superior to the former, blind determination of the occluder’s characteristics poses a challenge for exploiting its presence, and even with a known occluder, imaging performance will be limited by its size and shape.

There are many avenues for future research that build upon the work we have reported. Here we shall list just a few of the possibilities. First, because diffuse transmission (and, for the NLoS case, diffuse reflection) creates laser speckle, our assumption that we can measure the speckle-averaged, short-time average irradiance needs to examined. Toward that end, it is worth noting that Liu et al.’s experiments Liu2018 did not suffer any obvious ill effects of laser speckle. Second, it remains to be seen how occlusion-aided imaging with modulated light might benefit from synergy between the approaches we have examined. A third avenue to pursue is evaluating $\mathcal{P}$-field imaging of specular objects. Next, because Liu et al. Liu2018 used ps-duration pulsed illumination to obtain three-dimensional scene reconstructions—and such illumination violates our quasimonochromatic-light assumption—a fourth item on our plate would be to treat the pulsed case, including the value of synthesizing desirable input $\mathcal{P}$ fields. Fifth on our list is to extend our propagation primitives beyond the paraxial regime, i.e., to replace Fresnel diffraction with Rayleigh-Sommerfeld diffraction. Moreover, we need to address NLoS imaging explicitly, rather than its transmissive proxy, and include more than just three-bounce returns. It is also possible—and potentially interesting—to extend our TFSWD transport model to account for arbitrary linear transformations of the $E$ field of the type given by Eq. (48). Finally, the work we have presented could be fruitfully specialized to sinusoidal $E$-field modulation and wedded to the $\mathcal{P}$-field optics introduced and demonstrated in Reza et al. Reza2019 .

Funding

This work was supported by the DARPA REVEAL program under Contract HR0011-16-C-0030.

Appendix A Propagation Calculations

In this appendix we provide derivations for the TFSWD’s propagation primitives given earlier in Eqs. (51)–(54).

**Propagation through a diffuser:

**Consider propagation through one of our diffusers: assume that we know $W_{\mathcal{E}_{z}}(\bm{\rho}_{+},{\bf k},\omega_{+},\omega_{-})$ and we want to find $W_{\mathcal{E}^{\prime}_{z}}(\bm{\rho}_{+},{\bf k},\omega_{+},\omega_{-})$, where

[TABLE]

with

[TABLE]

In this case we immediately get

[TABLE]

Physically, the ${\bf k}$ dependence of the TFSWD carries the field’s spatial-frequency information, i.e., its directionality. The result we have just obtained shows that the diffuser has completely destroyed the directionality of $\mathcal{E}_{z}(\bm{\rho},\omega)$, because $W_{\mathcal{E}^{\prime}_{z}}(\bm{\rho}_{+},{\bf k},\omega_{+},\omega_{-})$ is independent of ${\bf k}$.

**Propagation through a deterministic occluder:

**Now consider propagation through a deterministic transmission mask. Here we want to find $W_{\mathcal{E}^{\prime}_{z}}(\bm{\rho}_{+},{\bf k},\omega_{+},\omega_{-})$ given $W_{\mathcal{E}_{z}}(\bm{\rho}_{+},{\bf k},\omega_{+},\omega_{-})$ and a deterministic $P(\bm{\rho})$, where

[TABLE]

For this case we have that

[TABLE]

[TABLE]

where

[TABLE]

is the spatial Wigner distribution of $P(\bm{\rho})$. In words, Eq. (91) shows that multiplying $\mathcal{E}_{z}(\bm{\rho},\omega)$ by a deterministic field-transmission mask implies that $W_{\mathcal{E}^{\prime}_{z}}(\bm{\rho}_{+},{\bf k},\omega_{+},\omega_{-})$ is obtained from a ${\bf k}$-space convolution of $W_{\mathcal{E}_{z}}(\bm{\rho}_{+},{\bf k},\omega_{+},\omega_{-})$ with the field-transmission mask’s spatial Wigner distribution. Moreover, Eq. (92), together with Eq. (50), immediately leads to

[TABLE]

as could have been directly obtained from Eq, (87) and the $\mathcal{P}$-field’s definition.

**Propagation through a specular-plus-diffuser mask:

**Combining the approaches for the diffuser and deterministic transmission mask allows us to model the propagation through a specular-plus-diffuser mask. We take such a mask to be a multiplicative random process $F(\bm{\rho}_{1})$ having nonzero mean $\langle F(\bm{\rho}_{1})\rangle\neq 0$, and covariance, $\langle\Delta F(\bm{\rho}_{+}+\bm{\rho}_{-}/2)\Delta F^{*}(\bm{\rho}_{+}-\bm{\rho}_{-}/2)\rangle\approx\lambda_{0}^{2}\mathcal{F}(\bm{\rho}_{+})\delta(\bm{\rho}_{-})$ where $0\leq\mathcal{F}(\bm{\rho}_{+})\leq 1$ and $\Delta F(\bm{\rho})\equiv F(\bm{\rho})-\langle F(\bm{\rho})\rangle$. The propagation analysis follows from combining the two previous analyses:

[TABLE]

From expanding $F(\bm{\rho})$ into a sum of its (deterministic) mean and zero-mean random portions, it follows that

[TABLE]

**Fresnel diffraction:

**Our final task is to find $W_{\mathcal{E}_{L}}(\bm{\rho}_{+},{\bf k},\omega_{+},\omega_{-})$ when

[TABLE]

i.e., for Fresnel diffraction over a distance $L$ footnote2 . This calculation turns out to be more complicated than its predecessors in this section. We start from

[TABLE]

Exploiting $\Delta\omega\ll\omega_{0}$, and making the coordinate transformation from $\bm{\rho}_{0}$ and $\bm{\rho}_{0}^{\prime}$ to $\bm{\rho}_{0_{+}}\equiv(\bm{\rho}_{0}+\bm{\rho}^{\prime}_{0})/2$ and $\bm{\rho}_{0_{-}}\equiv\bm{\rho}_{0}-\bm{\rho}_{0}^{\prime}$, we can reduce Eq. (101) to

[TABLE]

Rearranging terms allows us to put the $\bm{\rho}_{-}$ integral inside the $\bm{\rho}_{0_{+}}$ and $\bm{\rho}_{0_{-}}$ integrals, i.e.,

[TABLE]

Performing the $\bm{\rho}_{-}$ integral then yields

[TABLE]

which, after some terms cancel, gives

[TABLE]

The term

[TABLE]

in Eq. (106)’s integrand behaves like the impulse $\delta[\bm{\rho}_{0+}-\bm{\rho}_{+}+kcL/(\omega_{0}+\omega_{+})]$. This delta-function behavior follows because: (1) The term in question is a highly-oscillatory function outside of a narrow slow-oscillation region that is centered at $\bm{\rho}_{+}-kcL/(\omega_{0}+\omega_{+})$ with nominal width $\sqrt{\omega_{-}cL}/2(\omega_{0}+\omega_{+})$, and $\omega_{0}\gg\max|\omega_{+}|$ implies that it integrates to one. (2) The other $\bm{\rho}_{0_{+}}$-dependent terms in Eq. (106) are the oscillatory term, $\exp(i\omega_{-}|\bm{\rho}_{+}-\bm{\rho}_{0_{+}}|^{2}/2cL)$, which varies much more slowly than its predecessor, because $\omega_{0}\gg\max|\omega_{-}|$, and the Wigner distribution, whose $\bm{\rho}_{0_{+}}$ dependence can reasonably be assumed to be nearly constant over regions of diameter $\sqrt{\omega_{-}cL}/2(\omega_{0}+\omega_{+})$. So, using the delta-function approximation in Eq. (106), we get

[TABLE]

Finally, again making use $\omega_{0}\gg\omega_{+}$, we have

[TABLE]

As a consistency check on Eq. (108), let us use it to calculate $\mathcal{P}_{L}(\bm{\rho}_{+},\omega_{-})$ when $z=0$ illumination with TFSWD $W_{\mathcal{E}_{0}}(\bm{\rho}_{0_{+}},{\bf k},\omega_{+},\omega_{-})$ passes through the diffuser specified in Eq. (81) before undergoing Fresnel diffraction over a distance $L$. We then have that

[TABLE]

Using Eq. (86) now gives us

[TABLE]

Changing variables so that ${\bf k}=\omega_{0}(\bm{\rho}_{+}-\bm{\rho}_{0})/cL$ leaves us with

[TABLE]

which reduces to the result from Sec. II,

[TABLE]

by virtue of Eq. (50).