Paraxial Theory of Phasor-Field Imaging

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.


I. INTRODUCTION
Non-line-of-sight (NLoS) imaging, colloquially known as imaging around corners, is an important and growing area of research in the imaging community [1][2][3][4][5][6][7][8]. Reza et al. [9] recently introduced the phasor-field (P-field) representation for light transport that involves diffuse reflection (such as occurs in NLoS imaging) or diffuse transmission. Attempting to apply their light transport model to NLoS geometries that include intermediate occluding objects or non-Lambertian reflections will reveal that the P-field is an insufficient summary of the underlying field at the site of such features. Nevertheless, Liu et al. [10] used the P-field approach to propose and demonstrate that line-of-sight imaging techniques can be fruitfully applied, in a computational manner, to NLoS operation, even in the presence of intermediate occluders and non-Lambertian reflections. 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 P-field propagation to scenarios more general than purely Lambertian reflections. The success of Liu et al.'s experiments is impressive and 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 P-field-already provide such a transport model, but the aforementioned work has demonstrated the intuitive utility of the P-field approach. Consequently, we believe it is worthwhile to pursue propagation primitives that can establish the P-field input-output relation for the initial and final Lambertian reflections when occluders and non-Lambertian reflectors are present in the intervening space.
In this paper, we develop a set of propagation primitives that extend the P-field formalism to scenarios that go beyond what was considered in Ref. [9]. For convenience, we assume a transmissive geometry (without reflections) that is an unfolded proxy for occlusion-aided, three-bounce NLoS imaging [11,12] and use scalar-wave, paraxial optics although these restrictions are not essential. We begin by tracing light propagation through an example transmissive geometry wherein a natural definition for the P field presents itself. Continuing this analysis, we arrive at a paraxial P-field propagator analogous to that reported by Reza et al. [9]. Using this result, we analyze the performance of P-field imaging for unoccluded transmissive geometries. Next, moving beyond the P-field, we introduce the two-frequency spatial Wigner distribution and derive 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 P-field input-output relation for occlusion-aided, diffuse-object, transmissive imaging. With that analysis in hand, we compare the P-field point-spread function for diffuse-object imaging using modulated light in the absence of an occluder with those for diffuse-object imaging using unmodulated light that is aided by the presence of either a Gaussian-pinhole occluder or a Gaussian-pinspeck occluder.

II. P-FIELD PROPAGATION AND IMAGING
In this section we consider electromagnetic field propagation through a paraxial, transmissive geometry that serves as a surrogate for an around-the-corner imaging configuration.
As in [9], we define the P-field as the Fourier transform of the short-time average irradiance.
Using this definition, we derive a formula for paraxial propagation of the P-field, which we find to be similar to the traditional Fresnel-diffraction formula for the propagation of the electromagnetic field, as reported in [9]. We then apply this understanding of the P-field to the task of imaging through diffusers and analyze the associated performance. Figure 1 shows the transmissive geometry we shall address in this paper for P-field propagation within the paraxial regime, i.e., wherein Fresnel diffraction applies. Here, E 0 (ρ 0 , t)

A. Setup for Paraxial Propagation through Multiple Diffusers
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 = (x 0 , y 0 ), and time, t. This field has center frequency ω 0 and bandwidth ∆ω ω 0 , so that the optical-frequency field is Re[E 0 (ρ 0 , t)e −iωot ]. Its units are W/m 2 , making I 0 (ρ 0 , t) = |E 0 (ρ 0 , t)| 2 the short-time average irradiance [13] illuminating the z = 0 plane. It will be assumed, in all that follows, that ∆ω is such that available photodetectors 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 1 , and z = L 1 + L 2 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 (ρ 1 ), is to be imaged using the light that emerges from the z = L 1 + L 2 plane.
can fully resolve the time dependence of I 0 (ρ 0 , t). As soon will be seen, it will be valuable to employ the time-domain Fourier transform of E 0 (ρ 0 , t), viz., for use analyzing the Fig. 1 configuration.
After propagating through the thin diffuser h 0 (ρ 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 h 0 (ρ 0 )/c time delay. Because it is unreasonable to presume we can accurately account for this delay as a deterministic quantity, we shall suppress its average value-across an ensemble of statistically identical diffusers-and consider h 0 (ρ 0 ) to be a zero-mean, homogeneous, isotropic, Gaussian random function of ρ 0 , with covariance function K h (|∆ρ|) = h 0 (ρ 0 + ∆ρ)h 0 (ρ 0 ) , where angle brackets denote ensemble average. Moreover, in keeping with h 0 (ρ 0 )'s being a diffuser, we shall take its standard deviation, σ h = K h (0) to be much greater than the center wavelength, λ 0 = 2πc/ω 0 , and its coherence length ρ c -the transverse distance beyond which K h (|∆ρ|) vanishes-to be at most a few λ 0 . Furthermore-and this condition is essential to there being a useful P-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 E 1 (ρ 1 , t), the field illuminating the z = L 1 plane.
This illumination results in being the time-domain Fourier transform of E 1 (ρ 1 , t), the field that emerges at z = L 1 + , after propagation through a deterministic thin transmission screen with intensity transmission pattern T (ρ 1 ), and a thin diffuser, h 1 (ρ 1 ), that we will take to be statistically independent of, but identically distributed as, h 0 (ρ 0 ).
Paraxial propagation to z = L 1 + L 2 , now gives us and propagation through the thin diffuser at z = L 1 + L 2 results in being the time-domain Fourier transform of E 2 (ρ 2 , t), the field that emerges at z = (L 1 + L 2 ) + . We will assume that h 2 (ρ 2 ) is statistically independent of, but identically distributed as, h 0 (ρ 0 ) and h 1 (ρ 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 = L 1 diffuser represent a Lambertian-reflecting hidden wall with spatially-varying albedo T (ρ 1 ). The z = L 1 + L 2 diffuser represents a second Lambertian reflection at the visible wall, where statistical independence from the first visible-wall reflection can be ensured by the NLoS imaging sensor's viewing a different section of that wall than what the laser illuminates. The goal of threebounce 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 P-field propagator for the preceding transmission geometry.
Equation (15)-which coincides with the result of applying the Fresnel approximation to the Rayleigh-Sommerfeld P-field propagator from Ref. [9]-is our essential result for L 1 -m paraxial P-field propagation. 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 L 1 -m paraxial P-field propagation at frequency ω − governed by a modified version of the E-field's Fresnel-diffraction formula, viz., one in which the exponent's optical frequency in the E-field Fresnel formula is replaced by the P-field's modulation frequency and the E-field formula's ω 0 /i2πcL 1 factor is replaced by the P-field's 1/L 2 1 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 which has the following pleasing physical interpretation: Paraxial propagation of the shorttime average irradiance from the diffuser's output to the z = L 1 presumes that can be employed, and results in I 1 (ρ 1 , t) being governed by the paraxial form of geometric optics, viz., the differential contribution of I 0 (ρ 0 , t) to I 1 (ρ 1 , t) is time delayed by L 1 /c + |ρ 1 − ρ 0 | 2 /2cL 1 and attenuated by the inverse-square-law factor 1/L 2 1 . Paralleling the previous development, it is now easy to show that where the averaging brackets in Eq. (19) represent averaging over the h 0 (ρ 0 ) and the h 1 (ρ 1 ) ensembles. Combining this result with what we have already obtained for relating Before continuing, it is crucial to note the behavior of P 2 (ρ 2 , 0). From Eq. (21) we immediately find that indicating that there is no spatial information about T (ρ 1 ) available in P 2 (ρ 2 , 0). This behavior is a consequence of using the paraxial approximation. Going beyond the paraxialpropagation regime-to Rayleigh-Sommerfeld diffraction-will yield a P 2 (ρ 2 , 0) containing some spatial information about T (ρ 1 ), but the inverse problem for recovering T (ρ 1 ) from P 2 (ρ 2 , 0) will still be poorly conditioned in the Fig. 1 configuration. This behavior has been seen in Refs. [11] and [12]'s 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.
Equation (21) shows that the intensity transmission pattern, T (ρ 1 ), we wish to reconstruct is illuminated by P 1 (ρ 1 , ω − ), the P-field that results from propagation of the laser illumination's P 0 (ρ 0 , ω − ) from z = 0 to z = L 1 . After transmission through T (ρ 1 ) and the diffuser h 1 (ρ 1 ), P-field propagation from to z = L 1 + L 2 results in P 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 = L 1 + L 2 plane, and let E 2 (ρ 2 , ω) be its time-domain Fourier transform. After Fresnel propagation from z = L 1 + L 2 to z = L 1 + L 2 + L 3 , propagation through the diameter-D circular-pupil, focal-length-f , thin lens, and Fresnel propagation over an additional Fourier transform given by Performing the integration over ρ 3 results in where J 1 (·) 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 h 2 (ρ 2 ) makes which together with Eq. (25) yields and hence So, by measuring I im (ρ im , t) , i.e., the diffuser-averaged, short-time average, image-plane irradiance, we obtain a 1.
. From that irradiance image we can then compute a 1.22λ 0 /D-angularresolution image of P 2 (ρ 2 , ω − ) at any modulation frequency of interest.
For reconstructing T (ρ 1 ), let us suppose that the z = 0 illumination is a t 0 -s-duration, cosinusoidally-modulated, collimated Gaussian-beam laser field where ∆ωt 0 1, i.e., with P 0 t 0 /2 being the energy illuminating the z = 0 plane. This field's short-time average irradiance is then which leads to and hence because ∆ωt 0 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 ∼ 1 m, we have cL 1 /∆ωd 2 1 from which it follows that the spatial extent of P 1 (ρ 1 , ∆ω) will be In other words, the effect of the diffuser h 0 (ρ 0 ) is to ensure that a finite, but much larger than diameter-d, region of the z = L 1 plane is illuminated by the frequency-∆ω P-field.
To proceed further, assume we have generated the computed image, of P 2 (ρ 2 , ∆ω) from the I im (ρ im , t) measurement. We can computationally invert Eq. (20) to obtain a reconstruction of T (ρ 1 )P 1 (ρ 1 , ∆ω) and use our knowledge of P 1 (ρ 1 , ∆ω) to obtain a T (ρ 1 ) image. In particular, suppose we measure I im (ρ im , t) for |ρ im | ≤ d im /2, and then defineT (ρ 1 ) bỹ where D ≡ d im L 3 /L im . Neglecting noise, and assuming that the 1.22λ 0 /D angular resolution is sufficient to makeP for |ρ 2 | ≤ D /2, Eq. (34) leads tõ Thus, over the region in the z = L 1 plane wherein |P 1 (ρ 1 , ∆ω)| has an appreciable value, the P-field imager using cosinusoidal E-field modulation at frequency ∆ω/2 achieves a spatial resolution of 1. The P-field propagator derived above does not suffice to provide a forward model for the Fig. 3 setup because of the occluders and F (ρ 1 )'s specular component. The P-field itself lacks sufficient detail to recover all relevant portions of the underlying field's interaction with these features. These features, however, can be properly accounted by using the two-frequency, spatial Wigner distribution. The spatial Wigner distribution has long been transmission mask with field-transmission function F (ρ 1 ), whose associated intensity-transmission pattern is to be imaged using the light that emerges from the z = L 1 + L 2 plane. That imaging process is aided by the presence of occluders in the z = L 1 − L d and z = L 1 + L d planes, whose field-transmission functions are P (ρ d ) and P (ρ d ), respectively.
recognized as a useful tool in optics, see, e.g., [14-16], but we need its generalization to two-frequency operation, viz., for the z-plane Fourier-domain field E z (ρ, ω) this function is The z-plane P field can be found from that plane's two-frequency spatial Wigner distribution as follows: From this result we see that the two-frequency Wigner distribution allows us to realize the goal of analyzing occluded phasor-field imaging if we can: (1) propagate W Ez (ρ + , k, ω + , ω − ) through a z-plane field-transmission mask, whether that be a diffuser, deterministic occluder, or specular-plus-diffuser mask; and (2) propagate W Ez (ρ + , k, ω + , ω − ) through L m of Fresnel diffraction. All of these propagation calculations will be done Sec. III B. For convenience, we summarize these results below: Propagation through a diffuser: .

B. Propagation Calculations
Consider propagation through one of our diffusers: assume that we know W Ez (ρ + , k, ω + , ω − ) and we want to find with In this case we immediately get Physically, the k dependence of the two-frequency spatial Wigner distribution 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 E z (ρ, ω), because Now consider propagation through a deterministic transmission mask. Here we want to For this case we have that where is the spatial Wigner distribution of P (ρ). In words, Eq. (53) shows that multiplying E z (ρ, ω) by a deterministic field-transmission mask implies that W E z (ρ + , k, ω + , ω − ) is obtained from a k-space convolution of W Ez (ρ + , k, ω + , ω − ) with the field-transmission mask's spatial Wigner distribution. Moreover, Eq. (54) immediately leads to as could have been directly obtained from Eq, (49) and the P-field's definition.
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, . 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 Our final task in this section is to find i.e., for L m of Fresnel diffraction [17]. 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. (64) 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, gives The term in Eq. (69)'s integrand behaves like the impulse δ[ρ 0+ − ρ + + kcL/(ω 0 + ω + )]. This deltafunction 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 √ ω − cL/2(ω 0 + ω + ), and ω 0 ω + implies that it integrates to one. (2) The other ρ 0 + -dependent terms in Eq. (69) are the oscillatory term, exp(iω − |ρ + − ρ 0 + | 2 /2cL), which varies much more slowly than its predecessor, because ω 0 ω − , and the Wigner distribution, whose ρ 0 + dependence can reasonably be assumed to be nearly constant over regions of diameter √ ω − cL/2(ω 0 + ω + ). So, using the delta-function approximation in Eq. (69), we get Finally, again making use ω 0 ω + , we have As a consistency check on Eq. (71), let us use it to calculate P L (ρ + , ω − ) when z = 0 illumination with two-frequency spatial Wigner function W E 0 (ρ 0 + , k, ω + , ω − ) passes through the diffuser specified in Eq. (43) before undergoing L m of Fresnel diffraction. We then have that Using Eq. (48) now gives us Changing variables so that k = ω 0 (ρ + − ρ 0 )/cL leaves us with which reduces to the result from Sec. II, by virtue of Eq. (38).

C. Relations with Other Second-Order Field Quantities
From the preceding results it follows that In general, P z (ρ, ω − ) and I z (ρ, t) provide equivalent information about the underlying (possibly random) complex field envelope E z (ρ, t), but that information is insufficient to determine W Ez (ρ + , k, ω + , ω − ), whereas knowledge of W Ez (ρ + , k, ω + , ω − ) does determine P z (ρ + , ω − ) and I z (ρ + , t) . All three of these functions are second-order field quantities, thus it is of interest to relate them to two better-known second-order field quantities: the time-dependent specific intensity from small-angle approximation linear transport theory [18], and the space-time autocorrelation function, that is used in parabolic-approximation propagation theory through random media [19].
The time-dependent specific intensity 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 be found from knowledge of the time-dependent specific intensity.
The space-time autocorrelation function is equivalent to the two-frequency spatial Wigner distribution because we have that and To complete our set of relationships, we note that and P z (ρ + , ω − ) = dt ds I z (ρ + , s, t)e iω − t = dt Γ z (ρ + , ρ + , t, t)e iω − t .

D. Occlusion-Aided Imaging
In Sec. II we noted that, in the paraxial limit, unoccluded imaging configurations without modulated light are unconditioned with respect to reconstructing the target mask's albedo.
Moreover, we showed that the addition of modulation enabled reconstruction of the target mask's albedo at a resolution limited by the bandwidth of that modulation. What remains then is to examine the unmodulated and modulated cases for occluded geometries. For clarity and convenience, we will consider a simplified version of Fig. 3 in which the first occluder is absent and the screen at z = L 1 is purely diffuse. In the NLoS analogy, this corresponds to a geometry in which a single occluding object is encountered in the hidden space only on the light's return trip from the Lambertian hidden wall. Further convenience, without appreciable loss of generality, is afforded by our assuming that the laser light incident on the z = 0 plane is a +z-going plane wave of short-time average irradiance I 0 (t), and that the distances in Fig. 3 satisfy L 1 = L 2 = L, and L d = L/2 The two-frequency, spatial Wigner function 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 find At z = L 1 this Wigner distribution encounters a diffuse target mask, i.e., one whose fieldtransmission function F (ρ 1 ) has zero mean and covariance ∆F (ρ 1 )∆F * (ρ 2 ) = λ 2 0 F[(ρ 1 + ρ 2 )/2]δ(ρ 1 − ρ 2 ), which results in Fresnel propagation to z = 3L/2 now gives us and passage through the occluder in that plane leads to Another L/2 m of Fresnel propagation then gives from which we get Now, using and changing variables to k − = k − k and k + = (k + k )/2 we have 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. (98), 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 occlusionaided 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 0 (ρ 0 )e −iω 0 t ]. In this unmodulated case we can use the usual spatial Wigner distribution, i.e., of the z = 0-plane field, in lieu of the two-frequency spatial Wigner distribution. The propagation primitives given earlier for the two-frequency spatial Wigner distribution 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. (98) assuming that E 0 (ρ 0 ) = √ I 0 is a constant, we get where and we have used the evanescence cutoff to justify replacing dk I 0 /(2π) 2 with πI 0 /λ 2 0 . Equations (100) and (101) show that this unmodulated case offers no spatial information about F(ρ) in the absence of an occluder, i.e., we get G(ρ) = π/L 2 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. (100) and (101), 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. [11], 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 where k 0 ≡ ω 0 /c = 2π/λ 0 is the wave number at the optical frequency, and Ω ≡ 4k 0 ρ 2 0 /L is the Fresnel number for pinhole's propagation geometry. The spatial resolution of G ph (ρ) improves with decreasing ρ 0 when Ω > 1, and degrades with decreasing ρ 0 when Ω < 1.
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 G ps (ρ)/G ps (∞) versus ρ/ρ res (Ω) for ρ = (x, 0) and Ω = 0.1, 1, and 10, where ρ res (Ω) is the Gaussian pinhole's spatial resolution.

IV. DISCUSSION
In summary, we have presented a complete light transport model, in phasor-field terms, capable of describing propagation through a transmissive, paraxial geometry-including intermediate occluders and a specular-plus-diffuser mask-that serves as an unfolded proxy for occlusion-aided, three-bounce NLoS imaging. For imaging purely diffuse objects without intermediate occluders we phrased our analysis in terms of the P 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 derived propagation primitives for the two-frequency spatial Wigner distribution. 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.
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 P-field imaging of specular objects. Next, because Liu et al. [10] used ps-duration pulsed illumination to obtain three-dimensional scene reconstructionsand 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 P fields. Fifth on our list is to extend our propagation primitives beyond the paraxial regime, i.e., to replace Fresnel diffraction with Rayleigh-Sommerfeld diffraction. Finally, we need to address NLoS imaging explicitly, rather than its transmissive proxy, and include more than just three-bounce returns.