Abstract
The phasor field has been shown to be a valuable tool for non-line-of-sight imaging. We present a formal analysis of phasor-field imaging using paraxial wave optics. Then, we derive a set of propagation primitives—using the two-frequency, spatial Wigner distribution—that extend the purview of phasor-field imaging. We use these primitives to analyze a set of simple imaging scenarios involving occluded and unoccluded geometries with modulated and unmodulated light. These scenarios demonstrate how to apply the primitives in practice and reveal what kind of insights can be expected from them.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
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, E0(ρ0, t) is the baseband, complex-field envelope for a quasimonochromatic, scalar-wave, modulated laser field entering the z = 0 plane, expressed as a function of the transverse spatial coordinates, ρ0 = (x0, y0), and time, t. This field has center frequency ω0 and bandwidth Δω ≪ ω0, so that the optical-frequency field is Re[E0(ρ0, t)e−iωot]. Its units are , making I0(ρ0, t) = |E0(ρ0, t)|2 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 I0(ρ0, t). As soon will be seen, it will be valuable to employ the time-domain Fourier transform of E0(ρ0, t), viz. [15],
for use analyzing the Fig. 1 configuration.After propagating through the thin diffuser h0(ρ0), the Fourier-domain field at z = 0+ is
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 h0(ρ0)/c time delay. Because it is unreasonable to presume we can accurately account for this delay as a deterministic quantity, we shall suppress its average value—across an ensemble of statistically identical diffusers—and consider h0(ρ0) to be a zero-mean, homogeneous, isotropic, Gaussian random function of ρ0, with covariance function Kh(|Δρ|) = 〈h0(ρ0 + Δρ)h0(ρ0)〉, where angle brackets denote ensemble average. Moreover, in keeping with h0(ρ0)’s being a diffuser, we shall take its standard deviation, to be much greater than the center wavelength, λ0 = 2πc/ω0, and its coherence length ρc—the transverse distance beyond which Kh(|Δρ|) vanishes—to be at most a few λ0. 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
is the time-domain Fourier transform of E1(ρ1, t), the field illuminating the z = L1 plane. This illumination results in being the time-domain Fourier transform of E′1(ρ1, t), the field that emerges at z = L1+, after propagation through a deterministic thin transmission screen with intensity transmission pattern T(ρ1), and a thin diffuser, h1(ρ1), that we will take to be statistically independent of, but identically distributed as, h0(ρ0).Paraxial propagation to z = L1 + L2, now gives us
and propagation through the thin diffuser at z = L1 + L2 results in being the time-domain Fourier transform of E′2(ρ2, t), the field that emerges at z = (L1 + L2)+. We will assume that h2(ρ2) is statistically independent of, but identically distributed as, h0(ρ0) and h1(ρ1).Before proceeding further, let us briefly comment on how the Fig. 1 geometry relates to three-bounce NLoS active imaging. The z = 0 diffuser, which is illuminated by modulated laser light, represents a Lambertian-reflecting visible wall with a uniform albedo. The combination of the intensity transmission pattern T(ρ1) and the z = L1 diffuser represent a Lambertian-reflecting hidden wall with spatially-varying albedo T(ρ1). The z = L1 + L2 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(ρ1). 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 〈I1(ρ1, t)〉, where I1(ρ1, t) ≡ |E1(ρ1, t)|2 is the short-time average irradiance illuminating the z = L1 plane and angle brackets denote averaging over the statistics of h0(ρ0). Going to the temporal-frequency domain, we have that
where * denotes complex conjugate, ω+ ≡ (ω + ω′)/2, ω− ≡ ω − ω′, and we have introduced the 𝒫 field at the z = L1 plane as the Fourier transform of 〈I1(ρ1, t)〉. Next, employing Eqs. (2) and (3), we get where, as before, ω+ ≡ (ω + ω′)/2 and ω− ≡ ω − ω′. Because Δω ≪ ω0 and σh ≪ Δλ, the preceding result can be reduced toA standard result for Gaussian random functions gives us [16]
Then, because σh ≫ λ0 and ρc ∼ λ0 we can use an impulse approximation, viz., in Eq. (11) to obtain Here, the 𝒫 field at z = 0 is 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 L1. 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 L1 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 ω0/i2πcL1 factor is replaced by the 𝒫 field’s factor. By inverse Fourier transformation of Eq. (15), we see that irradiance propagation from the diffuser at z = 0 to the z = L1 plane is governed by
which has the following pleasing physical interpretation: Paraxial propagation of the short-time average irradiance from the diffuser’s output to the z = L1 presumes that can be employed, and results in 〈I1(ρ1, t)〉 being governed by the paraxial form of geometric optics, viz., the differential contribution of I0(ρ0, t) to 〈I1(ρ1, t)〉 is time delayed by L1/c+|ρ1−ρ0|2/2cL1 and attenuated by the inverse-square-law factor .Paralleling the previous development, it is now easy to show that
where the averaging brackets in Eq. (19) represent averaging over the h0(ρ0) and the h1(ρ1) ensembles. Combining this result with what we have already obtained for relating 𝒫1(ρ1, ω−) to 𝒫0(ρ0, ω−) we getBefore continuing, it is crucial to note the behavior of 𝒫2(ρ2, 0). From Eq. (21) we immediately find that
indicating that there is no spatial information about T(ρ1) available in 𝒫2(ρ2, 0). This behavior is a consequence of using the paraxial approximation. Going beyond the paraxial-propagation regime—to Rayleigh-Sommerfeld diffraction—will yield a 𝒫2(ρ2, 0) containing some spatial information about T(ρ1), but the inverse problem for recovering T(ρ1) from 𝒫2(ρ2, 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(ρ1) reconstruction in the paraxial regime 𝒫-field formalism
Equation (21) shows that the intensity transmission pattern, T(ρ1), we wish to reconstruct is illuminated by 𝒫1(ρ1, ω−), the 𝒫 field that results from propagation of the laser illumination’s 𝒫0(ρ0, ω−) from z = 0 to z = L1. After transmission through T(ρ1) and the diffuser h1(ρ1), 𝒫-field propagation from to z = L1 + L2 results in 𝒫2(ρ2, ω−), 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(ρ1).
Let E′2(ρ2, t) be the baseband, complex-field envelope emerging from the diffuser in the z = L1 + L2 plane, and let ℰ′2(ρ2, ω) be its time-domain Fourier transform. After Fresnel propagation from z = L1 + L2 to z = L1 + L2 + L3, propagation through the diameter-D circular-pupil, focal-length-f, thin lens, and Fresnel propagation over an additional Lim distance where 1/f = 1/L3 + 1/Lim, the resulting image-plane field Eim(ρ, t) has time-domain Fourier transform given by
Performing the integration over ρ3 results in where J1(·) is the first-order Bessel function of the first kind, and we have used πD/λ0 in lieu of (ω0 + ω)D/2c in the Airy pattern because Δω ≪ ω0.The presence of the diffuser h2(ρ2) makes
which together with Eq. (25) yields and hence So, by measuring 〈Iim(ρim, t)〉, i.e., the diffuser-averaged, short-time average, image-plane irradiance, we obtain a 1.22λ0/D-angular-resolution, image of 〈I2(ρ2, t − (L3 + Lim)/c − |ρ2|2/2cL3)〉. From that irradiance image we can then compute a 1.22λ0/D-angular-resolution image of 𝒫2(ρ2, ω−) at any modulation frequency of interest.For reconstructing T(ρ1), let us suppose that the z = 0 illumination is a duration t0, cosinusoidally-modulated, collimated Gaussian-beam laser field where Δωt0 ≫ 1, i.e.,
with P0t0/2 being the energy illuminating the z = 0 plane. This field’s short-time average irradiance is then which leads to and hence because Δωt0 ≫ 1. Although this expression can be evaluated analytically, we shall not bother. We just note that with Δω/2π ∼ 1 GHz, d ∼ 1 mm, and L1 ∼ 1 m, we have cL1/Δωd2 ≫ 1 from which it follows that the spatial extent of 𝒫1(ρ1, Δω) will be ∼cL1/Δωd ≫ d. In other words, the effect of the diffuser h0(ρ0) is to ensure that a finite, but much larger than diameter-d, region of the z = L1 plane is illuminated by the frequency-Δω 𝒫 field.To proceed further, assume we have generated the computed image,
of 𝒫2(ρ2, Δω) from the 〈Iim(ρim, t)〉 measurement. We can computationally invert Eq. (20) to obtain a reconstruction of T(ρ1)𝒫1(ρ1, Δω) and use our knowledge of 𝒫1(ρ1, Δω) to obtain a T(ρ1) image. In particular, suppose we measure 〈Iim(ρim, t)〉 for |ρim| ≤ dim/2, and then define T̃(ρ̃1) by where D′ ≡ dimL3/Lim. Neglecting noise, and assuming that the 1.22λ0/D angular resolution is sufficient to make for |ρ2| ≤ D′/2, Eq. (34) leads to Thus, over the region in the z = L1 plane wherein |𝒫1(ρ1, Δω)| has an appreciable value, the 𝒫-field imager using cosinusoidal E-field modulation at frequency Δω/2 achieves a spatial resolution of 1.22ΔλL2/D′, where: Δλ = 2πc/Δω; L2 is the distance from the transparency-containing plane to the plane visible to the sensor; and D′ = dimL3/Lim, with L3 being the distance from the plane visible to the sensor to the sensor’s entrance pupil, Lim being the distance from that entrance pupil to the image plane where irradiance measurements are made, and dim 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 = L1 − Ld and z = L1 + L′d planes, and the z = L1 plane contains a field-transmission mask F(ρ1) that has both specular and diffuse components. In the NLoS analogy, the two occluders represent objects in the hidden space—encountered by the light as it propagates towards and returns from the hidden wall, respectively—and the generalized field-transmission mask accounts for more general, non-Lambertian hidden walls. This configuration—if F(ρ1) is purely diffuse with a space-varying albedo that is to be imaged, i.e., equivalent to the stacked intensity-transmission mask and thin diffuser from Fig. 1—is our unfolded proxy for Xu et al.’s experiments [12].
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 = L1 − Ld, z = L1, or z = L1 + 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]:
In computer vision, this quantity is known as the 5D light field [18–20]. By replacing 2πs/λ0 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., which has long been recognized as a useful tool in optics, see, e.g., [21–23]. The short-time average irradiance is obtained from Iz(ρ+, s, t) by integrating out its directionality information, 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:
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): from which the time-dependent specific irradiance can be obtained viaThe merit of the TFSWD’s added generality can be seen by considering the space-time autocorrelation function,
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 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 where t+ ≡ (t1 + t2)/2, t− ≡ t1 − t2, andFor E-field propagation through an arbitrary linear transformation of the form
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, 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:
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
. - Propagation through a deterministic occluder:
With WP(ρ+, k) ≡ ∫d2ρ− P(ρ+ + ρ−/2)P*(ρ+ − ρ−/2)e−ik·ρ−, we have
- Propagation through a specular-plus-diffuser mask:
With F(ρ1) having nonzero mean 〈F(ρ1)〉 ≠ 0, and covariance, where 0 ≤ ℱ(ρ+) ≤ 1, we get
- Fresnel diffraction:
For Fresnel diffraction from the z = 0+ plane to the z = L1 − Ld plane, we get
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 = L1 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 I0(t), and that the distances in Fig. 3 satisfy L1 = L2 = L, and Ld = L/2
The TFSWD of the plane-wave laser light is easily shown to be
where After the diffuser in the z = 0 plane we get and after propagation to the z = L plane, we findAt z = L1 this Wigner distribution encounters a diffuse target mask, i.e., one whose field-transmission function F(ρ1) has zero mean and covariance , which results in
Fresnel propagation to z = 3L/2 now gives us and passage through the occluder in that plane leads to Fresnel propagation over another L/2 distance then gives from which we getNow, using
and changing variables to k− = k − k′ and k+ = (k + k′)/2 we have where we have suppressed the ρ+ argument of 𝒫0(ρ+, ω−) because that field has no such dependence for the plane-wave source we have assumed. We define a new function With this definition we have Changing variables again, ρ̃ = ρ+ − cLk+/ω0, we get our final result 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[E0(ρ0)e−iω0t]. In this unmodulated case we can use the usual spatial Wigner distribution, i.e.,
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 is a constant, we get where and we have used the evanescence cutoff, |k| ≤ 2π/λ0, to justify replacing ∫ d2k I0/(2π)2 with .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(ρ) = π/L2 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
and the Gaussian pinspeck, where ρ0 is the e−1/2-attenuation radius of the Gaussian functions. The Gaussian-pinhole camera can be analyzed with far less complication than our approach to obtaining Eqs. (70) and (71), but (after accounting for image inversion) its point-spread function (psf) Gph(ρ) is revealing. The Gaussian-pinspeck camera, on the other hand, is more relevant to the experiments of Xu et al. [12], but its psf Gph(ρ) 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
where k0 ≡ ω0/c = 2π/λ0 is the wave number at the optical frequency and is the Fresnel number for the pinhole’s propagation geometry. The spatial resolution of Gph(ρ) improves with decreasing ρ0 when Ω > 1, and degrades with decreasing ρ0 when Ω < 1. Thus the Gaussian pinhole’s resolution-optimized psf, is obtained when . The optimized psf’s spatial resolution—taken to be its e−π-attenuation radius—is then , 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 Gph(0) and its spatial resolution changing, i.e., we have that for the Gaussian pinhole.For the Gaussian pinspeck, we get
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 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 Gps(ρ)/Gps(∞) versus ρ/ρres(Ω) for ρ = (x, 0) and Ω = 0.1, 1, and 10, where ρres(Ω) is the Gaussian pinhole’s spatial resolution.Note that in the near-field region, wherein Ω ≫ 1, Eq. (77) reduces to the geometric optics result,
which is analogous to the geometric optics treatment used by Xu et al. [12] and Thrampoulidis et al. [13] for the hard-aperture, circular occluder4. 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 , where
with In this case we immediately get 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 is independent of k.Propagation through a deterministic occluder:
Now consider propagation through a deterministic transmission mask. Here we want to find given Wℰz (ρ+, k, ω+, ω−) and a deterministic P(ρ), where
For this case we have that where is the spatial Wigner distribution of P(ρ). In words, Eq. (91) shows that multiplying ℰz(ρ, ω) by a deterministic field-transmission mask implies that 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 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(ρ1) having nonzero mean 〈F(ρ1)〉 ≠ 0, and covariance, where 0 ≤ ℱ(ρ+) ≤ 1 and ΔF(ρ) ≡ F(ρ) − 〈F(ρ)〉. The propagation analysis follows from combining the two previous analyses:
From expanding F(ρ) into a sum of its (deterministic) mean and zero-mean random portions, it follows that Fresnel diffraction:Our final task is to find WℰL (ρ+, k, ω+, ω−) when
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 Exploiting Δω ≪ ω0, and making the coordinate transformation from ρ0 and ρ′0 to ρ0+ ≡ (ρ0 + ρ′0)/2 and ρ0− ≡ ρ0 − ρ′0, we can reduce Eq. (101) to Rearranging terms allows us to put the ρ− integral inside the ρ0+ and ρ0− integrals, i.e., Performing the ρ− integral then yields which, after some terms cancel, givesThe term
in Eq. (106)’s integrand behaves like the impulse δ[ρ0+ − ρ+ + kcL/(ω0 + ω+)]. 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/(ω0 + ω+) with nominal width , and ω0 ≫ max |ω+| implies that it integrates to one. (2) The other ρ0+-dependent terms in Eq. (106) are the oscillatory term, exp(iω−|ρ+ − ρ0+|2/2cL), which varies much more slowly than its predecessor, because ω0 ≫ max |ω−|, and the Wigner distribution, whose ρ0+ dependence can reasonably be assumed to be nearly constant over regions of diameter . So, using the delta-function approximation in Eq. (106), we get Finally, again making use ω0 ≫ ω+, we haveAs a consistency check on Eq. (108), let us use it to calculate 𝒫L(ρ+, ω−) when z = 0 illumination with TFSWD Wℰ0(ρ0+, k, ω+, ω−) passes through the diffuser specified in Eq. (81) before undergoing Fresnel diffraction over a distance L. We then have that
Using Eq. (86) now gives us Changing variables so that k = ω0(ρ+ − ρ0)/cL leaves us with which reduces to the result from Sec. 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 Ta satisfying ω0Ta ≫ 1 and ΔωTa ≪ 1.
15. In what follows, integrals without explicit limits are over the integration variable’s entire domain.
16. Because h0(ρ) is a zero-mean Gaussian process, its samples at ρ0 and ρ′0 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.