A Physical Light Transport Model for Non-Line-of-Sight Imaging Applications

In a non-line-of-sight (NLOS) imaging scenario, a target is hidden from the direct line of sight of the imaging system. NLOS 3D image recovery methods exploit the temporal information encoded in the photons that have scattered off a visible relay surface in the visible scene. A capture system illuminates a relay wall with short pulses and then retrieves the time-of-flight of photons that have scattered multiple times and interacted with the target. The goal of this paper is to provide the basis for a comprehensive mathematical framework for NLOS imaging that is directly derived from physical concepts and establishes a direct analogy between NLOS imaging systems and conventional line of sight (LOS) systems. We start by defining an irradiance Phasor Field (P-Field), an abstract quantity for the irradiance, akin to the complex envelope of the Electrical Field (E-Field). While the phase of the E-Field is lost in the reflection off a diffuse relay wall and thus destroys imaging information, phase of the P-Field remains a known and controllable quantity. We show that the P-Field propagates according to a propagator analogous to the Huygens-Fresnel propagator that propagates solutions to the EM wave equation. Using this formalism we can therefore describe NLOS imaging systems and and imaging phenomena with the same tools and processing methods that are available for LOS imaging. Simulated scenarios for different cases show that the irradiance Phasor Field is an accurate representation.


Introduction
In a LOS imaging scenario, such as the one depicted in Figure 1a, the goal is to reconstruct an image of a scene with direct LOS to the receiving camera and a transmitting light source (considered co-located in Figure 1a). After the interaction between the target and the signal emitted by the transmitter, the receiver is tasked with recovering the reflected signal from the target which is processed to create an image. In optical imaging this image formation process is typically performed by a lens.
In a NLOS imaging scenario, the goal is to reconstruct an image of a scene with a direct line of sight to neither the transmitter nor the receiver. This is achieved by illuminating a relay surface in the visible scene and collecting returned light from the object off a relay surface.
An example of NLOS optical experimental setup is shown in Figure 1b. The transmitter is a laser with nanosecond to picosecond long pulses. For the receiver, we consider an ultra-fast camera, such a single photon avalanche diode (SPAD) [1,2]. The setup also comprises a relay wall, with a diffusely reflecting ideally Lambertian surface. The light pulse generated by the laser travels to the relay wall arriving at point p and scatters from the wall in all directions; part of the photons reach the target, and a fraction of these travel back to the wall. The ultra-fast camera is focused on q and measures the light flux reflected at q as a function of time. The acquired data is a 7-dimensional space, because it is a function of the 3D coordinates of p and q, as well as a function of the time, t. The problem of reconstructing a 2D or 3D image of the hidden object is an inverse light transport problem. The light transport theory models the propagation of light through a scene (see, for example, [3]) and allows us to infer about the scene by analyzing the data captured by the camera [4].
Previous approaches have used ray optics and attempted to model the light propagation through the scene as a linear operator that can be inverted with a variety of well-studied inverse methods [5][6][7][8][9][10][11][12][13][14][15]. If posed in this way, the reconstruction problem is only approximately linear for very simplistic scenes. The ray optics model also poorly captures the underlying physical light transport processes. Nonlinear inverse methods for more complex scenes have been proposed, but the added level of complexity makes their application challenging. Model complexity and an inaccurate modeling of real light transport also make it challenging to conduct more fundamental studies. There are unanswered questions regarding null-spaces, attainable resolutions and contrast, how to deal with multiple refections in the hidden scene, the role of the bi-directional reflectance function (BRDF) of the surfaces in the scene, and relation between reconstructions and scene complexity.
In this paper, we introduce a new signal processing framework that ties NLOS TOF light transport to conventional physical light transport by introducing a new quantity, the Phasor Field (or P-Field) to describe an NLOS imaging system, such as one shown in Figure 1b, as a (a) (b) Fig. 1. A LOS scenario, shown in (a), is characterized by a target (a sphere) that is in direct path of both the transmitter (laser) and the receiver (camera). In a NLOS scenario, shown in (b), the target is hidden by an occluder from both the transmitter and the receiver.
conventional imaging system that is projected onto the relay wall. We show that the data collected by a NLOS imaging system is fundamentally equivalent to the data collected by a conventional LOS imaging system, such as a camera, ultrasound, or microwave imaging system that is placed at the relay wall and observing the hidden scene directly. The parameters of this proxy imaging system, including illumination strength, sensitivity, aperture shape and size, and wavelength are computed from the properties of the NLOS imaging system used. This approach, therefore, provides an intuitive understanding of the methodology and capabilities of NLOS imaging. Our method models the entirety of the physical light transport process as a linear, time invariant system. The motivations for this approach are as follows 1. make predictions about the best possible performance of a given NLOS imaging system, 2. translate insights into efficient imaging system design from conventional imaging systems, 3. address fundamental questions about the NLOS reconstruction process by re-formulating them as questions about conventional imaging processes, which have been extensively studied, 4. make use of the insight that NLOS TOF light transport as described by our model is fundamentally a linear, time invariant process to analyze the NLOS reconstruction process and develop optimal reconstruction methods.
The core idea of our contribution relies on defining the P-Field. This quantity is the complex envelope of optical irradiance, for which we define an amplitude and a phase, analogous to the complex envelope of the E-Field. A phasor representation of radiance was discussed in [16], where the authors propose a framework to analyze the light transport in correlation-based TOF ranging. We introduce a different phasor representation that is linked to physical light transport in two different ways: (1) A P-Field wave propagates very similarly an E-Field wave and (2) within our assumptions the P-Field propagator models the physical light transport process. This allows us to link the limitations of our reconstructions to fundamental physical limitations.

Imaging Approaches with E-Fields and P-Fields
The goal of this section is to provide the necessary background and theory that allows us to define a new abstract mathematical construct: the irradiance P-Field P. The wave equation for time-harmonic E-Fields or magnetic fields, known as the Helmholtz Equation is stated as and it describes the propagation of energy in the form of an electromagnetic wave in any linear, isotropic and homogeneous medium. Eq.(1) is satisfied by all orthogonal scalar components of the E-Field ì E. The Green's Function-based solution of the wave equation for the E-Field is given by If represents an E-Field wavelet contribution from each location (x , y ), then (2) is expressed as This is also known as the Huygens-Fresnel principle [17,18] and explains the transfer of fields from any generic location (x , y ) inside plane A to an individual location (x, y) on another plane Σ separated by a distance z from A. This is depicted in Fig.2. The quantity |r | is the absolute distance between single location (x , y ) in A and (x, y) in Σ, E 0 (x , x ) is the amplitude of the E-Field at location (x , y ), K is the wavenumber of the E-Field and K E is the E-Field coefficient of proportionality. In other words, the Huygens-Fresnel integral comprehensively describes the physical process of light transport rendering, the propagation of an E-Field distribution or a wavefront from a light source through a scene to a camera, inverse rendering, and image formation (which is the propagation a wavefront from the aperture of a camera back into the scene). The Green's function propagator e j K |r | |r | is most-commonly known as the normalized Huygens E-Field wavelet.
In this paper, we model a temporally modulated light source and camera as emitters and detectors for a virtual wave of P-Field fluctuations of optical carrier irradiance i.e., In (5), the E-Field is integrated over at least one full cycle of its oscillation period. We show that the propagation of the complex P-Field amplitude is analogous to E-Field propagation and is described by a propagator P which is analogous to the Huygens-Fresnel propagator E for E-Fields and the P-Field distribution P(x, y) in Σ is given by The quantity is a P-Field wavelet contribution from location (x , y ) in A with an amplitude P 0 (x , y ), K P is the P-Field proportionality coefficient and β is the associated P-Field wavenumber. β is expressed in terms of the P-Field wavelength λ P , the corresponding P-Field frequency ω P and the refractive index of the medium of propagation n between A and Σ as We consider E-Field on an aperture plane, A and the goal is to calculate the E-Field distribution on an observation plane, Σ, whose distance from A is set to r.
Equation (6) is sufficient for many NLOS imaging scenarios where light reflects off diffuse surfaces, phase information in the E-Field component is lost, and the resulting intensity variations from the E-Field contributions are small resulting in carrier irradiance that is approximately uniform in space and time. In such cases, where the planes considered in the model are reflective or transmissive diffusers, the P-Field propagator can be used on its own.
If the scene involves non-Lambertian surfaces or if partial coherence effects such as speckle are present, the carrier propagation makes a significant contribution to the overall result and has to be modeled if a coherent (i.e. holographic) detector is used. In this case we show that if the coherence length of the carrier is shorter than the P-Field wavelength, the full propagator for the modulated light wave can be written as the product of unmodulated carrier propagator and the P-Field propagator. This, in turn, allows us to express the magnitude |I(x, y)| of the irradiance distribution I(x, y) in Σ to be expressed as In the remainder of this paper, we derive the properties stated above and show some examples of the applications of this model. We leave the exploration of all the implications of the model and it's use in experimental NLOS reconstructions to future publications. We present a brief overview of imaging with the Huygens-Fresnel integral and its evolution into much simpler Fresnel and Fraunhofer integrals under certain approximations. Following this, we demonstrate the derivation of the phasor field integral under further constraints and show that P-Field imaging follows essentially the same propagator as E-Field imaging and therefore our existing knowledge and signal processing techniques in E-Field imaging are directly applicable to P-Field imaging.

E-Field Propagation through Huygens Integral
For the scenario depicted in Figure 2, an electromagnetic wave propagates from plane A (which could be a scene of finite dimensions or an aperture) toward plane Σ. From the Huygens principle stated in (A.4), the E-Field intensity at a single location (x, y) on the observation plane Σ is the sum of E-Field spherical wavelet contributions from every location on the aperture plane A. Thus, the E-Field intensity at location (x, y) for K E = 1/ jλ E is expressed as where t(x , y ) is the E-Field transmission function at A, χ is the obliquity factor as discussed in Appendix A while |r | is the distance between any two generic locations (x , y ) and (x, y) given by (A.2). The relationship between the wavelength of the optical carrier λ E and the E-Field wavenumber is stated in (A.3). For the Fresnel approximation in (A.6) and (A.7) and the Fraunhofer approximation in (A.9) and (A.10), we can derive the corresponding E-Field distributions E N F (x, y) in (A.8) and E F F (x, y) in (A.11) for observation plane Σ located in the near-field and far-field regions.
Depending on the location of Σ, the corresponding optical irradiance distribution I H (x, y), I N F (x, y) or I F F (x, y) can be calculated from (A.12). For instance, I H (x, y) is stated in (A.13) as where is the proportionality coefficient which signifies the dependence of I H (x, y) amplitude on λ E . A comprehensive comparison between E-Field and P-Field imaging is presented in Table 1. Imaging with P-Fields is discussed in the following section.

Phasor Field propagation
Consider an optical source of DC average irradiance amplitude α 0 with its irradiance modulated by a non-zero time-varying modulating bias signal S(r, t) of peak amplitude S 0 be represented by S(r, t) has been chosen as non-zero because S(r, t) modulates the optical irradiance which cannot be a negative quantity (because there is no notion of negative light). The resulting modulated optical intensity signal M(t) is expressed as For simplicity we consider S(r, t) a monochromatic signal with angular frequency ω P and the proposed method works identically for each frequency component. Any P-Field time signal, such as an impulse can be Fourier transformed and different frequencies treated independently. In contrast to E-Fields, P-Field signals must always have a DC component as there can be no negative intensities. The presence of this DC component does not affect the generality of what is derived here. We demonstrate, under certain conditions, the possibility to state a Huygens-like integral for P-Fields i.e., expressing the magnitude of irradiance at any location as the sum of all irradiance P-Field wavelet contributions from the aperture plane. This P-field integral would take the form where T(x , y ) is the P-Field transmission coefficient at the source plane A. We begin by deriving the propagation of the P-Field from a point source to a single point in space.

Point-to-Point Modulated Irradiance Transfer
Consider the observation plane Σ in Fig.(2) . The time-averaged irradiance of the modulated optical carrier is evaluated at each location (x, y) as the sum of E-Field contributions squared from all (x , y ) locations at plane A. The time-averaging at each point is over an integration window of an integer number of periods of the carrier wave T . Here, it is relevant to state that we use an integration time which is far greater than one time period of the optical field (E-Field) but much less than the time period of the modulating function (P-Field), i.e., 2π From standard Huygens and Fresnel formulations, the modulated irradiance contribution P(x, y) T from one location (x , y ) in A to a single location (x, y) at the observation plane Σ over one integration time window, will be the time-averaged irradiance of the modulated optical carrier. This would be expressed as the time-averaged irradiance of the unmodulated optical carrier (which oscillates rapidly compared to the integration time T ) with the magnitude of this time-averaged irradiance adjusted by the time-averaged instance S(r, t) T of the modulating function S(r, t).We call it an instance of S(r, t) because of the condition in (16) where an integration time window is much smaller than any significant change in S(r, t). Therefore, in the J th integration time window, we can express the J th instance of S(r, t) T as For the condition in (16), the J th time-averaged instance of the modulated irradiance is thus expressed as a product of the point-to-point optical irradiance contribution I[(x , y ) → (x, y)] from (11) and the instance of modulating function S J (r, t) T (refer to Appendix B). The amplitude of the J th instance of the time-averaged modulated irradiance I J (x, y) T at (x, y) is thus given by Here I[(x , y ) → (x, y)] simply describes the propagation of E-Field wavelets from a point (x , y ) to a point (x, y) and it is expressed as The J th instance of the time-averaged modulated irradiance I J (x, y) T is hence expressed from (18) and (19) as Because S J (r, t) T ≈ S(r, JT ) and a constant irradiance amplitude I 0 is calculated from the E-Field amplitude of the optical carrier E 0 using (A.12) as (20) can be expressed in terms of I 0 as The modulation function is retrieved when adequate number of instances of S(r, t) are obtained.
In this case as T 2π/ω E , S(r, JT ) → S(r, t)). The phasor equivalent form S(r) of the time-harmonic modulation function S(r, t) with a peak amplitude S 0 is given by Equation (22) can be expressed in the phasor form as Next, we extend this understanding of the transfer of irradiance between a single location (x , y ) in A to a single location (x, y) in Σ and derive expressions for irradiance transfer from all locations in A to any unique location (x, y) in Σ.

Cumulative Modulated Irradiance Contribution
To obtain the irradiance I(x, y) at any single observation location (x, y) as the sum of field contributions from all locations in A, we integrate (24) over the entire aperture. In this case, the base irradiance contribution of the optical carrier is given by the summation of all E-Field wavelets as is shown in (11). This results in the total irradiance I(x, y) at each location (x, y) in Σ to be expressed as the sum of the product I H (x, y) × S(r), i.e., For constant E-Field contribution E 0 from all locations in A (i.e., E 0 (x , y ) = E 0 ), we expand (25) while introducing a phasor proportionality coefficient K P to obtain (26) In (26), the product of integrals, denoting the irradiance contribution of the optical carrier (without the contribution of the modulating signal S(r)) is calculated from (A.13) as This product can be pulled out of the main outermost integral (as is also shown in Appendix B) in (26) if the coherence length of the carrier is shorter than the P-Field wavelength λ P to yield, (28) We also know that for the Fresnel and Fraunhofer approximations, 1/|r | 2 ≈ 1/z 2 , (28) could be expressed as (29) Therefore, from (A.13), we express (29) as a product of the E-Field Huygens contributions and the P-Field Huygens-like contributions from A as The 1/z term in (30) can be included within the definition of K P . For a constant magnitude . This implies that (30) can be expressed as In (30), P and E denote a P-Field and E-Field wavelet respectively as described in (6) and (4). The P-Field wavelets P, in this case, have a magnitude of P 0 = S 0 . Moreover, if only the magnitude of the sum of P-Field contributions are measured, then the magnitude of |I(x, y)| can be expressed as a new quantity I P (x, y) where and ∬ A Pdx dy = P(x, y) = P H (x, y).
P H (x, y) is a quantity which we call the hyper P-Field at location (x, y) and it is somewhat analogous to optical irradiance for E-Fields denoted by ì E. ì E * . The hyper P-Field is given by In the following section, we show that for plane A coinciding with a diffuse surface (which is typically the case in NLOS imaging scenario), the summations of E-Field contributions with random phase terms φ R (x , y ) leads to the a measurable P-Field distribution function P(x, y) and its magnitude I P (x, y) as was stated in (6). We show that in this case (32) can be simply stated in terms of the absolute sum of P-Field wavelets as In (35), T 0 = t 2 0 and the summation of random E-Field phase contributions σ R (x , y ) from A is stated as For the ultra-near-field case when |r | z, the Fresnel approximation is not valid, therefore a Huygens correction factor C H (x, y, |r |) has to be introduced in (29) to correctly scale the the P-Field summation in (35) and obtain the correct amplitude of I(x, y). Also the 1/z term in (35) cannot be included within K P as it remains within the integral. For this case and a plane A with roughness, (35) is expressed as We could also express (37) as The correction factor allows us to equate It has to be noted that for the case when |r | z, the error is only in the amplitude estimation and not the phase estimation through the summation of the E-Field and P-Field wavelets. Therefore, the correction factor C H (x, y, |r | only corrects for the error in the magnitude of I(x, y) in Σ without altering the distribution obtained through addition of wavelet phases. Hence with the 1/z term included within K P , (35) for ultra-near-field calculations is expressed as If C H (x, y, |r |) is also incorporated within K P , then (40) is simply expressed as For simplicity of understanding, let us consider evaluating the expression in (31) when every location (x , y, ) within a square aperture plane A) of dimensions W × W provides an equal amplitude contribution of the modulated optical carrier. The plane Σ is the observation plane set in the far-field thus satisfying the Fraunhofer approximation stated in (A.9) and (A.10). From (29), the irradiance distribution in Σ is a product of the P-Field and E-Field distributions. The resulting observation plane irradiance distribution I(x, y) and its magnitude I P (x, y) in (32) is a scaled product of the integrals ∬ (28) also yields a correct solution for near-field (Fresnel) scenarios with the adequate correction factor applied. This multiplication of the E-Field and P-Field distributions in the observation/detection plane Σ is illustrated in Fig.3. Fig.3a depicts the situation when optical irradiance is directly modulated. For the typical case when λ E λ P , the spatial Fourier Transform of the optical carrier E-Field drops to zero, forcing the product of the two distribution functions to zero at x = λ E z/W. Therefore, the resulting irradiance distribution is dominated by the null-to-null width of the Fourier Transform of the E-Field contributions (i.e., I H (x, y)).
On the contrary, the case of an ummodulated optical carrier is depicted in Fig.3b. When the optical irradiance is not modulated, as is the case in all applications involving a monochromatic unmodulated optical source, the modulation function S(r, t) can be construed as having a constant unity amplitude and an infinite modulation wavelength (λ P = ∞) which results in a constant flat spatial P-Field distribution of ∬ A Pdx dy in the observation plane Σ. The product of the two spatial Fourier Transforms is then equal to only the spatial Fourier Transform of the E-Fields (as is expected in the case when the optical carrier has not been modulated in the first place!). In Fig.(3c), the aperture plane A is considered to be 'rough'. In the next section, we show that in this scenario, the spatial E-Field distribution in Σ from ∬ A Edx dy 2 is flat with a constant amplitude over the spatial coordinates (x, y). Hence the product of E-Field and P-Field distributions take the overall shape of the P-Field distributions. Analogous to the proportionality coefficient K E in (11) for E-Fields, the P-Field proportionality coefficient K P signifies the dependence of the amplitude of the P-Field distribution |P(x, y)| in Σ on the P-Field wavelength λ P . Therefore we state As we incorporated a 1/z term within K P in (29) as well as the amplitude correction factor C H (x, y, |r |) in (41), we state K P as For the ultra near-field scenario, the correction factor has to be adjusted further to obtain an accurate P-Field estimate as was done in (39). For a fixed aperture size, the larger the magnitude of λ P , the wider the P-Field spatial spectrum (or spatial distribution) in Σ, which results in a smaller peak amplitude of the energy distribution in Σ. Due to K P ∝ 1/λ P 2 , the total energy of the sum of P-Field contributions is preserved with corresponding smaller amplitudes of |P(x, y)| for a wider P-Field spatial spectrum. Next, we discuss the case when plane A represents a diffuse surface or coincides with a diffuse surface and how the magnitude of the P-Field distribution in (3c) is obtained through the summation of incoherent E-Field contributions.

Imaging through Apertures with Roughness
In this section, we consider the modulated optical irradiance propagating from a 'rough' plane A to the observation plane Σ. We consider a plane A 'rough' if the phase of the wavefront emanating from A is altered randomly at various different locations within A. This is true for transmissive and reflective apertures within an imaging system. For the visible spectrum, a frosted glass or a rough lens would be examples of partially-transmissive rough apertures whereas a painted rough wall, which is of great interest to research in NLOS imaging, is one example of a partially-reflective rough aperture. Now let us assume that for a rough aperture, the E-Field transmission function t(x, y) is given by where t 0 is the constant E-Field amplitude transmissivity of the aperture and ∆φ R (x , y ) is a random phase variable denoting a random E-Field phase contribution from any location (x , y ) in A. This random phase randomizes the phases of the E-Field contributions φ K = K |r | by altering the E-Field phase by ∆φ K (x , y ) = ∆φ R (x , y ). With φ β = β|r | denoting the phase of the P-Fields and due to a fixed relationship between φ K and φ β , the phase of the P-Field contributions are also consequently altered by ∆φ β . Let us define a coefficient of aperture roughness γ to be equal to a few (M) wavelengths λ E , i.e., Here M denotes the worst roughness magnitude over the entire aperture. |∆φ K | in terms of |γ| is given by From (45) and (46) Moreover, from (A.3) and (8), which results in This implies that ∆φ β can be expressed in terms of ∆φ K as Consequently, we obtain a range of possible phase distortions ∆φ β to the P-Field contributions from A as For a very realistic case when λ P λ E , we know from (48) that β/K 1. Hence from (51), With the addition of random phase offsets ∆φ β (x , y ) and ∆φ K (x , y ) to φ β and φ K respectively, (29) can be expressed as Additionally, from (49), (53) can be expressed as For the case when (β K) and (52) is valid, we conclude that the roughness-induced phase fluctuations ∆φ β of the P-Field contributions are minimal (often negligible) compared to the phase fluctuations ∆φ K incurred to the E-Fields (the optical carrier). Also, when M ≥ 1, it can be concluded from (47), that |∆φ K | ≥ π, i.e., the resulting E-Field phase contributions (φ K + ∆φ K ) from all locations (x , y ) randomize completely for |∆φ K | ≥ π. Also for β/K 1 and an aperture roughness K/β ≥ M ≥ 1, the P-Field contributions from A do not randomize entirely as |∆φ β | ≤ π. For β/K ≈ 0, the E-Field and P-Field phase contributions are hence given by and Under the approximations in (55) and (56), (54) can be expressed as The factor L denotes a reduction in the number of photons arriving at each location in Σ due to Lambertian scattering [19] of photons. Therefore, we downscale the sum of photons arriving at each detector location by a scattering coefficient L where L ≤ 1. Also the summation of random E-Field phase contributions is denoted by σ R , i.e., where T 0 = t 2 0 . The magnitude |I(x, y)| of the irradiance distribution I(x, y) in (57) is I P (x, y) and it is expressed solely in terms of the summation of P-Field wavelet contributions i.e., As for the ultra-near field case with neither of the Fresnel and Fraunhofer approximations applicable, as was stated earlier, a correction factor C H (x, y, |r |) is simply introduced to account for the discrepancy in the amplitudes of I P (x, y). Equation (59) for the ultra near-field case is expressed as I P (x, y) in (60) can be approximated to With the introduction of the correction factor as was previously shown in (39), the amplitude error in (61) can be completely corrected, thereby allowing us to express (60) as The correction factor as well as the 1/z term can be incorporated within K P as well as was suggested previously. At this juncture, designating what belongs to each of the quantities P 0 and K P in (35) is entirely arbitrary. If K P = C H (x, y, |r |)/zλ 2 P as per the definition in (41) and (42), and if we define P 0 as then P(x, y) = I P (x, y) given by The possibility to detect the magnitude of P-Field distribution P(x, y) allows us to obtain a pure P-Field magnitude at each location (x, y) in Σ in the case of a rough aperture. The scenario for pure P-Field wavelet summation is depicted in Fig.(3c) as discussed in the previous section. Equation (64) can be expressed as the square root of the hyper P-Field distribution For analogous comparison, the irradiance distribution in Σ as a result of the summation of unmodulated E-Field contributions of equal amplitudes E 0 from a smooth, non-rough aperture plane A is expressed as Comparing (65) to (66), we deduce that imaging through P-Field summation in the case of a rough plane A is analogous to imaging through E-Field summation in the case of a non-rough plane A. This implies that P-Field imaging involving a partially transmissive rough aperture or In summary, we deduce that the product of P-Field and E-Field contributions is simply equivalent to the summation of all P-Field contributions from the aperture plane A when the aperture roughness is sufficient to randomize the phases of all E-Field contributions but insufficient to completely scramble and randomize the phases of P-Field contributions from A. In the context of (40), where we deduced that the eventual irradiance distribution in the observation plane is determined by the product of E-Field contributions and P-Field contributions from A, the situation of imaging through a 'rough' aperture is pictorially illustrated in Fig.3c. A summation of random E-Field phases leads to a uniform irradiance distribution in Σ (represented by a constant horizontal line in Fig.3c) and the summation of P-Field contributions from A entirely determine the eventually irradiance distribution in Σ after multiplication with a constant. The distribution of irradiance in Σ is consequently determined by the P-Field properties such as the wavelength λ P .
In the case of a rough aperture, the accuracy of the amplitude approximation in P(x, y) generally depends on the contrast between β and K. While the condition β/K 1 (i.e. ∆φ β ) holds true, the additional amplitude error in P(x, y) compared to |I(x, y)| depends on the accuracy of the standard assumption r ≈ z. For the case of an ultra near-field estimation of P(x, y) -when the amplitude estimation error is largest -the overall amplitude estimation error in P(x, y) can be entirely corrected by the appropriate correction factor C H (x, y, |r |). To avoid confusion regarding notations used throughout the manuscript, we summarize the notations used in Table.2.
We can also calculate the approximate mean percentage error η (for M × N total (x,y) locations in Σ) between the actual estimate of |I(x, y)| from (60) and the P-Field estimate I P (x, y) from (64) for a varying separation distance z between A and Σ. This average estimation error η for a given separation z between A and Σ is given by In (67), the subscripts i and j denote the i th x location and j th y location out of a total M × N locations in Σ and We can also modify the definition of η in (67) to replace the 1/z term in (67) with 1/|r | AV (x, y) where |r | AV (x, y) is the average distance between any location (x, y) in Σ and a fixed central mean location ( x , y ) in A. In this modified definition, η = η AV which is given by In (69), |r | AV (x, y) is given by In the following section, we present simulation results for the far-field, near-field and ultra near-field scenarios to demonstrate the efficacy and accuracy of P-Field imaging through rough apertures.

Phasor Field Simulations
In this section, we present simulation results to demonstrate the calculation of a P-Field distribution P(x, y) for a modulated optical carrier. In the simulations, P(x, y) is obtained through the propagation and summation of P-Field contributions from a rough aperture as per (64). We perform three simulations to estimate P(x, y) for an observation plane Σ located in the far-field, near-field and the ultra-near field regions respectively.
For each of the three simulations, we compare estimates from (64) to the distribution |I(x, y)| obtained through the Huygens integral in (54). Simulation were repeated for different magnitudes of aperture roughness coefficient γ and we show that the magnitude of P-Field distributions calculated through conventional Huygens integral approximately yield similar results with the summation of P-Fields provided that aperture roughness is large enough to scramble the optical phase of the spherical E-Field wavelets, i.e.,γ λ E but small enough compared to the P-Field wavelength, i.e., γ λ P . For all of the three simulations, we set the E-Field and P-Field wavelengths to λ E = 1µm and λ P = 10cm respectively. These wavelength values correspond to β/K = 10 −5 =⇒ K/β = 10 5 . Also, for each simulation, the constant peak amplitude of the optical carrier E-Fields at the aperture plane A was set to E 0 = 0.1V/m. From (A.12), this corresponds to a uniform E-Field-based carrier illumination irradiance of I 0 = 265mW/m 2 .
As is expected, we demonstrate a strong agreement between |I(x, y)| computed from (54) and (64) for a smaller roughness magnitude. Subsequently, we also show a deterioration in the estimate from the P-Field integral with an increase in aperture roughness. this is particularly evident for very high roughness values when γ → K/β. In other words, the difference between the P-Field amplitude estimates increases with increasing M and additional specular noise appears in the actual Huygens estimate as M is increased.
To clearly demonstrate larger errors in P-Field estimates with increasing roughness, in Fig. 6, we also plot for 'simulation 1' a 2-D P-Field distribution estimate I P (0, y) from (64) along the y-axis at x = 0. For comparison, we also plot the 'ground truth' Huygens integral-based estimates of P-Field distribution |I(x, y)| along (0, y) for different roughnesses M. We clearly observe a deterioration in the P-Field estimation when M → K/β. The mean percentage errors in the P-Field estimates H and N F compared to the Huygens and near-field Fresnel estimates are calculated respectively as and In (72), |I N F (x, y)| j and I P(j) (x, y) denote the j th instances of |I N F (x, y)| and I P (x, y) respectively for the j th unique location in Σ out of a total number J locations. From (71) and (72), we plot the mean errors e H and e N F for 'simulation 1' in Fig.(8) for a fixed ratio β/K = 1 × 10 −5 (with λ E = 1µm and λ P = 10cm). The roughness was varied such that 0 ≥ M ≥ 0.9 × 10 5 . This corresponds to a roughness variation 0 ≥ M ≥ 0.9K/β. As is expected, the estimation error for the P-Field irradiance distribution increases with increasing magnitude of |γ| from |γ| = 0 to |γ| = 4.50 × 10 4 (i.e.−4.5 × 10 4 ≤ γ ≤ +4.5 × 10 4 where |γ| = 4.5 × 10 4 =⇒ |∆φ β | = 0.9π. Additionally, to highlight the dependence of the P-Field approximation accuracy in (64) on the assumption β/K ≈ 0, we set the aperture roughness to a low value (M = 10) and reduce the value of β/K by keeping β constant to β = 20π (corresponding to λ P = 10cm) and changing the magnitude of K. For this fixed β value, K was altered such that β/K changed from 0.01 to 0.07. This was done by varying λ E in the range 1 × 10 −3 > λ E > 7 × 10 −3 (instead of of a fixed λ E = 1µm used earlier). This results in 0.7 ≥ M β/K ≥ 0.1. We plot the 2-D P-Field distribution estimates |I(0, y)| and I P (0, y) from the Huygens and P-Field integrals for these varying β/K values along the y-axis in Fig.(7). Results clearly show a good agreement between P-Field and Fresnel/Huygens estimates for low β/K values and, in accordance with (50), an increase in the estimation error for larger λ E (i.e., β/K) values.

'Simulation 2': Near-Field Imaging Scenario
In 'simulation 2', we repeat 'simulation 1' for a near-field scenario with the location of the observation plane Σ set to z = 10m. The dimensions of the aperture as well as the fixed values of λ E and λ P are set the same as for the the results in Fig.4 for 'simulation 1'. The x-y plane P-Field distributions I P (x, y) and |I(x, y)| for this scenario, computed through the integrals in (54) and (64) respectively, are plotted in Fig.5. The simulation was repeated for three roughness values of M = 0.1 × 10 5 , M = 0.3 × 10 5 and M = 0.5 × 10 5 . These values of M with β/K = 10 −5 correspond to |∆φ β | ≤ 0.1π, |∆φ β | ≤ 0.3π and |∆φ β | ≤ 0.5π respectively. As is the case with the far-field estimates in 'simulation 1', estimates using from the P-Field integral in (64) deteriorate with respect to the actual Huygens estimate in (54) with increasing M. In this case, an increasing M results in an increase of the P-Field specular noise.

'Simulation 3': Ultra Near-Field Imaging Scenario and its Comparison to Near-Field and Far-Field Scenarios
The purpose of performing 'simulation 3' is to demonstrate a larger amplitude error in the P-Field estimates of irradiance for the ultra near-field case (as was discussed in (40),) and a greater necessity to incorporate a correction factor C H (x, y, |r |) when the separation distance z is significantly reduced between planes A and Σ. 'Simulation 3' was performed for a 1m × 1m uniformly-illuminated square aperture and a varying distance 'z' between the aperture and observation planes A and Σ respectively. The dimension of Σ was set to 2m × 2m. The aperture was considered to be minimally rough, i.e., M ≈ 0. A comparison of P-Field distribution estimates between the actual and P-Field integrals in (60) and (61) for z = 1m, z = 10m, and z = 100m are plotted in Fig.9 in the x-z plane (i.e., along (x, 0, z) at y = 0). From Fig.9, we clearly see that even for an ultra near-field scenario, a P-Field distribution estimate from (61) only yields an amplitude estimation error but no phase estimation error when compared to the actual P-Field distribution calculated from (60) and consequently only a larger amplitude correction factor is required in (62) for the ultra near-field case in order to obtain the actual P-Field distribution in (60).

Application Examples
The Phasor Field virtual wave approach can be used to describe any time based lens-less imaging system, such as the ones presented in [20]. Our primary interest is however in its use for virtual camera projection in NLOS imaging as illustrated in Figure 1b. The objective is to capture an image of the scene via a diffuse reflection at a relay wall. The imaging system can project temporally modulated intensities on the relay wall and detect incoming intensities that strike the wall. The wall thus becomes the aperture of a holographic P-Field projection and detection system. After detection of the P-Field at all points on the aperture any conceivable imaging system can be realized through digital post processing. A simple imaging lens, for example would apply a position dependent phase delay to the signal followed by a summation of the fields on the camera pixels and a time integral over the absolute value to reconstruct a 2D image of a scene (cf. where P(R) is the P-Field of the on the image plane, R, P(q) is the P-Field at the relay wall at point q, Q is the set of relay wall points, and φ l is the phase shift applied by the computational lens which in a real lens is achieve by traveling through different amounts of glass. The details of this reconstruction operation and the design and demonstration of a virtual camera projection system are beyond the scope of this work. Here we present some insights about the properties of this virtual imaging system that can be derived from the phasor field formalism.

Light source and detector design
To project a light source onto the relay wall a collimated laser is usually used to create a small spot that can be treated as a P-Field point emitter. The spacing between these point emitters can be derived from phased array theory. Phased antenna arrays require an antenna spacing of λ/2 to be able to project and detect wavefronts of arbitrary shape without artifacts. Consequently the spacing of laser positions on the wall should be about half the P-Field wavelength λ P . The situation for the virtual detector positions is similar. However, an ideal imaging system would also have to maximize the total area of relay wall that is sampled, i.e. maximize the fill factor of the virtual camera. The detector array is thus composed of closely spaced area detectors of side length λ/2. While the local density of camera and projector spots should be λ/2, parts of the relay wall regions can be left out without affecting the image quality. This is analogous to blocking part of the aperture of a camera. It only lowers the signal strength but has very little effect on the image.

Resolution
It is important to realize that image quality is fundamentally only limited by the signal-to-noise Ratio (SNR). Resolution of the final image does not have to equal the resolution of the imaging system itself. Common methods of 'super resolution' work by changing the definition of resolution or by changing the properties of the scene. The resolution of an imaging system is typically defined as the size of the systems airy disc which is a measure of the highest spatial frequency that can be faithfully resolved by the imaging system. The distance from the central maximum to the first minimum of the airy disc in the far-field is given by the Rayleigh resolution limit where p E is the size of the resolved patch, r E is the distance between the patch and the aperture, λ E is the light wavelength, and d E is the diameter of the aperture. Applying the P-Field principle we can thus estimate the resolution of a virtual camera as where p P is the size of the resolved patch, r P is the distance between the patch and the aperture, λ P is the virtual wavelength and d P is the diameter of the aperture. The Rayleigh limit can be derived using the Rayleigh-Sommerfeld diffraction integral. To formally derive it for P-Fields, one simply has to replace the quantities in that derivation by the corresponding P-Field quantities.

Connection to LOS
The concept of P-Field imaging is useful for applications in NLOS imaging systems which image around a corner or a wall as in [6,13], as simple LOS imaging systems with rough apertures as is also shown in Fig. 11. NLOS systems generally consider objects to be hidden or occluded by virtue of their positions with respect to a the laser transmitter and the receiving sensor. In other words, there is no direct or straight path between the detector and the object as there are light-blocking components between the two. In such cases, a third (or possibly more) object such as a wall is used to enable light reflecting off of it to propagate from the object and reach the detector.
In a LOS scenario, as the one shown in Figure 12, we consider an imaging system, composed of a source and a camera and a point-like target, which we consider our object under investigation. The source emits a monochromatic wave that travels to the object, passing through a lens, whose goal is to refocus the electromagnetic on the object (we assume that the object is located on the focal spot of the lens). After the light interacts with the object, it travels to the camera, once again passing through the a lens which focuses the diffuse light into the receiver.
We demonstrate that a reflective structure such as a wall can be treated as a 'non-rough' P-Field aperture, similar to a P-Field lens, where despite the loss of spatial coherence of light upon reflection, the P-Field contributions are preserved enabling us to treat NLOS imaging as conventional LOS imaging but rather in the realm of P-Fields instead of E-Fields of the optical carrier. As the effect of the surface roughness is minimal to the P-Field phase as compared to the phase of the optical field, we can effectively treat any rough surface in a multi-bounce NLOS imager as a P-Field aperture of negligible roughness. Consequently, this allows us to describe a NLOS imaging system as a LOS P-Field imaging system.

Conclusions
In this paper, we introduced the concept of irradiance Phasor Fields, a quantity with an amplitude and a phase term. We show that P-Fields provide an ideal representation for imaging through apertures which exhibit a roughness in the order of the optical (E-Field) wavelength. This provides a means to model any NLOS imaging system, such as imagers that image around corners, as an LOS imaging system. The benefit of modeling with P-Fields is the inherent ability to use well-known techniques in LOS imaging to model any NLOS system by considering partially reflective and scattering rough surfaces as P-Field apertures and using existing knowledge of light transport in LOS imaging to model any NLOS system from there onwards.
We back our claims with in-depth near-field and far-field simulation results for uniformly illuminated rectangular and square apertures. We show that the P-Field approximation holds valid even for aperture roughnesses that far exceed the E-field wavelengths. We also demonstrate a reduced accuracy of this approximation with increasing roughness values.

A. Huygens, Fresnel and Fraunhofer Integrals
Let us consider a scenario as shown in Fig. 2, where we assume that an electromagnetic wave propagates from an aperture plane A to an observation plane Σ. The goal is to find the E-Field distribution and the resulting irradiance distribution over Σ. Given that (x , y ) defines a location ∈ A, a spherical scalar E-Field wavelet E(x , y , |r |) emanating from (x , y ) and observed at a location (x, y) in Σ, is expressed as In (A.1), |r | is the absolute distance between the aperture location (x , y ) and any generic location (x, y) given by and K is the wavenumber associated with each E-Field wavelength λ E (or corresponding angular frequency ω E ). K is expressed as As per the Huygens principle, the E-Field intensity at a single location (x, y) in Σ is the sum of scalar E-Field spherical wavelet contributions from every location on the aperture plane A. Thus, the E-Field intensity at location (x, y) is expressed as where we substitute A.1 in A.4a to obtain A.4b. t(x , y ) is the E-Field transmission function at A is the ratio between the transmitted and the incident E-Fields at the aperture and χ is the obliquity factor which represents the inclination between the planes A and Σ. If, for simplicity, A and Σ are considered parallel (i.e., χ ≈ 1) in (A.4) separated by a distance z λ E , then it is possible to approximate 1 |r | ≈ 1 z to express (A.4) as where we considered a uniform aperture field illumination E(x , y ) = E 0 , ∀(x , y ) ∈ A. For Moreover, the parabolic phase term, e x 2 +y 2 2z , quickly goes to zero when, in addition to the condition in (A.6), the following also holds true 2z k (x 2 + y 2 + x 2 + y 2 ). This results in a further simplification of (A.8) to yield the far-field Fraunhofer integral E F F (x, y) for E-fields which is given by (A. 13) In (A.12), ζ is the medium impedance defined in terms of the free-space impedance ζ 0 as ζ = ζ 0 µ r r where r is the relative permittivity and µ r is the relative permeability. Similarly, expressions for the Fresnel and Fraunhofer irradiance distributions I N F (x, y) and I F F (x, y) in Σ can be obtained from I N F (x, y) = |E N F (x, y)| 2 2ζ , (A.14) and I F F (x, y) = |E F F (x, y)| 2 2ζ (A. 15) respectively.

B. Separation of E-Field and P-Field integrals
Let us consider the product of the optical field squared function [E(r, t)] 2 and an amplitude modulating function S(r, t) which modulates the function E 2 . Let E(r, t) be defined as E(r, t) = E 0 sin(ω E t + φ E (r)), (B.1) where φ E (r) is the phase accumulated by the E-Field due to propagation from a location (x , y ) in A to a location (x, y) in Σ. φ E (r) = K |r |. This optical irradiance associated with the E-Field is modulated by a positive modulating function S(r, t) expressed as S(r, t) = S 0 [1 + sin(ω P t + φ P (r))].