Open Access
Add to BibSonomy
Abstract
Numerical propagation calculation is a fundamental research topic in optical engineering. The standard angular spectrum method (ASM) is accurate but time- and memory-consuming, especially for high-NA systems. In this work, we propose a fast and simple numerical propagation method, the resampling ASM (RS-ASM). Numerical propagation can be accelerated by combining a resampling technique with interpolation methods in the angular spectrum domain of a constrained object at the focal plane. RS-ASM has three main advantages: simple implementation, faster calculation than the standard ASM, and SNR enhancement. Here we validate RS-ASM using theory, simulation and experiment. Using the “bilinear” ASM with a proper resampling factor can result in a speed-up factor of up to 20x (for a transformation from the angular spectrum to the E field) and 4x (for a transformation from E field to the angular spectrum), together with a SNR improvement of approximately 2x. For an application example of Gerchberg-Saxton phase reconstruction, the “bilinear” RS-ASM can converge 2.6x faster than the standard ASM.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Numerical propagation is a preliminary topic of investigation for calculating the diffraction of wave fields in the optical sciences and has a wide range of applications, including digital holography (DH) [1,2], phase reconstruction techniques [3,4], coherent diffraction imaging (CDI) [5,6], computer-generated holography (CGH) [7,8], and diffraction tomography (DT) [9]. Two numerical propagation calculation methods are widely used: the paraxial Fresnel diffraction method [10] (FDM) and the angular spectrum method [11] (ASM). FDM is limited to a small numerical aperture (NA) and hence, a small spatial resolution. In contrast, ASM is a rigorous solution directly derived from the Helmholtz equation without the paraxial assumption and can be applied to a high-NA system [12,13]. However, the main drawback of ASM is the enormous calculation complexity resulting from having to perform a fast Fourier transform (FFT) twice, especially when the image size ($N^2$) is large. This calculation complexity is amplified in a high-NA system because of the small spatial pixel size $\Delta _x$ and large field size $L_x$ used to cover the full diffraction field at the distant plane, which leads to a large $N=L_x/\Delta _x$.
Several methods have been developed to speed up ASM for strongly diffractive cases, such as high-NA systems or long-distance propagation. The tile super-position propagation technique (TSP) [13,14] has been demonstrated for digital in-line holography. Scaled ASM using the nonuniform fast Fourier transform (NUFFT) [15,16] reduces the calculation waste of a FFT and performs well in a high-NA scheme. Propagation between spherical surfaces [17–19] can be accelerated by several orders over that achieved using the standard ASM between planes due to angle-dependent and thus self-rescaling pixels. Nonuniform ASM [20] has been proven to speed up and generalize propagation in a complex inhomogeneous medium. Band-extended ASM [21] can be used to calculate near-field and far-field diffraction with high accuracy and has a low computational complexity by NUFFT- based rearrangement of the sampling points in the spatial frequency domain. A novel sampling strategy with compact space bandwidth product representation [22,23] that minimizes the space bandwidth product for a given optical signal has been demonstrated to exhibit high performance and high accuracy. Moreover, a detailed frequency sampling strategy for numerical diffraction calculations [24] has been presented for accurate and efficient computation in different use cases.
In this work, we propose a fast and simple numerical propagation method, the resampling ASM (RS-ASM), to accelerate high-NA propagation and improve the SNR. Generally in a high-NA system, the wave field is well focused and constrained in a small region at the focal plane $z_0$, while spreading out with a large spatial extension at the diffraction plane $z_H$. The standard ASM requires that the same sampling interval $\Delta _x$ and image size $L_x$ be used between $z_0$ and $z_H$. To reach the desired spatial resolution $\delta _x=\lambda /NA$, the sampling rate should be at least twice the signal bandwidth, thus $\Delta _x<0.5\delta _x$[14]. Meanwhile, the large spatial distribution at $z_H$ requires a large $L_x$, resulting in an enormous total image size $N^2=(L_x/\Delta _x)^2$ (typically $N=10^3$ to $10^5$ in one dimension), which is highly challenging computationally. We can employ RS-ASM to accelerate the calculation at $z_0$, where the signal region is localized at the image center, surrounded by background noise. By means of a sampling conversion together with an interpolation method at the angular spectrum domain, we can significantly reduce the calculation complexity while enhancing the SNR. RS-ASM has three main advantages: Firstly, the method is simple and easy to apply in a standard ASM use case without changing any hardware. Secondly, the calculation time can be reduced as a function of the resampling factor $\beta$. For an application example, the convergence time for G-S phase reconstruction can be reduced from 128 s (the standard ASM) to 50 s (RS-ASM) with $N=8192$ and $\beta =32$ (corresponding to a 2.6x speed up). Thirdly, the SNR can be typically enhanced by a factor of $2$ due to noise suppression during interpolation. By comparison with other works, TSP [13,14] can have a faster speed but also a higher noise level than our RS-ASM due to the summation of all tiles. RS-ASM does not change the spatial sampling interval $\Delta _x$ as opposed to scaled ASM [15,16] so the images in the spatial domain are always with the highest resolution during propagation. For most situations, the application of RS-ASM is straightforward rather than the propagation of the spherical surfaces [17–19] which may need a spherical coordination conversion. Moreover, RS-ASM may be combined with others such as band-extended ASM [21] to further reduce the computation effort at the object plane and get additional acceleration. We demonstrate our method by theoretical proof, simulations of angular spectrum transformation and G-S phase reconstruction, and an experiment on cold-atom in-line digital holography.
2. Methods
Figure 1 shows a typical high-NA in-line digital holography system. The object is located at $z=z_0$ and illuminated by a known reference laser $E_r$ (a spherical wave is considered as an example). The object is excited and generates a coherent scattering wave $E_s$ behind $z_0$. The intensity distribution of the interference between $E_r$ and $E_s$(i.e., a hologram) is recorded by a CCD placed at the out-of-focus diffracted plane $z=z_H$. A high-NA system is required to image small spatial details of objects. The high-NA can be obtained by using a wide CCD sensor size (typically several $mm$) and a short imaging distance (typically several $mm$). The holographic reconstruction procedure can be divided into two steps. At first, the captured holograms and the known reference wave $E_r$ are used to retrieve the object wave $E_s$ at $z_H$ using in-line holographic methods [25–28]. The second step consists of numerically propagating $E_s$ from $z_H$ to $z_0$ to obtain the focused object image. The popular ASM can be applied for short distance propagation ($z<\frac {N(\Delta x)^2}{\lambda }$) [2] but usually requires a resizing factor $\beta$ to achieve the desired spatial resolution at both the $z_0$ and $z_H$ plane. For simplicity, we consider the same pixel size $\Delta x$ and total pixel numbers $N$ in both the x and y direction.
Fig. 1. Scheme of high-NA in-line digital holography. A laser beam with a known wavefront $E_r$ acts as a reference light. The object at $z_0$ is excited and generates a coherent scattering wave $E_s$. A CCD is placed at the diffraction plane $z_H$ to record the interferogram between $E_r$ and $E_s$ (hologram).
2.1 Standard ASM
The relation between an object wave $E_s(x,y)$ and its angular spectrum $A(m,n)$ at plane z is shown as:
where $m,n,x,y = 0,1, \ldots N-1$. Using standard ASM, we can calculate the propagation of angular spectrum: where $\lambda$ is the wavelength. $\Delta k=\frac {2\pi }{N\Delta x}$ is the pixel size of Fourier space. $L$ is propagation distance. The standard ASM procedure is described in Fig. 2: propagating from $z_0$ to $z_H$ is Fig. 2(a$\to$b$\to$c$\to$d) with $L=z_H-z_0$, and the reversed propagation is Fig. 2(d$\to$c$\to$b$\to$a) with $L=z_0-z_H$. ASM requires to perform 2D FFT/IFFT twice and multiplication once. In a high-NA system, the FFT calculation could be the most expensive part due to wide $E_s$ extension at diffraction plane $z_H$ and small feature size at focal plane $z_0$ so a large $N$ is required. Figure 2(g) shows with increasing $N$, the time cost of FFT$[N^2]$ can increase remarkably. This time-consuming problem restricts ASM applied in fast imaging of high-NA systems [13,14]. Fig. 2(g) also shows the high memory cost of a 2D complex matrix at large $N$. We perform all our calculations using MATLAB 2018b on PC platform Intel Core i7-6700K CPU @ 4.00GHz.Fig. 2. Schematic diagram of numeric propagation of standard ASM (a$\leftrightarrow$b$\leftrightarrow$c$\leftrightarrow$d) and proposed RS-ASM (a$\leftrightarrow$e$\leftrightarrow$f$\leftrightarrow$b$\leftrightarrow$c$\leftrightarrow$d). (a,d) The original object wave $E_s(x,y)$ at $z_0$ (a) and $z_H$ (d) with $N\times N$ size. (b,c) The original angular spectrum $A(m,n)$ at $z_0$ (b) and $z_H$ (c) with $N\times N$ size. (e) Cut object wave $E_s^D(x',y')$ at $z_0$ with $N'\times N'$ size (the white box in (a)) where $\beta =N/N'$ is the resampling factor. In Fig.2, we take $N=512$, $N'=64$ and $\beta =8$ as an example and zoom (e) and (f) by 5x for visualization. (f) Down-sampled angular spectrum $A^D(m,n)$ at $z_0$ with $N'\times N'$ size. (g) 2D FFT time cost and 2D complex matrix (double data type) memory occupation versus $N$.
2.2 Resampling ASM by the “pick-up” and “0-insert”
We developed RS-ASM to speed up ASM propagation for high-NA systems. Unlike $E_s(z_H)$ which has a large spatial extension (Fig. 2(d)), $E_s(z_0)$ is well focused and constrained within a small region (Fig. 2(a)), which leads to considerable calculation waste of the FFT due to the large blank background area. To avoid a direct FFT$[N^2]$ calculation (Fig. 2(a$\leftrightarrow$b)), we use the down-sampling/up-sampling method with the FFT$[N'^2]$ calculation where $N'$ is much smaller than $N$, as shown in Fig. 2(a$\leftrightarrow$ e$\leftrightarrow$f$\leftrightarrow$b).
To accelerate $A$ to $E_s$ conversion at $z_0$, we firstly down-sample $A(m,n)$ by a factor of $\beta = N/N'$ and use an intermediate variable $A^{(1)}(m,n)$ by “pick-up” of $A(m,n)$:
2.3 Resampling ASM with interpolation
However, the limitation of the “pick-up” RS-ASM is energy loss because only a few parts of $A(m,n)$ are sampled. This procedure can lead to SNR loss when noise is not negligible. To achieve both fast calculation and SNR enhancement, we combine RS-ASM with interpolation algorithms, such as “bilinear” and “bicubic” interpolation.
We can use “bilinear” interpolation ($f_{bilinear}$) [29] to replace $A^{(1)}$ in Eq. (5):
Fig. 3. Interpolation weighting kernels. (a) linear kernels. (b) cubic kernels. Blue: normal kernels; Yellow: rescaled kernels with $\beta =2$; Red: rescaled kernels with $\beta =4$.
The yellow(red) lines in Fig. 3(a,b) show the rescaled linear and cubic weighting kernels under $\beta =2$($4$). Compared to the normal kernel, the rescaled kernels have been stretched along x but reduced in amplitude by a factor of $\beta$ to keep the normalization. So the new sampled points of bilinear(bicubic) is a weighted average of its adjacent $2\beta \times 2\beta$ ($4\beta \times 4\beta$) pixels. The rescaled kernel can be seen as a strong low-pass filter (“smooth” effect in Ref. [30]) and hence it can enhance SNR for images mainly made of low-frequency components (simulation demonstration in Supplement 1 section 2). In this paper, we always use the name “bilinear” and “bicubic” to apply the normal kernels for up-sampling and the rescaled kernels for down-sampling.
As a summary, RS-ASM combines FFT calculation and the up/down-sampling procedure with various interpolation kernels to achieve $A\leftrightarrow E_s$ conversion. We will show in the following sections that for a large $\beta$, the time cost of FFT[$N'^2$] combined with the interpolation and resampling procedure can still be considerably lower than that for the full FFT[$N^2$]. In addition to calculation acceleration and obviously lower memory cost, RS-ASM can enhance the SNR in both the $A\to E_s$ and $E_s\to A$ transformations.
3. Simulations
To determine the time cost and SNR characteristics of the proposed RS-ASM, we first perform simulations of the transformation between $E_s$ and $A$ for the standard USAF target, where the results are shown in Fig. 4, Fig. 5, Fig. 6, and Fig. 7. The USAF pattern is $256\times 256$ in size and located at the center of a $8192\times 8192$ images surrounded by “0-padding”. Figure 4 shows the object wave field $E_s$ in $A \to E_s$ procedure implemented by different methods. $E_s$ is calculated by the standard ASM (Fig. 4(a)), “pick-up” RS-ASM (Fig. 4(b)), “bilinear” RS-ASM (Fig. 4(c)) and “bicubic” RS-ASM (Fig. 4(d)). The resampling factor $\beta$ is chosen as 4 in Fig. 4(b1,c1 and d1), 32 in Fig. 4(b2,c2 and d2) and 64 in Fig. 4(b3,c3 and d3). With $\beta =64$, we meet aliasing noises in Fig. 4(b3,c3 and d3), which originate from the cutting of FOV (yellow box in Fig. 4(a)). So the aliasing can not be fully removed even though the rescaled interpolation kernel (Eq. (13)) was used. Practically we need to evaluate the max spatial extension of $E_s$ and choose an appropriate $\beta$ ($32$ is suggested in this case). Hence we drop the $\beta =64$ group in the following analysis. We plot the time cost (Fig. 4(e)) and SNR (Fig. 4(f)) versus $\beta$ for the different RS-ASM methods and the standard ASM as the baseline. The simulated random noise with an amplitude of 20% of the peak signal is added. The SNR is defined by the ratio between the averaged signal in the patterned region and the RMS value in the nonpatterned region.
Fig. 4. Simulation of the object wave fields in $A\to E_s$ procedure using USAF target ($N=8192$). (a-d) $|E_s|$ calculated by different methods and $\beta$. (a) standard ASM. (b) “pick-up” RS-ASM. (c) “bilinear” RS-ASM. (d) “bicubic” RS-ASM. The re-sampling factor $\beta =4$ in (b1,c1,d1), $\beta =32$ in (b2,c2,d2) and $\beta =64$ in (b3,c3,d3). The green circles in (c3) and (d3) indicate the noticeable aliasing noises. (e) Time cost and (f) SNR versus $\beta$ of different methods. For comparison, all sub-figures in (a-d) have the same pixel size $\Delta x$ and are zoomed into the same size except $\beta =64$ as its original smaller size (half of $\beta =32$).
Fig. 5. Simulation of the angular spectrums in $A\to E_s$ procedure using USAF target ($N=8192$). $A$ is calculated by different methods and $\beta$, corresponding to Fig. 4. (a-d): amplitude image $|A|$. (e-h): phase image $\phi (A)$. First column: standard ASM. Second column: “pick-up” RS-ASM. Third column: “bilinear” RS-ASM. Fourth column: “bicubic” RS-ASM. The re-sampling factor $\beta =4$ in (b1,c1,d1,f1,g1,h1), $\beta =32$ in (b2,c2,d2,f2,g2,h2) and $\beta =64$ in (b3,c3,d3,f3,g3,h3). All sub-figures are rescaled for comparison, and the scale bar in each sub-figure is identical to $\frac {\pi }{32\Delta x}$.
Fig. 6. Simulation of the object wave fields in $E_s\to A$ procedure using USAF target($N=8192$). For a rigorous comparison, a final standard $IFFT[N^2]$ is applied to convert $A$ back to $E_s$. (a-d) $|E_s|$ calculated by different methods and adjusted into the same size: (a) standard ASM. (b) “0-insert” RS-ASM. (c) “bilinear” RS-ASM. (d) “bicubic” RS-ASM. The re-sampling factor $\beta =4$ in (b1,c1,d1) and $\beta =32$ in (b2,c2,d2). (e) Time cost and (f) SNR versus $\beta$ of different methods.
Fig. 7. Simulation of ASMs under different noise level. The SNR of $E_s\to A$ (a) and $A\to E_s$ (f) of standard ASM and proposed RS-ASMs versus various noise level from 10% to 60% of the peak signal. (b-e) and (g-j) show the standard ASM (b,g), “0-insert” RS-ASM (c), “pick-up” RS-ASM (h), “bilinear” RS-ASM (d,i), “bicubic” RS-ASM (e,j) under the noise level equal to 30%.
Figure 4(e) shows that compared to the standard ASM, all the three RS-ASM methods can speed up the $A\to E_s$ procedure. However, in Fig. 4(f), compared to the SNR of the standard ASM, only the “bilinear” and “bicubic” methods can boost the SNR, whereas the “pick-up” deteriorates the SNR. With increasing $\beta$, this acceleration and SNR difference become increasingly prominent. In the extreme case of $\beta =32$, the time cost is 1.22 s (for the standard ASM), 0.01 s (for the “pick-up” RS-ASM), 0.06 s (for the “bilinear” RS-ASM), 0.12 s (for the “bicubic” RS-ASM), and the corresponding SNR are 4.8, 0.8, 10.5 and 8.0. Overall, the “bilinear” RS-ASM exhibits the best performance because there is no significant energy loss during the down-sampling procedure. Noise is effectively suppressed for a sufficient FOV and the procedure is simple and less time-consuming than the “bicubic” method.
To check the interpolation performance in $A$, in Fig. 5, we show the angular spectrums corresponding to Fig. 4 at different methods and $\beta$. Figures 5(a),(e) show the true amplitude and phase distribution of $A$ by standard ASM, while others are from RS-ASMs at various $\beta$. The efficient sampling of $\beta =4$ and $32$ is guaranteed by both the amplitude and phase part of $A$. At $\beta =64$, the under-sampling effect is obvious which can lead to aliasing in $E_s$ (Fig. 4(b3,c3,d3)).
In Fig. 6, we also simulate the reverse procedure $E_s\to A$. For a rigorous comparison, a final standard IFFT[$N^2$] is applied to convert $A$ back to $E_s$. The results obtained using the standard ASM (Fig. 6(a)), “0-insert” RS-ASM (Fig. 6(b)), “bilinear” RS-ASM (Fig. 6(c)) and “bicubic” RS-ASM (Fig. 6(d)) are shown. Figures 6(e),(f) show that all the three RS-ASM methods have a lower time cost and a higher SNR than the standard ASM. The SNR never decreases because all the three up-sampling methods do not lose energy or information. With increasing $\beta$, the calculation speed increases, and the SNR is further enhanced. In the extreme case $\beta =32$, the time cost is 1.09 s (for the standard ASM), 0.24 s (for the “0-insert” RS-ASM), 0.26 s (for the “bilinear” RS-ASM) and 0.32 s (for the “bicubic” RS-ASM), and the corresponding SNR is 4.9, 8.1, 13.1 and 9.5. Once more, the “bilinear” RS-ASM exhibits the best overall performance.
To evaluate the robustness of the different methods against noise, we simulate both the $A\to E_s$ and $E_s\to A$ procedures under different levels of uniform random noise. Figures 7(a)-(e) show $E_s\to A$, and Figs. 7(f)-(j) show $A\to E_s$. Figure 7(a) and Fig. 7(f) show that the “bilinear” RS-ASM has the best SNR performance under all noise level because of a sufficient bandwidth (FOV) for resampling, efficient noise suppression and no significant energy loss. Considering the case of 30% noise as an example, the SNRs for the different ASMs are as follows: $E_s\to A$ (Fig. 7(b),(c),(d), and (e)): 3.6 (the standard ASM), 5.7 (the “0-insert” RS-ASM), 9.1 (the “bilinear” RS-ASM), 6.5 (the “bicubic” RS-ASM); $A\to E_s$ (Fig. 7(g),(h),(i), and (j)): 3.3 (the standard ASM), 0.9 (the “pick-up” RS-ASM), 7.1 (the “bilinear” RS-ASM), 5.4 (the “bicubic” RS-ASM). Therefore, the “bilinear” RS-ASM has the best SNR improvement (of approximately 2x) compared to the standard ASM. In addition to uniform random noises, we made more simulations to proof that the bilinear RS-ASM has the highest SNR improvement among various noise types (in Supplement 1 section 3).
However, it should be pointed out that besides the simpleness and high-speed implementation, the bilinear interpolation has its limitations such as variance (contrast) loss, pyramid or Mach bands artifacts [32]. It is suggested to try different interpolation kernels (bicubic, B-spline, Lanczos or other specially designed kernels) for various types of signals [33,34], or apply extra operations to improve the contrast, such as a “Histogram-Preserving Blending Operator” [35].
4. Experiments
For experimental demonstration, we built an inline holographic microscopy system for cold $^{39}K$ atoms, as shown in Fig. 8. About $10^3$ atoms are laser-cooled to a temperature of several tens of $\mu K$ and trapped in a vacuum cell. A focused laser with $\lambda =770$ nm serves as the reference wave $E_r$. The frequency of $E_r$ is tuned to near the $^{39}K$ D1 resonance of the $4^2S_{1/2}\,F=1$ to $4^2P_{1/2}\,F'=2$ transition with detuning $\Delta$ ($\Delta$ can vary from $-\Gamma$ to $\Gamma$, where $\Gamma =2\pi \times 6$ MHz is the natural linewidth of the excitation state $4^2P_{1/2}$.) to obtain a strong atomic coherent scattering wave $E_s$. A far red-detuned single optical dipole trap (ODT) laser with $\lambda '=780\ nm$ is applied as a harmonic trap to constrain the atoms. Due to the limited space of the vacuum cell, we use an imaging system with $NA=0.3$ to image both $E_r$ and $E_s$ in free space with $z_0$ as the atomic focal plane. A CCD with $1040\times 1392$ pixels each with dimensions of $6.45\times 6.45\ \mu m^2$ is mounted on a movable translation stage to record intensity images at different z positions. The diffraction limit is $\lambda /NA=2.6\ \mu m$. The atom sample has a Gaussian shape with $\sigma _x = 15\ \mu m$ and $\sigma _y = \sigma _z = 1\ \mu m$, which is near the diffraction limit. To achieve the desired imaging resolution, we split each original pixel into 36 subpixels by a factor of 6 in both the x and y directions before holographic reconstruction. So the total image size is $6240\times 6240$ (the non-signal region near the boundary is dropped and the whole image is cut into a square) and each pixel size has dimensions of $1.075\times 1.075\ \mu m^2$.
Fig. 8. Cold $^{39}K$ atoms inline holographic microscopy experiment. (a) Experiment scheme. Cold atoms are trapped in a vacuum cell by a single optical dipole trap (ODT) and illuminated by a focused laser $E_r$ with wavelength $\lambda = 770$ nm. An imaging system with $NA=0.3$ is applied to collect both the reference $E_r$ and the atomic coherent scattering light $E_s$ with atomic focal plane at $z_0$. A CCD is mounted on a movable translation stage to record intensity images (including holograms $|E_r+E_s|^2$) at different z positions. (b) The amplitude and (c) sin of phase difference (refer to spherical wave phase $\phi _0$) distribution of precisely reconstructed reference wave field $E_r$ at $z_0$. (d) Hologram difference image $H_{diff}=|E_s+E_r|^2-|E_r|^2$ at $z_H$. The amplitude $|E_s|$ (e) and phase shift $\phi _{shift}$ (f) images of holographically reconstructed atomic wave $E_s$ at $z_0$ ($\Delta =\Gamma$ as an example). The phase shift is defined as $\phi _{shift}=arg(\frac {E_r+E_s}{E_r})$. The unit “cnt” means the photon count of each pixel.
The accurate holographic reconstruction of $E_s$ has two steps: (1) To reconstruct the reference complex wave field $E_r$; (2) With the full knowledge of $E_r$, to reconstruct the atomic complex wave field $E_s$ from holograms $|E_s+E_r|^2$. In the first step, we turn off the single ODT so atoms are absent. Then we move CCD at different z-planes to take a series of intensity images of reference wave $I_r(z_i)=|E_r(z_i)|^2, (i=1,2,3\cdots )$. Then we carefully retrieve the phase distribution of $E_r$ using a modified Gerchberg-Saxton (G-S) [36,37] algorithm. Of course it would usually be acceptable to just assume $E_r$ as an ideal spherical wave ($E_{r0}\backsim \frac {e^{ikR}}{R}$) in conventional holography experiments [12]. However, in cold atom experiment a narrow laser line-width ($\backsim 1$ MHz) is required due to the narrow atomic natural line-width ($\Gamma = 2\pi \times 6$ MHz for $^{39}K$ $4^2P_{1/2}$), which leads to strong speckle noises and reduces the imaging quality [38]. The retrieved $E_r$ is shown in Fig. 8(b) (the amplitude) and Fig. 8(c) (sin of the phase difference from an ideal spherical wave phase). The “multi-ring patterns” in $E_r$ is an artificial mask placed near the laser focus to increase the complexity of $E_r$ and evaluate the phase reconstruction accuracy. In the second step, we put the CCD at fixed diffraction plane $z_H$ to take holograms. Along with the fully reconstructed $E_r$ from the first step, we can get:
In Eq. (14), $H_{diff}=|E_s+E_r|^2-|E_r|^2$ at $z_H$ is shown in Fig. 8(d). There are three terms at the right hand side: $E_s$ (the target term), twin image and DC image in order. We then apply an iterative algorithm [38] to remove the twin and DC images from Eq. (14). Finally we get the precisely reconstructed $E_s$ as shown in Fig. 8(e) (the amplitude) and Fig. 8(f) (the phase shift). (The holographic reconstruction details are in Supplement 1 section 4. And recently we have built another similar $^{87}Rb$ cold atom inline holographic microscopy setup in Ref. [39].) However, it should be noted that the RS-ASM can be applied to numerically propagate a complex wave field from any complex imaging method (not just holography).To evaluate RS-ASM performance on experimentally retrieved $E_s$, we apply standard ASM and RS-ASMs of $E_s$ at $z_0$. This cold atom system is a good test case because $E_s$ is well focused in a $52\times 6$ region (shown by the white dashed box in Fig. 9(a) and Fig. 11(a)) of a total $6240\times 6240$ matrix, which results in enormous calculation waste using the standard ASM. Figure 9 shows $A\to E_s$ procedure. Compared to the standard ASM (Fig. 9(a)), all the three RS-ASMs (Fig. 9(b),(c) and (d)) have lower time cost, and the “bilinear” and “bicubic” RS-ASMs also have better SNR. For the $\beta =312$ group in Fig. 9(b3,c3,d3), the aliasing noise with green circles shows the destructive interference between the atomic signal and its aliasing image. So we drop $\beta =312$ in the following analysis due to inefficient sampling. Figures 9(e),(f) show the time cost and SNR of all the ASMs with increasing $\beta$. In the extreme case $\beta =120$, the time cost is 0.84 s (for the standard ASM), 0.01 s (for the “pick-up” RS-ASM), 0.04 s (for the “bilinear” RS-ASM) and 0.09 s (for the “bicubic” RS-ASM), and the corresponding SNR is 13.9, 1.45, 21.0 and 17.2.
Fig. 9. Experiment of $^{39}K$ cold atom of the object wave fields in $A\to E_s$ procedure ($N=6240$, $\Delta =0$). (a-d) $|E_s|$ calculated by different methods and $\beta$. (a) standard ASM. (b) “pick-up” RS-ASM. (c) “bilinear” RS-ASM. (d) “bicubic” RS-ASM. The re-sampling factor $\beta =16$ in (b1,c1,d1), and $\beta =120$ in (b2,c2,d2) $\beta =312$ in (b3,c3,d3). The green circles in (c3) and (d3) indicate the noticeable aliasing noises. (e) Time cost and (f) SNR versus $\beta$ of different methods. For comparison, all sub-figures in (a-d) have the same pixel size $\Delta x$ and are zoomed into the same size except $\beta =312$ as its original smaller size than $\beta =120$.
Fig. 10. Experiment of $^{39}K$ cold atom of the angular spectrums in $A\to E_s$ procedure ($N=6240$, $\Delta =0$). $A$ is calculated by different methods and $\beta$, corresponding to Fig. 9. (a-d): amplitude image $|A|$. (e-h): phase image $\phi (A)$. First column: standard ASM. Second column: “pick-up” RS-ASM. Third column: “bilinear” RS-ASM. Fourth column: “bicubic” RS-ASM. The re-sampling factor $\beta =16$ in (b1,c1,d1,f1,g1,h1), $\beta =120$ in (b2,c2,d2,f2,g2,h2) and $\beta =312$ in (b3,c3,d3,f3,g3,h3). All sub-figures are rescaled for comparison, and the scale bar in each sub-figure is identical to $\frac {\pi }{4\Delta x}$.
Fig. 11. Experiment of $^{39}K$ cold atom $E_s\to A$ procedure ($N=6240$, $\Delta =0$). To check performance, a final standard $IFFT[N^2]$ is applied to convert $A$ back to $E_s$. (a-d) $|E_s|$ calculated by different methods and adjusted into the same size: (a) standard ASM. (b) “0-insert” RS-ASM. (c) “bilinear” RS-ASM. (d) “bicubic” RS-ASM. The re-sampling factor $\beta =16$ in (b1,c1,d1) and $\beta =120$ in (b2,c2,d2). (e) Time cost and (f) SNR versus $\beta$ of different methods.
To analyze the SNR enhancement and the alias in experimental data, in Fig. 10, we show the angular spectrums corresponding to Fig. 9 at different methods and $\beta$. Figures 10(a)(e) show the true amplitude and phase distribution of $A$ by standard ASM, while others are from RS-ASMs at various $\beta$. Similar to the simulation result in Fig. 5, the sampling of amplitude and phase part of $A$ using $\beta =16$ and $120$ is sufficient, while $\beta =312$ suffers from the under-sampling effect and leads to alias in $E_s$ (Fig. 9(b3,c3,d3)).
Figure 11 similarly shows the results for the $E_s\to A$ procedure. These results are slightly different from those obtained for $A\to E_s$: all the three RS-ASMs (Fig. 11(b),(c) and (d)) have lower time cost and higher SNR than the standard ASM (Fig. 11(a)) because all the three up-sampling methods never lose energy. Figures 11(e),(f) show the time cost and SNR of all ASMs with increasing $\beta$. In the extreme case $\beta =120$, the time cost is 0.53 s (for the standard ASM), 0.14 s (for the “0-insert” RS-ASM), 0.15 s (for the “bilinear” RS-ASM) and 0.19 s (for the “bicubic” RS-ASM), and the corresponding SNR is 14.6, 17.7, 21.3 and 17.1. Similar to the simulation result (Fig. 4 and Fig. 6), the experimental data show that the “bilinear” RS-ASM exhibits the best overall performance in terms of SNR improvement and calculation acceleration.
5. Application example
In the previous sections, we demonstrated that the “bilinear” RS-ASM produces the highest improvement in the SNR and accelerates both the $A\to E_s$ and $E_s\to A$ procedures at $z_0$ compared to the standard ASM. In this section, we present a practical example based on diffraction calculation, the G-S phase retrieval algorithm [36,37], to evaluate the acceleration performance of the “bilinear” RS-ASM.
To apply the G-S iteration, we use a simulated complex virus image $E_s(z_0)=|E_0|e^{i\phi _0}$ with an amplitude (Fig. 12(c)) and phase distribution (Fig. 12(d)) at $z_0=0$ and an original size $N^2=128\times 128$ and a pixel size $\Delta _x=1.075\ \mu m$. To cover the large spread at diffraction plane $z_H=9.7\ mm$ with $NA=0.4$, we expand the original images to $8192\times 8192$ by placing “0-padding” around the images. With only two amplitude images of $E_s$ at $z_0$ and $z_H$ and a random phase image as input, we can apply a G-S iterative algorithm [37] to reconstruct the correct phase after several iterations. We use the “bilinear” RS-ASM with $\beta =32$ to speed up numerical propagation in the G-S algorithm and compare the results with those obtained using the standard ASM. To describe the reconstruction speed, we define the convergence as the $E_s$ difference between the $i^{th}$ (marked as $E_s^{i}$) and $(i-1)^{th}$ (marked as $E_s^{i-1}$) iteration steps: $\sqrt {\frac {1}{N^2}\sum (|E_s^{i}|-|E_s^{i-1}|)^2}$ at $z_0$, where $\sum$ means adding all pixels in a square box in which the signal is fully covered. We also define the phase fidelity to quantify the similarity between the iterated $E_s^{i}$ and the correct one $E_c$: $\frac {\langle E_s^{i}|E_c \rangle }{\sqrt {\langle E_c|E_c \rangle \langle E_s^{i}|E_s^{i} \rangle }}$ at $z_H$, where $\langle E_1|E_2 \rangle = \sum |E_1 E_2|$ represents for the inner product between two amplitudes $|E_1|$ and $|E_2|$. This phase fidelity definition is based on the assumption that any phase imperfections of $E_s^i$ at $z_0$ would lead to the amplitude deviation from the correct answer $E_c$ at diffraction plane $z_H$. Figure 12(a) and Fig. 12(b) show the convergence and phase fidelity versus iteration time, respectively. The “bilinear” RS-ASM converges faster than the standard ASM and always has higher phase fidelity for the same iteration time. As a quantitative example, to achieve a convergence of $=0.8\times 10^{-4}$, the standard ASM takes 128 s, whereas the “bilinear” RS-ASM only takes 50 s. To obtain a phase fidelity of $=98.3\%$, the standard ASM takes 68 s, whereas the “bilinear” RS-ASM only takes 54 s. The final phase images retrieved by the standard ASM and “bilinear” RS-ASM are shown in Fig. 12(e) and Fig. 12(f), respectively (the full size $E_s$ amplitude image and angular spectrum images are shown in Supplement 1 sections 5 and 6). This G-S example proves that the proposed “bilinear” RS-ASM can accelerate diffracting-calculation-related algorithms.
Fig. 12. G-S phase retrieval algorithm accelerated by “bilinear” RS-ASM. (a) Convergence versus iteration time. (b) Phase fidelity versus iteration time. Insets in (a,b) are partial enlarged drawing. The correct amplitude (c) and phase (d) distribution of a simulated virus image at $z_0$. Retrieved phase image by standard ASM (e) and “bilinear” RS-ASM (f).
6. Conclusion
In this study, we have developed simple and efficient numerical propagation methods, RS-ASMs, which are especially suitable for high-NA systems requiring a large image size $N^2$. To avoid significant calculation waste at the focal plane of a constrained object, we combine a resampling technique with interpolation methods in the angular spectrum domain. In addition to being simple and convenient to implement directly, RS-ASM is faster, less memory cost, and has a higher SNR performance than the standard ASM. RS-ASMs are validated using theory, simulations of a USAF target and a G-S phase retrieval procedure, and a cold-atom in-line holography experiment. Using the “bilinear” RS-ASM approach with an appropriate $\beta$ can result in a speed-up factor of up to 20x ($A\to E_s$) and 4x ($E_s \to A$), together with an SNR enhancement of approximately 2x due to noise suppression of the interpolation process. Considering an application example, the “bilinear” RS-ASM has 2.6x faster convergence speed than the standard ASM for the G-S phase reconstruction algorithm. We expect this acceleration and SNR enhancement to be more efficient for larger NA or longer propagating distances. RS-ASM can be applied to a wide range of fields, including CDI, digital holographic microscopy, and real-time holography monitoring.
Acknowledgments
The author gratefully acknowledges Saijun Wu lab for providing experimental support during which this work is conceived. The author also thanks Yami Fang, Yuzhuo Wang, and Yizun He for experimental contribution leading to the Fig. 8-11 application example of this algorithm.
Disclosures
The author declares no conflicts of interest.
Data availability
No data were generated or analyzed in the presented research.
Supplemental document
See Supplement 1 for supporting content.
References
1. D. Gabor, “A new microscopic principle,” Nature 161(4098), 777–778 (1948). [CrossRef]
2. M. K. Kim, “Principles and techniques of digital holographic microscopy,” SPIE Reviews 1(1), 018005 (2010). [CrossRef]
3. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21(15), 2758–2769 (1982). [CrossRef]
4. J. R. Fienup, “Reconstruction of an object from the modulus of its fourier transform,” Opt. Lett. 3(1), 27–29 (1978). [CrossRef]
5. J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of x-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature 400(6742), 342–344 (1999). [CrossRef]
6. J. Miao, D. Sayre, and H. Chapman, “Phase retrieval from the magnitude of the fourier transforms of nonperiodic objects,” J. Opt. Soc. Am. A 15(6), 1662–1669 (1998). [CrossRef]
7. T. Nishitsuji, T. Shimobaba, T. Kakue, N. Masuda, and T. Ito, “Fast calculation of computer-generated hologram using the circular symmetry of zone plates,” Opt. Express 20(25), 27496–27502 (2012). [CrossRef]
8. J. Weng, T. Shimobaba, N. Okada, H. Nakayama, M. Oikawa, N. Masuda, and T. Ito, “Generation of real-time large computer generated hologram using wavefront recording method,” Opt. Express 20(4), 4018–4023 (2012). [CrossRef]
9. W. Choi, C. Fang-Yen, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Extended depth of focus in tomographic phase microscopy using a propagation algorithm,” Opt. Lett. 33(2), 171–173 (2008). [CrossRef]
10. D. Wang, J. Zhao, F. Zhang, G. Pedrini, and W. Osten, “High-fidelity numerical realization of multiple-step fresnel propagation for the reconstruction of digital holograms,” Appl. Opt. 47(19), D12–D20 (2008). [CrossRef]
11. J. W. Goodman, Introduction to fourier optics. 3rd, Roberts and Company Publishers3 (2005).
12. M. Kanka, R. Riesenberg, P. Petruck, and C. Graulig, “High resolution (NA= 0.8) in lensless in-line holographic microscopy with glass sample carriers,” Opt. Lett. 36(18), 3651–3653 (2011). [CrossRef]
13. M. Kanka, R. Riesenberg, and H. Kreuzer, “Reconstruction of high-resolution holographic microscopic images,” Opt. Lett. 34(8), 1162–1164 (2009). [CrossRef]
14. M. Kanka, A. Wuttig, C. Graulig, and R. Riesenberg, “Fast exact scalar propagation for an in-line holographic microscopy on the diffraction limit,” Opt. Lett. 35(2), 217–219 (2010). [CrossRef]
15. T. Shimobaba, K. Matsushima, T. Kakue, N. Masuda, and T. Ito, “Scaled angular spectrum method,” Opt. Lett. 37(19), 4128–4130 (2012). [CrossRef]
16. T. Shimobaba, T. Kakue, M. Oikawa, N. Okada, Y. Endo, R. Hirayama, and T. Ito, “Nonuniform sampled scalar diffraction calculation using nonuniform fast fourier transform,” Opt. Lett. 38(23), 5130–5133 (2013). [CrossRef]
17. Y. Takaki and H. Ohzu, “Fast numerical reconstruction technique for high-resolution hybrid holographic microscopy,” Appl. Opt. 38(11), 2204–2211 (1999). [CrossRef]
18. J. Garcia-Sucerquia, D. Alvarez-Palacio, and H. Kreuzer, “High resolution talbot self-imaging applied to structural characterization of self-assembled monolayers of microspheres,” Appl. Opt. 47(26), 4723–4728 (2008). [CrossRef]
19. H. Kreuzer, K. Nakamura, A. Wierzbicki, H.-W. Fink, and H. Schmid, “Theory of the point source electron microscope,” Ultramicroscopy 45(3-4), 381–403 (1992). [CrossRef]
20. R. Xu, M. Feng, Z. Chen, J. Yang, D. Han, J. Xie, and F. Song, “Non-uniform angular spectrum method in a complex medium based on iteration,” Opt. Lett. 47(8), 1972–1975 (2022). [CrossRef]
21. W. Zhang, H. Zhang, and G. Jin, “Band-extended angular spectrum method for accurate diffraction calculation in a wide propagation range,” Opt. Lett. 45(6), 1543–1546 (2020). [CrossRef]
22. T. Kozacki and K. Falaggis, “Angular spectrum-based wave-propagation method with compact space bandwidth for large propagation distances,” Opt. Lett. 40(14), 3420–3423 (2015). [CrossRef]
23. M. Chlipała, T. Kozacki, H.-J. Yeom, J. Martinez-Carranza, R. Kukołowicz, J. Kim, J.-H. Yang, J. H. Choi, J.-E. Pi, and C.-S. Hwang, “Wide angle holographic video projection display,” Opt. Lett. 46(19), 4956–4959 (2021). [CrossRef]
24. W. Zhang, H. Zhang, and G. Jin, “Frequency sampling strategy for numerical diffraction calculations,” Opt. Express 28(26), 39916–39932 (2020). [CrossRef]
25. G. Koren, F. Polack, and D. Joyeux, “Iterative algorithms for twin-image elimination in in-line holography using finite-support constraints,” J. Opt. Soc. Am. A 10(3), 423–433 (1993). [CrossRef]
26. G. Koren, D. Joyeux, and F. Polack, “Twin-image elimination in in-line holography of finite-support complex objects,” Opt. Lett. 16(24), 1979–1981 (1991). [CrossRef]
27. T. Latychevskaia and H.-W. Fink, “Solution to the twin image problem in holography,” Phys. Rev. Lett. 98(23), 233901 (2007). [CrossRef]
28. L. Rong, Y. Li, S. Liu, W. Xiao, F. Pan, and D. Wang, “Iterative solution to twin image problem in in-line digital holography,” Opt. Laser. Eng. 51(5), 553–559 (2013). [CrossRef]
29. B. M. Ratliff, C. F. LaCasse, and J. S. Tyo, “Interpolation strategies for reducing ifov artifacts in microgrid polarimeter imagery,” Opt. Express 17(11), 9112–9125 (2009). [CrossRef]
30. R. Keys, “Cubic convolution interpolation for digital image processing,” IEEE Trans. Acoust., Speech, Signal Process. 29(6), 1153–1160 (1981). [CrossRef]
31. D. Schumacher, “General filtered image rescaling,” in Graphics Gems III (IBM Version), (Elsevier, 1992), pp. 8–16.
32. Bart Wronski, Bilinear texture filtering – artifacts, alternatives, and frequency domain analysis, (2020). https://bartwronski.com/2020/04/14/bilinear-texture-filtering-artifacts-alternatives-and-frequency-domain-analysis/.
33. S. Fadnavis, “Image interpolation techniques in digital image processing: an overview,” Int. J. Eng. Res. Appl. 4, 70–73 (2014).
34. M. Hu and J. Tan, “Adaptive osculatory rational interpolation for image processing,” J. Comput. Appl. Math. 195(1-2), 46–53 (2006). [CrossRef]
35. E. Heitz and F. Neyret, “High-performance by-example noise using a histogram-preserving blending operator,” Proc. ACM Comput. Graph. Interact. Tech. 1(2), 1–25 (2018). [CrossRef]
36. R. W. Gerchberg, “Phase determination for image and diffraction plane pictures in the electron microscope,” Optik 34, 275 (1971).
37. R. W. Gerchberg, “A practical algorithm for the determination of plane from image and diffraction pictures,” Optik 35, 237–246 (1972).
38. J. Sobol and S. Wu, “Imaging cold atoms with shot-noise and diffraction limited holography,” New J. Phys. 16(9), 093064 (2014). [CrossRef]
39. Y. Wang, J. Zhao, X. Huang, L. Qiu, L. Ji, Y. Ma, Y. He, J. P. Sobol, and S. Wu, “Imaging moving atoms by holographically reconstructing the dragged slow light,” Phys. Rev. Appl. 18(1), 014065 (2022). [CrossRef]
