Incorporating compact, computing-empowered diffractive optics into a wide range of imagers and sensory systems has emerged into a demanding research field, with its promising applications in holography [1], lens design [2,3], optical computing [4], and computational imaging [5–7]. However, most reported solutions are demonstrated with a limited field of view, due to the lack of an accurate yet efficient off-axis diffraction modeling scheme in the design space. That said, although fundamentals of diffraction have been well established [8,9] with many algorithms developed for numerical implementation [10,11], modeling light propagation at non-paraxial and large-scale fields is still problematic, requiring a significantly large number of samples. Even off-the-shelf optics design software such as Zemax [12] cannot faithfully support the simulation in a computing resource-efficient manner.

To simultaneously quest for computational accuracy and efficiency, one often adopts the angular spectrum method (ASM) based on fast Fourier transform (FFT) for off-axis diffraction modeling [9]. Unfortunately, conventional FFT-based ASM needs non-trivial zero-padding to avoid aliasing [13], suffering from a severe trade-off between computational efficiency and accuracy. Representative solutions to alleviate this sampling problem include the shifted band-limited ASM [14], generalized transfer function method [15], and shifted band-extended ASM (shift-BEASM) [16]. However, the samples are still rigidly determined subject to the respective requirements of the naive incident field and transfer function, which overlooks the intrinsic characteristics of the incident field and the link between the angular spectrum and the transfer function, resulting in excessive oversampling or harmful undersampling. More specifically, these works assume the sampling of the transfer function is independent of the incident field, which is incomplete since the correct computation of the final diffractive field requires a proper sampling of the product of the angular spectrum and the transfer function, rather than that from independent derivations.

To address the above issues, herein, we make the following contributions: (1) establishing an effective modeling scheme that enables accurate, fast, and flexible wave propagation simulation with a large incident angle, termed least-sampling ASM (LS-ASM), which is achieved by identifying the optimal sampling strategy of ASM; (2) employing the shifting property of the Fourier transform to eliminate redundant samples induced by the off-axis configuration; (3) developing a universal approach to identifying both lower and upper sampling bounds that ensure correct results for all types of fields. The proposed method is implemented with a matrix-based discrete Fourier transform, enabling an arbitrary selection of sampling windows and intervals.

Unlike prior work, where the off-axis diffraction is considered with plane waves, we start with a more general and practical case. A spherical wave emitted from an off-axis point source at $({x_0},{y_0},{z_0})$ impinges the aperture plane with the incident angle defined by $\sin {\theta _x} = {x_0}/{r_0}$ and $\sin {\theta _y} = {y_0}/{r_0}$ (${r_0} = \sqrt {x_0^2 + y_0^2 + z_0^2}$ is the distance from the point source to the origin of the aperture plane $O$), yielding a diffractive field on the observation plane [Fig. 1(a)]. Here we assume a circular aperture with a diameter of $D$ for simplicity. Using ASM, the diffractive field at a distance $z$ can be calculated by [9]

 

Fig. 1. (a) Diagram of off-axis diffraction. (b) The LPC shifts the angular spectrum of the input field to the frequency center. $\delta \eta$ is the sampling interval along $\eta$ axis in the spatial domain.

Download Full Size | PDF

In the scope of this paper, we mainly discuss the sampling of the input field’s phase with the temporal assumption that the amplitude is a gradually varying function of $({\xi ,\eta})$ and the phase is smooth. Without loss of generality, here we employ an additional phase modulation to the incident field, which is denoted by ${\phi _{{\rm mod}}}$. Then the phase of the input field can be expressed as

Analyzing the effective angular spectrum of the incident wave reveals that it has asymmetric coordinates relative to the origin, primarily caused by a linear phase component. This asymmetry causes an increase in the sampling rate as the effective angular spectrum is located in the high-frequency area, according to the Nyquist sampling theorem [9]. To tackle this problem, we apply a linear phase compensation (LPC) to the input field, such that the center of the effective angular spectrum can be shifted to the origin, as determined through an exploration of the shifting property of the Fourier transform, i.e.,

The Nyquist sampling barrier of the input field is now lifted. In contrast, without LPC, the naive ${\phi _u}$ requires additional samples due to the band asymmetry to avoid aliasing. Figure 1(b) illustrates that sampling an off-axis spherical wave following the reduced requirement results in an aliased angular spectrum, while with LPC, we correctly calculate it.

Although the LPC reduces the sampling requirement of off-axis waves, they still need to be properly sampled according to the requirement of on-axis waves. Notably, while phase gradient analysis can well determine the on-axis wave’s bandwidth, the truncation by the finite-size aperture broadens the angular spectrum to the infinite. Here we seek a straightforward solution based on the optical insight that a positive thin lens’s back focal plane is the Fourier transform of the field at the front focal plane [9].

We apply a virtual lens to focus light in a limited region on the back focal plane, acting as a Fourier transform operator and identifying the effective bandwidth of the input field, i.e., the on-axis wave [Fig. 2(a)]. The angular spectrum of the input field can be viewed as the ideal one convolved with the Fraunhofer diffraction field of the aperture without modulation. Therefore, the effective bandwidth is the sum of the phase gradient term and the size of the Airy disk. Given the focal length of virtual lens ${d_{\rm v}}$, the effective width of the Airy disk is specified as $W = \frac{{41.2\lambda {d_{\rm v}}}}{D}$, which covers around 99% energy [8].

 

Fig. 2. (a) Fourier transform interpreted as a virtual thin lens model. The green and yellow areas denote spherical and plane waves, respectively. (b) Influence of oversampling factor. The inner squares indicate the sampling region of a spherical wave (left) and a plane wave (right) with oversampling factor $s$ of one (top) and 1.2 (bottom), respectively. A larger region is visualized than the sampling rate to illustrate the aliasing artifact.

Download Full Size | PDF

The virtual lens can also determine the effective bandwidth of a spherical wave under the paraxial assumption, but it may be inaccurate in cases where the assumption is invalid or when the input field has a complicated phase. In such cases, we adopt the phase gradient analysis, which determines bandwidth from maximum gradients (Supplement 1, Sections 1B and 1C). Finally, the sampling rate of the input field can be well determined as

After establishing an optimal sampling for the input field in the spatial domain, we examine the counterpart in the frequency domain. We demonstrate that manipulating the observation window (OW) can substantially decrease the sampling requirement in this domain. Specifically, we virtually shift the OW from $({{x_c},{y_c}})$ to $({{x_c} - \Delta x,{y_c} - \Delta y})$ and incorporate the shift as an additional linear phase in the transfer function, expressed as

Similar to the spatial domain, the frequency-domain sampling rate must be analyzed collectively, taking into account the angular spectrum $U$ of the input field and the shifted transfer function ${\hat H_z}$. While the angular spectrum is generally not expressible in a closed form, it can be approximately analyzed via the virtual lens and an upper-bound approach, which includes amplitude sampling as it cannot be assumed to be slowly varying. We also observe that ${\hat H_z}$ can be regarded as a converging spherical wave with a negative quadratic phase. Depending on the characteristics of the angular spectrum, the quadratic phase may counteract or accumulate to the corresponding one of $U$ when considered together, resulting in fewer or more samples than those of sampling $U$ and ${\hat H_z}$ separately, as illustrated in Fig. 3. The maximum gradient of ${\phi _U} + {\phi _{{{\hat H}_z}}}$, which is at the boundary, is much smaller than that of ${\phi _U}$ and ${\phi _{{{\hat H}_z}}}$. In addition, a broadening factor is included to account for the frequency cut-off, similar to the effect induced by the finite-sized aperture in the spatial domain. The detailed derivations of combined sampling are in Supplement 1, Section 2B.

 

Fig. 3. Effect of combined sampling in frequency domain. From left to right, the first column is the zoomed view from the right ones. Second to fourth columns are the phase of the angular spectrum of an off-axis converging spherical wave with LPC, the shifted transfer function, and the superposition, respectively.

Download Full Size | PDF

Determined by not only the combined sampling of $U$ and ${\hat H_z}$ but also the OW size $w$, the sampling rate in the frequency domain with $f \in \{{{f_X},{f_Y}} \}$ can be derived as

Adaptive sampling involves arbitrary samplings across the two domains, which cannot be handled by conventional FFT. To improve computational flexibility, algorithms such as chirp-Z transform (CZT) [17], Non-uniform FFT (NUFFT) [16], or matrix triple product (MTP) [18–20] can be used. However, CZT requires three conventional FFTs, and NUFFT often suffers from numerical precision issues. As such, we adopt MTP due to its numerical stability and intrinsic linear algebraic operation, which can be highly optimized with modern parallel computing techniques. To provide a detailed explanation of LS-ASM’s procedure, we summarize the overall computational flow in Supplement 1, Section 3.

We validate the proposed LS-ASM by comparing it with the state-of-the-art shift-BEASM [16] to compute the off-axis complex point spread function (PSF) for a thin lens camera at 500 nm wavelength with a focal length of 35 mm and $f$-number of 16. The PSF is formed by focusing an off-axis spherical wave that originates from a distance of ${z_0} = 1.7\;{\rm m} $ away from the aperture. The OW size is set to $D/8$ with $512 \times 512$ pixels. The oversampling factors in both domains are set as 1.5. All numerical simulations are run on a workstation with an AMD Ryzen Threadripper Pro 3955WX CPU, using $\texttt{Python 3.7, PyTorch 1.10}$ for the matrix product in LS-ASM, and the $\texttt{finufft}$ library [21] for NUFFT implementation in shift-BEASM.

Figure 4 presents visual examples at a 3° incident angle and quantitative results from 0° to 20°. LS-ASM follows the sampling strategy mentioned before, while shift-BEASM is implemented in two different cases to highlight the superiority of our approach. The results from the first Rayleigh–Sommerfeld (RS) integral are used as the ground truth [9]. We assess the accuracy of output fields using the signal-to-noise ratio (SNR) [16] and root-mean-square error (RMSE). For visualization purposes, a common linear phase term $k({x\sin {\theta _x} + y\sin {\theta _y}})$ is removed from the results, and the remainder is denoted as ${\hat \phi _{{u_z}}}(x,y)$.

 

Fig. 4. (a) Visual results of the amplitude (first row) and phase (second row) of the complex PSF ${u_z}(x,y)$, and the amplitude for the angular spectrum $U({f_X},{f_Y})$ (third row) of the input field at 3°. The amplitude of the PSF and the angular spectrum is normalized by respective maxima. In the third row, the red dots denote the center of the angular spectrum. (b) SNR and runtime w.r.t. incident angles. (c) Sampling number w.r.t. incident angles. More results with different incident angles and metrics can be found in Supplement 1, Section 4.

Download Full Size | PDF

In case 1, where the samples in shift-BEASM are set to achieve the same level of SNR as LS-ASM, LS-ASM consistently performs at a remarkably short runtime even at large angles (e.g., 20°), outperforming shift-BEASM by over $35 \times$. This is because the effective angular spectrum of shift-BEASM is decentered, requiring more samples to avoid aliasing, whereas LS-ASM always places the effective angular spectrum at the origin in the frequency domain. In addition, LS-ASM optimally specifies the sampling rate in the frequency domain, whereas shift-BEASM does not. It is worth noting that RS is in fact biased towards shift-BEASM as the unshifted transfer function is equivalent to the Fourier transform of the first RS impulse response [9]. Nevertheless, LS-ASM still exhibits $7 \times$ fewer spatial samples and $3 \times$ fewer frequency samples at 20°, resulting in less memory usage and enables larger-batch parallel computing, which further accelerates the computation.

In case 2, where the samples in shift-BEASM are set as the same rate as LS-ASM, the effective angular spectrum of the input field is partially cropped due to insufficient spatial samples, resulting in an incorrect result even at 5°. In contrast, LS-ASM enables efficient computations at much larger angles.

We further verify our method by showcasing the results at an ultra-large angle of 35°, the computation of which requires only 1.2 s. LS-ASM uses $17 \times$ and $3 \times$ fewer samples in each domain than shift-BEASM. Another example of a complex input field is presented with a cubic phase plate placed against the thin lens at 10°. The cubic phase is a special case that possesses an intrinsic frequency offset but without any linear phase terms. Our method can accurately handle such cases. In addition, a comparison with Zemax and an ablation study on the oversampling factor $s$ is also presented, from which we observe a balance between accuracy and runtime when $s = 1.5$ (Supplement 1, Section 4).

While our scheme has demonstrated flexibility and efficiency with assumptions of slowly varying amplitude and smooth phase, it applies also to non-smooth, i.e., suddenly varying, phase and amplitude. This can be achieved by applying Eq. (4) to local sub-apertures where the field is smooth. The largest incident angle that our method can handle is limited by the validity of the scalar diffraction theory and available computational resources. In addition, although a long propagation distance may increase the sampling burden, it can be alleviated by the semi-analytical Fourier transform method, which enables analytical computation of quadratic phase [22]. Relevant analysis and demonstrations are provided in Supplement 1, Section 5.

In conclusion, we have explored a flexible and efficient numerical wave propagation modeling scheme that is tailored to represent off-axis diffractive propagation. Compared to prior works, our approach provides the LS strategy for optical fields and is applicable to all kinds. We envision a wide range of applications where efficient and accurate modeling of wave propagation is essential, particularly where a large off-axis incident angle is expected.