1. INTRODUCTION
Numerical simulation of light behavior has greatly boosted the development of modern optics, including lens design [1–3], optical imaging [4–6], holographic display [7–12], laser manufacturing [13,14], optical trapping [15], and optical computing [16,17]. Wave propagation modeling plays a vital role in computational optics. Particularly, the simulation between non-parallel planes provides remarkable flexibility, yielding benefits in diverse fields. For instance, wave propagation through micro-lenses can be efficiently simulated by discretizing the curved interface into tangential plane elements [18]. In microscopy, one may calculate the point spread function on a tilted plane [19]. Moreover, it enables high-fidelity three-dimensional holographic display [20–22], and generates an optical focus array on a tilted plane to efficiently trap microparticles [23]. As a result, there has been significant interest in developing methods that can accurately and efficiently simulate wave propagation between arbitrarily tilted planes. When dealing with parallel planes, one may adopt the angular spectrum method aided with fast Fourier transform (FFT) to achieve both high computational accuracy and efficiency [24–26]. However, it becomes more challenging for non-parallel planes, as the (spatial) frequency coordinate is not directly linked with the desired spatial counterpart by Fourier transform, which restricts the flexibility of computational optics.
Recently, several solutions have been developed, which can be roughly categorized into three groups. The first category involves the combination of conventional FFT with angular spectrum interpolation [27–29]. This approach interpolates the angular spectrum on an arbitrarily tilted plane from that on a parallel plane and employs FFT to derive the field distribution, allowing for comparable efficiency to that of the parallel case. However, the transform between two uniform frequency domains introduces interpolation errors and fails at large intersection angles due to the additional Jacobian [27]. Alternatively, the second category involves the direct computation of the diffractive field with nonuniform FFT algorithms, eliminating the explicit interpolation [21,30]. While this approach partially guarantees computational accuracy, it often takes more than tenfold time compared to standard FFT. Notably, instead of directly resorting to diffraction theories, emerging methods use deep learning to approximate wave propagation. These methods have been employed to improve computational efficiency when dealing with freeform surfaces [10] and close the gap between simulations and experiments [11,12]. However, learning-based methods require significant time and data for training. Additionally, while they perform well on fixed configurations, they may still encounter errors when generalized to different setups without retraining.
In this work, we present a novel approach for wave propagation modeling between two arbitrary non-parallel planes. Our method involves a reorganization of the angular spectrum from a uniform frequency domain to a nonuniform one, which incorporates projection and selective merging. Subsequently, a matrix-based Fourier transform is applied to enable flexible but efficient computations; it applies to both scalar and vectorial diffraction modeling. Compared to existing methods, our algorithm enables a concurrent enhancement in accuracy and efficiency, even between two orthogonal planes. By tuning the merging threshold, we can also control the accuracy and efficiency as needed. In addition, we experimentally explore its application in computer-generated holography (CGH).
2. METHOD
Before introducing our method for wave propagation modeling between non-parallel planes, it is beneficial to revisit the parallel case. Here we consider two theories, the scalar angular spectrum method (SASM) [25] and the vectorial one (VASM) [31–33]. SASM can be used when the polarization distribution of the beam is not complex or the numerical aperture (NA) of the system is not too large. However, it becomes insufficient when the beam’s polarization varies rapidly on the source plane, such as cylindrical vector beams [34]. On the other hand, a high NA, corresponding to a large diffraction angle and bandwidth, necessitates the use of VASM because polarization significantly affects the diffractive field in this case [35]. For instance, in typical holographic display configurations, the beam usually exhibits simple linear polarization, and the system’s NA is low [22]. Under these conditions, SASM can yield accurate results. However, when the polarization becomes complex or the system’s NA increases a lot [36–38], we must rely on VASM.
Theoretically, SASM can be generalized to VASM by treating each polarization component independently. This extension is straightforward in low-NA cases. However, it becomes challenging for high-NA systems because directly extending SASM requires a large number of samples to discretize a highly converging spherical wavefront. Consequently, a more popular VASM derived from the far-field approximation is often used [33]. In this approach, the converging spherical wavefront is analytically removed, requiring only one Fourier transform. Therefore, we adopt this more popular version of VASM to demonstrate our method on vectorial diffraction modeling. While SASM and VASM differ in the treatment of light polarization and are preferred for optical systems with different NAs, they share the same fundamental interpretation of diffraction—the propagation of the angular spectrum. Notably, the transformation from the source field to the frequency domain is often different in the two methods. The diagrams of scalar and vectorial diffraction are illustrated in Supplement 1, Figs. S1 and S2. Herein, we only present the main formulas relevant to our method; the detailed derivations can be found in Supplement 1, Sections 1 and 2.
Regarding SASM, we assume that the field on the source plane is s?s(?,?); then the diffractive field ? on the parallel observation plane (?,?) at a distance ? can be calculated by [25,27,28]
where ? and ?−1 represent the Fourier transform and its inverse, respectively. ?0 is the angular spectrum in the frequency domain rT?r=(??,??,??)T (T represents transpose). Notably, the coordinate symbols, such as ? and ?, represent the discrete samples instead of their continuous counterparts in the real world.
VASM, unlike the scalar case, often applies to calculate the tightly focused field of a high-NA objective lens [33,39], where the source plane is located at the entrance pupil of the objective lens. Given the source field s?s, the diffractive field ? on the parallel plane with a distance of ? from the focal plane can be formulated as [24,29,31–33,39,40]
where ? and ? are a coefficient and transformation matrix, respectively, which link the source field with the angular spectrum ?0. Notably, we define the source field in frequency coordinates (??,??) instead of the spatial ones (?,?), following the relation (?,?)=(−???/?,−???/?) since the focused field can be completely determined by the NA of the objective lens and the refractive index ? in focal space. Here, ?=2??/?0 is the wavenumber in medium with the vacuum wavelength ?0, and ? represents the focal length of the objective lens.
To extend these theories to arbitrary non-parallel planes, additional steps are required to account for the relationship between the spatial and frequency coordinates. The diffraction modeling, as illustrated in Fig. 1(a), comprises three steps. First, the field defined on the source plane (?,?) undergoes a transformation into the frequency domain r?r, yielding the angular spectrum. Next, the angular spectrum is projected onto a new coordinate system tT?t=(??,??,??)T, depending on the intersection angle (?,?) between the two non-parallel planes. Herein, ? and ? denote the polar and azimuthal angles, respectively, of the spherical coordinate system. By merging angular spectrums located at similar ?? or ??, a condensed angular spectrum distribution is obtained. Finally, the diffractive field on the observation plane (?,?) is computed as the Fourier transform of the angular spectrum.
Fig. 1. Wave propagation between two non-parallel planes. (a) Coordinate transformation, where the ?? plane can be viewed as the ?? plane that undergoes two successive rotations around the ?-axis by an angle of ?, and then ?-axis by an angle of ?. (b) Principles of our method.
The projection between r?r and t?t follows tr?t=??r, where ? is the projective matrix [Supplement 1, Eq. (S4)]. Here we reasonably assume that the source (?,?), observation (?,?), and (??,??) planes are all uniformly sampled. As ??=?2−??2−??2, ?? is nonuniformly sampled, indicating that the angular spectrum matrix in t?t space is nonuniformly and sparsely distributed when the source and observation planes are not parallel [Fig. 1(b)]. Such a matrix retains the potential to be further compressed to enhance computational efficiency.
As a result, we consider merging angular spectrums located at similar ?? or ?? to achieve a compact angular spectrum distribution. Specifically, our merging strategy consists of three steps. First, identify similar ?? or ?? based on a predefined threshold. Next, find an equivalent ?? or ?? for these similar ?? or ??. Finally, superpose the corresponding angular spectrums. As illustrated in Fig. 1(b), these angular spectrums with similar ?? or ?? in the yellow and red boxes are superposed, respectively, yielding only one element at one ?? or ?? position. This process can be interpreted by approximating angular spectrums, i.e., plane waves, in similar directions to one certain direction, leading to a much tighter angular spectrum matrix. As for the selection of the merging threshold, we assume that the number of samples after merging in the frequency domain approaches, for example, ? samples for ??. Specifically, we first sort ?? and determine an initial threshold according to the ?-th largest gap between adjacent samples. Then we iteratively adjust the threshold to ensure that the number of samples after merging is close to ?.
Supposing we have identified similar ?? as {??(1),??(2),?,??(?)} along with their corresponding angular spectra {?(1),?(2),?,?(?)}, where ? is a positive integer, the equivalent ?? is determined by the weighted average
The equivalent ?? can be determined similarly. This merging strategy, together with the aforementioned projection, is termed angular spectrum rearrangement. Notably, by tuning the merging threshold, i.e., the maximum interval between two ?? or ?? allowed for merging, we can control the samples in the frequency domain, thereby flexibly balancing the computational accuracy and efficiency.
Although the angular spectrum matrix has been condensed, it remains nonuniform. To calculate the diffractive field on (?,?) plane from the nonuniform angular spectrum in t?t space, we employ the matrix triple product (MTP) as a flexible but efficient two-dimensional Fourier transform [41,42]. This allows handling the case of nonuniform sampling without interpolation and enables arbitrary selection of the sampling intervals and regions without zero-padding. This operation can be described by
where T?(?,??)=exp?(???T?) and T?(?,??)=exp?(??T??). ??, ??, ?, and ? are the coordinates, represented by row vectors, in the frequency and spatial domains, respectively. Aided by the MTP, our approach can accurately but efficiently calculate the diffractive field from the nonuniformly distributed angular spectrum.
3. RESULTS
To establish the efficacy of our method, we compared it with a control (Ctrl) method that combined the FFT with angular spectrum interpolation [28,29]. In addition, we replaced the vanilla FFT in the Ctrl method with a generalized variant known as the chirp-Z transform for a fair comparison, which also overcomes the sampling constraint between the spatial and frequency domains [24,40,41]. The ground truth (GT) is derived from naive point-by-point integration [41]. The light wavelength utilized in our simulations is 785 nm. All simulations were executed by MATLAB on a personal computer equipped with an Intel i5 10400F CPU and evaluated with the timeit function.
A. Scalar Diffraction
We first demonstrated our method on scalar diffraction. We assumed the incident beam passed through a thin lens with a focal length of 225 mm, and calculated the diffractive field at the focus. The incident field was set to be slightly random to avoid any computational ambiguity and then superposed with a vortex phase to produce some noticeable features in the diffractive field [Fig. 2(a)]. The dimensions of the source and observation planes were mmmm6.39mm×6.39mm and mmmm0.25mm×0.25mm, respectively.
Fig. 2. Simulation results of scalar diffraction. (a) Amplitude and phase distribution on the source plane. (b) Calculated intensity and phase of the diffractive fields of GT (left column), our method (middle column), and the Ctrl method (right column). Note that a common linear phase term, induced by the intersection angle, is removed for visualization. (c) Field deviations, computed by the absolute value of the direct subtraction, between the two methods and GT. (d) Mean error under different incident fields of the two methods as a function of intersection angles. The error function is defined by Supplement 1, Eq. (S17). (e) Volumetric intensity distribution of the diffractive field, where the observation plane is colored in blue. (f) Time consumption as a function of errors. Note that the errors in (d) and (f) are shown in logarithmic values, since our method is far more accurate than the Ctrl one.
We investigated the diffractive field on the observation plane with a randomly selected angle pair (?,?)=(50°,30°) [Fig. 2(b)]. Both the source and observation planes were sampled with 256×256 points. The position of the observation plane and the volumetric intensity distribution obtained through the integration method are depicted in Fig. 2(e) as a reference. By adjusting the merging threshold, our method demonstrated an almost equivalent time consumption to the Ctrl method (about 0.06 s). Remarkably, as depicted in Fig. 2(c), the deviation of our method is approximately one-tenth that of the Ctrl method. In contrast to the nearly uniform distribution of the error induced by interpolation in the Ctrl method, the deviation in our method is characterized by the minimum at the center, escalating towards the periphery. This arises from the merging process within the angular spectrum rearrangement. For merging multiple ?? or ?? into a single ?? or ??, it can be interpreted that the multi-beam interference along approximate directions is equivalent to a plane wave. Consequently, within a confined field of view (FOV), the interference patterns can be well emulated by a uniform beam of light, leading to the best accuracy around the FOV center. Further elucidation of this phenomenon can be found in Supplement 1, Section 3.
We further quantified the computational accuracy at different angles. As illustrated in Fig. 2(d), our approach consistently achieves superior accuracy when maintaining similar computational efficiency of the two methods to that in Fig. 2(b). In particular, when ? and ? are set to 0°, 90°, or 180°, our method exhibits the maximum improvement, which is attributed to reduced merging. At these angles, a tremendous number of duplicated elements of ?? and ?? are generated, leading to a tighter angular spectrum matrix after projection. Consequently, with comparable computational time to other angles, this matrix closely approximates the GT with less merging, resulting in higher accuracy. In contrast, apart from these angles, the duplicated elements are much fewer, resulting in a slightly increased error with more merging, but still 102−104 times more accuracy than that of the Ctrl method. The maximum accuracy enhancement, about 109 times compared to the Ctrl method, is attained at (?,?)=(90°,0°) and (90°, 180°), with errors of 3.20×10−10 and 0.46 for our method and Ctrl, respectively. In contrast, the Ctrl method exhibits greater computational errors, especially when ?=90°, induced not only by interpolation but also by Jacobian determinant (Supplement 1, Section 1C).
Furthermore, our method enables a flexible balance between accuracy and efficiency. As demonstrated in Fig. 2(f), the time consumption of our method diminishes progressively as the error increases with two different angles. As the error increases from 0 to 0.01, the time consumption drops from 0.970 to 0.034 s for (?,?)=(50°,30°) and from 0.062 to 0.029 s for (?,?)=(90°,90°), respectively. This indicates that our method can achieve more accurate and efficient calculation results compared to the Ctrl method by trading accuracy for efficiency, and vice versa. Moreover, even in the slowest case, the time consumption of our method is less than one second, while nearly eliminating errors. In contrast, the Ctrl method fails to reduce errors at various angles due to interpolation, rendering further accuracy or efficiency improvement unattainable. We also provide results of other intersection angles and that of an additional representative source field with a Gaussian amplitude (Supplement 1, Section 4). Moreover, we present an example where the Ctrl method suffers from severe numerical instability due to the Jacobian determinant (Supplement 1, Section 4C). These results indicate that our method consistently outperforms the Ctrl method.
B. Vectorial Diffraction
We further demonstrated our method on vectorial diffraction, which is more effective in calculating the tightly focused field formed by a high-NA objective lens (e.g., NANA=1.35 and ?=1.40). Similarly, the incident field on the entrance pupil of the objective lens includes slightly random polarization, amplitude, and phase, to represent a general situation [Fig. 3(a)].
Fig. 3. Simulation results of vectorial diffraction. (a) Polarization, amplitude, and phase distribution on the pupil plane. The pseudo color for polarization visualization denotes ellipticity. (b) Calculated intensity and phase of the ?-, ?-, ?-polarized component, and the total intensity, of the diffractive fields of GT (top row), our method (middle row), and the Ctrl method (bottom row). (c) Volumetric intensity distribution of the diffractive field, where the observation plane is colored in blue. (d) Field deviations between the two methods and GT of the ?-, ?-, ?-polarized components and the total intensity. (e) Mean error under different pupils of the two methods as a function of intersection angles. (f) Time consumption as a function of errors.
The optical field on the observation plane was calculated at (?,?)=(130°,30°) [Fig. 3(b)], with the location in the volumetric intensity distribution depicted in Fig. 3(c). Both methods are calculated with approximately the same time consumption (about 0.02 s). The samples on the pupil and the observation planes (µmµm3.10µm×3.10µm) were 128×128 and 100×100, respectively. Similar to the scalar one, more accurate results can be obtained in the center of the FOV due to the equivalence of merging, as demonstrated in Fig. 3(d). Additionally, our method yields 10−104 times more accurate results for most cases, achieving a maximum accuracy improvement of 8.35×104 times at (?,?)=(36°,0°) and (144°, 180°), with errors of 1.52×10−4 and 12.72 for our approach and the Ctrl one, respectively [Fig. 3(e)]. Notably, the error of the Ctrl approach is unacceptable, deriving from the interpolation and the Jacobian (Supplement 1, Section 2B). Besides, it features a symmetric distribution as shown in Fig. 3(e), which arises from the fact that the two observation planes coincide when their ? values are complementary and ? are 0° and 180°, respectively.
For vectorial diffraction modeling, the balance between accuracy and efficiency is consistent with that in the scalar case [Fig. 3(f)]. As the error increases from 0 to 0.1, the time consumption decreases from 0.149 to 0.008 s for (?,?)=(130°,30°) and from 0.021 to 0.007 s for (?,?)=(90°,90°), respectively. Consequently, our method can enhance both accuracy and efficiency through adjustment of the merging threshold compared to the Ctrl one as in the scalar case. In addition, our method has the potential for acceleration by recording the merging process, which remains consistent when system parameters remain unchanged, such as the intersection angle, NA, refractive index, and wavelength. We also provide further quantification on other intersection angles and another representative demonstration, in which the incident field carries a left-handed circular polarization, a Gaussian amplitude, and a vortex phase (Supplement 1, Section 5). In addition, as mentioned earlier, the Ctrl method may suffer from severe numerical instability due to the Jacobian, while our method does not have this drawback. This fact is further supported by additional results (Supplement 1, Section 5C).
C. Application in CGH
After theoretical validation, we employ our approach for CGH, which is a representative application of diffraction theory. We targeted to realize phase-only Fourier holography on a tilted plane, as presented in Fig. 4(a). Herein, we simplified our analysis by neglecting the polarization effects. The target image, shaped like a heart, was set on the (?,?)=(35°,0°) observation plane with the size of µmµm7.68µm×7.68µm. The hologram was generated by the weighted Gerchberg–Saxton algorithm [42,43], where our method and the Ctrl one were employed for both forward and backward wave propagation modeling. The samples of the hologram and observation planes were 432×432 and 256×256 pixels, respectively.
Fig. 4. CGH on a tilted plane. (a) Simplified schematic of the optical setup. The black arrow represents the incident direction of the laser. The bottom right patterns are the holograms calculated from our method and the Ctrl one, respectively; L1–L2, lenses. (b) Generated patterns utilizing our method (left column) and the Ctrl one (right column) in the simulation. (c) Experimental results. (d) Volumetric maximum intensity projection, revealing the highest intensity structure within a volume, of diffractive field calculated from our method (left column) and the Ctrl one (right column) in the experiment.
We first validated our method through simulation, as illustrated in Fig. 4(b). The algorithms converged after 100 iterations. The intensity distribution generated by our method matches well with the target. However, that obtained from the Ctrl method exhibits poor contrast.
In the experimental validation, the two holograms were loaded on a phase-only spatial light modulator (SLM), which was conjugated with the entrance pupil of an objective lens (1.35 NA) using a 4? system. As shown in Fig. 4(d), we first acquired a volumetric intensity by laser scanning imaging of a single gold bead with 150 nm diameter [44]. The intensity distribution on the tilted plane was then obtained by linear interpolation. Our method yields a uniformly shaped heart distribution [Fig. 4(c)], agreeing well with the simulation. In contrast, the result from the Ctrl method remarkably deviates from the target, and even fails to generate the outline.
4. DISCUSSION AND CONCLUSION
Angular spectrum rearrangement relies on merging the angular spectrum at similar ?? and ?? coordinates. The error introduced by merging is influenced by the angular spectrum distribution in the original sparse (??,??) coordinates and their sparsity. Moreover, the distribution of the angular spectrum is dependent on the source field, while the sparsity is influenced by the intersection angles ? and ? (Supplement 1, Section 6A). Consequently, it is challenging to determine an optimal threshold range that consistently yields better results in all cases. Nevertheless, as previously mentioned, we can set the threshold by ensuring the number of samples after merging matches that of the source field. This strategy generally ensures accurate results (Supplement 1, Section 6B), although slight fluctuations may occur in different cases. In addition, the balance between accuracy and efficiency can be further adjusted by tuning the threshold. A larger threshold can result in a more compact angular spectrum matrix, enabling faster calculations at the expense of accuracy, and vice versa.
Besides the merging threshold, we also discuss the determination of the effective frequency coordinates after merging. Although we previously adopted the weighted average of the angular spectrum’s amplitude, alternative weighting methods are applicable and may exhibit different performances. Therefore, we compared our weighting method with two other methods where the weights are evaluated by intensity and median, respectively. The results indicate that our weighting method and the others introduce similar errors (Supplement 1, Section 7), with ours showing slightly lower errors. This is reasonable because the merging operation, as mentioned in the Methods section, can be regarded as approximating multiple-beam interference with one plane wave, where the interference pattern can be dominated by the beam with the maximum amplitude.
In conclusion, our study demonstrates accurate and efficient diffraction modeling between arbitrary non-parallel planes, applicable to both scalar and vectorial diffraction models. This is accomplished by compressing the sparse angular spectrum matrix in a physics-informed manner coupled with the MTP algorithm. We calculated the diffractive field with different intersection angles, consistently proving that our approach leads to higher accuracy than that of existing methods. Importantly, our method enables a flexible balance between accuracy and efficiency by adjusting the merging threshold. We can decrease/increase the threshold to further improve the accuracy/efficiency at the cost of each other. Based on this insight, our method can be much faster than the state-of-the-art ones, with acceptable accuracy. Furthermore, with the potential for leveraging advanced parallel computing technologies, our method provides access to further acceleration. Besides numerical verifications, we experimentally demonstrated our method in CGH on a randomly selected tilted plane. In contrast to existing methods, our approach can well generate the target and agrees with the simulation. In addition, we also note that our method can serve as a cornerstone of the rising learning-based approaches to further improve the computational accuracy and efficiency. This work represents an important step towards high-performance computational optics.
