Paraxial theory of phasor-field imaging

Justin Dove; Jeffrey H. Shapiro

doi:10.1364/OE.27.018016

1. 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. [1] 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. [2] 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. [3] who applied similar techniques with success to data collected by relatively inexpensive photonic-mixer-device (PMD) time-of-flight sensors. Buttafava et al. [4] 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. [5] 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. [6] 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. [7] 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. [8] 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. [9] recently introduced the phasor-field (𝒫-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 𝒫 field is an insufficient representation of the underlying field at the site of such features. Nevertheless, Liu et al. [10] used the 𝒫-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 𝒫-field propagation to scenarios more general than purely Lambertian reflections. Very recently, Reza et al. [11] reported an elegant series of experiments that verify the 𝒫 field’s legitimacy. These experiments clearly demonstrate the 𝒫 field’s wave-like properties, which offer the possibility of NLoS imaging without the need for computational reconstructions by using a 𝒫-field lens instead.

The success of Liu et al.’s experiments is impressive, and Reza et al.’s 𝒫-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 𝒫 field—already provide such a transport model, but the aforementioned works have demonstrated the intuitive utility of the 𝒫-field approach. Consequently, we believe it is worthwhile to pursue propagation primitives that can readily establish the 𝒫-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 𝒫-field formalism to scenarios that go beyond what was considered in [9] 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 [12,13] and use scalar-wave, paraxial optics although these restrictions are not essential. In Sec. 2 we present our own development and analysis of the 𝒫-field notion. We begin by tracing light propagation through an example transmissive geometry wherein a natural definition for the 𝒫 field presents itself. Continuing this analysis, we arrive at a paraxial 𝒫-field propagator analogous to that reported by Reza et al. [9]. Using this result, we analyze the performance of 𝒫-field imaging for unoccluded transmissive geometries. Next, moving beyond the 𝒫 field, in Sec. 3 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 𝒫-field input-output relation for occlusion-aided, diffuse-object, transmissive imaging. With that analysis in hand, we compare the 𝒫-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. 4 we summarize our results and consider directions for further research.

2. 𝒫-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. [9], we define the 𝒫 field as the Fourier transform of the short-time average irradiance. Using this definition, we derive a formula for paraxial propagation of the 𝒫 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 [9]. We then apply this understanding of the 𝒫 field to the task of imaging through diffusers and analyze the associated performance.

2.1. Setup for paraxial propagation through multiple diffusers

Figure 1 shows the transmissive geometry we shall address in this paper for 𝒫-field propagation within the paraxial regime, i.e., wherein Fresnel diffraction applies. Here, E₀(ρ₀, 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, ρ₀ = (x₀, y₀), and time, t. This field has center frequency ω₀ and bandwidth Δω ≪ ω₀, so that the optical-frequency field is Re[E₀(ρ₀, t)e^−iω_ot]. Its units are $\sqrt{W / m^{2}}$ , making I₀(ρ₀, t) = |E₀(ρ₀, t)|² the short-time average irradiance [14] illuminating the z = 0 plane. It will be assumed, in all that follows, that Δω is such that available photodetectors can fully resolve the time dependence of I₀(ρ₀, t). As soon will be seen, it will be valuable to employ the time-domain Fourier transform of E₀(ρ₀, t), viz. [15],

ℰ_{0} (ρ_{0}, ω) \equiv \int d t E_{0} (ρ_{0}, t) e^{i ω t},

for use analyzing the Fig. 1 configuration.

Fig. 1 Unfolded geometry for three-bounce NLoS active imaging. Scalar, paraxial diffraction theory is assumed, with {E_k(ρ_k, t) : 0 ≤ k ≤ 2} being the baseband complex-field envelopes illuminating the z = 0, z = L₁, and z = L₁ + L₂ planes, respectively, written as functions of the transverse spatial coordinates, {ρ_k = (x_k, y_k) : 0 ≤ k ≤ 2}, in those planes and time, t. The blue rectangles represent thin transmissive diffusers, and the black line represents a thin transmission screen whose intensity transmission pattern, T(ρ₁), is to be imaged using the light that emerges from the z = L₁ + L₂ plane.

Download Full Size | PDF

After propagating through the thin diffuser h₀(ρ₀), the Fourier-domain field at z = 0₊ is

ℰ_{0}^{'} (ρ_{0}, ω) = ℰ_{0} (ρ_{0}, ω) exp [i (ω_{0} + ω) h_{0} (ρ_{0}) / c],

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₀(ρ₀)/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₀(ρ₀) to be a zero-mean, homogeneous, isotropic, Gaussian random function of ρ₀, with covariance function K_h(|Δρ|) = 〈h₀(ρ₀ + Δρ)h₀(ρ₀)〉, where angle brackets denote ensemble average. Moreover, in keeping with h₀(ρ₀)’s being a diffuser, we shall take its standard deviation,

σ_{h} = \sqrt{K_{h} (0)}

to be much greater than the center wavelength, λ₀ = 2πc/ω₀, and its coherence length ρ_c—the transverse distance beyond which K_h(|Δρ|) vanishes—to be at most a few λ₀. Furthermore—and this condition is essential to there being a useful 𝒫-field propagator—we shall assume that σ_h is much smaller than the wavelength of the modulation bandwidth, Δλ = 2πc/Δω.

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

ℰ_{1} (ρ_{1}, ω) = \int d^{2} ρ_{0} ℰ_{0}^{'} (ρ_{0}, ω) \frac{exp [i (ω_{0} + ω) L_{1} / c + i (ω_{0} + ω) {| (ρ_{1} - ρ_{0}) |}^{2} / 2 c L_{1}] (ω_{0} + ω)}{i 2 π c L_{1}},

is the time-domain Fourier transform of E₁(ρ₁, t), the field illuminating the z = L₁ plane. This illumination results in

ℰ_{1}^{'} (ρ_{1}, ω) = ℰ_{1}^{'} (ρ_{1}, ω) \sqrt{T (ρ_{1})} exp [i (ω_{0} + ω) h_{1} (ρ_{1}) / c],

being the time-domain Fourier transform of E′₁(ρ₁, t), the field that emerges at z = L_1₊, after propagation through a deterministic thin transmission screen with intensity transmission pattern T(ρ₁), and a thin diffuser, h₁(ρ₁), that we will take to be statistically independent of, but identically distributed as, h₀(ρ₀).

Paraxial propagation to z = L₁ + L₂, now gives us

ℰ_{2} (ρ_{2}, ω) = \int d^{2} ρ_{1} ℰ_{1}^{'} (ρ_{1}, ω) \frac{exp [i (ω_{0} + ω) L_{2} / c + i (ω_{0} + ω) {| (ρ_{2} - ρ_{1}) |}^{2} / 2 c L_{2}] (ω_{0} + ω)}{i 2 π c L_{2}},

and propagation through the thin diffuser at z = L₁ + L₂ results in

ℰ_{2}^{'} (ρ_{2}, ω) = ℰ_{2} (ρ_{2}, ω) exp [i (ω_{0} + ω) h_{2} (ρ_{2}) / c],

being the time-domain Fourier transform of E′₂(ρ₂, t), the field that emerges at z = (L₁ + L₂)₊. We will assume that h₂(ρ₂) is statistically independent of, but identically distributed as, h₀(ρ₀) and h₁(ρ₁).

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(ρ₁) and the z = L₁ diffuser represent a Lambertian-reflecting hidden wall with spatially-varying albedo T(ρ₁). The z = L₁ + L₂ 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(ρ₁). In the next section, we will derive the 𝒫-field propagator for the preceding transmission geometry.

2.2. 𝒫-field propagator in the paraxial regime

To start our derivation, consider 〈I₁(ρ₁, t)〉, where I₁(ρ₁, t) ≡ |E₁(ρ₁, t)|² is the short-time average irradiance illuminating the z = L₁ plane and angle brackets denote averaging over the statistics of h₀(ρ₀). Going to the temporal-frequency domain, we have that

〈 I_{1} (ρ_{1}, t) 〉 = \int \frac{d ω}{2 π} \int \frac{d ω^{'}}{2 π} 〈 ℰ_{1} (ρ_{1}, ω) ℰ_{1}^{*} (ρ_{1}, ω^{'}) 〉 e^{- i (ω - ω^{'}) t}

= \int \frac{d ω_{-}}{2 π} [\int \frac{d ω_{+}}{2 π} 〈 ℰ_{1} (ρ_{1}, ω_{+} + ω_{-} / 2) ℰ_{1}^{*} (ρ_{1}, ω_{+} - ω_{-} / 2) 〉] e^{- i ω_{-} t}

= \int \frac{d ω_{-}}{2 π} 𝒫_{1} (ρ_{1}, ω_{-}) e^{- i ω_{-} t},

where ^* denotes complex conjugate, ω₊ ≡ (ω + ω′)/2, ω₋ ≡ ω − ω′, and we have introduced the 𝒫 field at the z = L₁ plane as the Fourier transform of 〈I₁(ρ₁, t)〉. Next, employing Eqs. (2) and (3), we get

\begin{array}{l} 𝒫_{1} (ρ_{1}, ω_{-}) = \int d^{2} ρ_{0} \int d^{2} ρ_{0}^{'} \int \frac{d ω_{+}}{2 π} ℰ_{0} (ρ_{0}, ω) ℰ_{0}^{*} (ρ_{0}^{'}, ω^{'}) 〈 e^{i [(ω_{0} + ω) h_{0} (ρ_{0}) - (ω_{0} + ω^{'}) h_{0} (ρ_{0}^{'})] / c} 〉 \\ \times (ω_{0} + ω) (ω_{0} + ω^{'}) e^{i (ω - ω^{'}) L_{1} / c + i [(ω_{0} + ω) {| ρ_{1} - ρ_{0} |}^{2} - (ω_{0} + ω^{'}) {| ρ_{1} - ρ_{0}^{'} |}^{2}] / 2 c L_{1}} / {(2 π c L_{1})}^{2}, \end{array}

where, as before, ω₊ ≡ (ω + ω′)/2 and ω₋ ≡ ω − ω′. Because Δω ≪ ω₀ and σ_h ≪ Δλ, the preceding result can be reduced to

\begin{array}{l} 𝒫_{1} (ρ_{1}, ω_{-}) & = \int d^{2} ρ_{0} \int d^{2} ρ_{0}^{'} \int \frac{d ω_{+}}{2 π} ℰ_{0} (ρ_{0}, ω) ℰ_{0}^{*} (ρ_{0}^{'}, ω^{'}) 〈 e^{i ω_{0} [h_{0} (ρ_{0}) - h_{0} (ρ_{0}^{'})] / c} 〉 ω_{0}^{2} / {(2 π c L_{1})}^{2} \\ \times e^{i (ω - ω^{'}) L_{1} / c + i [(ω_{0} + ω) {| ρ_{1} - ρ_{0} |}^{2} - (ω_{0} + ω^{'}) {| ρ_{1} - ρ_{0}^{'} |}^{2}] / 2 c L_{1}} . \end{array}

A standard result for Gaussian random functions gives us [16]

〈 e^{i ω_{0} [h_{0} (ρ_{0}) - h_{0} (ρ_{0}^{'})] / c} 〉 = exp {- ω_{0}^{2} [σ_{h}^{2} - K_{h} (| ρ_{0} - ρ_{0}^{'} |)] / c^{2}} .

Then, because σ_h ≫ λ₀ and ρ_c ∼ λ₀ we can use an impulse approximation, viz.,

〈 e^{i ω_{0} [h_{0} (ρ_{0}) - h_{0} (ρ_{0}^{'})] / c} 〉 \approx λ_{0}^{2} δ (ρ_{0} - ρ_{0}^{'}),

in Eq. (11) to obtain

𝒫_{1} (ρ_{1}, ω_{-}) = \int d^{2} ρ_{0} \int \frac{d ω_{+}}{2 π} ℰ_{0} (ρ_{0}, ω) ℰ_{0}^{*} (ρ_{0}, ω^{'}) e^{i (ω - ω^{'}) L_{1} / c + i (ω - ω^{'}) {| ρ_{1} - ρ_{0} |}^{2} / 2 c L_{1}} / L_{1}^{2} .

= \int d^{2} ρ_{0} 𝒫 (ρ_{0}, ω_{-}) e^{i ω_{-} L_{1} / c + i ω_{-} {| ρ_{1} - ρ_{0} |}^{2} / 2 c L_{1}} / L_{1}^{2} .

Here, the 𝒫 field at z = 0 is

𝒫_{0} (ρ_{0}, ω_{-}) = \int \frac{d ω_{+}}{2 π} ℰ_{0} (ρ_{0}, ω_{+} + ω_{-} / 2) ℰ_{0}^{*} (ρ_{0}, ω_{+} - ω_{-} / 2),

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 𝒫-field propagator [9]—is our essential result for paraxial 𝒫-field propagation over a distance L₁. 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 𝒫-field propagation at frequency ω₋ over a distance L₁ that is governed by a modified version of the ℰ-field’s Fresnel-diffraction formula, viz., one in which the exponent’s optical frequency in the ℰ-field Fresnel formula is replaced by the 𝒫 field’s modulation frequency and the ℰ-field formula’s ω₀/i2πcL₁ factor is replaced by the 𝒫 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₁ plane is governed by

〈 I_{1} (ρ_{1}, t) 〉 = \int d^{2} ρ_{0} I_{0} (ρ_{0}, t - L_{1} / c - {| ρ_{1} - ρ_{0} |}^{2} / 2 c L_{1}) / L_{1}^{2},

which has the following pleasing physical interpretation: Paraxial propagation of the short-time average irradiance from the diffuser’s output to the z = L₁ presumes that

\frac{exp [i ω \sqrt{L_{1}^{2} + {| ρ_{1} - ρ_{0} |}^{2} / c}]}{\sqrt{L_{1}^{2} + {| ρ_{1} - ρ_{0} |}^{2}}} \approx \frac{exp (i ω L_{1} / c + i ω {| ρ_{1} - ρ_{0} |}^{2} / 2 c L_{1})}{L_{1}}, for | ω | \leq Δ ω

can be employed, and results in 〈I₁(ρ₁, t)〉 being governed by the paraxial form of geometric optics, viz., the differential contribution of I₀(ρ₀, t) to 〈I₁(ρ₁, t)〉 is time delayed by L₁/c+|ρ₁−ρ₀|²/2cL₁ and attenuated by the inverse-square-law factor

1 / L_{1}^{2}

.

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

𝒫_{2} (ρ_{2}, ω_{-}) \equiv \int \frac{d ω_{+}}{2 π} 〈 ℰ_{2} (ρ_{2}, ω_{+} + ω_{-} / 2) ℰ_{2}^{*} (ρ_{2}, ω_{+} - ω_{-} / 2) 〉

= \int d^{2} ρ_{1} 𝒫_{1} (ρ_{1}, ω_{-}) T (ρ_{1}) exp (i ω_{-} L_{2} / c + i ω_{-} {| ρ_{2} - ρ_{1} |}^{2} / 2 c L_{2}) / L_{2}^{2},

where the averaging brackets in Eq. (19) represent averaging over the h₀(ρ₀) and the h₁(ρ₁) ensembles. Combining this result with what we have already obtained for relating 𝒫₁(ρ₁, ω₋) to 𝒫₀(ρ₀, ω₋) we get

\begin{array}{l} 𝒫_{2} (ρ_{2}, ω_{-}) & = \int d^{2} ρ_{1} (\int d^{2} ρ_{0} 𝒫_{0} (ρ_{0}, ω_{-}) exp (i ω_{-} L_{1} / c + i ω_{-} {| ρ_{1} - ρ_{0} |}^{2} / 2 c L_{1}) / L_{1}^{2}) \\ \times T (ρ_{1}) exp (i ω_{-} L_{2} / c + i ω_{-} {| ρ_{2} - ρ_{1} |}^{2} / 2 c L_{2}) / L_{2}^{2} . \end{array}

Before continuing, it is crucial to note the behavior of 𝒫₂(ρ₂, 0). From Eq. (21) we immediately find that

𝒫_{2} (ρ_{2}, 0) = \int d^{2} ρ_{1} T (ρ_{1}) \int d^{2} ρ_{0} 𝒫_{0} (ρ_{0}, 0) / {(L_{1} L_{2})}^{2},

indicating that there is no spatial information about T(ρ₁) available in 𝒫₂(ρ₂, 0). This behavior is a consequence of using the paraxial approximation. Going beyond the paraxial-propagation regime—to Rayleigh-Sommerfeld diffraction—will yield a 𝒫₂(ρ₂, 0) containing some spatial information about T(ρ₁), but the inverse problem for recovering T(ρ₁) from 𝒫₂(ρ₂, 0) will still be poorly conditioned in the Fig. 1 configuration. This behavior has been seen by Xu et al. [12] and Thrampoulidis et al. [13] 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.

2.3. T(ρ₁) reconstruction in the paraxial regime 𝒫-field formalism

Equation (21) shows that the intensity transmission pattern, T(ρ₁), we wish to reconstruct is illuminated by 𝒫₁(ρ₁, ω₋), the 𝒫 field that results from propagation of the laser illumination’s 𝒫₀(ρ₀, ω₋) from z = 0 to z = L₁. After transmission through T(ρ₁) and the diffuser h₁(ρ₁), 𝒫-field propagation from to z = L₁ + L₂ results in 𝒫₂(ρ₂, ω₋), 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(ρ₁).

Fig. 2 Thin-lens imaging setup. A focal-length f thin lens casts an inverted image of the intensity pattern that emerges from the diffuser at z = L₁ + L₂. The image is located in the plane—shown as a black dashed line—a distance L_im behind the lens, where 1/f = 1/L₃ + 1/L_im.

Download Full Size | PDF

Let E′₂(ρ₂, t) be the baseband, complex-field envelope emerging from the diffuser in the z = L₁ + L₂ plane, and let ℰ′₂(ρ₂, ω) be its time-domain Fourier transform. After Fresnel propagation from z = L₁ + L₂ to z = L₁ + L₂ + L₃, propagation through the diameter-D circular-pupil, focal-length-f, thin lens, and Fresnel propagation over an additional L_im distance where 1/f = 1/L₃ + 1/L_im, the resulting image-plane field E_im(ρ, t) has time-domain Fourier transform given by

\begin{array}{l} ℰ_{im} (ρ_{im}, ω) & = \int_{| ρ_{3} | \leq D / 2} d^{2} ρ_{3} \frac{e^{i (ω_{0} + ω) L_{im} / c + i (ω_{0} + ω) {| ρ_{im} - ρ_{3} |}^{2} / 2 c L_{im} - i (ω_{0} + ω) {| ρ_{3} |}^{2} / 2 c f}}{i λ_{0} L_{im}} \\ \times \int d^{2} ρ_{2} ℰ_{2}^{'} (ρ_{2}, ω) \frac{e^{i (ω_{0} + ω) L_{3} / c + i (ω_{0} + ω) {| ρ_{3} - ρ_{2} |}^{2} / 2 c L_{3}}}{i λ_{0} L_{3}} \end{array}

\begin{array}{l} = e^{i (ω_{0} + ω) {| ρ_{im} |}^{2} / 2 c L_{im}} \int d^{2} ρ_{2} ℰ_{2}^{'} (ρ_{2}, ω) \frac{e^{i (ω_{0} + ω) (L_{3} + L_{im}) / c + i (ω_{0} + ω) {| ρ_{2} |}^{2} / 2 c L_{3}}}{i λ_{0} L_{3}} \\ \times \int_{| ρ_{3} | \leq D / 2} d^{2} ρ_{3} \frac{e^{- i (ω + ω_{0}) ρ_{3} \cdot (ρ_{2} / L_{3} + ρ_{im} / L_{im}) / c}}{i λ_{0} L_{im}} . \end{array}

Performing the integration over ρ₃ results in

\begin{array}{l} ℰ_{im} (ρ_{im}, ω) = e^{i (ω_{0} + ω) {| ρ_{im} |}^{2} / 2 c L_{im}} \\ \times \int d^{2} ρ_{2} ℰ_{2}^{'} (ρ_{2}, ω) \frac{e^{i (ω_{0} + ω) (L_{3} + L_{im}) / c + i (ω_{0} + ω) {| ρ_{2} |}^{2} / 2 c L_{3}}}{- λ_{0}^{2} L_{3} L_{im}} \frac{π D^{2}}{4} \frac{J_{1} (\frac{π D}{λ_{0}} | \frac{ρ_{2}}{L_{3}} + \frac{ρ_{im}}{L_{im}} |)}{\frac{π D}{2 λ_{0}} | \frac{ρ_{2}}{L_{3}} + \frac{ρ_{im}}{L_{im}} |}, \end{array}

where J₁(·) is the first-order Bessel function of the first kind, and we have used πD/λ₀ in lieu of (ω₀ + ω)D/2c in the Airy pattern because Δω ≪ ω₀.

The presence of the diffuser h₂(ρ₂) makes

〈 ℰ_{2}^{'} (ρ_{2}, ω) ℰ_{2}^{' *} (ρ_{2}^{'}, ω^{'}) 〉 \approx λ_{0}^{2} 〈 ℰ_{2} (ρ_{2}, ω) ℰ_{2}^{*} (ρ_{2}, ω^{'}) 〉 δ (ρ_{2} - ρ_{2}^{'}),

which together with Eq. (25) yields

\begin{array}{l} 𝒫_{im} (ρ_{im}, ω_{-}) = \int d^{2} ρ_{2} 𝒫_{2} (ρ_{2}, ω_{-}) \\ \times e^{i ω_{-} (L_{3} + L_{im}) / c + i ω_{-} {| ρ_{2} |}^{2} / 2 c L_{3} + i ω_{-} {| ρ_{im} |}^{2} / 2 c L_{im}} {[\frac{π D^{2}}{4 λ_{0} L_{3} L_{im}} \frac{J_{1} (\frac{π D}{λ_{0}} | \frac{ρ_{2}}{L_{3}} + \frac{ρ_{im}}{L_{im}} |)}{\frac{π D}{2 λ_{0}} | \frac{ρ_{2}}{L_{3}} + \frac{ρ_{im}}{L_{im}} |}]}^{2} . \end{array}

and hence

\begin{array}{l} 〈 I_{im} (ρ_{im}, t) 〉 & = \int d^{2} ρ_{2} 〈 I_{2} (ρ_{2}, t - (L_{3} + L_{im}) / c - {| ρ_{2} |}^{2} / 2 c L_{3} - {| ρ_{im} |}^{2} / 2 c L_{im}) 〉 \\ \times {[\frac{π D^{2}}{4 λ_{0} L_{3} L_{im}} \frac{J_{1} (\frac{π D}{λ_{0}} | \frac{ρ_{2}}{L_{3}} + \frac{ρ_{im}}{L_{im}} |)}{\frac{π D}{2 λ_{0}} | \frac{ρ_{2}}{L_{3}} + \frac{ρ_{im}}{L_{im}} |}]}^{2} . \end{array}

So, by measuring 〈I_im(ρ_im, t)〉, i.e., the diffuser-averaged, short-time average, image-plane irradiance, we obtain a 1.22λ₀/D-angular-resolution, image of 〈I₂(ρ₂, t − (L₃ + L_im)/c − |ρ₂|²/2cL₃)〉. From that irradiance image we can then compute a 1.22λ₀/D-angular-resolution image of 𝒫₂(ρ₂, ω₋) at any modulation frequency of interest.

For reconstructing T(ρ₁), let us suppose that the z = 0 illumination is a duration t₀, cosinusoidally-modulated, collimated Gaussian-beam laser field where Δωt₀ ≫ 1, i.e.,

E_{0} (ρ_{0}, t) = {\begin{array}{l} \sqrt{\frac{8 P_{0}}{π d^{2}}} e^{- 4 {| ρ_{0} |}^{2} / d^{2}} cos (Δ ω t / 2), & for | t | \leq t_{0} / 2, \\ 0, & otherwise, \end{array}

with P₀t₀/2 being the energy illuminating the z = 0 plane. This field’s short-time average irradiance is then

I_{0} (ρ_{0}, t) = {\begin{array}{l} \frac{8 P_{0}}{π d^{2}} e^{- 8 | ρ_{0}^{2} / d^{2}} {cos}^{2} (Δ ω t / 2) = \frac{4 P_{0}}{π d^{2}} e^{- 8 | ρ_{0}^{2} / d^{2}} [1 + cos (Δ ω t)], & for | t | \leq t_{0} / 2, \\ 0, & otherwise, \end{array}

which leads to

𝒫_{0} (ρ_{0}, ω_{-}) = \frac{8 P_{0} t_{0}}{π d^{2}} e^{- 8 | ρ_{0} |^{2} / d^{2}} [\frac{sin (ω_{-} t_{0} / 2)}{ω_{-} t_{0} / 2} + \frac{sin [(ω_{-} + Δ ω) t_{0} / 2]}{(ω_{-} + Δ ω) t_{0}} + \frac{sin [(ω_{-} - Δ ω) t_{0} / 2]}{(ω_{-} - Δ ω) t_{0}}],

and hence

𝒫_{1} (ρ_{1}, Δ ω) \approx \int d^{2} ρ_{0} \frac{4 P_{0} t_{0}}{π d^{2}} e^{- 8 | ρ_{0} |^{2} / d^{2}} \frac{exp (i Δ ω L_{1} / c + i Δ ω | ρ_{1} - ρ_{0} |^{2} / 2 c L_{1})}{L_{1}^{2}},

because Δωt₀ ≫ 1. Although this expression can be evaluated analytically, we shall not bother. We just note that with Δω/2π ∼ 1 GHz, d ∼ 1 mm, and L₁ ∼ 1 m, we have cL₁/Δωd² ≫ 1 from which it follows that the spatial extent of 𝒫₁(ρ₁, Δω) will be ∼cL₁/Δωd ≫ d. In other words, the effect of the diffuser h₀(ρ₀) is to ensure that a finite, but much larger than diameter-d, region of the z = L₁ plane is illuminated by the frequency-Δω 𝒫 field.

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

{\tilde{𝒫}}_{2} (ρ_{2}, Δ ω) \equiv {(L_{im} / L_{3})}^{2} 𝒫_{im} (- ρ_{2} L_{im} / L_{3}, Δ ω) e^{- i Δ ω (L_{3} + L_{im}) / c - i Δ ω | ρ_{2} |^{2} / 2 c L_{3} - i ω_{-} | ρ_{im} |^{2} / 2 c L_{im}},

of 𝒫₂(ρ₂, Δω) from the 〈I_im(ρ_im, t)〉 measurement. We can computationally invert Eq. (20) to obtain a reconstruction of T(ρ₁)𝒫₁(ρ₁, Δω) and use our knowledge of 𝒫₁(ρ₁, Δω) to obtain a T(ρ₁) image. In particular, suppose we measure 〈I_im(ρ_im, t)〉 for |ρ_im| ≤ d_im/2, and then define T̃(ρ̃₁) by

\tilde{T} ({\tilde{ρ}}_{1}) | 𝒫_{1} ({\tilde{ρ}}_{1}, Δ ω) | = | \int_{| ρ_{2} | \leq D' / 2} d^{2} ρ_{2} {\tilde{𝒫}}_{2} (ρ_{2}, Δ ω) \frac{e^{- i Δ ω | ρ_{2} |^{2} / 2 c L_{2} + i Δ ω ρ_{2} \cdot {\tilde{ρ}}_{1} / c L_{2}}}{Δ λ^{2}} |,

where D′ ≡ d_imL₃/L_im. Neglecting noise, and assuming that the 1.22λ₀/D angular resolution is sufficient to make

{\tilde{𝒫}}_{2} (ρ_{2}, Δ ω) \approx 𝒫_{2} (ρ_{2}, Δ ω),

for |ρ₂| ≤ D′/2, Eq. (34) leads to

\begin{array}{l} \tilde{T} ({\tilde{ρ}}_{1}) | 𝒫_{1} ({\tilde{ρ}}_{1}, Δ ω) | & = | \int d^{2} ρ_{1} 𝒫_{1} (ρ_{1}, Δ ω) T (ρ_{1}) e^{i Δ ω | ρ_{1} |^{2} / 2 c L_{2}} \\ \times \frac{π}{4} {(\frac{D^{'}}{Δ λ L_{2}})}^{2} \frac{J_{1} (π D^{'} | {\tilde{ρ}}_{1} - ρ_{1} | / Δ λ L_{2})}{π D^{'} | {\tilde{ρ}}_{1} - ρ_{1} | / 2 Δ λ L_{2}} | . \end{array}

Thus, over the region in the z = L₁ plane wherein |𝒫₁(ρ₁, Δω)| has an appreciable value, the 𝒫-field imager using cosinusoidal E-field modulation at frequency Δω/2 achieves a spatial resolution of 1.22ΔλL₂/D′, where: Δλ = 2πc/Δω; L₂ is the distance from the transparency-containing plane to the plane visible to the sensor; and D′ = d_imL₃/L_im, with L₃ being the distance from the plane visible to the sensor to the sensor’s entrance pupil, L_im being the distance from that entrance pupil to the image plane where irradiance measurements are made, and d_im being the diameter of the image-plane region over which those measurements are made.

3. 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 𝒫 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. [12] and Thrampoulidis et al. [13].

3.1. Setup for paraxial propagation through multiple diffusers with occlusion

Figure 3 shows a generalized setup for transmissive 𝒫-field imaging. Here, two occluders, having field-transmission functions P(ρ_d) and P′(ρ′_d), have been introduced in the z = L₁ − L_d and z = L₁ + L′_d planes, and the z = L₁ plane contains a field-transmission mask F(ρ₁) 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(ρ₁) 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 [12].

Fig. 3 Unfolded geometry for three-bounce, occlusion-aided NLoS active imaging. Scalar, paraxial diffraction theory is assumed, with E₀(ρ₀, t) being the baseband complex-field envelope illuminating the z = 0 plane and E′₂(ρ₂, t) being the baseband complex-field envelope emerging from the z = L₁ + L₂ plane. These fields are written as functions of their transverse spatial coordinates, {ρ_k = (x_k, y_k) : k = 0, 2}, in their respective planes and time, t. The blue rectangles represent thin transmissive diffusers, and the black line at z = L₁ represents a thin specular-plus-diffuser transmission mask with field-transmission function F(ρ₁), whose associated intensity-transmission pattern is to be imaged using the light that emerges from the z = L₁ + L₂ plane. That imaging process is aided by the presence of occluders in the z = L₁ − L_d and z = L₁ + L′_d planes, whose field-transmission functions are P(ρ_d) and P′(ρ′_d), respectively.

Download Full Size | PDF

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 𝒫 field—given the short-time average irradiance, or its associated 𝒫 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 𝒫-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 𝒫-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₁ − L_d, z = L₁, or z = L₁ + L′_d in Fig. 3, we are interested in planes that do not contain pure diffusers, the 𝒫 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 𝒫 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 𝒫 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 𝒫-field input-output relation built up from propagating the 𝒫 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 [17]:

I_{z} (ρ_{+}, s, t) \equiv \int \frac{d^{2} ρ_{-}}{λ_{0}^{2}} 〈 E_{z} (ρ_{+} + ρ_{-} / 2, t) E_{z}^{*} (ρ_{+} - ρ_{-} / 2, t) 〉 e^{- i 2 π s \cdot ρ_{-} / λ_{0}} .

In computer vision, this quantity is known as the 5D light field [18–20]. By replacing 2πs/λ₀ with 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.,

W (ρ_{+}, k) \equiv \int d^{2} ρ_{-} E_{z} (ρ_{+} + ρ_{-} / 2) E_{z}^{*} (ρ_{+} - ρ_{-} / 2) e^{- i k \cdot ρ_{-}},

which has long been recognized as a useful tool in optics, see, e.g., [21–23]. The short-time average irradiance is obtained from I_z(ρ₊, s, t) by integrating out its directionality information,

〈 I_{z} (ρ_{+}, t) 〉 = \int d^{2} s I_{z} (ρ_{+}, s, t),

and the 𝒫 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:

\begin{array}{l} \begin{array}{l} I_{z} (ρ_{+}, s, t) \\ = \int \frac{d ω}{2 π} \int \frac{d ω^{'}}{2 π} \int \frac{d^{2} ρ_{-}}{λ_{0}^{2}} 〈 ℰ_{z} (ρ_{+} + ρ_{-} / 2, ω) ℰ_{z}^{*} (ρ_{+} - ρ_{-} / 2, ω^{'}) 〉 e^{- i 2 π s \cdot ρ_{-} / λ_{0}} e^{- i (ω - ω^{'}) t} \end{array} \end{array}

= \int \frac{d ω_{-}}{2 π} [\int \frac{d ω_{+}}{2 π} (\int \frac{d^{2} ρ_{-}}{λ_{0}^{2}} 〈 ℰ_{z} (ρ_{+} + ρ_{-} / 2, ω) ℰ_{z}^{*} (ρ_{+} - ρ_{-} / 2, ω^{'}) 〉 e^{- 2 π s \cdot ρ_{-} / λ_{0}})] e^{- i ω_{-} t},

where ω₊ ≡ (ω + ω′)/2 and ω₋ ≡ ω − ω′ as we employed in Sec. 2. 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. 2 analysis, this quantity is the directionality-augmented analog of the 𝒫 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 𝒫 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):

W_{ℰ_{z}} (ρ_{+}, k, ω_{+}, ω_{-}) \equiv \int d^{2} ρ_{-} 〈 ℰ_{z} (ρ_{+} + ρ_{-} / 2, ω_{+} + ω_{-} / 2) ℰ_{z}^{*} (ρ_{+} - ρ_{-} / 2, ω_{+} - ω_{-} / 2) 〉 e^{- i k \cdot ρ_{-}},

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

I_{z} (ρ_{+}, s, t) = \frac{1}{λ_{0}^{2}} \int \frac{d ω_{-}}{2 π} \int \frac{d ω_{+}}{2 π} W_{ℰ_{z}} (ρ_{+}, 2 π s / λ_{0}, ω_{+}, ω_{-}) e^{- i ω_{-} t} .

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

Γ_{z} (ρ_{1}, ρ_{2}, t_{1}, t_{2}) \equiv 〈 E_{z} (ρ_{1}, t_{1}) E_{z}^{*} (ρ_{2}, t_{2}) 〉,

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

I_{z} (ρ_{+}, s, t) = \int \frac{d^{2} ρ_{-}}{λ_{0}^{2}} Γ_{z} (ρ_{+} + ρ_{-} / 2, ρ_{+} - ρ_{-} / 2, t, t) e^{- i 2 π s \cdot ρ_{-} / λ_{0}},

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

\begin{array}{l} W_{ℰ_{z}} (ρ_{+}, k, ω_{+}, ω_{-}) & = \int d^{2} ρ_{-} \int d t_{1} \int d t_{2} Γ_{z} (ρ_{+} + ρ_{-} / 2, ρ_{+} - ρ_{-} / 2, t_{+} + t_{-} / 2, t_{+} - t_{-} / 2) \\ \times e^{i (ω_{+} t_{-} + ω_{-} t_{+} - k \cdot ρ_{-})}, \end{array}

where t₊ ≡ (t₁ + t₂)/2, t₋ ≡ t₁ − t₂, and

\begin{array}{l} Γ_{z} (ρ_{+} + ρ_{-} / 2, ρ_{+} - ρ_{-} / 2, t_{+} + t_{-} / 2, t_{+} - t_{-} / 2) & = \int \frac{d^{2} k}{{(2 π)}^{2}} \int \frac{d ω_{+}}{2 π} \int \frac{d ω_{-}}{2 π} W_{ℰ_{z}} (ρ_{+}, k, ω_{+}, ω_{-}) \\ \times e^{- i (ω_{+} t_{-} + ω_{-} t_{+} - k \cdot ρ_{-})} . \end{array}

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

E_{z}^{'} (ρ^{'}, t) = \int d τ \int d^{2} ρ E_{z} (ρ, τ) h (ρ^{'}, ρ; t, τ),

the input’s space-time autocorrelation function suffices to determine the output’s space-time autocorrelation function, and hence the output-plane 𝒫 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 𝒫 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,

I_{z} (ρ_{+}, s, ω_{+}, t) \equiv \frac{1}{λ_{0}^{2}} \int \frac{d ω_{-}}{2 π} W (ρ_{+}, 2 π s / λ_{0}, ω_{+}, ω_{-}) e^{- i ω_{-} t},

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

The z-plane 𝒫 field can be found from that plane’s TFSWD as follows:

𝒫_{z} (ρ_{+}, ω_{-}) = \int \frac{d ω_{+}}{2 π} \int \frac{d^{2} k}{{(2 π)}^{2}} W_{ℰ_{z}} (ρ_{+}, k, ω_{+}, ω_{-}) .

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_{ℰ_z} (ρ₊, k, ω₊, ω₋) through a z-plane field-transmission mask, whether that be a diffuser, deterministic occluder, or specular-plus-diffuser mask; and (2) propagate W_{ℰ_z} (ρ₊, k, ω₊, ω₋) 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
$W_{ℰ_{0}^{'}} (ρ_{+}, k, ω_{+}, ω_{-}) = λ_{0}^{2} \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} W_{ℰ_{0}} (ρ_{+}, k^{'}, ω_{+}, ω_{-}) .$ .
Propagation through a deterministic occluder:
With W_P(ρ₊, k) ≡ ∫d²ρ₋ P(ρ₊ + ρ₋/2)P^*(ρ₊ − ρ₋/2)e^−ik·ρ₋, we have
$W_{ℰ_{L_{1} - L_{d}}^{'}} (ρ_{+}, k, ω_{+}, ω_{-}) = \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} W_{ℰ_{L_{1} - L_{d}}} (ρ_{+}, k^{'}, ω_{+}, ω_{-}) W_{P} (ρ_{+}, k - k^{'}) .$
Propagation through a specular-plus-diffuser mask:
With F(ρ₁) having nonzero mean 〈F(ρ₁)〉 ≠ 0, and covariance, $〈 Δ F (ρ_{+} + ρ_{-} / 2) Δ F^{*} (ρ_{+} - ρ_{-} / 2) 〉 \approx λ_{0}^{2} ℱ (ρ_{+}) δ (ρ_{-})$ where 0 ≤ ℱ(ρ₊) ≤ 1, we get
$\begin{array}{l} W_{ℰ_{L_{1}}^{'}} (ρ_{+}, k, ω_{+}, ω_{-}) = & \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} W_{ℰ_{L_{1}}} (ρ_{+}, k^{'}, ω_{+}, ω_{-}) W_{〈 F 〉} (ρ_{+}, k - k^{'}) \\ + λ_{0}^{2} ℱ (ρ_{+}) \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} W_{ℰ_{L_{1}}} (ρ_{+}, k^{'}, ω_{+}, ω_{-}) . \end{array}$
Fresnel diffraction:
For Fresnel diffraction from the z = 0₊ plane to the z = L₁ − L_d plane, we get
$W_{ℰ_{L_{1} - L_{d}}} (ρ_{+}, k, ω_{+}, ω_{-}) = W_{ℰ_{0}^{'}} (ρ_{+} - c (L_{1} - L_{d}) k / ω_{0}, k, ω_{+}, ω_{-}) e^{i [ω_{-} (L_{1} - L_{d}) / c] (1 + c^{2} {| k |}^{2} / 2 ω_{0}^{2})} .$

3.2. Occlusion-aided imaging

In Sec. 2 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₁ 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₀(t), and that the distances in Fig. 3 satisfy L₁ = L₂ = L, and L_d = L/2

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

W_{ℰ_{0}} (ρ_{+}, k, ω_{+}, ω_{-}) = W_{in} (ω_{+}, ω_{-}) {(2 π / λ_{0})}^{2} δ (k),

where

W_{in} (ω_{+}, ω_{-}) = λ_{0}^{2} \int d t \sqrt{I_{0} (t)} e^{i (ω_{+} + ω_{-} / 2) t} \int d u \sqrt{I_{0} (u)} e^{- i (ω_{+} - ω_{-} / 2) u} .

After the diffuser in the z = 0 plane we get

W_{ℰ_{0}^{'}} (ρ_{+}, k, ω_{+}, ω_{-}) = W_{in} (ω_{+}, ω_{-}),

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

W_{ℰ_{L}} (ρ_{+}, k, ω_{+}, ω_{-}) = W_{in} (ω_{+}, ω_{-}) e^{i (ω_{-} L / c) (1 + c^{2} {| k |}^{2} / 2 ω_{0}^{2})} .

At z = L₁ this Wigner distribution encounters a diffuse target mask, i.e., one whose field-transmission function F(ρ₁) has zero mean and covariance $〈 Δ F (ρ_{1}) Δ F^{*} (ρ_{2}) 〉 \approx λ_{0}^{2} ℱ [(ρ_{1} + ρ_{2}) / 2] δ (ρ_{1} - ρ_{2})$ , which results in

W_{ℰ_{L}^{'}} (ρ_{+}, k, ω_{+}, ω_{-}) = ℱ (ρ_{+}) W_{in} (ω_{+}, ω_{-}) e^{i ω_{-} L / c} 2 π i c / ω_{-} L .

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

W_{ℰ_{3 L / 2}} (ρ_{+}, k, ω_{+}, ω_{-}) = ℱ (ρ_{+} - c L k / 2 ω_{0}) W_{in} (ω_{+}, ω_{-}) e^{i ω_{-} 3 L / 2 c} e^{i ω_{-} c L {| k |}^{2} / 4 ω_{0}^{2}} 2 π i c / ω_{-} L,

and passage through the occluder in that plane leads to

\begin{array}{l} W_{ℰ_{3 L / 2}^{'}} (ρ_{+}, k, ω_{+}, ω_{-}) = W_{in} (ω_{+}, ω_{-}) \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} ℱ (ρ_{+} - c L k^{'} / 2 ω_{0}) e^{i ω_{-} 3 L / 2 c} e^{i ω_{-} c L {| k^{'} |}^{2} / 4 ω_{0}^{2}} \\ \times W_{P} (ρ_{+}, k - k^{'}) 2 π i c / ω_{-} L . \end{array}

Fresnel propagation over another L/2 distance then gives

\begin{array}{l} W_{ℰ_{2 L}} (ρ_{+} k, ω_{+}, ω_{-}) = W_{in} (ω_{+}, ω_{-}) \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} ℱ (ρ_{+} - c L (k^{'} + k) / 2 ω_{0}) e^{i ω_{-} 2 L / c} \\ \times e^{i ω_{-} c L ({| k |}^{2} + {| k^{'} |}^{2}) / 4 ω_{0}^{2}} W_{P} (ρ_{+} - c L k / 2 ω_{0}, k - k^{'}) 2 π i c / ω_{-} L, \end{array}

from which we get

\begin{array}{l} 𝒫_{2 L} (ρ_{+}, ω_{-}) = \int \frac{d ω_{+}}{2 π} W_{in} (ω_{+}, ω_{-}) \int \frac{d^{2} k}{{(2 π)}^{2}} \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} ℱ (ρ_{+} - c L (k^{'} + k) / 2 ω_{0}) e^{i ω_{-} 2 L / c} \\ \times e^{i ω_{-} c L ({| k |}^{2} + {| k^{'} |}^{2}) / 4 ω_{0}^{2}} W_{P} (ρ_{+} - c L k / 2 ω_{0}, k - k^{'}) 2 π i c / ω_{-} L . \end{array}

Now, using

𝒫_{0} (ρ_{+}, ω_{-}) = \int \frac{d ω_{+}}{2 π} \int \frac{d^{2} k}{{(2 π)}^{2}} W_{ℰ_{0}} (ρ_{+}, k, ω_{+}, ω_{-}) = \int d t I_{0} (t) e^{i ω_{-} t},

and changing variables to k₋ = k − k′ and k₊ = (k + k′)/2 we have

\begin{array}{l} 𝒫_{2 L} (ρ_{+}, ω_{-}) & = λ_{0}^{2} 𝒫_{0} (ω_{-}) e^{i ω_{-} 2 L / c} \int \frac{d^{2} k_{+}}{{(2 π)}^{2}} \int \frac{d^{2} k_{-}}{{(2 π)}^{2}} ℱ (ρ_{+} - c L k_{+} / ω_{0}) \\ \times e^{i ω_{-} c L (2 | k_{+} |^{2} + | k_{-} |^{2} / 2) / 4 ω_{0}^{2}} W_{P} (ρ_{+} - c L (k_{+} / 2 + k_{-} / 4) / ω_{0}, k_{-}) 2 π i c / ω_{-} L, \end{array}

where we have suppressed the ρ₊ argument of 𝒫₀(ρ₊, ω₋) because that field has no such dependence for the plane-wave source we have assumed. We define a new function

G (ρ, ω_{-}) = \int \frac{d^{2} k_{-}}{{(2 π)}^{2}} e^{i ω_{-} c L | k_{-} |^{2} / 8 ω_{0}^{2}} W_{P} (- ρ / 2 - c L k_{-} / 4 ω_{0}, k_{-}) 2 π i c / ω_{-} L .

With this definition we have

\begin{array}{l} 𝒫_{2 L} (ρ_{+}, ω_{-}) & = λ_{0}^{2} 𝒫_{0} (ω_{-}) e^{i ω_{-} 2 L / c} \int \frac{d^{2} k_{+}}{{(2 π)}^{2}} ℱ (ρ_{+} - c L k_{+} / ω_{0}) \\ \times G (- 2 ρ_{+} + c L k_{+} / ω_{0}, ω_{-}) e^{i ω_{-} c L | k_{+} |^{2} / 2 ω_{0}^{2}} . \end{array}

Changing variables again, ρ̃ = ρ₊ − cLk₊/ω₀, we get our final result

𝒫_{2 L} (ρ_{+}, ω_{-}) = 𝒫_{0} (ω_{-}) e^{i ω_{-} 2 L / c} \int d^{2} \tilde{ρ} ℱ (\tilde{ρ}) G (- ρ_{+} - \tilde{ρ}, ω_{-}) \frac{e^{i ω_{-} | ρ_{+} \tilde{ρ} |^{2} / 2 c L}}{L^{2}} .

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 Re[E₀(ρ₀)e^−iω₀t]. In this unmodulated case we can use the usual spatial Wigner distribution, i.e.,

W_{E_{0}} (ρ_{+}, k) \equiv \int d^{2} ρ_{-} E_{0} (ρ_{+} + ρ_{-} / 2) E_{0}^{*} (ρ_{+} - ρ_{-} / 2) e^{- i k \cdot ρ_{-}},

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 ω₋ = 0 in the Fresnel-diffraction primitive. Paralleling the development that led to Eq. (68) assuming that

E_{0} (ρ_{0}) = \sqrt{I_{0}}

is a constant, we get

I_{2 L} (ρ_{+}) \equiv 〈 | E_{2 L} (ρ_{+}) |^{2} 〉 = I_{0} \int d^{2} \tilde{ρ} ℱ (\tilde{ρ}) G (- ρ_{+} - \tilde{ρ}),

where

G (ρ) \equiv \frac{π}{L^{2}} \int \frac{d^{2} k_{-}}{{(2 π)}^{2}} W_{P} (- ρ / 2 - c L k_{-} / 4 ω_{0}, k_{-}),

and we have used the evanescence cutoff, |k| ≤ 2π/λ₀, to justify replacing ∫ d²k I₀/(2π)² with

π I_{0} / λ_{0}^{2}

.

Equations (70) and (71) show that this unmodulated case offers no spatial information about ℱ(ρ) in the absence of an occluder, i.e., we get G(ρ) = π/L² when P(ρ) = 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

P_{ph} (ρ) = e^{- | ρ |^{2} / 2 ρ_{0}^{2}},

and the Gaussian pinspeck,

P_{ps} (ρ) = 1 - e^{- | ρ |^{2} / 2 ρ_{0}^{2}},

where ρ₀ 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_ph(ρ) is revealing. The Gaussian-pinspeck camera, on the other hand, is more relevant to the experiments of Xu et al. [12], but its psf G_ph(ρ) 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

G_{ph} (ρ) = \frac{π Ω^{2}}{L^{2} (1 + Ω^{2})} exp [- \frac{Ω^{2}}{1 + Ω^{2}} \frac{| ρ |^{2}}{4 ρ_{0}^{2}}],

where k₀ ≡ ω₀/c = 2π/λ₀ is the wave number at the optical frequency and

Ω \equiv 4 k_{0} ρ_{0}^{2} / L

is the Fresnel number for the pinhole’s propagation geometry. The spatial resolution of G_ph(ρ) improves with decreasing ρ₀ when Ω > 1, and degrades with decreasing ρ₀ when Ω < 1. Thus the Gaussian pinhole’s resolution-optimized psf,

G_{ph}^{opt} (ρ) = \frac{π exp (- π | ρ |^{2} / λ_{0} L)}{2 L^{2}},

is obtained when

ρ_{0} = \sqrt{L / 4 k_{0}} = \sqrt{λ_{0} L / 8 π}

. The optimized psf’s spatial resolution—taken to be its e^−π-attenuation radius—is then

\sqrt{λ_{0} L}

, which is far superior to the 1.22ΔλL/D′ for the unoccluded, modulated case governed by Eq. (36). For example, with λ = 1 μm and L = 1 m the optimum spatial resolution of occlusion-aided unmodulated imaging is 1 mm, while that of unoccluded modulated imaging, with Δλ = 3 cm (Δω/2π = 10 GHz) and D′ = 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 Ω, with only its overall amplitude G_ph(0) and its spatial resolution

ρ_{res} (Ω) \equiv \sqrt{4 π (1 + Ω^{2})} ρ_{0} / Ω

changing, i.e., we have that

G_{ph} (ρ) / G_{ph} (0) = exp [- π | ρ^{2} | ρ_{res}^{2} (Ω)],

for the Gaussian pinhole.

For the Gaussian pinspeck, we get

G_{ps} (ρ) = \frac{π}{L^{2}} {| 1 - \frac{Ω}{\sqrt{1 + Ω^{2}}} exp [- \frac{Ω}{1 + Ω^{2}} \frac{| ρ |^{2}}{8 ρ_{0}^{2}} (Ω - i) - i {tan}^{- 1} (1 / Ω)] |}^{2} .

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 Ω = 1, in which case we get

G_{ps}^{opt} (ρ) = \frac{π}{L^{2}} {| 1 - \frac{exp (- π | ρ |^{2} (1 - i) / 2 λ_{0} L - i π / 4)}{\sqrt{2}} |}^{2} .

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_ps(ρ)/G_ps(∞) versus ρ/ρ_res(Ω) for ρ = (x, 0) and Ω = 0.1, 1, and 10, where ρ_res(Ω) is the Gaussian pinhole’s spatial resolution.

Fig. 4 Plots of G_ps(ρ)/G_ps(∞) for the Gaussian pinspeck versus ρ/ρ_res(Ω) for ρ = (x, 0) and Ω = 0.1, 1, and 10.

Download Full Size | PDF

Note that in the near-field region, wherein Ω ≫ 1, Eq. (77) reduces to the geometric optics result,

G_{ps} (ρ) / G_{ps} (\infty) = {[1 - exp (- | ρ |^{2} / 8 ρ_{0}^{2})]}^{2},

which is analogous to the geometric optics treatment used by Xu et al. [12] and Thrampoulidis et al. [13] for the hard-aperture, circular occluder

P (ρ) = circ (2 ρ / d) \equiv {\begin{array}{l} 1, & for | ρ | \leq d / 2 \\ 0, & otherwise . \end{array}

4. 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 𝒫 field and provided a straightforward derivation of its behavior, analogous to that reported by Reza et al. [9]. 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 [10] 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 𝒫-field imaging of specular objects. Next, because Liu et al. [10] 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 𝒫 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 𝒫-field optics introduced and demonstrated in Reza et al. [11].

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_{ℰ_z} (ρ₊, k, ω₊, ω₋) and we want to find $W_{ℰ_{z}^{'}} (ρ_{+}, k, ω_{+}, ω_{-})$ , where

ℰ_{z}^{'} (ρ, ω) = ℰ_{z} (ρ, ω) e^{i (ω_{0} + ω) h_{z} (ρ) / c} \approx ℰ_{z} (ρ, ω) e^{i ω_{0} h_{z} (ρ) / c},

with

〈 e^{i ω_{0} [h_{z} (ρ) - h_{k} (ρ')] / c} 〉 \approx λ_{0}^{2} δ (ρ - ρ') .

In this case we immediately get

W_{ℰ_{z}^{'}} (ρ_{+}, k, ω_{+}, ω_{-}) = \int d^{2} ρ_{-} 〈 ℰ_{z}^{'} (ρ_{+} + ρ_{-} / 2, ω_{+} + ω_{-} / 2) ℰ_{z}^{' *} (ρ_{+} - ρ_{-} / 2, ω_{+} - ω_{-} / 2) 〉 e^{- i k \cdot ρ_{-}}

\begin{array}{l} = \int d^{2} ρ_{-} 〈 ℰ_{z} (ρ_{+} + ρ_{-} / 2, ω_{+} + ω_{-} / 2) ℰ_{z}^{*} (ρ_{+} - ρ_{-} / 2, ω_{+} - ω_{-} / 2) 〉 \\ \times 〈 e^{i ω_{0} [h_{z} (ρ_{+} + ρ_{-} / 2) - h_{k} (ρ_{+} - ρ_{-} / 2)] / c} 〉 e^{- i k \cdot ρ_{-}} \end{array}

= λ_{0}^{2} 〈 ℰ_{z} (ρ_{+}, ω_{+} + ω_{-} / 2) ℰ_{z}^{*} (ρ_{+}, ω_{+} - ω_{-} / 2) 〉

= λ_{0}^{2} \int \frac{d^{2} k'}{{(2 π)}^{2}} W_{ℰ_{z}} (ρ_{+}, k', ω_{+}, ω_{-}) .

Physically, the 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 ℰ_z(ρ, ω), because

W_{ℰ_{z}^{'}} (ρ_{+}, k, ω_{+}, ω_{-})

is independent of k.

Propagation through a deterministic occluder:

Now consider propagation through a deterministic transmission mask. Here we want to find $W_{ℰ_{z}^{'}} (ρ_{+}, k, ω_{+}, ω_{-})$ given W_{ℰ_z} (ρ₊, k, ω₊, ω₋) and a deterministic P(ρ), where

ℰ_{z}^{'} (ρ, ω) = ℰ_{z} (ρ, ω) P (ρ) .

For this case we have that

W_{ℰ_{z}^{'}} (ρ_{+}, k, ω_{+}, ω_{-}) = \int d^{2} ρ_{-} 〈 ℰ_{z}^{'} (ρ_{+} + ρ_{-} / 2, ω_{+} + ω_{-} / 2) ℰ_{z}^{' *} (ρ_{+} - ρ_{-} / 2, ω_{+} - ω_{-} / 2) 〉 e^{- i k \cdot ρ_{-}}

\begin{array}{l} = \int d^{2} ρ_{-} 〈 ℰ_{z} (ρ_{+} + ρ_{-} / 2, ω_{+} + ω_{-} / 2) ℰ_{z}^{*} (ρ_{+} - ρ_{-} / 2, ω_{+} - ω_{-} / 2) 〉 \\ \times P (ρ_{+} + ρ_{-} / 2) P^{*} (ρ_{+} - ρ_{-} / 2) e^{- i k \cdot ρ_{-}} \end{array}

= \int \frac{d^{2} k'}{{(2 π)}^{2}} W_{ℰ_{z}} (ρ_{+}, k', ω_{+}, ω_{-}) \int d^{2} ρ_{-} P (ρ_{+} + ρ_{-} / 2) P^{*} (ρ_{+} - ρ_{-} / 2) e^{- i (k - k') \cdot ρ_{-}}

= \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} W_{ℰ_{z}} (ρ_{+}, k^{'}, ω_{+}, ω_{-}) W_{P} (ρ_{+}, k - k^{'}),

where

W_{P} (ρ_{+}, k) \equiv \int d^{2} ρ_{-} P (ρ_{+} + ρ_{-} / 2) P^{*} (ρ_{+} - ρ_{-} / 2) e^{- i k \cdot ρ_{-}}

is the spatial Wigner distribution of P(ρ). In words, Eq. (91) shows that multiplying ℰ_z(ρ, ω) by a deterministic field-transmission mask implies that

W_{ℰ_{z}^{'}} (ρ_{+}, k, ω_{+}, ω_{-})

is obtained from a k-space convolution of W_{ℰ_z} (ρ₊, k, ω₊, ω₋) with the field-transmission mask’s spatial Wigner distribution. Moreover, Eq. (92), together with Eq. (50), immediately leads to

𝒫_{z} (ρ_{+}, ω_{-}) = \int \frac{d ω_{+}}{2 π} \int \frac{d^{2} k}{{(2 π)}^{2}} W_{ℰ_{z}} (ρ_{+}, k, ω_{+}, ω_{-})

= \int \frac{d ω_{+}}{2 π} \int \frac{d^{2} k}{{(2 π)}^{2}} \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} W_{ℰ_{z}} (ρ_{+}, k^{'}, ω_{+}, ω_{-}) W_{P} (ρ_{+}, k - k^{'})

= \int \frac{d ω_{+}}{2 π} \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} W_{ℰ_{z}} (ρ_{+}, k^{'}, ω_{+}, ω_{-}) {| P (ρ_{+}) |}^{2} = 𝒫_{z} (ρ_{+}, ω_{-}) {| P (ρ_{+}) |}^{2},

as could have been directly obtained from Eq, (87) and the 𝒫-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(ρ₁) having nonzero mean 〈F(ρ₁)〉 ≠ 0, and covariance, $〈 Δ F (ρ_{+} + ρ_{-} / 2) Δ F^{*} (ρ_{+} - ρ_{-} / 2) 〉 \approx λ_{0}^{2} ℱ (ρ_{+}) δ (ρ_{-})$ where 0 ≤ ℱ(ρ₊) ≤ 1 and ΔF(ρ) ≡ F(ρ) − 〈F(ρ)〉. The propagation analysis follows from combining the two previous analyses:

\begin{array}{l} W_{ℰ_{z}^{'}} (ρ_{+}, k, ω_{+}, ω_{-}) \\ = \int d^{2} ρ_{-} 〈 ℰ_{z}^{'} (ρ_{+} + ρ_{-} / 2, ω_{+} + ω_{-} / 2) ℰ_{z}^{' *} (ρ_{+} - ρ_{-} / 2, ω_{+} - ω_{-} / 2) 〉 e^{- i k \cdot ρ_{-}} \end{array}

\begin{array}{l} = \int d^{2} ρ_{-} 〈 ℰ_{z} (ρ_{+} + ρ_{-} / 2, ω_{+} + ω_{-} / 2) ℰ_{z}^{*} (ρ_{+} - ρ_{-} / 2, ω_{+}, - ω_{-} / 2) 〉 \\ \times 〈 F (ρ_{+} + ρ_{-} / 2) F^{*} (ρ_{+} - ρ_{-} / 2) 〉 e^{- i k \cdot ρ_{-}} . \end{array}

From expanding F(ρ) into a sum of its (deterministic) mean and zero-mean random portions, it follows that

\begin{array}{l} W_{ℰ_{z}^{'}} (ρ_{+}, k, ω_{+}, ω_{-}) = \int d^{2} ρ_{-} 〈 ℰ_{z} (ρ_{+} + ρ_{-} / 2, ω_{+} + ω_{-} / 2) ℰ_{z}^{*} (ρ_{+} - ρ_{-} / 2, ω_{+} - ω_{-} / 2) 〉 \\ \times (〈 F (ρ_{+} + ρ_{-} / 2) 〉 〈 F^{*} (ρ_{+} - ρ_{-} / 2) 〉 + 〈 Δ F (ρ_{+} + ρ_{-} / 2) Δ F^{*} (ρ_{+} - ρ_{-} / 2) 〉) e^{- i k \cdot ρ_{-}} \end{array}

= \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} W_{ℰ_{L_{1}}} (ρ_{+}, k^{'}, ω_{+}, ω_{-}) W_{〈 F 〉} (ρ_{+}, k - k^{'}) + λ_{0}^{2} ℱ (ρ_{+}) \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} W_{ℰ_{L_{1}}} (ρ_{+}, k^{'}, ω_{+}, ω_{-}) .

Fresnel diffraction:

Our final task is to find W_{ℰ_L} (ρ₊, k, ω₊, ω₋) when

ℰ_{L} (ρ_{L}, ω) = \int d^{2} ρ_{0} ℰ_{0} (ρ_{0}, ω) \frac{(ω_{0} + ω) e^{i (ω_{0} + ω) (L / c + {| ρ_{L} - ρ_{0} |}^{2} / 2 c L)}}{i 2 π c L},

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

\begin{array}{l} W_{ℰ_{L}} (ρ_{+}, k, ω_{+}, ω_{-}) & = \int d^{2} ρ_{-} \int d^{2} ρ_{0} \int d^{2} ρ_{0}^{'} 〈 ℰ_{0} (ρ_{0}, ω_{+} + ω_{-} / 2) ℰ_{0}^{*} (ρ_{0}^{'}, ω_{+} - ω_{-} / 2) 〉 \\ \times e^{i ω_{-} L / c} e^{- i k \cdot ρ_{-}} \frac{(ω_{0} + ω_{+} + ω_{-} / 2) e^{i (ω_{0} + ω_{+} + ω_{-} / 2) {| ρ_{+} + ρ_{-} / 2 - ρ_{0} |}^{2} / 2 c L}}{i 2 π c L} \\ \times \frac{(ω_{0} + ω_{+} - ω_{-} / 2) e^{- i (ω_{0} + ω_{+} - ω_{-} / 2) {| ρ_{+} - ρ_{-} / 2 - ρ_{0}^{'} |}^{2} / 2 c L}}{- i 2 π c L} . \end{array}

Exploiting Δω ≪ ω₀, and making the coordinate transformation from ρ₀ and ρ′₀ to ρ_0₊ ≡ (ρ₀ + ρ′₀)/2 and ρ_0₋ ≡ ρ₀ − ρ′₀, we can reduce Eq. (101) to

\begin{array}{l} W_{ℰ_{L}} (ρ_{+}, k, ω_{+}, ω_{-}) \\ = \int d^{2} ρ_{-} \int d^{2} ρ_{0_{+}} \int d^{2} ρ_{0_{-}} 〈 ℰ (ρ_{0_{+}} + ρ_{0_{-}} / 2, ω_{+} + ω_{-} / 2) ℰ_{0}^{*} (ρ_{0_{+}} - ρ_{0_{-}} / 2, ω_{+} - ω_{-} / 2) 〉 \\ \times \frac{e^{i ω_{-} L / c}}{{(λ_{0} L)}^{2}} e^{i (ω_{0} + ω_{+}) + (ρ_{+} - ρ_{0_{+}}) \cdot (ρ_{-} - ρ_{0 -}) / c L} e^{i ω_{-} ({| ρ_{+} - ρ_{0_{+}} |}^{2} + {| ρ_{-} - ρ_{0 -} |}^{2} / 4) / 2 c L} e^{- i k \cdot ρ_{-}} . \end{array}

Rearranging terms allows us to put the ρ₋ integral inside the ρ_0₊ and ρ_0₋ integrals, i.e.,

\begin{array}{l} W_{ℰ_{L}} (ρ_{+}, k, ω_{+}, ω_{-}) \\ = \int d^{2} ρ_{0_{+}} \int d^{2} ρ_{0_{-}} 〈 ℰ_{0} (ρ_{0_{+}} + ρ_{0_{-}} / 2, ω_{+} + ω_{-} / 2) ℰ_{0}^{*} (ρ_{0_{+}} - ρ_{0_{-}} / 2, ω_{+} - ω_{-} / 2) 〉 \frac{e^{i ω_{-} L / c}}{{(λ_{0} L)}^{2}} \\ \times e^{- i (ω_{0} + ω_{+}) (ρ_{+} - ρ_{0_{+}}) \cdot ρ_{0 -} / c L} e^{i ω_{-} ({| ρ_{+} - ρ_{0_{+}} |}^{2} / 2 c L + {| ρ_{0_{-}} |}^{2} / 8 c L)} \\ \times \int d^{2} ρ_{-} e^{i ω_{-} {| ρ_{-} |}^{2} / 8 c L} e^{- i [k - (ω_{0} + ω_{+}) (ρ_{+} - ρ_{0_{+}}) / c L + ω_{-} ρ_{0_{-}} / 4 c L] \cdot ρ_{-}} . \end{array}

Performing the ρ₋ integral then yields

\begin{array}{l} W_{ℰ_{L}} (ρ_{+}, k, ω_{+}, ω_{-}) \\ = \int d^{2} ρ_{0_{+}} \int d^{2} ρ_{0_{-}} 〈 ℰ_{0} (ρ_{0_{+}} + ρ_{0_{-}} / 2, ω_{+} + ω_{-} / 2) ℰ_{0}^{*} (ρ_{0_{+}} - ρ_{0_{-}} / 2, ω_{+} - ω_{-} / 2) 〉 \frac{e^{i ω_{-} L / c}}{{(λ_{0} L)}^{2}} \\ \times e^{- i (ω_{0} + ω_{+}) (ρ_{+} - ρ_{0_{+}}) \cdot ρ_{0 -} / c L} e^{i ω_{-} {| ρ_{+} - ρ_{0_{+}} |}^{2} / 2 c L} e^{i ω_{-} {| ρ_{0_{-}} |}^{2} / 8 c L} (i 8 π c L / ω_{-}) \\ \times e^{- 2 i c L {| k - (ω_{0} + ω_{+}) (ρ_{+} - ρ_{0_{+}}) / c L + ω_{-} ρ_{0_{-}} / 4 c L |}^{2} / ω_{-}}, \end{array}

which, after some terms cancel, gives

\begin{array}{l} W_{ℰ_{L}} (ρ_{+}, k, ω_{+}, ω_{-}) \\ = \int d^{2} ρ_{0_{+}} \int d^{2} ρ_{0_{-}} 〈 ℰ_{0} (ρ_{0_{+}} + ρ_{0_{-}} / 2, ω_{+} + ω_{-} / 2) ℰ_{0}^{*} (ρ_{0_{+}} - ρ_{0_{-}} / 2, ω_{+} - ω_{-} / 2) 〉 \\ \times \frac{e^{i ω_{-} L / c}}{{(λ_{0} L)}^{2}} e^{i ω_{-} {| ρ_{+} - ρ_{0_{+}} |}^{2} / 2 c L} e^{- 2 i c L {| k - (ω_{0} + ω_{+}) (ρ_{+} - ρ_{0_{+}}) / c L |}^{2} / ω_{-}} e^{- i k \cdot ρ_{0_{-}}} (i 8 π c L / ω_{-}) \end{array}

\begin{array}{l} = \int d^{2} ρ_{0_{+}} W_{ℰ_{0}} (ρ_{0_{+}}, k, ω_{+}, ω_{-}) \frac{e^{i ω_{-} L / c}}{{(λ_{0} L)}^{2}} e^{i ω_{-} {| ρ_{+} - ρ_{0_{+}} |}^{2} / 2 c L} \\ \times e^{- 2 i c L {| k - (ω_{0} + ω_{+}) (ρ_{+} - ρ_{0_{+}}) / c L |}^{2} / ω_{-}} (i 8 π c L / ω_{-}) . \end{array}

The term

e^{- 2 i c L {| k - (ω_{0} + ω_{+}) (ρ_{+} - ρ_{0_{+}}) / c L |}^{2} / ω_{-}} i 8 π c L / ω_{-} {(λ_{0} L)}^{2}

in Eq. (106)’s integrand behaves like the impulse δ[ρ₀₊ − ρ₊ + kcL/(ω₀ + ω₊)]. 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 ρ₊ − kcL/(ω₀ + ω₊) with nominal width

\sqrt{ω_{-} c L} / 2 (ω_{0} + ω_{+})

, and ω₀ ≫ max |ω₊| implies that it integrates to one. (2) The other ρ_0₊-dependent terms in Eq. (106) are the oscillatory term, exp(iω₋|ρ₊ − ρ_0₊|²/2cL), which varies much more slowly than its predecessor, because ω₀ ≫ max |ω₋|, and the Wigner distribution, whose ρ_0₊ dependence can reasonably be assumed to be nearly constant over regions of diameter

\sqrt{ω_{-} c L} / 2 (ω_{0} + ω_{+})

. So, using the delta-function approximation in Eq. (106), we get

W_{ℰ_{L}} (ρ_{+}, k, ω_{+}, ω_{-}) = W_{ℰ_{0}^{'}} (ρ_{+} - c L k / (ω_{0} + ω_{+}), k, ω_{+}, ω_{-}) e^{i (ω_{-} L / c) (1 + c^{2} {| k |}^{2} / 2 {(ω_{0} + ω_{+})}^{2})} .

Finally, again making use ω₀ ≫ ω₊, we have

W_{ℰ_{L}} (ρ_{+}, k, ω_{+}, ω_{-}) = W_{ℰ_{0}^{'}} (ρ_{+} - c L k / ω_{0}, k, ω_{+}, ω_{-}) e^{i (ω_{-} L / c) (1 + c^{2} {| k |}^{2} / 2 ω_{0}^{2})} .

As a consistency check on Eq. (108), let us use it to calculate 𝒫_L(ρ₊, ω₋) when z = 0 illumination with TFSWD W_ℰ₀(ρ_0₊, k, ω₊, ω₋) passes through the diffuser specified in Eq. (81) before undergoing Fresnel diffraction over a distance L. We then have that

𝒫_{L} (ρ_{+}, ω_{-}) = \int \frac{d ω_{+}}{2 π} \int \frac{d^{2} k}{{(2 π)}^{2}} W_{ℰ_{0}^{'}} (ρ_{+} - c L k / ω_{0}, k, ω_{+}, ω_{-}) e^{i (ω_{-} L / c) (1 + c^{2} {| k |}^{2} / 2 ω_{0}^{2})} .

Using Eq. (86) now gives us

𝒫_{L} (ρ_{+}, ω_{-}) = λ_{0}^{2} \int \frac{d ω_{+}}{2 π} \int \frac{d^{2} k}{{(2 π)}^{2}} \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} W_{ℰ_{0}^{'}} (ρ_{+} - c L k / ω_{0}, k^{'}, ω_{+}, ω_{-}) e^{i (ω_{-} L / c) (1 + c^{2} {| k |}^{2} / 2 ω_{0}^{2})} .

Changing variables so that k = ω₀(ρ₊ − ρ₀)/cL leaves us with

𝒫_{L} (ρ_{+}, ω_{-}) = \int \frac{d ω_{+}}{2 π} \int d^{2} ρ_{0} \int \frac{d^{2} k^{'}}{{(2 π)}^{2}} W_{ℰ_{0}} (ρ_{0}, k^{'}, ω_{+}, ω_{-}) \frac{e^{i (ω_{-} L / c) (1 + {| ρ_{+} - ρ_{0} |}^{2} / 2 L^{2})}}{L^{2}} .

which reduces to the result from Sec. 2,

𝒫_{L} (ρ_{+}, ω_{-}) = \int d^{2} ρ_{0} 𝒫_{0} (ρ_{0}, ω_{-}) \frac{e^{i ω_{-} L / c} e^{i ω_{-} {| ρ_{+} - ρ_{0} |}^{2} / 2 c L}}{L^{2}},

by virtue of Eq. (50).

Funding

DARPA REVEAL program (HR0011-16-C-0030).

References

1. A. Kirmani, T. Hutchison, J. Davis, and R. Raskar, “Looking around the corner using ultrafast transient imaging,” Int. J. Comput. Vision 95, 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, 745 (2012). [CrossRef] [PubMed]

3. F. Heide, L. Xiao, W. Heidrich, and M. B. Hullin, “Diffuse mirrors: 3D reconstruction from diffuse indirect illumination using inexpensive time-of-flight sensors,” in Proc. IEEE Conf. Comput. Vis. Pattern Recog., pp. 3222–3229 (2014).

4. M. Buttafava, J. Zeman, A. Tosi, K. Eliceiri, and A. Velten, “Non-line-of-sight imaging using a time-gated single photon avalanche diode,” Opt. Express 23, 20997–21011 (2015). [CrossRef] [PubMed]

5. G. Gariepy, F. Tonolini, R. Henderson, J. Leach, and D. Faccio, “Detection and tracking of moving objects hidden from view,” Nat. Photonics 10, 23–27 (2015). [CrossRef]

6. A. Kadambi, H. Zhao, B. Shi, and R. Raskar, “Occluded imaging with time-of-flight sensors,” ACM Trans. Graph. 35, 1–12 (2016). [CrossRef]

7. J. Klein, M. Laurenzis, and M. Hullin, “Transient imaging for real-time tracking around a corner,” Proc. SPIE 9988, 998802 (2016). [CrossRef]

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

9. 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].

10. X. Liu, I. Guillén, M. La Manna, J. H. Nam, S. A. Reza, T. H. Le, D. Gutierrez, A. Jarabo, and A. Velten, “Virtual wave optics for non-line-of-sight imaging,” arXiv:1810.07535 [cs.CV].

11. S. A. Reza, M. La Manna, S. Bauer, and A. Velten, “Wave-like properties of phasor fields: experimental demonstrations,” arXiv:190401565 [physics.optics].

12. 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, 9945 (2018). [CrossRef] [PubMed]

13. 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. Imag. 4, 419 (2018). [CrossRef]

14. The short-time average z-plane irradiance is the instantaneous irradiance averaged over a time T_a satisfying ω₀T_a ≫ 1 and ΔωT_a ≪ 1.

15. In what follows, integrals without explicit limits are over the integration variable’s entire domain.

16. Because h₀(ρ) is a zero-mean Gaussian process, its samples at ρ₀ and ρ′₀ are zero-mean jointly Gaussian random variables whose joint characteristic function is as given in Eq. (12).

17. A. Ishimaru, Wave Propagation and Scattering in Random Media, Vol. 1: Single Scattering and Transport Theory (Academic, New York, 1978).

18. A. Gershun, “The light field,” J. Math. Phys. 18, 51–151 (1939). [CrossRef]

19. E. H. Adelson and J. R. Bergen, “The plenoptic function and the elements of early vision,” in M. S. Landy and J. A. Movshon, eds., Computational Models of Visual Processing, (MIT Press, 1991), pp. 3–20.

20. M. Levoy and P. Hanrahan, “Light field rendering,” in Proc. SIGGRAPH (ACM, New York, NY, USA, 1996), pp. 31–42.

21. A. Walther, “Radiometry and coherence,” J. Opt. Soc. Am. 58, 1256 (1968). [CrossRef]

22. M. J. Bastiaans, “Wigner distribution and its application to first-order optics,” J. Opt. Soc. Am. 69, 1710–1716 (1980). [CrossRef]

23. M. A. Alonso, “Wigner functions in optics: describing beams as ray bundles and pulses as particles,” Adv. Opt. Photon. 3, 272–365 (2011). [CrossRef]

24. A. Ishimaru, Wave Propagation and Scattering in Random Media, Vol. 2: Multiple Scattering, Turbulence, Rough Surfaces, and Remote Sensing (Academic, New York, 1978).

25. For notational convenience, we have assumed that the diffraction takes place between the z = 0 and z = L planes, but the result we obtain will apply for +z-going Fresnel diffraction over a distance L starting from an arbitrary z plane.

Paraxial theory of phasor-field imaging

Abstract

1. Introduction

2. 𝒫-field propagation and imaging

2.1. Setup for paraxial propagation through multiple diffusers

2.2. 𝒫-field propagator in the paraxial regime

2.3. T(ρ₁) reconstruction in the paraxial regime 𝒫-field formalism

3. Two-frequency spatial wigner distribution and occlusion-aided imaging

3.1. Setup for paraxial propagation through multiple diffusers with occlusion

3.2. Occlusion-aided imaging

4. Discussion

A. Propagation calculations

Funding

References

Cited By

Figures (4)

Equations (113)

Optics Express

Abstract

1. Introduction

2. 𝒫-field propagation and imaging

2.1. Setup for paraxial propagation through multiple diffusers

2.2. 𝒫-field propagator in the paraxial regime

2.3. T(ρ1) reconstruction in the paraxial regime 𝒫-field formalism

3. Two-frequency spatial wigner distribution and occlusion-aided imaging

3.1. Setup for paraxial propagation through multiple diffusers with occlusion

3.2. Occlusion-aided imaging

4. Discussion

A. Propagation calculations

Funding

References

Cited By

Figures (4)

Equations (113)

Optics Express

2.3. T(ρ₁) reconstruction in the paraxial regime 𝒫-field formalism