1. Introduction

As Joseph W. Goodman says, “The phenomenon known as diffraction plays a role of utmost importance in the branches of physics and engineering that deal with wave propagation [1].” The wave equation describes the propagation of light through most media, and leads us to regard optical imaging as a specific mapping of object light into image light distributions. Advanced computational techniques along with the developed digital devices have promoted the contemporary computational imaging [2], which demands precise and high-speed numerical wave propagation that can work over a wide range of, distances, especially in situations that require numerous wave propagations, or where the target fields consist of a large amount of data, such as three-dimensional (3D) imaging, which requires volume diffraction [3,4] or multiple-slice beam propagation methods [5–8], computer-generated holograms (CGH) [9–11] for 3D [12] AR/VR [13] display and imaging with inverse problem solving [4,7,14].

Fast-Fourier-transform (FFT) convolution approach is widely used to implement the diffraction formulas for its speed. However, it induces sampling issues that affect the numerical accuracy of computing diffraction [15,16]. A lot of efforts have been made to achieve accurate propagation in a wide range of distances [17], such as intermediate plane techniques [18], band-limit angular spectrum method (ASM) to avoid aliasing [19], padding to expand the band-width [17,20,21], adaptive-sampling to utilize the space-bandwidth product while reducing the computational complexity [22], non-uniform FFT to achieve fast calculating [22,23]. Most methods focus on sampling issues [24–26].

This paper improves the Rayleigh-Sommerfeld diffraction computation without concentrating on the sampling issue. Instead, by Taylor expanding the exponential term in the diffraction convolution kernel, the angular spectrum diffraction formula can be approximated such that a similar numerical accuracy is attainable with fewer bits formats, e.g., single-precision floating-point data type. This cuts off the required computational memory by half compared to the usual double-precision floating-point data types, and significantly reducing run-time consequentially.

2. Theory

We consider free-space propagation of a fully coherent monochromatic scalar wave field, in the spatial domain. Let $\lambda$ be the wavelength of the scalar field, and $k = 2\pi /\lambda$ be the wavenumber. In a 3D Cartesian coordinate system $(x, y, z)$, we consider wave propagation along the $z$-axis, and the transverse 2D plane $\boldsymbol {r}_{\perp } = (x, y)$ is of our interest. $|\cdot |$ denotes vector lengths in $\ell _2$. We are interested in how an unbounded initial field $u(\boldsymbol {r}_{\perp }, 0)$ evolves into a new field $u(\boldsymbol {r}_{\perp }, z)$ [1], after a propagation distance of $z$:

Except for the sampling issues in the numerical diffraction implementation, we found that accuracy depends on floating-point precision. Because the visible light’s wavelength is hundreds of nano-meter and the sensor pixel pitch is at the micro-meter scale, the terms $\sqrt {1-\lambda ^{2}|\boldsymbol {\rho }_{\perp }|^{2}}$ in Eq. (4) and $\exp (\text {j} k \boldsymbol {r})$ in Eq. (3) cannot be accurately represented by a single-precision floating point number, i.e., the conventional numerical diffraction requires double-precision representations. This limitation is demonstrated in Fig. 1. This limitation presents an obstacle in particular for implementations on graphics processing units (GPUs), which have both limited memory resources, and are significantly faster at processing single vs. double precision numbers. Single precision arithmetic is commonly used in deep learning frameworks [29] and physics-based neural network training [3], but naïvely using single precision for diffraction simulations leads to inaccurate results due to numerical errors in the phase terms.

 

Fig. 1. (a) Computed PSFs and (b) CTFs from Eq. (3) and Eq. (4) using different computational approaches. The Naïve way computes the kernels by brute force computing the square root, resulting in large numerical errors. In contrast, Taylor expanding the square root leads to smaller numerical errors, and further converging to the ground truth with more Taylor terms, here demonstrated by the first two terms.

Download Full Size | PDF

In the following, we show the proposed Hybrid Taylor Rayleigh-Sommerfeld diffraction (HTRSD) method. With a proper Taylor expansion and hybrid FFT-DI and ASM, we can achieve precise numerical wave propagation with a single-precision float number. We also show the time and memory efficiency. For convenience, we name ASM-single and ASM-double as the conventional band-limited ASM [19] method work in single- and double-precision float format, respectively.

For the PSF in Eq. (3), direct sampling may lead to severe numerical rounding errors caused by $\exp (\text {j} k \boldsymbol {r})$, in a general setting when $k \gg 1/z$, we can approximate $r$ using Taylor series, that is valid when $|\boldsymbol {r}_{\perp }| < z$:

As mentioned in previous paragraph, the double-precision float number requirement leads to significant time and memory consumption, thus limiting the numerical diffraction’s applications. Here we show that with the Taylor expansion, we can eliminate the double-precision float number requirement in the Rayleigh-Sommerfeld diffraction. Figure 1 demonstrates the limitations when computing the PSF and CTF naïvely and with the Taylor expansion in single-precision float number. Figure 1(a) shows the real parts of $\exp (\text {j} (k\boldsymbol {r} - kz) )$ with respect to $|\boldsymbol {r}_{\perp }|/z$, and Fig. 1(b) shows the real parts of $\exp \small (\text {j} kz (\sqrt {1 - \lambda ^{2} |\boldsymbol {\rho }_{\perp }|^{2}} - 1) \small )$ with respect to $\lambda |\boldsymbol {\rho }_{\perp }|$. The difference (red lines) between the naïve and the Taylor expansion show that numerical precision could decrease when using single floating points for computing the PSF and CTF via Eq. (3) and Eq. (4). A naïve approach suffers from severe numerical rounding errors (first rows). In contrast, using Taylor series and by approximating the phase terms, converges to the ground truths computed using double precision accuracy (second rows). Approximation accuracy further increases when adding more terms (third rows). As for the IEEE floating point standard, for single precision the smallest spacing is $1.1921 \times 10^{-7}$ and for double precision it is up to $2.2204 \times 10^{-16}$ [30]. The accuracy could not beyond this.

Since the ASM is not suitable for long distance propagation [26,31], combine the ASM and FFT-DI with the conditions of the Taylor expansion, the finial diffraction formula is:

3. Numerical verification

Following the band-limited ASM, in the implementation, we pad the original field with zeros for two times and crop the field to the original size after propagation, as shown in Fig. 2. We implemented the code with PyTorch 1.10 [29]. For generalization, we set the input field as a volume of size $N_x \times N_y \times N_z$. The CTF or PSF are the same size as the input field. This can be easily applied to both 2D and 3D wave propagation by setting $N_z=1$ or $N_z > 1$. We set the Taylor expansion order as $n=10$. The following numerical experiments were conducted on a workstation with two Intel Xeon E5-2690 2.60GHz CPUs.

 

Fig. 2. Schematic for Computing optical wave free-space propagation. To avoid boundary artifacts, the object plane is padded, and the detection plane is cropped after propagation.

Download Full Size | PDF

The proposed method is verified on circular and square apertures with analytical expressions. For a circular aperture of radius $a$, at different propagation distance $z$, under normal uniform plane-wave illumination [1], the on-axis scalar fields on the optical axis can be expressed exactly [32,33]:

 

Fig. 3. Magnitude and phase profile vs. propagation distance.

Download Full Size | PDF

The on-axis field can not reflect the performance of the whole field propagation since we usually detect the 2D fields perpendicular to the optical axis in practice. Here we perform the comparison of circular and rectangle apertures. In the following, we set the pixel size of the original field as $N_x = N_y = 255$, the width and length are $L_x = L_y = {0.5}\;\textrm {m}$, $a=0.125$, and $\lambda ={500}\;\textrm {nm}$ (same parameters as in Chapter 5 of Ref. [31]). We measured the diffracted fields at $z=\{1000,2000,4000,20000\}\textrm {m}$ respectively, which correspond to Fresnel number $N_f =\frac {a^{2}}{\lambda z} = \{31.25,15.625,7.813,1.562\}$. This setting ensures that the measured diffraction images are within the Fresnel region. Figure 4 shows the comparison of circular aperture diffraction. From both the 2D images in Fig. 4(a) and the plots of the crossed center line in Fig. 4(b), we can observe that ASM-single method doesn’t work at all, and on the contrary, ASM-double and the HTRSD-single work well. The diffraction images of a rectangle aperture in Fig. 5 shows the similar results as in the circular aperture diffraction case. In this case, we also provide the plot of the RS calculation with numerical integral (RSI) method in Fig. 5(b). The proposed HTRSD method shows a comparable results to both the RSI and ASM with double precision. From both the tests on circular and rectangle apertures, we can say that the proposed HTRSD-single can achieve comparable results as the conventional ASM.

 

Fig. 4. Intensity (a) and profile (b) comparison on a circular aperture.

Download Full Size | PDF

 

Fig. 5. Intensity (a) and (b) profile comparison on a rectangle aperture.

Download Full Size | PDF

The previous paragraphs show that the proposed HTRSD can achieve even better accuracy in a single-precision float number format than the conventional band-limited ASM method with a double-precision float number. Here we show the time and memory cost comparison. To compare the time cost, we propagate a 2D complex field of pixel size $1024 \times 1024$ for various times with the ASM in double float precision and the proposed method in single float precision. Figure 6(a) shows the time cost with respect to the propagation numbers. It shows that the time cost is linearly proportional to the propagation number in both cases, while the conventional method costs 1.5 times over the proposed method. This is not surprising since both GPU and CPU can compute higher number operations if numbers have less precision. Notably, calculation of the Fourier transform of the PSF in HTRSD-single was included in the comparison. Even so, the speed of HTRSD-single is still much faster than the conventional ASM in double precision. For the memory costs, the target object is a volume with the same lateral size as in the speed test but with varying depth numbers. Similar to the time cost, the memory cost is also linearly proportional to the target volume size, while the conventional method costs double memory to the proposed method. This is because the double float is a 64-bit point number, and the single float is 32-bit, so a single float number uses half of the memory as the same size double float number uses. It should be noted that the HTRSD might perform more accurately if double float precision is used, however this may come at the expense of speed and memory usage.

 

Fig. 6. Time (a) and memory (b) costs of the conventional ASM and proposed HTRSD.

Download Full Size | PDF

4. Conclusion

We present a numerical diffraction technique, HTRSD, that can achieve more accurate Rayleigh-Sommerfeld diffraction at a much higher speed but consume dramatically less computational resources as well as memory. This could promote the development of contemporary computational optics, especially computational imaging that includes optimization-based inverse imaging and extensive volume diffraction calculation in both optical imaging and displays. This work could also help computational imaging techniques that employ GPUs due to the power of GPU computing on single-precision floating points [34].