Scalable angular spectrum propagation

Rainer Heintzmann; Rainer Heintzmann; Lars Loetgering; Lars Loetgering; Lars Loetgering; Lars Loetgering; Felix Wechsler; Felix Wechsler

doi:10.1364/OPTICA.497809

1. INTRODUCTION

Accurate and performant wave propagation algorithms are a key ingredient to various active fields of research. Diverse applications, ranging from digital holography, ptychography, wavefront sensing, and synthetic aperture radar to scatterometry and crystallography [1–6], employ lens-less experimental detection geometries. In such configurations, the measured signal is typically a diffraction pattern (or a sequence of such). Phase-sensitive detection schemes [7] or phase retrieval algorithms [8–10] in conjunction with suitable wave propagation models [11–14] are commonly used to refocus the reconstructed field yielding an image of the specimen of interest or the wavefront at a plane of choice. On the other side of the application spectrum are lens-based technologies. Although in such systems a specimen can directly be imaged, a growing number of applications use data-driven techniques to enhance the recorded data. Recent examples include Fourier ptychography, optical diffraction tomography, structured illumination microscopy, and quantitative phase imaging [15–20]. Thus, a common theme among modern computational imaging (CI) architectures, both lens-less and lens-based, is the need to formulate the underlying forward model of the system at hand to extract information from the measured data. This implies the need for both accurate and performant wave optical modeling, as CI systems often use iterative reconstruction algorithms where the propagation step has to be repeated many times.

A suitable starting point for numerical wave propagation is the so-called angular spectrum (AS) method [11]. In its simplest form, the AS enables computation of a diffracted wave field between two parallel planes by a three-step process [Fig. 1(a)]: (1) Fourier transform the field given in the source plane, (2) multiply by a free-space transfer function, and (3) inverse Fourier transform to obtain the field distribution in the destination plane. Thus the AS can be computed at a cost proportional to two fast Fourier transform operations. Importantly, the AS is an exact solution to the Helmholtz equation [21], which governs wave propagation of an electric field component in homogeneous media. It thus outperforms approximate models such as Fresnel and Fraunhofer diffraction in terms of numerical precision. However, while computationally simple and precise, the AS approach faces several complications. A drawback is the fixed pixel pitch of the angular spectrum propagator [Fig. 1(a)]. In practical scenarios, the measured intensity is sampled at the native pixel pitch of the detector, which is often much coarser than the diffraction-limited resolution with respect to the destination amplitude distribution set by the detection NA in a given measurement geometry. Naïve application of the AS method can thus inhibit the exploitation of the significantly reduced sampling requirements of a given imaging architecture. To propagate rays up to the sampling limit, the AS method usually needs extensive zero-padding, which slows down the processing enourmously.

Fig. 1. Operator diagrams propagating a source field of source pixel pitch ${\Delta _s}$ over a distance $z$ to a destination field of pixel pitch ${\Delta _d}$. (a) Standard AS propagation and (b) proposed SAS propagation algorithm jointly applying a pre-compensation (${{ H}_{\rm{AS}}}{H}_{\rm{Fr}}^*$) term followed by a single-step Fresnel propagation. [i] FFT denotes the [inverse] fast Fourier transform, $\lambda$ the wavelength, and $N$ the number of pixels per dimension in source or destination plane. The black, red, and blue colors distinguish among the source, Fourier, and destination planes, respectively. See Section 2 for the definitions of the transfer functions ($H$) and Fresnel single-step quadratic phase exponentials ${{Q}_{1,2}}$.

Download Full Size | PDF

However, in the non-paraxial regime, a large variety of field propagation methods directly based on the AS have been proposed that allow for pixel rescaling or truncation of the field of view in the destination plane. These methods typically make use of three classes of numerical algorithms: the nonuniform fast Fourier transform (NUFFT), the matrix triple product (MTP), and the chirp Z transformation (CZT). Shimobaba et al. used the NUFFT to allow for scaling of the destination pixel pitch by replacing the inverse Fourier transform step in the AS propagation algorithm with a NUFFT [22]. Zhang et al. proposed the so-called band-extended AS method, which used two successive NUFFTs to relax sampling requirements on the free-space transfer function and to enable scalable pixel pitch in the destination plane [23]. Another method to calculate off-axis diffraction is termed shifted AS, where the computational cost can be reduced by considering only a finite observation window away from the optical axis by means of the CZT [24]. Wei et al. proposed a wave propagation algorithm using the MTP for computing off-axis diffraction [25]. Yu et al. presented the windowed AS, a propagation method based on two successive CZTs, which allows for pixel rescaling in the destination plane [26]. If the convolution operations involved are preallocated in form of transfer functions, the latter method comes at the computational cost of four FFTs. High numerical aperture vectorial focusing methods with zoom-in capability have been proposed, which, similar to off-axis diffraction, utilize the truncation ability of the CZT [27]. Hillenbrand et al. [28,29] presented a tile-based scheme that is combined with the CZT. Both CZTs and MTPs are useful when the destination field is truncated, as is often the case for the computation of point spread functions (PSFs), focused fields, and off-axis diffraction. However, when the destination field has a similar number of pixels as the source field, CZT- and MTP-based methods have a higher computational complexity as compared to FFTs [30]. While the FFT and NUFFT theoretically have the same computational complexity, the latter requires an additional convolution step and several exponential evaluations, which gives rise to higher prefactors in the computational cost as compared to the standard FFT [31]. A variety of approximate methods have been proposed. The Fresnel scaling theorem (FST) [32] states that a diffraction pattern collected under divergent illumination is equivalent to a diffraction pattern collected under plane wave illumination at shorter propagation distance and smaller effective pixel pitch—an effective tool for pixel downsizing. In addition, scalable two-step Fresnel propagators [33–35] and Fresnel propagators with zoom-in capability based on the CZT [36] have been proposed that allow for variable scaling of the pixel pitch between source and destination planes. Unfortunately, the validity of these propagation models is limited to the paraxial diffraction regime.

Another practical difficulty in the numerical evaluation of the AS is the possibility of aliasing. For large propagation distances, it is increasingly challenging to appropriately sample the transfer function underlying the AS. Sampling problems arise by waves propagating beyond the limit of the (circularly connected) computational array. This statement is equivalent to an undersampling of the associated phase modification of the propagator. A sampled representation of such wave fields is restricted to finite pixel pitch and field of view, as ultimately limited by the available computing memory. Thus, because the standard AS method preserves the pixel pitch in the source and destination planes, wave fields propagating over large distances eventually surpass the spatial boundaries of the computation grid, thereby causing aliasing artifacts. This problem is particularly severe for high spatial frequencies in a wave field, which propagate under large angles with respect to the optical axis. Matsushima and Shimobaba [37] proposed a bandlimited angular spectrum method that mitigates this problem by suppressing those spatial frequencies in the Fourier domain that surpass the lateral limits of the sampling domain. While this technique reduces aliasing issues, it tends to shrink the angular spectrum over large propagation distances while oversampling the field in the destination plane, resulting in the need for a more efficient sampling scheme. Moreover, the method is not flexible in terms of pixel scaling. Asoubar et al. [38] proposed a scalable angular spectrum propagation method requiring three FFT operations, where quadratic phase factors in both the free-space transfer function and the field to be propagated are analytically treated by means of the so-called semi-analytic Fourier transform. However, the authors did not derive bounds on the propagation distances covered by the method nor derived padding schemes for the convolutions involved to be guaranteed to avoid wrap-around artifacts. Since they apparently determine the required amount of zero-padding directly from the sampling requirements, comparably large arrays need to be Fourier-transformed. In addition, the method requires a fitting procedure that estimates quadratic phases present in the signal to be propagated (see also follow-up work [39]).

Here we propose a scalable angular spectrum (SAS) propagation algorithm, based on three uniformly spaced FFTs [Fig. 1(b), details presented below]. This method features numerical precision and applicability at high numerical apertures. The SAS method generalizes the standard AS method based on two FFTs by allowing for pixel scaling at a computational cost proportional to only one additional FFT as compared to the standard AS propagation method. We derive sampling conditions for the validity of SAS propagation, including an effective cutoff spatial frequency region, from which practical magnification limits arise. We show that the SAS method outperforms approximate models and the numerically padded AS in terms of numerical accuracy at high NA and computational complexity, respectively. We thus pave a precise, yet efficient path towards performant, high NA wave propagation.

2. THEORY

A. AS and Fresnel Diffraction

Free-space wave propagation methods aim at solving the homogeneous Helmholtz equation

(1)$${\nabla ^2}\psi + {k^2}\psi = 0,$$

subject to boundary conditions, either in an exact or in an approximate sense. In Eq. (1), $\psi (x,y,z)$ is a scalar component of the electromagnetic field at location $({x,y,z})$ in three-dimensional space. If the underlying boundary conditions of the problem provide knowledge of the field in a two-dimensional plane, it is straightforward to compute the field at a parallel plane separated by a propagation distance $z$ by means of AS propagation. Mathematically, AS propagation can be formulated as

(2)$${\psi _z} = {{\cal F}^{- 1}}\left[{{H_{{\rm AS}}}{\cal F}[{{\psi _0}} ]} \right],$$

where ${\psi _z}$ and ${\psi _0}$ are shorthand notation for $\psi (x,y,z)$ and $\psi (x,y,0)$, respectively, and ${\cal F}$ denotes two-dimensional Fourier transformation mapping from spatial coordinates $({x,y})$ to spatial frequency coordinates $({{f_x},{f_y}}) = ({\frac{{{k_x}}}{{2\pi}},\frac{{{k_y}}}{{2\pi}}})$. The AS transfer function is given by [11,21]

(3)$${H_{{\rm AS}}}({{f_x},{f_y}} ) = \exp \!\left({i2\pi z\sqrt {\frac{1}{{{\lambda ^2}}} - f_x^2 - f_y^2}} \right),$$

where $\lambda$ denotes wavelength and $z$ the propagation distance along the optical axis. We obtain an approximate form of the AS transfer function by Taylor-expanding the latter into

(4)$${H_{{\rm Fr}}}({{f_x},{f_y}} ) \approx \exp \!\left({i2\pi z\left[{\frac{1}{\lambda} - \frac{{\lambda f_x^2}}{2} - \frac{{\lambda f_y^2}}{2}} \right]} \right),$$

which we refer to as the Fresnel transfer function. The Fresnel diffraction formalism can be rewritten in real space ([11,21]):

(5)$$\begin{split}{\psi _z} &\approx {{\cal F}^{- 1}}\big[{{\cal F}[{{\psi _0}} ]{H_{{\rm Fr}}}} \big]\\ & = \frac{{\exp ({ikz} )}}{{i\lambda z}}\int \psi ({{x^\prime},{y^\prime},0} ) \\&\quad\cdot \exp \!\left[{i\frac{\pi}{{\lambda z}}\left({{{({x - {x^\prime}} )}^2} + {{({y - {y^\prime}} )}^2}} \right)} \right] {\rm d}{x^\prime} {\rm d}{y^\prime}\;\;\;{\rm (CV \text{-} FR)}\\ & = {Q_2}({x,y} )\int {Q_1}({{x^\prime},{y^\prime}} )\psi ({{x^\prime},{y^\prime},0} )\\&\quad\cdot\exp \!\left[{- i\frac{{2\pi}}{{\lambda z}}({x{x^\prime} + y{y^\prime}} )} \right] {\rm d}{x^\prime} {\rm d}{y^\prime}\;\;\;{\rm (SFT \text{-} FR)},\end{split}$$

where the quadratic phase exponentials are defined as

(6)$${Q_1}\big({{x^\prime},{y^\prime}} \big) = \exp \!\left[{i\frac{k}{{2z}}\big({{x^{\prime 2}} + {y^{\prime 2}}} \big)} \right]$$

and

(7)$${Q_2}({x,y} ) = \frac{{\exp ({ikz} )}}{{i\lambda z}}\exp \!\left[{i\frac{k}{{2z}}({{x^2} + {y^2}} )} \right]$$

in source (primed) and destination (unprimed) coordinates, respectively. Here we refer to the second line in Eq. (5) as the convolution-based Fresnel integral (CV-FR), while the third line is referred as the single-Fourier-transform Fresnel integral (SFT-FR). Formally, the CV-FR and SFT-FR methods are equivalent. However, numerically, they are different: the CV-FR method practically involves two Fourier transforms (assuming ${H_{{\rm Fr}}}$ is analytically computed in the Fourier domain), while the SFT-FR method only requires a single Fourier transformation leading to scaled destination coordinates. While for the AS and CV-FR methods, assuming that no cropping or padding in Fourier space is performed, we have equal destination and source pixel pitches,

(8)$${\Delta _d} = {\Delta _s},$$

the SFT-FR method rescales the destination plane pixel pitch to

(9)$${\Delta _d} = \frac{{\lambda z}}{{N \cdot {\Delta _s}}},$$

where for simplicity we assume a square input grid, with isotropic pixel pitch and $N$ samples per dimension [12,13]. Extensions to non-isotropic pixel pitches and number per samples are obtained by applying Eq. (9) to each lateral spatial dimension independently. In Eq. (9), we can alternatively solve for the source pixel pitch if the detector pixel pitch is given. It is this scaling property of the SFT-FR method that is desirable in many practical computational imaging systems. Even if the ${Q_2}$ phase exponential is massively undersampled at the destination pixel pitch, this often does not matter as long as the intensity is the measured quantity. ${Q_2}$ can then entirely be neglected. In fact, ${Q_2}$ is always undersampled for a magnification $M = {\Delta _d}/{\Delta _s} \gt 1$. For clarity, Table 1 provides an overview of the methods referred to throughout this work, including abbreviations, the computational complexity, and the respective pixel pitch.

Table 1. Overview of Commonly Used Scalar Wave Propagation Methods as Compared to the Proposed SAS Method^a

View Table | View all tables in this article

B. Pre-compensation and the Scalable Angular Spectrum Algorithm

The key idea of our proposed method of wave propagation is to pre-compensate for the approximation error that each spatial frequency of the source field experiences upon Fresnel propagation. Single-step Fresnel propagation scales its pixel size upon propagation, a feature that we wish to retain here. Yet, we want to achieve the numerical precision of the more accurate angular spectrum (AS) propagator, which commonly keeps the pixel size fixed.

The action of both angular spectrum and Fresnel propagation is conveniently seen in Fourier space. Each spatial frequency emanating from the object undergoes a different phase modulation upon free-space propagation. The propagators are described by $\exp (iz{k_z})$ with the propagation distance $z$ and the axial spatial frequency ${k_z}$. For AS propagation, ${k_z}$ resides on a sphere (${k_{z,{\rm sphere}}}$), whereas Fresnel propagation approximates ${k_z}$ by a parabola (Fig. 2). Thus the phase error incurred under the Fresnel approximation is given $z\Delta {k_z}$ with $\Delta {k_z} = ({k_{z,{\rm parabola}}} - {k_{z,{\rm sphere}}})$.

Fig. 2. Graphical depiction of the error of Fresnel versus angular spectrum (AS) propagation. The dashed $k$-sphere describes all possible propagating waves solving the Helmholtz equation. Propagation along $z$ requires a phase change of $z \cdot {k_z}$. Yet Fresnel propagation approximates this sphere by a Fresnel-parabola using a wrong ${k_z}$. This phase error to pre-compensate before Fresnel propagation is thus given by $z \cdot \Delta {k_z}$.

Download Full Size | PDF

We may use an approximate wave propagation model and still get the analytically correct result if we succeed in pre-compensating for errors introduced in the approximate model. In the previous subsection, we observed a formal equivalence between the CV-FR and SFT-FR methods. Since the parabolic part of the pre-compensation represents the CV-FR kernel in Fourier space, the formal equivalence suggests that applying the aforementioned pre-compensation kernel should be a valid approach. SAS propagation is thus obtained by spatially filtering the original signal with the difference phase $\Delta {k_z}$,

(10)$$\begin{split}{H_{\rm{AS}}}H_{\rm{Fr}}^* &= \exp ({i\Delta {k_z}z} ) \\&= \exp \!\left({i\frac{{2\pi}}{\lambda}z\left[{\sqrt {1 - {{({\lambda {f_x}} )}^2} - {{({\lambda {f_y}} )}^2}} }\right.}\right.\\&\quad- \left.{\left.{\left({1 - \frac{{{{({\lambda {f_x}} )}^2}}}{2} - \frac{{{{({\lambda {f_y}} )}^2}}}{2}} \right)} \right]} \right),\end{split}$$

and subsequently subjecting it to a single-step Fresnel (SFT-FR) propagation [see Fig. 1(b)]. We refer to the combination of (1) the pre-compensation step [Eq. (10)] and (2) the SFT-FR operation as the scalable angular spectrum (SAS) algorithm. The result is a numerically correct wave field propagated in the sense of the AS. However, unlike the AS, the propagated wave field resulting from the SAS algorithm exhibits a scaled pixel size, thanks to the presence of the SFT-FR propagator in step 2.

C. SAS Propagation in Operator Form

In this subsection, we give an alternative viewpoint on the pre-compensation method described in the last section. Both of these viewpoints are formally equivalent. The overall propagation is described by the following series of operations (read from right to left):

(11)$${\psi _z} = \underbrace {{Q_2}{\cal F}{Q_1}}_{{\rm SFT} \text{-} {\rm FR}}\underbrace {{{\cal F}^{- 1}}{H_{{\rm AS}}}H_{{\rm Fr}}^*{\cal F}}_{{\rm pre} \text{-} {\rm compensation}}{\psi _0}.$$

The pre-compenstation assures numerical correctness, whereas the single-step Fresnel transform (SFT-FR) assures appropriate scaling of the pixel pitch.

Interestingly, by introducing an identity as a pair of additional Fourier transformations, we can reinterpret the SAS method:

(12)$${\psi _z} = \underbrace {{Q_2}{\cal F}{Q_1}}_{{\rm SFT} \text{-} {\rm FR} ({+ z} )}\underbrace {{{\cal F}^{- 1}}{H_{{\rm AS}}}{\cal F}}_{{\rm AS} ({+ z} )}\underbrace {{{\cal F}^{- 1}}H_{{\rm Fr}}^*{\cal F}}_{{\rm CV}\text{-}{\rm FR} ({- z} )}{\psi _0},$$

where the complex conjugation in the CV-FR transfer function accounts for a back-propagation and ${Q_1}$ for the SFT-FR results in forward propagation by $z$. In this interpretation, we have a cascade of three operations: (1) a convolution-based Fresnel transform backwards (${\rm CV}\text{-}{\rm FR}(- { z})$) followed by (2) an angular spectrum propagation forwards (${\rm AS}(+ { z})$) before applying (3) the final single-step Fresnel (${\rm SFT}\text{-}{\rm FR}(+ { z})$) forward transform taking care of the pixel-pitch scaling. Since the two Fresnel-transform steps are formally equivalent and cancel another, we are left with an angular spectrum propagation supporting a scaled pixel pitch.

The quadratic phase exponential ${Q_2}$ can be omitted if only the intensity at the destination plane is of interest. In the simplified form of Eq. (11), we may interpret the SAS in an alternative way: the conjugate CV-FR transfer function pre-compensates for the numerical error introduced by the SFT-FR method. As compared to the SFT-FR method, the SAS algorithm requires two additional FFTs and a multiplication with a transfer function (the product ${H_{{\rm AS}}}H_{{\rm Fr}}^*$ can be pre-allocated). This reduction in speed is a minor price to pay for the gain in accuracy in many high NA applications, for which Fresnel diffraction is too imprecise. In addition, compared to standard AS, SAS propagation requires only a single additional FFT to provide pixel scaling. Thus SAS propagation combines the best of both worlds between the SFT-FR and AS methods at only moderately increased computational overhead.

D. Implementation Details

A pseudocode implementation is sketched in Algorithm 1. An important detail is the use of padding, which is required to avoid circular convolution artifacts in the pre-compensation step, which is essentially a convolution. The input field with size $(N \times N)$ is padded with zeros to obtain an array with size $(2N \times 2N)$. The complex-valued kernels ${Q_2},{Q_1},{H_{{\rm AS}}},H_{\rm{FR}}^*$ are fully evaluated on a grid with $(2N \times 2N)$. At the end of our algorithm, we return the cropped center region of size ($N \times N)$ only. This region is free of circular convolution artifacts. The total magnification of the SAS depends on the magnification of the Fresnel transform, the latter of which is given by

(13)$${M_{{\rm SFT}\text{-}{\rm FR}}} = \frac{{{\Delta _d}}}{{{\Delta _s}}} = \frac{{\lambda \cdot z \cdot N}}{{{L^2}}},$$

where we used the source pixel pitch ${\Delta _s} = L/N$ and the destination pixel pitch ${\Delta _d}$ as given by Eq. (9). To warrant proper compensation of both Fresnel-propagation steps, we Fresnel-propagate not the initial field of size $(N \times N)$ and side length $L$ but instead the two-fold zero-padded array. By inserting $L \to {L_p} = 2 \cdot L$ and $N \to {N_p} = 2 \cdot N$ into Eq. (13), we obtain the magnification of the SAS:

(14)$${M_{{\rm SAS}}} = \frac{{{\Delta _d}}}{{{\Delta _s}}} = \frac{{\lambda \cdot z \cdot N}}{{2 \cdot {L^2}}}.$$

Algorithm 1 mentions a ${z_{{\rm limit}}}$ and a bandlimit window $W$, which we discuss in the next subsection.

Algorithm 1. SAS implementation

View Table | View all tables in this article

E. Bandlimit Derivation

The first two Fourier transformations constitute a convolution of a twice zero-padded amplitude area with a full-sized kernel. Even though this convolution can cause wrap-around artifacts in the calculation array, the central region of interest is guaranteed to be free of such effects, independent of the kernel with any pixel still being able to be affected by any other in the central region via the kernel. Considering the final one-step propagation, the precompensation moves local spatial amplitude information each to a different starting position. Information that has been wrapped fully around close to the opposite boarder is then only affected by the application of the purely divergent ${Q_1}$ term away from the final region of interest. Interestingly, the precompensation kernel transfers a wavelet with a particular k-vector into the opposite direction, which leads to this wavelet being propagated in the final step away from the central region of interest. As long as ${Q_1}$ is not undersampled by itself, we do not expect such wrapped information to influence the destination area.

Previously, Matsushima and Shimobaba derived an analytical expression for a spatial frequency cutoff region for standard AS propagation, calculated on a regularly sampled array. The method is referred to herein as bandlimited AS propagation [37]. By applying the concept of a local signal frequency, the authors derived conditions under which ${H_{{\rm AS}}}$ is Nyquist sampled, resulting in a two-dimensional frequency bandpass filter support over which the AS transfer function is retained and outside of which it is set to zero. The aforementioned work showed that in the case of standard AS propagation, the spatial frequency support is the intersection of two elliptical regions. We follow a similar approach here with the difference that we aim at Nyquist sampling the argument of ${H_{{\rm AS}}}H_{{\rm Fr}}^*$ instead of ${H_{{\rm AS}}}$.

We start by extracting the phase of the transfer function of the SAS [Eq. (10)] as

(15)$$\begin{split}{\phi _{2D}} &= \arg \big({{H_{\rm{AS}}}H_{\rm{Fr}}^*} \big) \\&= \frac{{2\pi}}{\lambda}z\left[{\sqrt {1 - {{({\lambda {f_x}} )}^2} - {{({\lambda {f_y}} )}^2}} - \left({1 - \frac{{{{\left({\lambda {f_x}} \right)}^2}}}{2} - \frac{{{{({\lambda {f_y}} )}^2}}}{2}} \right)} \right]\!.\end{split}$$

The corresponding local signal frequencies in the $x$ and $y$ directions are given by

(16)$${\Omega _x} = \left| {\frac{1}{{2\pi}}\frac{{\partial {\phi _{2D}}}}{{\partial {f_x}}}} \right| = \left| {z\left[{\frac{{\lambda {f_x}}}{{\sqrt {1 - {{({\lambda {f_x}} )}^2} - {{({\lambda {f_y}} )}^2}}}} - \lambda {f_x}} \right]} \right|,$$

(17)$${\Omega _y} = \left| {\frac{1}{{2\pi}}\frac{{\partial {\phi _{2D}}}}{{\partial {f_y}}}} \right| = \left| {z\left[{\frac{{\lambda {f_y}}}{{\sqrt {1 - {{\left({\lambda {f_x}} \right)}^2} - {{({\lambda {f_y}} )}^2}}}} - \lambda {f_y}} \right]} \right|.$$

To avoid undersampling, we invoke the Shannon–Nyquist criterion on both local signal frequencies:

(18)$$\Delta {f_x} \le \frac{1}{{2{\Omega _x}}},$$

(19)$$\Delta {f_y} \le \frac{1}{{2{\Omega _y}}}.$$

Using $\Delta {f_{x,y}} = 1/{L_{p,x,y}}$ (where ${L_p}$ is the padded field size), we get

(20)$$\left| {\frac{{\lambda {f_x}}}{{\sqrt {1 - {{\left({\lambda {f_x}} \right)}^2} - {{({\lambda {f_y}} )}^2}}}} - \lambda {f_x}} \right| \le \frac{{{L_{p,x}}}}{{2z}},$$

(21)$$\left| {\frac{{\lambda {f_y}}}{{\sqrt {1 - {{\left({\lambda {f_x}} \right)}^2} - {{({\lambda {f_y}} )}^2}}}} - \lambda {f_y}} \right| \le \frac{{{L_{p,y}}}}{{2z}}.$$

We can manipulate this expression in such a way that it results in a set of four inequalities. These expressions, however, are significantly less compact and yield little practical insight, as they result in inequalities involving fourth-order polynomials in the spatial frequencies ${f_x}$ and ${f_y}$. We thus stick to the form of the conditions given by Ineqs. (20) and (21), in which they can be used to define binary masks in computer implementations of the SAS algorithm. Figure 3(b) shows the resulting bandlimit for the two-dimensional analysis discussed in this subsection. We observe that the spatial frequency support takes on the form of a square deformed into an oval shape. The inset in Fig. 3(c) highlights that the Nyquist limit is precisely reached at the spatial frequency cutoff: a phase shift of $2\pi$ is exactly represented by two pixels along the diagonal direction, resulting in a checkerboard appearance of the depicted phase.

F. Derivation of the Propagation Distance Limit

The SAS exhibits a vignetting artifact for large propagation distances, originating from a conflict between the aforementioned bandlimit for the pre-compensation transfer function ${H_{{\rm AS}}}H_{{\rm Fr}}^*$ and the magnification. The bandlimit, given by Ineqs. (20) and (21), becomes increasingly restrictive with larger propagation distances, effectively clipping parts of the angular spectrum that exhibit large angles with respect to the optical axis. Yet the magnification of the SAS, given by Eq. (14), increases proportionally to the propagation distance. Thus, there exists a propagation distance ${z_{{\rm limit}}}$ where the size of the magnified destination FoV is larger than the size of the bandlimited FoV restricted by the pre-compensation term. We note that we distinguish between field of view (FoV), which is the full size of a field including padding, and the region of interest (RoI), which is only half of the size of the FoV, that is, padding is excluded.

The vignetting effect is illustrated in Fig. 4, where we simulated a diffraction pattern from a randomized diffuser under variable magnification. The choice of a diffuser in this example is of no particular interest here other than the diffraction pattern having a wide range of spatial frequencies so that the vignetting effect can be observed clearly. In Fig. 4(a), we see a partially developed speckle field with no signs of vignetting, while in Fig. 4(b), larger magnification was chosen leading to vignetting. It is the goal of this subsection to derive the maximum propagation distance ${z_{{\rm limit}}}$ where the vignetting effect does not occur. We show below that this limit is conveniently expressed in terms of the dimensionless parameter $R = \Delta x/\lambda$, which is the source pixel pitch divided by the wavelength.

Fig. 3. SAS transfer function. (a) No bandlimit imposed. (b) Exact two-dimensional SAS bandlimit. The inset in (c) in the bottom right corner of the transfer function support shows a checkerboard with two pixels along $x$ or $y$ per $2\pi$ phase shift along a diagonal direction, indicating that the Nyquist limit is precisely reached. In each panel, hue corresponds to the phase $\phi = \arg ({{H_{{\rm AS}}}H_{{\rm Fr}}^*})$, and brightness represents the modulus of the complex-valued transfer function in (a). In (b), (c), it is one (bright) within and zero outside the cutoff frequency region, and the dashed yellow circle represents the evanescent wave cutoff with radius $1/\lambda$.

Download Full Size | PDF

Fig. 4. Diffraction from a diffuser illustrating the vignetting effect. (a) No vignetting occurs (${R} = {0.5},\;{ M} = {5}$). (b) At higher magnification the vignetting effect restricts the field inside the destination region of interest (${R} = {0.5},\;{M} = {15}$). In the main text, we describe how the allowable magnification relates to the dimensionless parameter $R = \Delta x/\lambda$, which is the source pixel pitch divided by the wavelength.

Download Full Size | PDF

We start by reviewing two inequalities first derived by Matsushima and Shimobaba in the context of a bandlimited AS with equal source and destination pixel pitch [37]:

(22)$$\frac{{f_x^2}}{{f_{{\rm limit}}^2}} + {({\lambda {f_y}} )^2} \leqslant 1,$$

(23)$${({\lambda {f_x}} )^2} + \frac{{f_y^2}}{{f_{{\rm limit}}^2}} \leqslant 1,$$

where

(24)$${f_{{\rm limit}}} = \frac{{{L_p}}}{{\lambda \sqrt {L_p^2 + 4{z^2}}}}$$

is the the maximum spatial frequency that can propagate in either the $x$ or $y$ direction. In two-dimensional propagation geometries, Ineqs. (22) and (23) constrain the angular spectrum of the bandlimited AS propagator to the intersection of two elliptical regions. Inequalities (20) and (21) derived in this work for the SAS propagator can be regarded as an analog to Ineqs. (22) and (23) derived for the AS propagator. The same inequalities as given by Matsushima and Shimobaba can be derived, when one considers the maximum diffraction angle that arrives at the edges of the destination RoI. In fact, spatial frequencies beyond the destination RoI edge would result in circular convolution artifacts—which is the very phenomenon the bandlimited AS seeks to effectively avoid, so such frequencies are suppressed. Thus the bandlimit in Eq. (24) has two interpretations (simplified here in only one dimension): (1) it is the maximum spatial frequency that can be Nyquist-sampled in the Fourier domain, and (2) it determines the angle of the geometric ray that would connect the optical axis and the edge of the destination FoV. It is the latter property that we will use for our purposes below.

The next step in the derivation of the propagation distance limit of the SAS is to combine Ineqs. (20) and (21) with Ineqs. (22) and (23). For simplicity, we regard only the corners of the destination RoI, where the vignetting effect is most restricting. At the corners, the maximum angle for which Ineqs. (22) and (23) hold is satisfied with equality, that is,

(25)$$\frac{{f_x^2}}{{f_{{\rm limit}}^2}} + {\lambda ^2}f_x^2 = 1$$

or, equivalently,

(26)$${f_x} = \frac{1}{{\sqrt {\frac{1}{{f_{{\rm limit}}^2}} + {\lambda ^2}}}}.$$

We only aim that half of the destination FoV is free of vignetting since we remove the zero-padding at the very end of the algorithm. Hence we can make further progress by setting ${L_p} \to {M_{{\rm SAS}}}L$. Then Eq. (24) reads

(27)$${f_{{\rm limit}}} = \frac{1}{{\lambda \sqrt {1 + 4{{\left({\frac{z}{{{M_{{\rm SAS}}}L}}} \right)}^2}}}},$$

and Eq. (26) becomes

(28)$$\lambda {f_x} = \frac{1}{{\sqrt {4{{\left({\frac{z}{{{M_{{\rm SAS}}}L}}} \right)}^2} + 2}}} \equiv s,$$

where we defined the variable $s$ as shorthand notation for the direction sine $\lambda {f_x}$. Defining the pixel pitch relative to the wavelength $R \equiv \Delta x/\lambda$, we can express the SAS magnification as ${M_{{\rm SAS}}} = z/({2 \cdot R \cdot L})$. Equation (28) can then be rewritten in a compact form:

Fig. 5. SAS distance propagation limit. The valid (green) range of propagation distances relative to the lateral field size of a square shaped field of view is shown in dependence of pixel pitch relative to the wavelength. The yellow area indicates the minimum SAS propagation distance. Since the pixel pitch defines the sampling, there is a corresponding numerical aperture ${{\rm NA}_{{\rm pixel}}}$, as given by the top axis. A pixel pitch of $\lambda /2$ means that a Nyquist sampled object grating would lead to a 90° diffracted beam of ${{\rm NA}_{{\rm pixel}}} = 1.0$. The green and purple marks correspond to Fig. 4 and the two black marks (${M} = {4}$ and ${M} = {8}$) to Figs. 7 and 8, respectively.

Download Full Size | PDF

(29)$$s = \lambda {f_x} = \frac{1}{{\sqrt {16{R^2} + 2}}}.$$

This expression can now be substituted into Ineq. (20), where we again set ${f_x} = {f_y}$, to get

(30)$$\left| {\frac{s}{{\sqrt {1 - 2{s^2}}}} - s} \right| \le \frac{{{L_p}}}{{2z}}.$$

This result can be formulated equivalently in terms of the maximum propagation distance:

(31)$$z \le {z_{{\rm limit}}} = \frac{L}{{\left| {\frac{1}{{4R}} - \frac{1}{{\sqrt {16{R^2} + 2}}}} \right|}},$$

where we used that a binning of two is assumed, ${L_p} = 2L$. Equation (31) is the main result of this subsection. We note that the maximum propagation distance results in a maximum possible SAS magnification:

(32)$${M_{{\rm SAS}}} \le \frac{{\lambda {z_{{\rm limit}}}}}{{2L\Delta x}}.$$

In Fig. 5, we illustrate the SAS propagation distance limit together with the maximum allowed magnification. The green region displays propagation geometries for which no vignetting occurs. Along the green line, we show the corresponding maximum magnification for a given pixel pitch. Each axis is scaled into dimensionless parameters, yielding parameters independent of wavelength and field of view. Since the last step of SAS propagation includes a single-step Fresnel transformation, there is also a minimum propagation distance as defined by the ${Q_1}$ term. In the case of two-fold padding, this can be shown to correspond to the yellow area: $\frac{z}{L} = 4\frac{{\Delta x}}{\lambda}$, which corresponds to a magnification of one. Note that all of this area is covered by bandlimited AS propagation [37].

3. NUMERICAL TESTS

In this section, we compare the wave propagation methods summarized in Table 1 by means of two canonical problems: diffraction by circular and rectangular apertures. The electrical input fields are shown in Figs. 6(a) and 6(b). The amplitude is encoded as brightness and the phase as color. The general situation is sketched in Fig. 6(c).

Fig. 6. Numerical examples on which we test the methods. (a), (b) Electrical input fields. (c) Schematic sketch of the considered situation. For space reasons, we do not show the final uncropped field of size $M \cdot {L_p}$. (a) Circular aperture; (b) rectangular aperture.

Download Full Size | PDF

A. Two Cases

In the first example, a circular aperture is illuminated by two oblique waves. The scalar electrical field immediately behind the aperture is

(33)$$\begin{split}{U_ \circ}({x,y} ) &= {\rm circ}\left({\frac{x}{{{D_ \circ}}},\frac{y}{{{D_ \circ}}}} \right) \cdot \left[{\exp \!\left({i\frac{{2\pi}}{\lambda}y\sin ({{\alpha _ \circ}} )} \right) }\right.\\&\quad+ \left.{\exp \!\left({i\frac{{2\pi}}{\lambda}x\sin ({- {\alpha _ \circ}} )} \right)} \right],\end{split}$$

where the circ-function is defined as

(34)$${\rm circ}({x,y} ) = \left\{{\begin{array}{ll}1,&{{x^2} + {y^2} \le {{\left({\frac{1}{2}} \right)}^2}}\\0,&{{\rm else}}\end{array}} \right.,$$

and the diameter is ${D_ \circ} = {L_ \circ}/8$ and ${\alpha _ \circ}={ 45^ \circ}$.

Fig. 7. Display of the intensity of ${U_ \circ}$ propagated over the distance ${z_ \circ}$ with (a) AS (computing time 1.4 s), (b) SFT-FR (0.04 s), and (c) SAS (0.11 s). The images show the absolute square of the resulting field.

Download Full Size | PDF

Fig. 8. Display of the intensity of ${U_\square}$ propagated over the distance ${z_\square}$ with (a) AS (computing time 6.1 s), (b) SFT-FR (0.04 s), and (c) SAS (0.11 s). The images show the absolute square of the resulting field. The SFT-FR is clearly not able to account for the expected distortion, whereas the AS and SAS produce nearly identical results.

Download Full Size | PDF

In the second example, a square-shaped aperture is illuminated at oblique incidence. The corresponding electric field immediately behind the aperture is

(35)$${U_\square} = {\rm rect}\left({\frac{x}{{{D_\square}}}} \right){\rm rect}\left({\frac{x}{{{D_\square}}}} \right)\exp \!\left({i\frac{{2\pi}}{\lambda}y\sin ({{\alpha _\square}} )} \right),$$

where the rect function is defined as

(36)$${\rm rect}(x ) = \left\{{\begin{array}{ll}1,&{\left| x \right| \le 1/2}\\0,&{{\rm else} .}\end{array}} \right.$$

The rectangular support has a side length of ${D_\square} = {L_\square}/16$, and the angle of incidence with respect to the optical axis is ${\alpha _\square}={ 20^ \circ}$.

In the examples, ${L_ \circ} = 64\;{\unicode{x00B5}{\rm m}},{L_\square} = 128\;{\unicode{x00B5}{\rm m}}$, is the side length of the array, $\lambda = 500\;{\rm nm}$ is wavelength, and the number of sample points per dimension is ${N_ \circ} = {N_\square} = 512$. The propagation distance from the source to the destination plane is ${z_ \circ} = 128\;{\unicode{x00B5}{\rm m}},{z_\square} = 1000\;{\unicode{x00B5}{\rm m}}$. This corresponds to a pixel magnification of ${M_ \circ} = {\Delta _d}/{\Delta _s} = 4$ for ${U_ \circ}$ and ${M_\square} = 8$ for ${U_\square}$, when using SAS algorithms to compute the field in the destination plane. The problems were designed for the investigation of variable magnification and aperture geometries. In addition, the oblique illumination aims at highlighting typical truncation or taper-off effects that can occur in bandlimited wave propagation methods towards high diffraction angles.

B. Results

The results of the numerical simulations (computed with an AMD Ryzen 9 5900X CPU) are presented in Figs. 7 and 8. For both examples, we display the absolute square of the field after propagation. We observe that the SFT-FR is clearly unable to reproduce the correct results of the padded AS and the SAS.

In Fig. 8, we display the intensity of the final field. The most notable difference between SAS and SFT-FR is that the latter is unable to describe the final field curvature appropriately. It should be noted that for the AS implementation, we had to zero-pad the initial field with zeros to prevent wrap-around artifacts. For the rectangular case, this resulted in an array that has 64 times more pixels than the initial one. Hence, the run time is much worse than the run time of the SAS method. SFT-FR and SAS perform almost equally fast; we used a resampling step in the post processing to achieve the same resolution and field size at the end for comparison. We also compare computation times spent to calculate the final fields. From Fig. 8, it can be seen that the precise AS method takes roughly 50 times longer than the SAS method, which is due to the aforementioned zero-padding. To compare the simulations quantitatively to each other, we calculated

(37)$$\sigma : = \frac{{\sum\nolimits_{x,y} {|{\Psi _z}(x,y) - {\Psi _{z,{\rm ref}}}(x,y{{)|}^2}}}}{{\sum\nolimits_{x,y} {|{\Psi _{z,{\rm ref}}}(x,y{{)|}^2}}}}$$

of the propagated amplitudes as described in [38]. The difference for the quadratic aperture is ${\sigma _{{\rm SAS},\square}} \approx 0.03\%$ and for the circular aperture ${\sigma _{{\rm SAS}, \circ}} \approx 1.3\%$, whereas the single-step Fresnel transformation yielded ${\sigma _{{\rm SFT}\text{-}{\rm FR},\square}} \approx 103\%$ and ${\sigma _{{\rm SFT}\text{-}{\rm FR}, \circ}} \approx 302\%$, respectively. The larger error for the SAS in the circular case is caused by waves that leave the FoV. Hence, most of the error is accumulated because of boundary handling artifacts.

4. DISCUSSION

The SAS method is based on pre-compensating for errors introduced by single-step Fresnel propagation by multiplying with the frequency-dependent pre-compensation factor in Fourier space. As described in Section 2.C, the precompensation may be viewed as the cascade of three operations: (1) CV-FR (${-}z$), (2) AS (${+}z$), and (3) SFT-FR (${+}z$). Since operation (1) pre-compensates for operation (3), we may introduce an additional degree of freedom into the SAS algorithm by using a modified cascade of operations: (1) CV-FR (${-}{z_1}$), (2) AS (${+}{z}$), and (3) SFT-FR (${+}{z_1}$).

This modified pre-compensation has two advantages: (1) the user can (within limits) freely choose the destination pixel size, which led to the name “scalable angular spectrum” rather than “scaled angular spectrum” method; (2) a smart choice of a (negative) pre-propagation distance can reduce the maximal slope of the pre-compensation term, extending the clipping-free propagation limits of the SAS method, as analyzed above, further. Further details will be described in future work.

It is worth noting that for the pre-compensation term, we use the analytical formulation in Fourier space, whereas this can also be interpreted as a convolution with a kernel in real space. Even though these two interpretations are formally equivalent, there are subtle differences particularly if the transfer function in Fourier space supersedes the frequency support governed by the fast Fourier transform operation. This amounts to undersampling in the source field, which has to be avoided.

The SAS method was presented in scalar terms. However, it is obvious that the considerations apply equally well to each of the three electric field components individually. We thus obtain the vectorial field propagation SAS method by propagating each complex-valued field component individually and reassembling the field vector in the destination plane. Alternatively, we can also propagate two field components (e.g., ${E_x}$ and ${E_y}$) and compute the third component (${E_z}$) from consistency considerations with respect to Maxwell’s equations. See [40] for details.

SAS propagation is particularly useful for propagation distances causing a magnification larger than one. Thus distances below ${z_{M = 1}} = 2RL$ should normally be treated by standard AS propagation with zero-padding to an array of twice the source size. However, the adjoint operator of the SAS algorithm is well defined and can be regarded as backwards propagation with supersampling of the field. Even though the SAS algorithm has many advantages as described, it also has its limits as seen by the white area in Fig. 5. Since these areas are essentially caused by the need to avoid vignetting effects, they can be covered by zero-padding beyond the factor of two as suggested above. While this leads to additional computational overhead, such approaches are still expected to be computationally much faster than standard AS propagation under similar conditions. A detailed investigation of these limits and also a detailed comparison to tile-based propagation [28,29] may be of interest in future studies.

Overall the SAS propagation method, as introduced here, should have many applications in fields involving inverse problems, phase recovery, physical optics, and optical modeling. Naturally, the wave fields are not restricted to electrical fields. Of course, applications in related fields dealing with wave propagation such as acoustics can also be expected to be plentiful.

APPENDIX A: COMPARISON TO SEMI-ANALYTICAL FOURIER TRANSFORM

To directly compare with the algorithm described in [38], we simulated (using an iterative Fourier-transformation approach) a synthetic hologram ($462 \times 462 \;{\rm pixels}$ covering $554.6\;{\unicode{x00B5}{\rm m}} \,\times\def\LDeqbreak{} 554.6\;{\unicode{x00B5}{\rm m}}$) generating $5 \times 5$ beams with a destination-plane separation distance of 0.5 mm in $x$ and $y$ directions, illuminated by a coherent Gaussian beam (532 nm of 100 µm beam waist diameter and propagated it by 10 mm as stated). Comparing the SAS result ($924 \times 924\;{\rm pixels}$ after two-fold padding as described) to AS propagation (padded to $18501 \times 18501\;{\rm pixels}$) yielded an agreement to better than 0.0013% using the above norm for comparison of the propagated amplitudes. As seen from the logarithmic plot in Fig. 3, the parameters of [38] are within the SAS bandlimit by a large margin, and the required number of pixels is well below the stated value of $2053 \times 2053\;{\rm pixels}$ for the semi-analytical SPW (Taylor) or $1252 \times 1252$ for the semi-analytical SPW (Avoort) methods [38].

Funding

Deutsche Forschungsgemeinschaft (316213987, TR 1278 Teilprojekt C04).

Acknowledgment

RH acknowledges financial support by the DFG. All simulations were performed in Julia v1.8.5 [41] and made use of FFTW.jl [42]. We thank Colin J.R. Sheppard for his suggestions and discussion and the reviewers for extensive comments and useful suggestions!

Disclosures

RH: no conflict of interest. LL: currently at Carl Zeiss Microscopy (E). FW: no conflict of interest.

Data availability

A Julia and Python+PyTorch implementation of the SAS algorithm can be found on our GitHub page [43].

REFERENCES

1. J. R. Fienup, “Phase retrieval algorithms: a personal tour,” Appl. Opt. 52, 45–56 (2013). [CrossRef]

2. U. Schnars, C. Falldorf, J. Watson, et al., “Digital holography and wavefront sensing,” in Digital Holography (Springer, 2015).

3. J. Rodenburg and A. Maiden, “Ptychography,” in Springer Handbook of Microscopy (Springer, 2019), p. 2.

4. J. Miao, T. Ishikawa, I. K. Robinson, et al., “Beyond crystallography: Diffractive imaging using coherent x-ray light sources,” Science 348, 530–535 (2015). [CrossRef]

5. D. E. Sands, Introduction to Crystallography (Courier, 1993).

6. B. Wang, M. Tanksalvala, Z. Zhang, et al., “Coherent Fourier scatterometry using orbital angular momentum beams for defect detection,” Opt. Express 29, 3342–3358 (2021). [CrossRef]

7. I. Yamaguchi, “Phase-shifting digital holography,” in Digital Holography and Three-Dimensional Display (Springer, 2006), pp. 145–171.

8. C. Falldorf, C. von Kopylow, and R. B. Bergmann, “Wave field sensing by means of computational shear interferometry,” J. Opt. Soc. Am. A 30, 1905–1912 (2013). [CrossRef]

9. M. Du, L. Loetgering, K. Eikema, et al., “Measuring laser beam quality, wavefronts, and lens aberrations using ptychography,” Opt. Express 28, 5022–5034 (2020). [CrossRef]

10. C. Zuo, J. Li, J. Sun, et al., “Transport of intensity equation: a tutorial,” Opt. Lasers Eng. 135, 106187 (2020). [CrossRef]

11. J. W. Goodman, Introduction to Fourier Optics, 4th ed. (WH Freeman, 2017).

12. D. G. Voelz, Computational Fourier Optics: A MATLAB Tutorial (SPIE, 2011).

13. J. Schmidt, Numerical Simulation of Optical Wave Propagation With Examples in MATLAB, 1st ed. (SPIE, 2010).

14. T.-C. Poon and J.-P. Liu, Introduction to Modern Digital Holography: With MATLAB (Cambridge University, 2014).

15. P. C. Konda, L. Loetgering, K. C. Zhou, et al., “Fourier ptychography: current applications and future promises,” Opt. Express 28, 9603–9630 (2020). [CrossRef]

16. G. Zheng, C. Shen, S. Jiang, et al., “Concept, implementations and applications of Fourier ptychography,” Nat. Rev. Phys. 3, 207–223 (2021). [CrossRef]

17. W. Choi, C. Fang-Yen, K. Badizadegan, et al., “Tomographic phase microscopy,” Nat. Methods 4, 717–719 (2007). [CrossRef]

18. J. M. Soto, J. A. Rodrigo, and T. Alieva, “Optical diffraction tomography with fully and partially coherent illumination in high numerical aperture label-free microscopy [Invited],” Appl. Opt. 57, A205–A214 (2018). [CrossRef]

19. K. Prakash, B. Diederich, S. Reichelt, et al., “Super-resolution structured illumination microscopy: past, present and future,” Philos. Trans. R. Soc. A 379, 20200143 (2021). [CrossRef]

20. Y. Park, C. Depeursinge, and G. Popescu, “Quantitative phase imaging in biomedicine,” Nat. Photonics 12, 578–589 (2018). [CrossRef]

21. H. Stark, Applications of Optical Fourier Transforms (Academic, 1982).

22. T. Shimobaba, K. Matsushima, T. Kakue, et al., “Scaled angular spectrum method,” Opt. Lett. 37, 4128–4130 (2012). [CrossRef]

23. W. Zhang, H. Zhang, and G. Jin, “Band-extended angular spectrum method for accurate diffraction calculation in a wide propagation range,” Opt. Lett. 45, 1543–1546 (2020). [CrossRef]

24. K. Matsushima, “Shifted angular spectrum method for off-axis numerical propagation,” Opt. Express 18, 18453–18463 (2010). [CrossRef]

25. H. Wei, X. Liu, X. Hao, et al., “Modeling off-axis diffraction with the least-sampling angular spectrum method,” Optica 10, 959–962 (2023). [CrossRef]

26. X. Yu, T. Xiahui, Q. Y. Xiong, et al., “Wide-window angular spectrum method for diffraction propagation in far and near field,” Opt. Lett. 37, 4943–4945 (2012). [CrossRef]

27. M. Leutenegger, R. Rao, R. A. Leitgeb, et al., “Fast focus field calculations,” Opt. Express 14, 11277–11291 (2006). [CrossRef]

28. M. Hillenbrand, A. Hoffmann, D. P. Kelly, et al., “Fast nonparaxial scalar focal field calculations,” J. Opt. Soc. Am. A 31, 1206–1214 (2014). [CrossRef]

29. M. Hillenbrand, D. P. Kelly, and S. Sinzinger, “Numerical solution of nonparaxial scalar diffraction integrals for focused fields,” J. Opt. Soc. Am. A 31, 1832–1841 (2014). [CrossRef]

30. A. S. Jurling, M. D. Bergkoetter, and J. R. Fienup, “Techniques for arbitrary sampling in two-dimensional Fourier transforms,” J. Opt. Soc. Am. A 35, 1784–1796 (2018). [CrossRef]

31. L. Greengard and J.-Y. Lee, “Accelerating the nonuniform fast Fourier transform,” SIAM Rev. 46, 443–454 (2004). [CrossRef]

32. D. Paganin, Coherent X-Ray Optics, 1st ed. (Oxford University, 2006).

33. S. Coy, “Choosing mesh spacings and mesh dimensions for wave optics simulation,” Proc. SPIE 5894, 589405 (2005). [CrossRef]

34. C. Rydberg and J. Bengtsson, “Efficient numerical representation of the optical field for the propagation of partially coherent radiation with a specified spatial and temporal coherence function,” J. Opt. Soc. Am. A 23, 1616–1625 (2006). [CrossRef]

35. D. G. Voelz and M. C. Roggemann, “Digital simulation of scalar optical diffraction: revisiting chirp function sampling criteria and consequences,” Appl. Opt. 48, 6132–6142 (2009). [CrossRef]

36. X. Deng, B. Bihari, J. Gan, et al., “Fast algorithm for chirp transforms with zooming-in ability and its applications,” J. Opt. Soc. Am. A 17, 762–771 (2000). [CrossRef]

37. K. Matsushima and T. Shimobaba, “Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields,” Opt. Express 17, 19662–19673 (2009). [CrossRef]

38. D. Asoubar, S. Zhang, F. Wyrowski, et al., “Efficient semi-analytical propagation techniques for electromagnetic fields,” J. Opt. Soc. Am. A 31, 591–602 (2014). [CrossRef]

39. Z. Wang, S. Zhang, O. Baladron-Zorita, et al., “Application of the semi-analytical Fourier transform to electromagnetic modeling,” Opt. Express 27, 15335–15350 (2019). [CrossRef]

40. F. Wyrowski and M. Kuhn, “Introduction to field tracing,” J. Mod. Opt. 58, 449–466 (2011). [CrossRef]

41. J. Bezanson, A. Edelman, S. Karpinski, et al., “Julia: a fresh approach to numerical computing,” SIAM Rev. 59, 65–98 (2017). [CrossRef]

42. M. Frigo and S. Johnson, “The design and implementation of FFTW3,” Proc. IEEE 93, 216–231 (2005). [CrossRef]

43. F. Wechsler and R. Heintzmann, “Implementations of the scalable angular spectrum method,” GitHub (2023) https://www.github.com/bionanoimaging/Scalable-Angular-Spectrum-Method-SAS.

Method	Abbreviation	# FFTs	Pixel Pitch
Angular spectrum	AS	2 FFTs	$Δ_{d} = Δ_{s}$
Convolution-based Fresnel	CV-FR	2 FFTs	$Δ_{d} = Δ_{s}$
Single-step Fresnel	SFT-FR	1 FFT	$Δ_{d} = \frac{λ z}{N \cdot Δ_{s}}$
Scalable angular spectrum (this work)	SAS	3 FFTs	$Δ_{d} = \frac{λ z}{2 N \cdot Δ_{s}}$

Method	Abbreviation	# FFTs	Pixel Pitch
Angular spectrum	AS	2 FFTs	$Δ_{d} = Δ_{s}$
Convolution-based Fresnel	CV-FR	2 FFTs	$Δ_{d} = Δ_{s}$
Single-step Fresnel	SFT-FR	1 FFT	$Δ_{d} = \frac{λ z}{N \cdot Δ_{s}}$
Scalable angular spectrum (this work)	SAS	3 FFTs	$Δ_{d} = \frac{λ z}{2 N \cdot Δ_{s}}$

Scalable angular spectrum propagation

Abstract

1. INTRODUCTION

2. THEORY

A. AS and Fresnel Diffraction

B. Pre-compensation and the Scalable Angular Spectrum Algorithm

C. SAS Propagation in Operator Form

D. Implementation Details

E. Bandlimit Derivation

F. Derivation of the Propagation Distance Limit

3. NUMERICAL TESTS

A. Two Cases

B. Results

4. DISCUSSION

APPENDIX A: COMPARISON TO SEMI-ANALYTICAL FOURIER TRANSFORM

Funding

Acknowledgment

Disclosures

Data availability

REFERENCES

Data availability

Cited By

Figures (8)

Tables (2)

Equations (37)

Optica