1. Bibliographic Information

• Title: Modeling off-axis diffraction with the least-sampling angular spectrum method
• Authors: Haoyu Wei, Xin Liu, Xiang Hao, Edmund Y. Lam, and Yifan Peng.
  • Affiliations: The University of Hong Kong, Hong Kong SAR, China; Zhejiang University, Hangzhou, China. The authors have expertise in electrical engineering, optical science, and computational imaging.
• Journal/Conference: Optica. This is a highly prestigious, peer-reviewed journal from the Optica Publishing Group (formerly OSA). It is known for publishing high-impact, significant advances in the field of optics and photonics.
• Publication Year: 2023 (Published July 19, 2023).
• Abstract: The paper addresses the challenge of accurately and efficiently simulating off-axis diffraction, a critical task for designing advanced computational optics. Existing methods are hampered by rigid sampling schemes. The authors propose a universal least-sampling angular spectrum method (LS-ASM). The core ideas are: (1) using the Fourier transform's shifting property to convert off-axis problems into more efficient quasi-on-axis ones, and (2) linking the angular spectrum to the transfer function to adaptively determine the minimum necessary sampling rates. Implemented with a flexible matrix-based Fourier transform, their method is demonstrated on coded-aperture imaging systems. The key findings include a remarkable ~36× speedup over the state-of-the-art at a 20° incident angle and the ability to compute ultra-large angles like 35° in seconds on a standard computer.
• Original Source Link: The paper is available at /files/papers/68e7b3ff066bd688e1396b16/paper.pdf and was published under the Optica Open Access Publishing Agreement.

2. Executive Summary

• Background & Motivation (Why):

  • Core Problem: Simulating how light diffracts when it comes from an off-axis source (i.e., at a large angle) is computationally very expensive. Traditional methods require an extremely large number of sampling points to avoid errors (aliasing), making the design and optimization of large-scale computational optics (like advanced lenses in VR/AR or cameras) slow and impractical.
  • Gaps in Prior Work: Previous methods, even advanced ones, determine the sampling requirements for the input light field and the propagation physics (transfer function) independently. This is suboptimal because it doesn't consider how these two components interact, leading to massive oversampling (wasting resources) or undersampling (producing incorrect results). They also struggle with the inherent asymmetry of off-axis problems.
  • Innovation: This paper introduces a new way to think about the sampling problem. Instead of treating the components separately, it analyzes their product to find the true, minimal sampling rate needed. It also cleverly simplifies the problem by mathematically "rotating" the off-axis wave to look like an on-axis one, which is much easier to simulate.

• Main Contributions / Findings (What):

  • A Novel Algorithm (LS-ASM): The authors developed the Least-Sampling Angular Spectrum Method (LS-ASM), an adaptive and efficient algorithm for modeling large-angle off-axis diffraction.
  • Optimal Sampling Strategy: The core contribution is a universal method to determine the absolute minimum sampling requirements for both the spatial (input plane) and frequency (propagation) domains. This is achieved by:
    1. Linear Phase Compensation (LPC): Converting the off-axis problem to a quasi-on-axis one, drastically reducing the required samples.
    2. Combined Frequency Analysis: Analyzing the product of the light's angular spectrum and the propagation transfer function to find the true minimum sampling rate, rather than wastefully sampling each part at its maximum rate.
  • Demonstrated Performance: The paper shows that LS-ASM achieves a speedup of approximately 36 times over the previous state-of-the-art method (shift-BEASM) for a 20° angle while maintaining high accuracy. It also demonstrates the feasibility of simulating extremely large angles (35°) within seconds, a task previously considered computationally prohibitive on standard hardware.

3. Prerequisite Knowledge & Related Work

• Foundational Concepts:

  • Diffraction: The phenomenon where light waves bend as they pass around an obstacle or through an aperture (an opening). It's a fundamental concept in wave optics that explains why light doesn't always travel in perfectly straight lines.
  • Angular Spectrum Method (ASM): A powerful technique for simulating wave propagation. It works by decomposing a complex light field into a sum of simple plane waves, each traveling in a different direction. The method then mathematically propagates each plane wave and sums them back up at the destination plane to get the final diffracted field. This process is efficiently calculated using the Fast Fourier Transform (FFT).
  • Fourier Transform: A mathematical tool that decomposes a signal (like a light field) into its constituent frequencies. In ASM, the Fourier transform of the light field at the aperture plane gives its angular spectrum (the collection of plane waves).
  • Nyquist Sampling Theorem: A fundamental principle in signal processing stating that to accurately represent a signal, one must sample it at a rate at least twice its highest frequency component. In diffraction simulation, violating this theorem leads to aliasing, an artifact where high-frequency information incorrectly appears as low-frequency information, corrupting the result.
  • Transfer Function (in ASM): A mathematical function that describes how each plane wave in the angular spectrum changes (specifically, its phase shifts) as it propagates from one plane to another.
  • Point Spread Function (PSF): The diffraction pattern produced by an optical system when viewing a single point source of light. It characterizes the system's imaging quality; a smaller, sharper PSF is better.

• Previous Works:

  • Conventional FFT-based ASM: The standard approach. Its major drawback is that it requires significant zero-padding (adding a large border of zeros around the data) to prevent aliasing, especially for off-axis or long-distance propagation. This drastically increases memory usage and computation time.
  • Shifted Band-Limited ASM & Generalized Transfer Function Method: These were early attempts to improve efficiency by "shifting" the calculation window to better align with the off-axis energy, but their sampling requirements remained rigid and often excessive.
  • Shifted Band-Extended ASM (shift-BEASM): The state-of-the-art method before this paper. It improved upon previous methods but still suffered from a key limitation: it determined the sampling rates for the input field and the transfer function independently. This leads to oversampling because it doesn't account for the possibility that the combined signal might have a lower overall frequency content.

• Differentiation: The proposed LS-ASM is fundamentally different from shift-BEASM and others in two ways:

  1. Quasi-on-axis Conversion: LS-ASM uses Linear Phase Compensation (LPC) to mathematically center the off-axis wave's spectrum around the origin (zero frequency). Previous methods work with the off-center spectrum, forcing them to use a much larger sampling window to capture it without aliasing.
  2. Combined Sampling Analysis: LS-ASM's key insight is to analyze the product of the input field's angular spectrum and the propagation transfer function. This allows it to find the true, minimal ("least") sampling rate required for the combined signal, whereas previous methods conservatively over-sample each component, leading to wasted computation.

4. Methodology (Core Technology & Implementation)

The core of LS-ASM is a multi-step process to identify and apply the minimum possible sampling rates in both the spatial and frequency domains.

• Principles: The method is built on the intuition that the sampling rate should be determined by the actual frequency content of the signal being computed at each step, not by a rigid, worst-case scenario. By manipulating the problem mathematically, the maximum frequency content can be minimized, thus minimizing the required samples.

• Steps & Procedures:

  1. Converting Off-axis to Quasi-on-axis via Linear Phase Compensation (LPC)

  The paper starts with the standard ASM formula for calculating the diffracted field $u_z(x, y)$ at a distance $z$: $$u _ { z } ( x , y ) = \mathcal { F } ^ { - 1 } [ U ( f _ { X } , f _ { Y } ) H _ { z } ( f _ { X } , f _ { Y } ) ]$$

  • $u_z(x, y)$: The complex light field at the observation plane.

  • $\mathcal{F}^{-1}$: The inverse Fourier transform.

  • $U(f_X, f_Y)$: The angular spectrum of the input field, which is the Fourier transform $\mathcal{F}$ of the field right after the aperture, $u(\xi, \eta)$.

  • $H_z(f_X, f_Y)$: The ASM transfer function.

    An off-axis wave has a phase that includes a large linear component, which shifts its angular spectrum $U$ away from the frequency origin $(0, 0)$. This is shown in Figure 1(b) (left panel, "without LPC"). According to the Nyquist theorem, a larger frequency range requires a higher sampling rate.

  To solve this, LS-ASM applies a compensating linear phase $\phi_{\mathrm{comp}}$ to the input field's phase $\phi_u$: $$\phi _ { \mathrm { on - a x i s } } = \phi _ { u } + \phi _ { \mathrm { c o m p } }$$

  • $\phi_{\mathrm{on-axis}}$: The new phase of the "quasi-on-axis" wave.

  • $\phi_u$: The original phase of the input field.

  • $\phi_{\mathrm{comp}}(\xi, \eta) = -2\pi(f_{Xc}\xi + f_{Yc}\eta)$: The compensating phase, which is a linear ramp perfectly opposing the inherent linear phase of the off-axis wave. $(f_{Xc}, f_{Yc})$ is the center of the original angular spectrum.

    This operation, leveraging the shifting property of the Fourier transform, effectively moves the angular spectrum back to the center of the frequency domain, as shown in Figure 1(b) (right panel, "with LPC"). The resulting "quasi-on-axis" wave has a much lower bandwidth and thus requires far fewer samples.

    

    2. Determining the Spatial Sampling Rate

  Even for an on-axis wave, its spectrum is broadened by the finite aperture. To find the effective bandwidth and thus the required sampling rate, the paper uses a "virtual thin lens" analogy, shown in Figure 2(a). This optical insight helps determine the extent of the spectral broadening. The final sampling rate for the spatial domain $(\xi, \eta)$ is given by: $$\tau _ { \mathrm { s p a t } } ^ { ( p ) } = s f _ { \mathrm { b } } ^ { ( p ) } = s \left( \frac { 1 } { \pi } \bigg | \frac { \partial } { \partial \boldsymbol { \rho } } \phi _ { \mathrm { o n - a x i s } } \bigg | _ { \operatorname* { m a x } } + \frac { 4 1 . 2 } { D } \right)$$

  • $\tau_{\mathrm{spat}}^{(p)}$: The required number of samples per unit length in the spatial domain.

  • $s \geq 1$: A heuristic oversampling factor (e.g., 1.2 or 1.5) to add a safety margin and reduce aliasing artifacts, as illustrated in Figure 2(b).

  • $|\frac{\partial}{\partial \boldsymbol{\rho}} \phi_{\mathrm{on-axis}}|_{\max}$: The maximum rate of change (gradient) of the quasi-on-axis phase, which corresponds to its highest spatial frequency.

  • $41.2/D$: A term accounting for the spectral broadening caused by the aperture of diameter $D$.

    

    3. Determining the Frequency Sampling Rate (Combined Analysis)

  This is the paper's most significant technical novelty. Instead of sampling the angular spectrum $U$ and the transfer function $\hat{H}_z$ independently at their respective maximum rates, LS-ASM analyzes the phase of their product, $U\hat{H}_z$.

  As shown in Figure 3, the angular spectrum $U$ (after LPC) can be viewed as a diverging spherical wave, and the transfer function $\hat{H}_z$ can be viewed as a converging spherical wave. When multiplied, their quadratic phase terms can partially cancel each other out. The phase of the product (rightmost panel) is much smoother and varies more slowly than the phases of the individual components (left and center panels). A slower-varying phase means a lower local frequency, which in turn requires a lower sampling rate.

  该图像为一组图表，展示了不同函数图像及其局部放大细节。图中从左到右依次为函数φ_u(f_x,f_y)、函数φ_Ĥ(f_x',f_y')及两者的叠加φ_u(f_x,f_y)+φ_Ĥ(f_x',f_y')，右侧配有色条，表示其取值范围为0到2π，反映了相位分布的变化特征。图中局部放大框突出展示了细节差异。

  The final sampling rate in the frequency domain is adaptively determined by the following formula, which considers this combined effect as well as the desired size of the observation window (w): $$\tau _ { \mathrm { f r e q } } ^ { ( f ) } = s \ \operatorname* { m a x } \left\{ \operatorname* { m i n } \left\{ 2 D + \tau _ { \gamma \hat { H } _ { z } } ^ { ( f ) } , D + \tau _ { \hat { H } _ { z } } ^ { ( f ) } \right\} + \frac { 4 1 . 2 } { \tau _ { \mathrm { s p a t } } ^ { ( p ) } } , 2 | q _ { c } ^ { \prime } | + w \right\}$$

  • $\tau_{\mathrm{freq}}^{(f)}$: The required number of samples per unit frequency.
  • $\tau_{\gamma \hat{H}_z}^{(f)}$: The sampling rate determined by the phase of the product $U\hat{H}_z$.
  • $\tau_{\hat{H}_z}^{(f)}$: The sampling rate determined by the phase of $\hat{H}_z$ alone.
  • The formula intelligently chooses the minimum necessary rate while also ensuring the final observation window of size $w$ is adequately sampled.

  4. Implementation with Matrix-based Fourier Transform

  Conventional FFT algorithms are rigid; they require uniform sampling and the number of samples must often be a power of two. Since LS-ASM calculates arbitrary, adaptive sampling rates, it requires a more flexible Fourier transform implementation. The authors use the Matrix Triple Product (MTP), which represents the Discrete Fourier Transform as a matrix multiplication. While potentially slower for uniform grids, MTP is perfectly suited for arbitrary sampling grids and can be highly accelerated on modern hardware like GPUs.

5. Experimental Setup

• Datasets: The study is based on numerical simulations, not experimental datasets. The primary task is to compute the off-axis complex point spread function (PSF) for a thin lens camera system with the following parameters:

  • Wavelength ($\lambda$): 500 nm
  • Focal length: 35 mm
  • f-number: 16
  • Source distance ($z_0$): 1.7 m

• Evaluation Metrics:

  • Signal-to-Noise Ratio (SNR):
    1. Conceptual Definition: SNR measures the quality of the computed result by comparing the power of the ground truth signal to the power of the error (the difference between the computed result and the ground truth). A higher SNR value (in decibels, dB) indicates a more accurate result.
    2. Mathematical Formula: $$\mathrm{SNR} = 10 \log_{10} \left( \frac{\sum |u_{gt}|^2}{\sum |u_{computed} - u_{gt}|^2} \right)$$
    3. Symbol Explanation:
       • $u_{gt}$: The ground truth complex field, obtained from the highly accurate but slow Rayleigh-Sommerfeld (RS) integral.
       • $u_{computed}$: The complex field computed by the method being evaluated (e.g., LS-ASM or shift-BEASM).
       • $\sum |\cdot|^2$: The sum of squared magnitudes (power) over all pixels.
  • Root-Mean-Square Error (RMSE):
    1. Conceptual Definition: RMSE measures the average magnitude of the error between the computed values and the ground truth values. A lower RMSE indicates a more accurate result.
    2. Mathematical Formula: $$\mathrm{RMSE} = \sqrt{\frac{1}{N} \sum_{i=1}^{N} |u_{computed, i} - u_{gt, i}|^2}$$
    3. Symbol Explanation:
       • $N$: The total number of pixels.
       • $u_{computed, i}$ and $u_{gt, i}$: The complex values of the $i$-th pixel for the computed and ground truth fields, respectively.

• Baselines:

  • Shift-BEASM: The primary baseline, representing the state-of-the-art before this work. The authors test it in two scenarios:
    • Case 1: shift-BEASM is given enough samples to achieve the same high accuracy (SNR) as LS-ASM. This allows for a fair comparison of computational time.
    • Case 2: shift-BEASM is forced to use the same (low) number of samples as LS-ASM. This is done to show that shift-BEASM fails (produces aliasing) under such efficient sampling conditions.
  • Rayleigh-Sommerfeld (RS) Integral: Not a baseline for speed, but used to generate the ground truth for accuracy comparisons. It is considered one of the most accurate methods for scalar diffraction but is extremely slow.

6. Results & Analysis

The results compellingly demonstrate the superiority of LS-ASM in both accuracy and efficiency.

该图像包含三部分：(a)为三种方法（Shift-BEASM、LS-ASM和RS）在两种案例下的衍射场幅度、相位和频谱分布对比；(b)为不同入射角θ下信噪比(SNR)和计算时间的曲线图，展示Shift-BEASM在20°时显著提升速度；(c)为空间采样点数和频率采样点数随入射角变化的折线图，显示LS-ASM和Shift-BEASM采样效率差异。整体体现了提出方法在高角度衍射模拟中的高效性和准确性。

• Core Results (Analysis of Figure 4):

  • Visual Quality (Figure 4a):
    • The amplitude and phase of the PSF computed by LS-ASM are visually identical to the RS ground truth.
    • Shift-BEASM (Case 1) is also accurate but, as shown later, is much slower. Its angular spectrum is visibly off-center.
    • Shift-BEASM (Case 2), using the same low sample count as LS-ASM, fails catastrophically. The angular spectrum is severely cropped (aliased), leading to a completely incorrect PSF. This proves that LS-ASM's sampling strategy is not just efficient but essential for correctness at low sample counts.
  • Quantitative Performance (Figure 4b):
    • Time: The green dashed lines show computation time. LS-ASM's runtime remains flat and extremely low (under 0.2 seconds) even as the incident angle increases to 20°. In contrast, the time for shift-BEASM (Case 1) explodes, reaching over 4 seconds at 20°. At 20°, LS-ASM is ~36 times faster than shift-BEASM for the same level of accuracy.
    • Accuracy (SNR): The blue lines show SNR. LS-ASM and shift-BEASM (Case 1) both maintain a high, stable SNR (~50 dB), proving their accuracy. shift-BEASM (Case 2) has a very low SNR that drops further with angle, confirming its failure.
  • Sampling Efficiency (Figure 4c):
    • This plot explains why LS-ASM is so much faster. The number of spatial samples ($N_s$) for shift-BEASM (green line) increases linearly and dramatically with the incident angle. For LS-ASM (blue line), $N_s$ remains constant and low because the LPC step removes the angle-dependent sampling burden. This is the key to its efficiency. The number of frequency samples ($N_f$) also shows an advantage for LS-ASM.

• Additional Findings:

  • Ultra-large Angle: The method successfully computes diffraction at an extreme angle of 35° in just 1.2 seconds, a feat that would be impractical with shift-BEASM.
  • Complex Fields: The paper also shows LS-ASM works correctly for more complex inputs, such as when a cubic phase plate is added to the lens, demonstrating its versatility.

7. Conclusion & Reflections

• Conclusion Summary: The paper successfully introduces and validates the Least-Sampling Angular Spectrum Method (LS-ASM), a novel and highly efficient algorithm for simulating off-axis diffraction. By converting the problem to a quasi-on-axis one and performing a combined analysis of the field and transfer function, LS-ASM establishes an optimal, adaptive sampling strategy. This leads to massive improvements in computational speed and memory usage without sacrificing accuracy, enabling the simulation of large-angle diffraction phenomena that were previously computationally prohibitive.

• Limitations & Future Work:

  • Non-smooth Fields: The method's core assumptions are a slowly varying amplitude and smooth phase. For fields with sharp discontinuities (e.g., a binary diffractive optical element), the authors suggest the method can be applied to local sub-apertures where the field is smooth.
  • Scalar Diffraction Limit: Like ASM in general, the method is based on scalar diffraction theory, which may become less accurate at extremely large angles or for features smaller than the wavelength.
  • Long Propagation Distance: Very long propagation can still increase the sampling burden. The authors suggest this can be mitigated by combining their method with semi-analytical Fourier transform techniques that handle quadratic phase components analytically.

• Personal Insights & Critique:

  • Impact: This work is a significant practical breakthrough for the field of computational optics. The ability to quickly and accurately model off-axis performance is a bottleneck in designing next-generation optical systems like metalenses, compact cameras, and AR/VR headsets, all of which require a wide field of view. LS-ASM directly addresses this bottleneck.
  • Methodological Elegance: The approach is both clever and rigorous. The insight to analyze the product $U\hat{H}_z$ rather than the components separately is a fundamental shift in perspective that yields enormous practical benefits. The use of the Fourier shift theorem via LPC is an elegant way to simplify a complex problem.
  • Implementation Choice: The choice of the matrix-based Fourier transform (MTP) is crucial for practicality. It unshackles the algorithm from the rigid constraints of FFT, allowing the theoretically-derived "least" sampling rates to be directly implemented. This highlights the synergy between theoretical algorithm design and hardware-aware implementation.
  • Open Questions: While highly effective, it would be interesting to see a more detailed analysis of the trade-off with the oversampling factor $s$. The paper finds $s=1.5$ to be a good balance, but a more formal guideline for choosing this hyperparameter for different types of fields could be valuable for practitioners. The application to full-color simulations (requiring runs at multiple wavelengths) would also be a compelling demonstration.