Is the Rayleigh-Sommerfeld diffraction always an exact reference for high speed diffraction algorithms?

In several areas of optics and photonics the behavior of the electromagnetic waves has to be calculated with the scalar theory of diffraction by computational methods. Many of these high-speed diffraction algorithms based on a fast-Fourier-transformation are approximations of the Rayleigh-Sommerfeld-diffraction (RSD) theory. In this article a novel sampling condition for the well-sampling of the Riemann integral of the RSD is demonstrated, the fundamental restrictions due to this condition are discussed, it will be demonstrated that the restrictions are completely removed by a sampling below the Abbe resolution limit and a very general unified approach for applying the RSD outside its sampling domain is given. © 2017 Optical Society of America under the terms of the OSA Open Access Publishing Agreement OCIS codes: (050.1960) Diffraction theory; (110.1758) Computational imaging; (110.1650) Coherence imaging; (090.1995) Digital holography. References and links 1. Sommerfeld, Optics, Lectures on Theoretical Physics, Vol. IV, New York (Academic, 1954). 2. J. W. Goodman, Introduction to Fourier Optics, 4th Ed., (W. H. Freeman, 2017). 3. M. Born and E. Wolf, Principles of Optics, New York (Cambridge University, 2005). 4. Y. M. Engelberg and S. Ruschin, “Fast method for physical optics propagation of high-numerical-aperture beams,” J. Opt. Soc. Am. A 21(11), 2135–2145 (2004). 5. V. Nascov and P. C. Logofătu, “Fast Computation algorithm for the Rayleigh-Sommerfeld diffraction formula using a type of scaled convolution,” Appl. Opt. 48(22), 4310–4319 (2009). 6. F. Shen and A. Wang, “Fast-Fourier-transform based numerical integration method for the Rayleigh-Sommerfeld diffraction formula,” Appl. Opt. 45(6), 1102–1110 (2006). 7. F. Zhang, G. Pedrini, and W. Osten, “Reconstruction algorithm for high-numerical-aperture holograms with diffraction-limited resolution,” Opt. Lett. 31(11), 1633–1635 (2006). 8. H. Scheffers, “Vereinfachte Ableitung der Formeln für die Fraunhoferschen Beugungserscheinungen,” Ann. Phys. 434, 211 (1942). 9. E. Lalor, “Conditions for the Sampling of the Angular Spectrum of Plane Waves,” J. Opt. Soc. Am. 58(9), 1235– 1237 (1968). 10. K. Matsushima, “Shifted angular spectrum method for off-axis numerical propagation,” Opt. Express 18(17), 18453–18463 (2010). 11. A. Ritter, “Modified shifted angular spectrum method for numerical propagation at reduced spatial sampling rates,” Opt. Express 22(21), 26265–26276 (2014). 12. T. C. Poon and J. P. Liu, Introduction to Modern Digital Holography with MATLAB (Cambridge University, 2014). 13. P. Picart and J. C. Li, Digital Holography (Wiley, 2012). 14. M.-K. Kim, Digital holographic Microscopy (Springer, 2011). 15. E. Wolf and E. W. Marchand, “Comparison of the Kirchhoff and the Rayleigh–Sommerfeld Theories of Diffraction at an Aperture,” J. Opt. Soc. Am. 54(5), 587–594 (1964). 16. M. W. Farn and J. W. Goodman, “Comparison of Rayleigh-Sommerfeld and Fresnel solutions for axial points,” J. Opt. Soc. Am. A 7(5), 948–950 (1990). 17. R. J. Marks II, The Joy of Fourier (Baylor University, 2006). 18. H. Gross, Handbook of Optical Systems, Volume 1: Fundamentals of Technical Optics (Wiley-VCH, 2005). 19. E. Kreyszig, Introductory Functional Analysis with Applications (John Wiley & Sons, 1989). Vol. 25, No. 24 | 27 Nov 2017 | OPTICS EXPRESS 30229 #306576 https://doi.org/10.1364/OE.25.030229 Journal © 2017 Received 7 Sep 2017; revised 7 Nov 2017; accepted 10 Nov 2017; published 17 Nov 2017


Introduction
The scalar theory of diffraction can be used in many different applications like wave propagation, digital holography, holographic microscopy, diffraction imaging, biomedical imaging and diffractive optics. One exact method for the scalar diffraction calculation is the Rayleigh-Sommerfeld diffraction (RSD). In contrast to approximations such as Fresnel or Fraunhofer diffraction, the RSD gives an exact solution for the output field of a given input field [1][2][3]. However, to the best of our knowledge there is no general analytical solution for the calculation of the RSD. Therefore, numerical methods have to be used. All numerical simulations are in principle based on a sampling of the analogue continuous field. Thus, with usual computational power, only high-speed algorithms make the utilization of the diffraction theory possible. These high-speed algorithms use approximations of the RSD integral such as a quadratic phase [2,4] or a frequency-cut in convolution RSD [5][6][7] and in the angular spectrum method (ASM) [6,[8][9][10][11]. Whereupon the latter is only valid for small propagation distances [2,12,13]. Another restriction in the ASM is the identical sampling spacing in the input and output plane [14]. In order to solve this problem some methods are developed. One proposed method uses interpolation [5], however, this causes a new numerical error dependent of the interpolation method. In another one, the object extension is limited [7]. Contrary to the Kirchhoff solution of the diffraction problem [2,3,15], the mathematical solution of the RSD is not inconsistent when the observation point is close to the diffracting screen. Thus, the accuracy of the high-speed alternatives, even for very short propagation distances, could be verified by referencing them with exact solutions given by the RSD. Since there exist only a few analytical solutions [16], they have to be compared with a discretized RSD. In this work, the RSD is treated as a Riemann integral, which has to be discretized. However, here we will demonstrate that the use of the discretized RSD can cause enormous aliasing. We will give the boundary conditions for the usage of the discretized RSD and show possibilities to completely remove the aliasing. Since in this paper we are presenting a method for the computer-based calculation of the Rayleigh-Sommerfeld diffraction integral without aliasing, throughout the paper we assume the object and all other images as already discretized. Figure 1 shows a typical setup for the calculation of the diffraction problem. The diffraction of a discretized input plane 1, generated by a coherent illumination of an analogue object, is reproduced by a discretized sensor like a CCD camera in plane 2 or 3.

Sampling condition for the Rayleigh-Sommerfeld-diffraction
According to the Rayleigh-Sommerfeld diffraction integral, the field distribution in an output plane 2 in the distance 12 z , parallel to the input plane is [1][2][3]:  is the field distribution in the input (object) plane 1 and k is the wave number.
The two vectors 1 r  , 2 r  are position vectors in plane 1 or 2, respectively and Ω is the integration area. Please note that, since all computer algorithms can only work with discretized data, although in real world an object would be analogue, for the computational calculation the input in plane 1 has to be discretized. Sometimes the object can be defined analytically by a plane wave or a finite chirp function. In such cases the requirement for the numerical treatment is discretizing all terms of the diffraction method and the object itself considering the Nyquist criterion [2]. However, in most cases, the image is discretized by a CCD-Camera. Consequently, for a numerical treatment of Eq. (1), only the propagation dependent harmonic term 2 1 exp( ) ik R r −   has to be sampled, according to the well-known Nyquist sampling criterion [17]. The sampling frequency must at least correspond to twice the highest frequency contained in the harmonic term. If we consider the phase of the propagation term as: ( ) with the transversal Cartesian coordinates in the input ( ) This spatial frequency is a monotonically increasing function of 1  x and 2 fp x are the distances of the farthest point from the center in x -direction in the relevant computational plane 1 or 2 with non-negligible amplitude value (see Fig. 1). Due to zero values of the amplitude around the boundary, the maximum width can be smaller than the real width of the computational domain. Thus, it follows for the maximum spatial frequency ,max The sampling frequency s f is related to the sampling spacing via Thus, it follows for the sampling according to the Nyquist criterion ( ) Therefore, the sampling condition for the numerical treatment of the RSD can be written as: 1 .
fp fp Consequently, for a given sampling spacing 1 x δ , 1 y δ in the object plane and a maximum size 1 fp x , 1 fp y of the object and its image 2 fp x , 2 fp y in the output computational plane, there is a critical minimum propagation distance 1 z allowed, in which the output plane can be calculated by: An analog derivation results in a similar condition for the y -direction: Thus, the condition for the minimum propagation distance (critical distance) is: The critical distance is the minimum distance in which an output field (image) can be numerically calculated by RSD without violating the Nyquist criterion for the interplay of the sampling conditions of input and output planes and the distance. The reconstruction of the object from the diffracted image is only possible, if the Nyquist theorem is fulfilled in the reverse direction too. This results in equations analogous to Eqs. (7)-(9) with reversed index 1 and 2.
Thus, the total critical distance c z for a forward and reverse transformation of the field has to be at least ( ) 1 2

Aliasing if
In this subsection, the aliasing due to the violation of the aforementioned sampling condition is shown. However, the discussion of aliasing is more appropriate and required in methods dealing with the reduction of the aliasing or methods, which investigate the domains with strong or weak aliasing. In these cases the aliasing can provide invaluable information about the advantage and capability of the developed methods. However, in section 4 a method for preventing the occurrence of aliasing in RSD is shown. Thus, aliasing is not discussed indepth, but is used to compare the RSD method suffering from aliasing with our aliasing-free proposed approach.
To show the aliasing, we have used an amplitude object (constant phase) as described in Figs. 2(a) and 2(b) and calculated its diffraction pattern at a distance c z z < . For the sake of simplicity but without loss of generality, the phase distributions are compared indirectly. Thus, the magnitude and real part of the complex amplitude will be discussed throughout the paper. As can be seen in Fig. 2, for this object both values are identical because it has a zero phase. The color-bar shows the normalized amplitude and real part values. The simulated distance between the object and the image plane is 12 20μm z = , which according to Eqs. (7)-(9), is much smaller than the critical distance 534μm c z = . The calculated amplitude and real part of the diffracted field 12 u at a distance 12 z can be seen in Figs. 3(a) and 3(b) respectively. The error due to the violation of the sampling condition can be seen by a reconstruction of the original object from the diffracted image 12 u . Thus, the reverse RSD (for 12 z z = − ) has to be applied. Here the term "reverse" instead of "inverse" transform will be used in order to avoid confusion.
The result can be seen in Figs. 3(c) and 3(d). A similar pattern like in the original object occurs but, with very strong aliasing. Especially the nonzero values of the field in areas, which originally exhibit zero values, lead to completely wrong results. The correlation coefficient r between 1 u and 121 u , i.e. the field after transforming from position 1 to 2 and back to 1 is 0.72 r = and 0.75 r = for the magnitude of the complex field and its real part, respectively. The obviously strong deviation to the input is a consequence of the insufficiency of the propagation distance is satisfied for the numerical calculation of the RSD, there will be no loss of information due to the transformation. For an object with a fixed sampling spacing the critical propagation distance is fixed. Accordingly, if the propagation distance is longer than c z , a correct treatment of the computational RSD should be achieved for the given sampling spacing.
In Figs. 4(a) and 4(b) the diffracted image is numerically calculated by the RSD under the assumption that the propagation distance 13 730μm z = satisfies the sampling condition. To confirm the correctness of the field 13 u in the plane 3 as a necessary condition, the reverse transform is considered to reconstruct the input object in the input plane, as reported in Figs. 4(c) and 4(d). The correlation coefficient is 0.97 r = for both, magnitude and real part of the complex amplitude. Since a small part of the whole information will be lost by the spatial limitation of the computational plane (The numerical aperture is smaller than one).
Comparing the correlation coefficients for the reconstructed object in Fig. 3 and Fig. 4 shows a more than 20% improvement. Therefore, at the distance 13 z almost the whole information of the object plane is preserved. However, in some cases the RSD might be applied as a reference for different algorithms at a propagation distance outside of the sampling domain. Additionally, a good object reconstruction by a forward and reverse calculation is just a necessary, but not a sufficient condition for testing the sampling of a diffraction algorithm. It does not necessarily mean, that the output corresponds to the expected result of the diffraction theory. In other words, a combination of an arbitrary propagation operator (physically or nonphysically) with its inverse always results in an identity operator, and consequently the reconstruction of the input is expected automatically. Thus, in the next subsection a sampling condition, which always satisfies the sampling condition and which can be used as a reference for the diffraction will be presented and in section 4 a general procedure, which makes the RSD a feasible method for arbitrary propagation distances will be shown.

Sampling spacing below the Abbe resolution limit for fine structures larger than the Abbe limit
According to inequality 6, the left side and the second term on the right side are always positive whereas, the first term on the right side can change its sign. For a sampling lower than half of the wavelength 1 / 2 x δ λ < , it will become negative and consequently the inequality will be fulfilled for all propagation distances. z . Thus, the sampling condition is always satisfied, if the sampling spacing of the harmonic term is smaller than the Abbe resolution limit. It should be emphasized that this condition only holds for the harmonic term. As will be shown in subsection 4, this does not contradict the Abbe resolution limit.
According to the discussion above, substructures smaller than the Abbe limit could be resolved and consequently be reconstructed by the Nyquist criterion. However, there are different meanings of "reconstruction" in respect to diffraction (restricted due to the Abbe limit) and in respect to the sampling theory, restricted by the Nyquist criterion. In the context of the sampling theory, a direct reconstruction of the original field after sampling will be possible, whereas for the diffraction theory the reconstruction of the original field is indirect, since it takes place after a propagation over the distance z . Thus, the sampled data is exposed to the diffraction effect and consequently restricted by the Abbe limit. Eventually, a sampling below the wavelength does not lead to the breaking of the Abbe rule for our approach, but it enables the calculation of a diffracted image without violating the RSD sampling condition. , except that for this simulation the object was sampled with a sampling spacing below the Abbe limit. However, the structures in the object are still larger than the Abbe limit. The reconstructed object is shown in Figs. 5(c) and 5(d). The correlation coefficient between the reconstructed and input object is 0.997 r = for both the magnitude and the real part of the complex amplitude. The minor loss of information is just due to the limited aperture. Thus, the sampling of the harmonic term with a sampling spacing smaller than the Abbe limit can be used as a reference for the evaluation of the quality of the numerical calculations.

Sampling spacing below the Abbe resolution limit for fine structures smaller than the Abbe limit
To investigate the effect of the sampling below the Abbe limit for under Abbe limit structures, the input area, the output area and the sampling spacing have been rescaled (10 and 20 times smaller than the object in Fig. 2), so that both the object's fine structures and the sampling spacing are below the Abbe limit. In Figs. 6(a)-(d) the rescaled object and the reconstructed object are shown respectively. 10.

F N =
As can be seen, due to the violation of the Abbe resolution limit, the fine structures of the object in Fig. 6 cannot be resolved anymore. The calculated correlation coefficients are 0.88 r = and 0.89 r = for the magnitude and the real part of the complex amplitude respectively.  field 12 u at the same distance. This field 12 u was calculated for a sampling spacing below the Abbe limit and can be used as a reference, as discussed in section 3-3. In Fig. 9 the reconstruction of the object 13231 u by the use of the operator 132  for the forward and 231  for the backward propagation is presented. The correlation coefficients for the magnitude and the real part of the complex amplitude are 0.97 r = . Again, the capability of the proposed approach can be seen.

Conclusion
In this paper the numerical treatment of the Rayleigh-Sommerfeld diffraction in its Riemannian integral form was investigated in detail. A sampling condition for the numerical calculation was derived. As have been shown, for a fixed sampling spacing, widths of the input and output and wavelength the allowed propagation distance is restricted to a minimum value. However, the restriction can be completely removed if the sampling spacing (not the structure in the object) is lower than the Abbe limit. As have been shown, this results in the maximum obtainable information in the output plane under the consideration of the limited computational domain, and was therefore used as a reference. Moreover, a very general approach for the calculation of the output field for arbitrary propagation distances was presented. This operator is based on a combination of forward and reverse RSD transforms and leads to very high correlation coefficients of 0.97 r = . An about 30% improvement in the magnitude of the complex amplitude and about 45% for the real part confirms the reliability of the new operator. A comparison of the results of the below Abbe limit sampling with the results of the composed operator is an additional verification of both methods and the consistency of the theoretically derived sampling condition for the RSD. The developed approach can be used as a reference for the testing of high-speed algorithms and other methods, which are based on the approximation of the exact scalar diffraction theory.
One of the most important advantages of this method over other methods is that it doesn't need any changes in the sampling spacing 1 x δ and 2 x δ to prevent aliasing.