1. Introduction

Non-line-of-sight (NLoS) imaging is a growing field of research concerned with the task of generating accurate reconstructions of scenes whose observation can only be accomplished by means of intervening diffuse scattering events, e.g., by penetrating diffuse transmissive scattering media like ground glass or fog, or reflecting off diffuse surfaces like conventional walls. Traditional approaches to this problem, which largely depend on time-of-flight (ToF) information and computational backprojection, have been facilitated by pulsed illumination and time-resolved detection [1,2]. More recent approaches have leveraged closed-form inversion techniques for ToF data [3,4], or the presence of occluders that obviate the need for ToF information [5,6]. In all of these schemes, the task of forming an image from collected data is computational, often requiring significant overhead.

Phasor-field (${{\mathcal {P}}}$-field) imaging is a new approach to NLoS imaging that exploits the wave-like propagation behavior of the oscillating envelope of amplitude-modulated, spatially-incoherent light [7–9]. It relies on the physical correlates of diffuse phase disruption often being large compared to the wavelength of the optical carrier field but insignificant relative to the much longer wavelength of radio-frequency (or even microwave) amplitude modulation. As a result, walls that diffusely scatter light, owing to their roughness at the optical-wavelength scale, appear smooth to the ${{\mathcal {P}}}$ field, and thin transmissive diffusers appear transparent to the ${{\mathcal {P}}}$ field. The upshot is that traditional wave-optical imaging techniques, e.g., lenses, can be applied to the ${{\mathcal {P}}}$ field despite the presence of these optically disruptive elements. To date, however, this capability has only been applied computationally to form NLoS images from ToF data [10,11].

Recently, Reza et al. [12] reported experiments verifying the ${{\mathcal {P}}}$ field’s physical wave-like properties. In one of them, they showed that a diffuse, concave reflector focuses the ${{\mathcal {P}}}$ field, even though it scatters the optical carrier. In another, they showed that an ordinary Fresnel lens focuses the optical carrier but not the ${{\mathcal {P}}}$ field. Reza et al.’s demonstrations led us to wonder how ${{\mathcal {P}}}$-field propagation might be controlled by an ordinary lens—as opposed to the engineered diffuse reflector they used—and, in particular, is there a ${{\mathcal {P}}}$-field physical optics configuration that can enable seeing around corners in real time without the nontrivial computational burden of prior NLoS techniques. In this paper, we answer those questions through a rigorous mathematical demonstration that ordinary lenses can indeed focus or project the ${{\mathcal {P}}}$ field through intervening diffusers and that they can image a scene hidden by such diffusers.

2. ${{\mathcal {P}}}$-field setup for computational imaging

We begin by reviewing our previously developed framework for computational ${{\mathcal {P}}}$-field imaging [8]. We assume paraxial [13], scalar-wave optics wherein an optical carrier at frequency $\omega _0$ is modulated by a baseband complex field envelope $E_z(\boldsymbol {\rho }_z,t)$ of bandwidth ${{\Delta \omega }} \ll \omega _0$ to produce an optical field $U_z(\boldsymbol {\rho }_z,t) = \textrm {Re}[E_z(\boldsymbol {\rho }_z,t) e^{-i\omega _0 t}]$, where $\boldsymbol {\rho }_z$ is the two-dimensional transverse spatial coordinate in the plane denoted by $z$. The amplitude modulation is characterized by the short-time-average (STA) irradiance, given by $I_z(\boldsymbol {\rho }_z,t) = |E_z(\boldsymbol {\rho }_z,t)|^{2}$, which can be measured by direct photodetection assuming detectors with sufficient bandwidth. The ${{\mathcal {P}}}$ field is defined to be the temporal Fourier transform of the STA irradiance, averaged over any diffusers present in the scenario, whose surface fluctuations are treated statistically:

We assume a transmissive geometry, shown in Fig. 1, that is a proxy for three-bounce reflective NLoS imaging. The frequency-domain complex envelope, ${{\mathcal {E}}}^{\prime}_z(\boldsymbol {\rho }_z,\omega )$, that emerges from the $z$-plane diffuser when that diffuser is illuminated by ${{\mathcal {E}}}_z(\boldsymbol {\rho }_z,\omega )$ is given by

 

Fig. 1. Unfolded geometry for three-bounce NLoS active imaging. The blue rectangles represent thin transmissive diffusers, which serve as analogs of diffuse reflections, and the black line represents a thin transmissivity mask whose intensity transmission pattern, $T(\boldsymbol {\rho }_1)$, represents the albedo pattern of a diffuse target in the hidden scene. Here, and elsewhere, we simplify subscripts involving $z$ where the meaning is clear from context.

Download Full Size | PDF

Equations (8) and (9) immediately imply a backprojection procedure to image the hidden albedo pattern, $T(\boldsymbol {\rho }_1)$. In particular, because

The central purpose of this paper is to show that the foregoing computational approach has a ${{\mathcal {P}}}$-field physical optics replacement. First, in Sec. 3, we prove that a convex lens illuminated by a plane-wave complex field envelope can cast a focused ${{\mathcal {P}}}$ field onto a plane that is hidden by an intervening diffuser which broadly scatters the optical carrier. Next, in Sec. 4, we establish a useful primitive for ${{\mathcal {P}}}$-field propagation from an initial diffuser through a convex lens and then to an output plane. The utility of that primitive is demonstrated in Sec. 5, where we use it to show that a convex lens can project an arbitrary ${{\mathcal {P}}}$ field from an initial diffuser to an output plane despite there being another diffuser located between the lens and the output plane. The primitive’s utility is also seen in Sec. 6, where we employ it to show that a convex lens enables a diffuser’s albedo pattern to be imaged onto a detector plane despite that diffuser’s being hidden by another diffuser.

3. Plane-wave ${{\mathcal {P}}}$-field focusing

The focusing capability of a convex lens is the natural starting point for that optical element’s use in conventional physical optics. Thus we begin our development of ${{\mathcal {P}}}$-field physical optics by showing how such a lens can focus the ${{\mathcal {P}}}$ field. Consider the configuration shown in Fig. 2, in which an infinite plane wave with complex field envelope

 

Fig. 2. Geometry for focusing a plane wave through an intervening diffuser onto a target plane containing a diffuser with an albedo pattern $T(\boldsymbol {\rho }_1)$. The dashed red arrow is a unit vector representing the plane wave’s propagation direction, with ${\textbf {s}}$ being its transverse component. The lens has a Gaussian field-transmission pupil $e^{-|\boldsymbol {\rho }_{\textrm{in}}|^{2}/2D^{2}}$.

Download Full Size | PDF

The temporal-frequency-domain complex field envelope at the first diffuser is given by

From the perspective of ${{\mathcal {P}}}$-field imaging, focusing enables us to raster scan the target as if the initial diffuser were not there. That said, although the lens focuses the ${{\mathcal {P}}}$ field, it does not focus the optical power, which is still spread out by the diffuser, as demonstrated by Reza et al. [12] for their diffuse, concave reflector. Because the ${{\mathcal {P}}}$ field is ultimately supported by the optical field, its peak strength, like the STA irradiance’s, is still subject to inverse-square law falloff, even in the presence of the lens. That Eq. (20) suggests otherwise, i.e., that increasing $D$ can offset the inverse-square law attenuation, is because we have assumed infinite-plane-wave illumination for which the power passing through the lens is proportional to $D^{2}$. Correcting for this scaling, it is clear that the peak ${{\mathcal {P}}}$ field has an inverse-square law falloff relative to the input power, regardless of the pupil diameter, i.e., regardless of how tightly the ${{\mathcal {P}}}$ field is confined in the target plane.

4. Lens primitive for ${{\mathcal {P}}}$-field projection and hidden-plane imaging

In this section we derive a ${{\mathcal {P}}}$-field propagation primitive for post-diffuser ${{\mathcal {P}}}$-field Fresnel diffraction through an intervening focal-length-$f$ thin lens that has a Gaussian field-transmission pupil $e^{-|\boldsymbol {\rho}_b |^{2}/2D^{2}}$, as depicted in Fig. 3. This primitive will prove useful for both the ${{\mathcal {P}}}$-field projection and ${{\mathcal {P}}}$-field imaging cases that follow. We have that

 

Fig. 3. Geometry for post-diffuser ${{\mathcal {P}}}$-field Fresnel propagation with an intervening lens. An optical field with complex envelope $E_a(\boldsymbol {\rho }_a,t)$ illuminates a diffuser $h_a(\boldsymbol {\rho }_a)$ located a distance $L_b$ in front of a convex lens with focal length $f$ and a Gaussian field-transmission pupil $e^{-|\boldsymbol {\rho}_b |^{2}/2D^{2}}$. This section’s objective is to relate the ${{\mathcal {P}}}$ field in the plane $L_c\,$m behind the lens to the ${{\mathcal {P}}}$ field entering the diffuser.

Download Full Size | PDF

5. ${{\mathcal {P}}}$-field projection

With the Eq. (26) primitive in hand, we now turn our attention to the task of projecting an arbitrary ${{\mathcal {P}}}$-field pattern from an input diffuser to a hidden target plane. This arrangement will allow structured-illumination paradigms for line-of-sight imaging with coherent light to be translated into structured-illumination NLoS imaging with ${{\mathcal {P}}}$ fields.

Consider an augmentation of the Fig. 1 scenario in which that figure’s initial diffuser is preceded by an instance of the lens primitive depicted in Fig. 3, as shown in Fig. 4. Modifying the Fig. 3 scenario’s placeholder notation, we will label the transverse coordinate of Fig. 4’s input plane as $\boldsymbol {\rho }_{\textrm{in}}$, its first distance as $L_{\textrm{in}_1}$, and its second distance as $L_{\textrm{in}_2}$. The lens primitive’s output-plane transverse coordinate remains as $\boldsymbol {\rho }_0$, leading into the same notation as Fig. 1 for the rest of that figure’s geometry. We take the lens to be configured to project the input ${{\mathcal {P}}}$ field onto the $z=L_1$ plane by choosing its focal length to obey $1/f = 1/L_{\textrm{in}_1} + 1/(L_{\textrm{in}_2}+L_1)$. For this configuration we have $\Pi (f,L_{\textrm{in}_1},L_{\textrm{in}_2}) = L_{\textrm{in}_2}(L_1+L_{\textrm{in}_2})/L_1$, and so the lens primitive gives us

 

Fig. 4. Geometry for projecting an arbitrary ${{\mathcal {P}}}$-field pattern from an input diffuser [23] to a hidden target plane. An optical field with complex envelope $E_{\mathrm{in}}(\boldsymbol {\rho }_{\mathrm{in}},t)$ illuminates an instance of the Fig. 3 geometry in which $E_{\textrm{in}}(\boldsymbol {\rho }_{\textrm{in}},t)$ takes the place of $E_a(\boldsymbol {\rho }_a,t)$, and $L_{\textrm{in}_1}$, $L_{\textrm{in}_2}$ take the places of $L_b$, $L_c$, respectively. The lens, whose focal length satisfies $1/f = 1/L_{\textrm{in}_1} + 1/(L_{\textrm{in}_2}+L_1)$, then projects the ${{\mathcal {P}}}$ field emerging from the initial diffuser onto the $z=L_1$ plane with magnification/minification factor $M = L_{\mathrm{prj}}/L_{\textrm{in}_1}$, where $L_{\textrm{prj}} \equiv L_{\textrm{in}_2} + L_1$.

Download Full Size | PDF

Ignoring inessential phase and scaling terms, this is a projected copy of the input ${{\mathcal {P}}}$ field,

6. ${{\mathcal {P}}}$-field imaging of a hidden plane

It turns out that the same kind of result just exhibited for ${{\mathcal {P}}}$-field projection can be obtained for physical-optics ${{\mathcal {P}}}$-field imaging of a hidden target plane. To do so, we place our lens primitive from Fig. 3 behind Fig. 1’s final diffuser so that the transverse coordinate of the first plane of the primitive is $\boldsymbol {\rho }_2$, as depicted in Fig. 5. We denote the lens primitive’s two distances as $L_{\textrm{im}_1}$ and $L_{\textrm{im}_2}$, and we denote the transverse coordinate of the primitive’s final plane as $\boldsymbol {\rho }_{\textrm{im}}$ with the intention in mind that a multipixel detector array capable of measuring the ${{\mathcal {P}}}$ field will be present there. The hidden plane is imaged onto this hypothetical detector by taking the focal length of the lens to obey $1/f=1/(L_2+L_{\textrm{im}_1})+1/L_{\textrm{im}_2}$. For this configuration we have $\Pi (f,L_{\textrm{im}_1},L_{\textrm{im}_2}) = L_{\textrm{im}_1}(L_2+L_{\textrm{im}_1})/L_2$, and so we get

 

Fig. 5. Geometry for physical-optics ${{\mathcal {P}}}$-field imaging of a hidden target plane. The complex field envelope $E_{1}(\boldsymbol {\rho }_{1},t)$ is the hidden-plane illumination in Fig. 1 that results when $E_0(\boldsymbol {\rho }_0,t)$ illuminates a diffuser in the $z=0$ plane. The rest of the setup shown here is an instance of the lens primitive from Fig. 3 in which $E_{2}(\boldsymbol {\rho }_{2},t)$ takes the place of $E_a(\boldsymbol {\rho }_a,t)$ and $L_{\textrm{in}_1}$, $L_{\textrm{in}_2}$ take the places of $L_b$, $L_c$, respectively. The lens, whose focal length satisfies $1/f = 1/(L_2+L_{\textrm{im}_1}) + 1/L_{\textrm{im}_2}$ then casts an image of ${{\mathcal {P}}}_1(\boldsymbol {\rho }_1,t)T(\boldsymbol {\rho }_1)$ onto the $z= L_1+L_2+L_{\textrm{im}_1}+L_{\textrm{im}_2}$ plane with minification/magnification factor $M = L_{\textrm{out}}/L_{\textrm{im}_2}$, where $L_{\textrm{out}} \equiv L_2 + L_{\textrm{im}_1}.$

Download Full Size | PDF

7. Summary and discussion

In this paper, we have mathematically investigated the possibility for NLoS imaging tasks to be carried out in the ${{\mathcal {P}}}$-field framework using physical optics in lieu of the more involved computational techniques employed to date. Using an unfolded proxy for three-bounce NLoS imaging and assuming paraxial propagation, we first demonstrated that plane-wave ${{\mathcal {P}}}$-field illumination can be focused onto a small target-plane region despite the presence of an intervening diffuser. Moreover, this ${{\mathcal {P}}}$-field focusing occurs even though the optical carrier, and thus the power, is broadly scattered. Next, we developed a useful primitive for ${{\mathcal {P}}}$-field propagation in scenarios involving lenses. We then applied that primitive to ${{\mathcal {P}}}$-field projection and ${{\mathcal {P}}}$-field imaging of a hidden plane’s albedo pattern. We found that arbitrary ${{\mathcal {P}}}$-field patterns can be projected, through an intervening diffuser, onto a target plane. Likewise, we showed that ${{\mathcal {P}}}$-field imaging of a target plane’s albedo pattern can be accomplished despite the presence of another intervening diffuser between that plane and the sensor. In these last two cases, the mathematical assumptions made in our analysis rely on rather large lenses. However, the qualifying lens size can be reduced by increasing the ${{\mathcal {P}}}$-field frequency. For all three cases—${{\mathcal {P}}}$-field focusing, projection, and imaging—we found that the spatial resolution was at the ${{\mathcal {P}}}$-field frequency’s diffraction limit, further encouraging the pursuit of higher ${{\mathcal {P}}}$-field frequencies.

Ultimately, the highest usable ${{\mathcal {P}}}$-field frequency for physical-optics imaging will be set by the bandwidth of available detectors, which is not likely to exceed 100 GHz. Despite this detector limitation, we allege that higher ${{\mathcal {P}}}$-field frequencies, perhaps up to 1 THz, can be achieved using Willomitzer et al.’s synthetic-wavelength holography [20], wherein coherently-detected outputs from sequentially illuminating an NLoS system with unmodulated light at two optical frequencies can be correlated and filtered to produce a synthetic ${{\mathcal {P}}}$ field at the difference frequency, without the need for ultrafast detectors. This technique offers a computational approach to ${{\mathcal {P}}}$-field generation and detection that stands in contrast to this paper’s physical-optics paradigm. However, the required computation is trivial in comparison to what is typically used, at present, for extracting images from NLoS datasets. So, even for this synthetic approach, the conclusion of this paper stands: there is real potential for NLoS imaging to be accomplished with traditional optical techniques in the ${{\mathcal {P}}}$-field framework using real, physical optics, obviating the need for the nontrivial computational inversion schemes that have been applied to date. Nevertheless, much more needs to be done to bring ${{\mathcal {P}}}$-field physical optics NLoS imaging to fruition. For one thing, many NLoS geometries will violate the validity conditions for paraxial propagation, thus requiring use of Rayleigh–Sommerfeld diffraction for accurate scene reconstructions from ${{\mathcal {P}}}$-field measurements. We have already taken initial steps [15,24] in this direction, in the transmissive geometry, but additional work is needed to convert them to the three-bounce NLoS scenario. For another, our transmissive geometry addresses 2D imaging, i.e., recovering the albedo pattern present on a hidden planar diffuser, but NLoS imaging almost invariably involves reconstructing a hidden 3D scene. Consequently, for NLoS imagers depth-of-focus and depth-of-eld considerations come into play. We expect that these issues can be addressed in our ${{\mathcal {P}}}$-field framework, but that is a subject for future work.

Appendix A: synthetic-wavelength holography

In this appendix, we will develop Willomitzer et al.’s synthetic-wavelength holography [20] in our ${{\mathcal {P}}}$-field framework. Within the limitations of current laser and detector technology, we will show that synthetic-wavelength holography (SWH) offers a practical method to access higher ${{\mathcal {P}}}$-field frequencies than are presently attainable with modulated illumination and direct detection.

Figure 6(a) shows the physical setup for SWH NLoS imaging. A tunable continuous-wave laser sequentially illuminates a spot on a visible wall at frequencies $\omega _0$ and $\omega _1 = \omega _0 - \Delta \omega$. Diffuse reflection from the visible wall then illuminates an object in the hidden space and diffuse reflections therefrom return to a large swath of the visible wall. Light reflected from a visible-wall patch that is disjoint from the laser-illuminated spot is imaged onto a heterodyne-detection array. This array gets its local-oscillator field from acousto-optic modulation (AOM) or electro-optic modulation (EOM) of the transmitter laser that produces an intermediate-frequency (IF) offset $\omega _{\textrm{IF}}$. The array uses a lock-in camera (or equivalent) that outputs the pixel-wise quadrature components of the IF signals. These pixel-wise quadratures are subsequently processed to generate a computational image of the hidden object.

 

Fig. 6. Synthetic wavelength holography (SWH). (a) Setup for NLoS SWH imaging as performed by Willomitzer et al. [20]. A tunable continuous-wave laser sequentially illuminates a spot on the visible wall at frequencies $\omega _0$ and $\omega _1 = \omega _0-\Delta \omega$. Diffuse reflection from that wall illuminates a hidden object. Diffuse reflection from that object illuminates the visible wall and a portion of that wall—disjoint from the original illumination spot—is imaged onto a heterodyne detection array. That array’s local oscillator is obtained by acousto-optic modulation (AOM) or electro-optic modulation (EOM) that produces an intermediate-frequency shift $\omega _{\textrm{IF}}$. The array uses, e.g., a lock-in camera that outputs the pixel-wise quadrature components of the IF signals that are used to generate a computational image of the hidden object. (b) Measurement and signal-processing chain for SWH imaging in our transmissive-geometry proxy for three-bounce NLoS imaging. Ideal heterodyne detection is assumed, with a continuum array whose signal-to-noise ratio is high enough that its outputs can be taken to be the complex field envelopes, ${{\mathcal {E}}}_{\textrm{out},\omega _0}(\boldsymbol {\rho })$ and ${{\mathcal {E}}}_{\textrm{out},\omega _1}(\boldsymbol {\rho })$, of the imaged fields. Multiplying ${{\mathcal {E}}}_{\textrm{out},\omega _0}(\boldsymbol {\rho })$ by the complex conjugate of ${{\mathcal {E}}}_{\textrm{out},\omega _1}(\boldsymbol {\rho })$ and passing that product through a spatial low-pass filter (LPF) then yields $\hat {{{\mathcal {P}}}}_{\textrm{out},\Delta \omega }(\boldsymbol {\rho })$, i.e., an estimate of the desired ${{\mathcal {P}}}$ field.

Download Full Size | PDF

Willomitzer et al. drew their inspiration from Gabor’s Nobel-prize winning work on holography, and their analysis/synthesis approach to SWH does not rely on statistics. Our objective, in this appendix, is to flesh out the ${{\mathcal {P}}}$-field SWH measurement and signal-processing chain shown in Fig. 6(b) to show how our statistical framework reproduces the basic concept from their paper in our transmissive-geometry proxy for NLoS imaging.

To begin, consider the ${{\mathcal {P}}}$ field associated with single-frequency irradiance modulation at angular frequency ${{\Delta \omega }}$:

(36)$${{\mathcal{P}}}(\boldsymbol{\rho},\omega_-) = {{\mathcal{P}}}^{*}_{{{\Delta\omega}}}(\boldsymbol{\rho})2\pi\delta(\omega_-+{{\Delta\omega}}) + {{\mathcal{P}}}_{0}(\boldsymbol{\rho})2\pi\delta(\omega_-) + {{\mathcal{P}}}_{{{\Delta\omega}}}(\boldsymbol{\rho})2\pi\delta(\omega_--{{\Delta\omega}}).$$

One way to generate this ${{\mathcal {P}}}$ field is by coherently summing optical fields at two angular frequencies, $\omega _0$ and $\omega _1 \equiv \omega _0 - {{\Delta \omega }}$, in which case the frequency-domain complex field envelope associated with such a signal is

(37)$${{\mathcal{E}}}(\boldsymbol{\rho},\omega) = {{\mathcal{E}}}_{\omega_0}(\boldsymbol{\rho})2\pi\delta(\omega) + {{\mathcal{E}}}_{\omega_1}(\boldsymbol{\rho})2\pi\delta(\omega+{{\Delta\omega}}).$$

The ${{\mathcal {P}}}$ field associated with this complex field envelope is given by

(38)$${{\mathcal{P}}}(\boldsymbol{\rho},\omega_-) = \int\! \frac{\textrm{d}\mathrm{\omega}_+}{2 \pi} \langle {{{\mathcal{E}}}(\boldsymbol{\rho},\omega_++\omega_-/2){{\mathcal{E}}}^{*}(\boldsymbol{\rho},\omega_+-\omega_-/2)} \rangle$$

(39)$$\begin{aligned}&= \langle {{{\mathcal{E}}}^{*}_{\omega_0}(\boldsymbol{\rho}){{\mathcal{E}}}_{\omega_1}(\boldsymbol{\rho})} \rangle2\pi\delta(\omega_-+{{\Delta\omega}}) + (\langle {\left| {{{\mathcal{E}}}_{\omega_0}(\boldsymbol{\rho})} \right|^{2}} \rangle+ \langle {\left| {{{\mathcal{E}}}_{\omega_1}(\boldsymbol{\rho})} \right|^{2}} \rangle) 2\pi\delta(\omega_-)\\ &+ \langle {{{\mathcal{E}}}_{\omega_0}(\boldsymbol{\rho}){{\mathcal{E}}}^{*}_{\omega_1}(\boldsymbol{\rho})} \rangle2\pi\delta(\omega_--{{\Delta\omega}}). \end{aligned}$$

This is a single-frequency ${{\mathcal {P}}}$ field at frequency ${{\Delta \omega }}$. In particular, we have that

(40a)$${{\mathcal{P}}}_{0}(\boldsymbol{\rho}) = \langle {\left| {{{\mathcal{E}}}_{\omega_0}(\boldsymbol{\rho})} \right|^{2}} \rangle+ \langle {\left| {{{\mathcal{E}}}_{\omega_1}(\boldsymbol{\rho})} \right|^{2}} \rangle$$

(40b)$${{\mathcal{P}}}_{{\Delta\omega}}(\boldsymbol{\rho}) = \langle {{{\mathcal{E}}}_{\omega_0}(\boldsymbol{\rho}){{\mathcal{E}}}^{*}_{\omega_1}(\boldsymbol{\rho})} \rangle.$$

For systems that are linear and time invariant with respect to the optical field—as is the case for typical NLoS scenarios—the frequency structure of the ${{\mathcal {P}}}$ field and its underlying complex field envelope are unaffected by that system. Consequently, propagating them through that system reduces to relating the input quantities to the output quantities at each frequency. It follows from Eq. (40) that the input-output relations for the complex field envelopes at $\omega _0$ and $\omega _1$ suffice to characterize the input-output behavior of the ${{\mathcal {P}}}$ field. More importantly, for linear time-invariant systems, the input-output relations that take inputs ${{\mathcal {E}}}_{\textrm{in},\omega _0}(\boldsymbol {\rho }_{\textrm{in}})$ and ${{\mathcal {E}}}_{\textrm{in},\omega _1}(\boldsymbol {\rho }_{\textrm{in}})$ and yield outputs ${{\mathcal {E}}}_{\textrm{out},\omega _0}(\boldsymbol {\rho }_{\textrm{out}})$ and ${{\mathcal {E}}}_{\textrm{out},\omega _1}(\boldsymbol {\rho }_{\textrm{out}})$ are the same as the complex field envelopes’ input-output relations for illuminating the system with unmodulated light at each frequency individually. The implication is that, provided we can accurately measure the output complex field envelopes and accomplish the required diffuser averaging, phasor-field imaging tasks can be carried out by sequentially illuminating the system with unmodulated inputs, meaning we are not burdened by needing direct detectors that are sufficiently fast to capture the ${{\mathcal {P}}}$-field modulation frequency. We will turn now to each of these concerns where, for convenience, we will use ${\boldsymbol {\rho }}$ instead of ${\boldsymbol {\rho }_{\textrm{out}}}$ for the output plane's coordinate vector.

To measure the output complex field envelope ${{\mathcal {E}}}_{\textrm{out},\omega _n}(\boldsymbol {\rho })$, for $n=0,1$, we will use balanced heterodyne detection [25,26]. The signal field is mixed on a 50–50 beam splitter with a strong, plane-wave local oscillator that is detuned from $\omega _n$ by an intermediate frequency $\omega _{\textrm{IF}} \ll \Delta \omega$ and has phase $\theta _{\textrm{LO}}$. The light emerging from the beam-splitter’s output arms are detected, at high spatial resolution, with detector arrays whose temporal bandwidths, including those of their post-detection electronics, exceed $\omega _{\textrm{IF}}$. For simplicity, we will treat these arrays as performing ideal continuum photodetection, i.e., they have unlimited spatial resolution. We will also ignore the effects of noise and assume the arrays have unity quantum efficiency. In effect, we take the arrays to accurately detect the STA irradiances at the beam-splitter outputs—denoted $I_+(\boldsymbol {\rho },t)$ and $I_-(\boldsymbol {\rho },t)$—at full spatial resolution. The difference of these outputs carries the desired complex field envelope at the intermediate frequency:

(41)$$I_+(\boldsymbol{\rho},t)-I_-(\boldsymbol{\rho},t) \propto \textrm{Re}\!\left[ {{{\mathcal{E}}}_{\textrm{out},\omega_n}(\boldsymbol{\rho}) e^{-i(\omega_{\textrm{IF}} t - \theta_{\textrm{ LO}})}} \right],$$

when the illumination frequency is $\omega _n$, for $n=0,1$. Thus the quadratures of that signal can be extracted by standard communication electronics and we obtain ${{\mathcal {E}}}_{\textrm{out},\omega _n}(\boldsymbol {\rho })$.

Now we can computationally form ${{\mathcal {E}}}_{\textrm{out},\omega _0}(\boldsymbol {\rho }){{\mathcal {E}}}^{*}_{\textrm{out},\omega _1}(\boldsymbol {\rho })$, but we still need to perform diffuser averaging. It is inadvisable to use arrays whose detector elements average over many speckles, because that destroys the signal of interest, i.e.,

(42)$$\left\langle {I_+(\boldsymbol{\rho},t)-I_-(\boldsymbol{\rho},t)} \right\rangle \propto \textrm{Re}\!\left[ {\langle {{{\mathcal{E}}}_{\textrm{out},\omega_n}(\boldsymbol{\rho})} \rangle e^{-i(\omega_{\textrm{IF}} t - \theta_{\textrm{ LO}})}} \right] = 0, \hbox{for} \;n=0,1.$$

Hence, it is critical that we use detector arrays whose spatial resolution is high enough to resolve the speckle pattern. In practice, we want no more than a few speckle cells to fall on each pixel in our detector arrays to avoid the associated signal attenuation. So, to perform diffuser averaging, we propose low-pass spatial filtering ${{\mathcal {E}}}_{\textrm{out},\omega _0}(\boldsymbol {\rho }){{\mathcal {E}}}^{*}_{\textrm{out},\omega _1}(\boldsymbol {\rho })$ with a Gaussian kernel to generate

(43)$$\hat{{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}(\boldsymbol{\rho}) \equiv \frac{1}{2\pi R^{2}} \int\! \rm{d}^{2}\boldsymbol{\tilde{\rho}}\; e^{-\left| {\boldsymbol{\rho}-\boldsymbol{\tilde{\rho}}} \right|^{2}/2R^{2}}{{\mathcal{E}}}_{\textrm{out},\omega_0}(\boldsymbol{\tilde{\rho}}){{\mathcal{E}}}^{*}_{\textrm{out},\omega_1}(\boldsymbol{\tilde{\rho}}).$$

Provided that our Gaussian kernel resolves the spatial features in the ${{\mathcal {P}}}$ field, viz., $R < \Delta \lambda \equiv 2\pi c / {{\Delta \omega }}$, then $\hat {{\mathcal {P}}}_{\textrm{out},{{\Delta \omega }}}(\boldsymbol {\rho })$ is approximately unbiased, i.e.,

(44)$$\langle {\hat{{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}(\boldsymbol{\rho})} \rangle = \frac{1}{2\pi R^{2}} \int\!{\textrm{d}}^{2}\boldsymbol{\tilde{\rho}}\; e^{-\left| {\boldsymbol{\rho}-\boldsymbol{\tilde{\rho}}} \right|^{2}/2R^{2}}\langle {{{\mathcal{E}}}_{\textrm{out},\omega_0}(\boldsymbol{\tilde{\rho}}){{\mathcal{E}}}^{*}_{\textrm{out},\omega_1}(\boldsymbol{\tilde{\rho}})} \rangle$$

(45)$$= \frac{1}{2\pi R^{2}} \int\! \rm{d}^{2}\boldsymbol{\tilde{\rho}}\; e^{-\left| {\boldsymbol{\rho}-\boldsymbol{\tilde{\rho}}} \right|^{2}/2R^{2}}{{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}(\boldsymbol{\tilde{\rho}})$$

(46)$$\approx {{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}(\boldsymbol{\rho}).$$

What remains then is to assess the variance of $\hat {{\mathcal {P}}}_{\textrm{out},{{\Delta \omega }}}(\boldsymbol {\rho })$. To do so, we appeal to a speckle analysis developed elsewhere [14,15]. We assume the unfolded three-bounce geometry from Fig. 1, with monochromatic (frequency $\omega _n$) illumination of power $P$ with spatial profile ${{\mathcal {E}}}_0(\boldsymbol {\rho }_0) = \sqrt {8P/\pi d_0^{2}} \exp (-4|\boldsymbol {\rho }_0|^{2}/d_0^{2})$. For simplicity, the target at the second bounce is replaced by a Gaussian field-transmissivity pupil, $\exp (-4|\boldsymbol {\rho }_1|^{2}/d_1^{2})$. Another such pupil, $\exp (-4|\boldsymbol {\rho }_2|^{2}/d_2^{2})$, is placed at the third-bounce location to characterize the finite size of the visible wall, and the distances propagated after each bounce are all $L$. We shall also assume that the Fresnel number products, $d_0^{2} d_1^{2} / \lambda _0^{2} L^{2}$ and $d_1^{2} d_2^{2} / \lambda _0^{2} L^{2}$, are large enough that the third-bounce speckle is well-approximated as single-bounce speckle from the final diffuser [14]. This condition leads to a third-bounce complex field envelope that is Gaussian distributed. Gaussian moment factoring then gives us

(47)$$\textrm{var}(\hat{{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}(\boldsymbol{\rho})) = \langle {|\hat{{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}(\boldsymbol{\rho})|^{2}} \rangle - |\langle {\hat{{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}(\boldsymbol{\rho})} \rangle|^{2}$$

(48)$$\begin{aligned}= &\frac{1}{4\pi^{2} R^{4}}\int\!\textrm{d}^{2}{\boldsymbol{\rho}_{\textrm{out}}}\, \int\!\textrm{d}^{2}{\boldsymbol{\tilde{\rho}}_{\textrm{out}}}\, e^{-(\left| {\boldsymbol{\rho}-\boldsymbol{\rho}_{\textrm{out}}} \right|^{2} + \left| {\boldsymbol{\rho}-\boldsymbol{\tilde{\rho}}_{\textrm{out}}} \right|^{2})/2R^{2}} \\ &\langle {{{\mathcal{E}}}_{\textrm{out},\omega_0}(\boldsymbol{\rho}_{\textrm{out}}){{\mathcal{E}}}^{*}_{\textrm{out},\omega_0}(\boldsymbol{\tilde{\rho}}_{\textrm{out}})} \rangle \langle {{{\mathcal{E}}}^{*}_{\textrm{out},\omega_1}(\boldsymbol{\rho}_{\textrm{out}}){{\mathcal{E}}}_{\textrm{out},\omega_1}(\boldsymbol{\tilde{\rho}}_{\textrm{out}})} \rangle, \end{aligned}$$

where $\mathcal {E}_{\textrm{out},\omega _n}(\boldsymbol {\rho }_{\textrm{out}})$—the spatial profile in the plane that lies a distance $L$ beyond Fig. 1’s output plane—satisfies

(49)$${{\mathcal{E}}}_{\textrm{out},\omega_n}(\boldsymbol{\rho}_{\textrm{out}}) = \frac{\omega_n}{i2\pi c L} e^{i\omega_n L/c} \int\!\textrm{d}^{2}{\boldsymbol{\rho}_2}\, {{\mathcal{E}}}_{2,\omega_n}(\boldsymbol{\rho}_2) e^{i\omega_n\left| {\boldsymbol{\rho}_{\textrm{out}}-\boldsymbol{\rho}_2} \right|^{2} / 2 c L} e^{i \omega_n h_2(\boldsymbol{\rho}_2) / c} e^{-4\left| {\boldsymbol{\rho}_2} \right|^{2}/d_2^{2}},$$

and similar relations hold for the previous bounces. Accordingly, we find that

(50)$$\langle {{{\mathcal{E}}}_{\textrm{out},\omega_n}(\boldsymbol{\rho}_{\textrm{out}}){{\mathcal{E}}}^{*}_{\textrm{out},\omega_n}(\boldsymbol{\tilde{\rho}}_{\textrm{out}})} \rangle = \langle {I_{\textrm{out}}} \rangle e^{-\omega_n^{2} d_2^{2} \left| {\boldsymbol{\rho}_{\textrm{out}}-\boldsymbol{\tilde{\rho}}_{\textrm{out}}} \right|^{2} / 32 c^{2} L^{2}} e^{i\omega_n(\left| {\boldsymbol{\rho}_{\textrm{out}}} \right|^{2}-\left| {\boldsymbol{\tilde{\rho}}_{\textrm{out}}} \right|^{2}) / 2 c L},$$

with $\langle {I_{\textrm{out}}} \rangle = P \pi ^{2} d_1^{2} d_2^{2} / 64 L^{6}$. This in turn implies

(51)$$\begin{aligned}\textrm{var}(\hat{{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}(\boldsymbol{\rho})) = &\frac{\langle {I_{\textrm{out}}} \rangle^{2}}{4\pi^{2}R^{4}} \int\!\textrm{d}^{2}{\boldsymbol{\rho}_{\textrm{out}}}\, \int\!\textrm{d}^{2}{\boldsymbol{\tilde{\rho}}_{\textrm{out}}}\, e^{-(\left| {\boldsymbol{\rho}-\boldsymbol{\rho}_{\textrm{out}}} \right|^{2} + \left| {\boldsymbol{\rho}-\boldsymbol{\tilde{\rho}}_{\textrm{out}}} \right|^{2})/2R^{2}}\\ &\times e^{-(\omega_0^{2}+(\omega_0 + {{\Delta\omega}})^{2}) d_2^{2}\left| {\boldsymbol{\rho}_{\textrm{out}}-\boldsymbol{\tilde{\rho}}_{\textrm{out}}} \right|^{2} / 32 c^{2} L^{2}} e^{-i{{\Delta\omega}}(\left| {\boldsymbol{\rho}_{\textrm{out}}} \right|^{2}-\left| {\boldsymbol{\tilde{\rho}}_{\textrm{out}}} \right|^{2}) / 2 c L} \end{aligned}$$

(52)$$ = \frac{\langle {I_{\textrm{out}}} \rangle^{2} \exp\left( {-{{\Delta\omega}}^{2} R^{2} \left| {\boldsymbol{\rho}} \right|^{2} / c^{2} L^{2} Z(R)} \right)}{Z(R)}, $$

where

(53)$$Z(R)\equiv 1 + \frac{{{\Delta\omega}}^{2} R^{4}}{c^{2} L^{2}} + \frac{d_2^{2} R^{2}}{8 c^{2} L^{2}} (\omega_0^{2}+\omega_1^{2}).$$

Equation (52) takes its maximum value at $\boldsymbol {\rho }={{\mathbf 0}}$, so we limit our attention to that bound. It is not hard to see that $|\langle {\hat {{\mathcal {P}}}_{\textrm{out},{{\Delta \omega }}}(\boldsymbol {\rho })} \rangle |^{2} = \langle {I_3} \rangle ^{2}$, and so the ratio of the squared mean to the variance is bounded below by

(54)$$\frac{|\langle {\hat{{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}(\boldsymbol{\rho})} \rangle|^{2}}{\textrm{var}(\hat{{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}(\boldsymbol{\rho}))} \geq \frac{|\langle {\hat{{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}(\boldsymbol{\rho})} \rangle|^{2}}{\textrm{var}(\hat{{\mathcal{P}}}_{\textrm{out},{{\Delta\omega}}}({\mathbf 0}))} = Z(R).$$

It is clear that the variance vanishes as $R$ grows without bound. However, we would like the variance to be significantly less than the squared mean for a small enough value of $R$ that the ${{\mathcal {P}}}$ field’s spatial detail is preserved. Taking reasonable parameter values $\lambda _0=532\,\textrm {nm}$, $L=1\,\textrm {m}$, $d_2=2\,\textrm {m}$, we ambitiously take the difference frequency to be ${{\Delta \omega }} / 2\pi = 1\,\textrm {THz}$ so that $\Delta \lambda =300\,{\mathrm{\mu} \textrm {m}}$. Even setting taking our Gaussian kernel to have $R=1\,{\mathrm{\mu} \textrm {m}}$, we find $Z(1\,{\mathrm{\mu} \textrm {m}})\approx 141$, implying that the squared mean greatly exceeds the variance, and thus $\hat {{\mathcal {P}}}_{\textrm{out},{{\Delta \omega }}}(\boldsymbol {\rho })\approx {{\mathcal {P}}}_{\textrm{out},{{\Delta \omega }}}(\boldsymbol {\rho })$. Of course, we have traded the need for exceedingly good time resolution in detection for exceedingly good spatial resolution, as sub-micron resolution is required to avoid speckle averaging at the detector arrays. A possible approach to realizing such arrays is to combine magnifying optics with lock-in cameras [27], as used by Willomitzer et al. [20] in their proofs-of-concept NLoS SWH experiments.

One more point remains to be dealt with, i.e., reconciling Willomitzer et al.’s nonstatistical approach to SWH imaging with the statistics-based ${{\mathcal {P}}}$-field version of SWH imaging that we have just described. It turns out this will not be hard to do. Our spatial low-pass filter in Fig. 6(b) averages out speckle so that $\hat {{{\mathcal {P}}}}_{\textrm{out},{{\Delta \omega }}}(\boldsymbol {\rho })$ converges to ${{\mathcal {P}}}_{\textrm{out},{{\Delta \omega }}}(\boldsymbol {\rho })$ without unduly compromising spatial resolution. Compare that to the spatial low-pass filtering that Willomitzer et al. report doing—in the discussion below Eq. (35) of their paper’s supplementary material—to suppress “parasitic interference”. Their parasitic interference is our speckle, making their (nonstatistical) latent synthetic-wavelength hologram very much like our diffuser-averaged $\hat {{{\mathcal {P}}}}$ field.

Funding

Defense Advanced Research Projects Agency (HR0011-16-C-0030).

Acknowledgments

The authors acknowledge fruitful interactions with members of the DARPA REVEAL teams from the University of Wisconsin and Southern Methodist University. They also thank Dr. Ravi Athale for organizing a valuable workshop on phasor-field imaging, and Dr. Jeremy Teichman for sharing an early version of his analysis of speckle effects in ${{\mathcal {P}}}$-field imaging.

Disclosures

The authors declare no conflicts of interest.

References

1. A. Kirmani, T. Hutchison, J. Davis, and R. Raskar, “Looking around the corner using ultrafast transient imaging,” Int. J. Comput. Vis. 95(1), 13–28 (2011). [CrossRef]  

2. A. Velten, T. Willwacher, O. Gupta, A. Veeraraghavan, M. G. Bawendi, and R. Raskar, “Recovering three-dimensional shape around a corner using ultrafast time-of-flight imaging,” Nat. Commun. 3(1), 745 (2012). [CrossRef]  

3. M. O’Toole, D. B. Lindell, and G. Wetzstein, “Confocal non-line-of-sight imaging based on the light-cone transform,” Nature 555(7696), 338–341 (2018). [CrossRef]  

4. D. B. Lindell, G. Wetzstein, and M. O’Toole, “Wave-based non-line-of-sight imaging using fast f-k migration,” ACM Trans. Graph. 38(4), 1–13 (2019). [CrossRef]  

5. F. Xu, G. Shulkind, C. Thrampoulidis, J. H. Shapiro, A. Torralba, F. N. C. Wong, and G. W. Wornell, “Revealing hidden scenes by photon-efficient occlusion-based opportunistic active imaging,” Opt. Express 26(8), 9945–9962 (2018). [CrossRef]  

6. C. Thrampoulidis, G. Shulkind, F. Xu, W. T. Freeman, J. H. Shapiro, A. Torralba, F. N. C. Wong, and G. W. Wornell, “Exploiting occlusion in non-line-of-sight active imaging,” IEEE Trans. Comput. Imaging 4(3), 419–431 (2018). [CrossRef]  

7. S. A. Reza, M. La Manna, and A. Velten, “A physical light transport model for non-line-of-sight imaging applications,” arXiv:1802.1823 [physics.optics].

8. J. Dove and J. H. Shapiro, “Paraxial theory of phasor-field imaging,” Opt. Express 27(13), 18016–18037 (2019). [CrossRef]  

9. J. A. Teichman, “Phasor field waves: a mathematical treatment,” Opt. Express 27(20), 27500–27506 (2019). [CrossRef]  

10. X. Liu, I. Guillén, M. La Manna, J. H. Nam, S. A. Reza, T. H. Le, A. Jarabo, D. Gutierrez, and A. Velten, “Non-line-of-sight imaging using phasor-field virtual wave optics,” Nature 572(7771), 620–623 (2019). [CrossRef]  

11. X. Liu, S. Bauer, and A. Velten, “Phasor field diffraction based reconstruction for fast non-line-of-sight imaging systems,” Nat. Commun. 11(1), 1645 (2020). [CrossRef]  

12. S. A. Reza, M. La Manna, S. Bauer, and A. Velten, “Phasor field waves: experimental demonstrations of wave-like properties,” Opt. Express 27(22), 32587–32608 (2019). [CrossRef]  

13. Paraxial means that Fresnel diffraction can be used in lieu of the more exact Rayleigh–Sommerfeld diffraction. A sufficient condition at wavelength λ for this to be so is that the longitudinal distance L between the input and output planes and the diameters d₀ and d_L of the transverse regions of interest within those planes satisfy $\bar {d}^{4}/8\lambda L^{3} \ll 1$, where $\bar {d} \equiv (d_0+d_L)/2$. An alternative condition is that the angular spectrum of the input field be confined to a paraxial region about the optical axis.

14. J. Dove and J. H. Shapiro, “Speckled speckled speckle,” Opt. Express, submitted, arXiv:2005.10165 [physics.optics].

15. J. Dove, “Theory of phasor-field imaging,” Ph.D. thesis, Massachusetts Institute of Technology (2020).

16. As shown in Ref. [14], this will be possible when the image-plane detector’s pixels are appreciably smaller than the spatial detail in 〈I₂(ρ₂, t)〉 and appreciably larger than the speckle size in that plane.

17. The focusing capability of a circular-pupil lens is limited by its diameter D. Our assumption of Gaussian pupil—here and in all that follows—enables us to obtain closed-form results while still capturing the essential physics of a finite pupil. In essence, this Gaussian pupil acts as an apodization that suppresses the Fresnel and Fraunhofer diffraction rings which result when a plane wave diffracts through a diameter-D circular pupil while preserving the transition from the near field to the far field.

18. Physically, $|\textrm {Re}[\mathcal {A}(\omega )]| \ll |\textrm {Im}[\mathcal {A}(\omega )]|$ is a near-field condition, i.e., it implies that D is much greater than the wavelength-λ₀ diffraction scale at propagation distance |1/L_in − 1/f|⁻¹ for a diameter-D collimated beam.

19. The assumption that $|\textrm {Re}[\mathcal {B}(\omega )]| \ll |\textrm {Im}[\mathcal {B}(\omega )]|$ is a near-field condition similar to the one explained in Ref. [18].

20. F. Willomitzer, P. V. Rangarajan, F. Li, M. M. Balaji, M. P. Christensen, and O. Cossairt, “Synthetic wavelength holography: an extension of Gabor’s holographic principle to imaging with scattered wavefronts,” arXiv:1912.11438 [physics.optics].

21. M. Saxena, G. Eluru, and S. S. Gorthi, “Structured illumination microscopy,” Adv. Opt. Photonics 7(2), 241–275 (2015). [CrossRef]  

22. P. Sen, B. Chen, G. Garg, S. R. Marschner, M. Horowitz, M. Levoy, and H. P. A. Lensch, “Dual photography,” ACM Trans. Graph. 24(3), 745–755 (2005). [CrossRef]  

23. This input diffuser is not part of the Fig. 1 scenario, i.e., it is not the transmissive analog of any of the walls in three- bounce NLoS imaging. Instead, it is a user-supplied diffuser whose purpose is to create the initial spatial incoherence that permits ${{\mathcal {P}}}$-field physical-optics projection over a distance $L_{\textrm{in}_1}\;+ \;L_{\textrm{in}_2}\;+\;L_1$ using an ordinary lens whose focal length satisfies the lens-maker's formula $1/f = 1/L_{\textrm{in}_1} + 1/(L_{\textrm{in}_2}+L_1)$.

24. J. Dove and J. H. Shapiro, “Nonparaxial phasor-field propagation,” Opt. Express, submitted, arXiv:2006.13775 [physics.optics].

25. J. H. Shapiro, “The Quantum Theory of Optical Communications,” IEEE J. Sel. Top. Quantum Electron. 15(6), 1547–1569 (2009). [CrossRef]  

26. Willomitzer et al.'s SWH implementation uses an unbalanced configuration in which only one output arm from the final beam splitter is detected, as shown in Fig. 6(a). We prefer to analyze the balanced case, in which both arms are detected, because that arrangement suppresses any excess noise that may be present on the local-oscillator beams without changing the essential physics of the measurement.

27. S. Foix, G. Alenya, and C. Torras, “Lock-in time-of-flight (ToF) cameras: A survey,” IEEE Sens. J. 11(9), 1917–1926 (2011). [CrossRef]