1. 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–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–7] and in the angular spectrum method (ASM) [6,8–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.

2. Sampling condition for the Rayleigh-Sommerfeld-diffraction

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.

 

Fig. 1 Coherent imaging of an object in the input plane (1). The second and third planes refer to the diffracted images of the object at the distance $z_{12}$ and$z_{13}$ from the object plane, respectively.

Download Full Size | PDF

According to the Rayleigh-Sommerfeld diffraction integral, the field distribution in an output plane 2 in the distance$z_{12}$, parallel to the input plane is [1–3]:

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 $\exp (ik\left|{\overrightarrow{R}}_2-{\overrightarrow{r}}_1\right|)$ 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:

This spatial frequency is a monotonically increasing function of $\left|x_1-x_2\right|$. To get its maximum value the conditions $\left|y_1-y_2\right|=0$ and $\left|x_1-x_2\right|=x_{1fp}+x_{2fp}$ must be fulfilled. Here $x_{1fp}$ and$x_{2fp}$ 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 $f_{x,\max }$ in $x$-direction:

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 $z_c$ for a forward and reverse transformation of the field has to be at least$z_c=\max \left(z_{1c},z_{2c}\right)$, which is the proposed sampling condition for a numerical treatment of the RSD.

According to the sampling condition, the correct calculation of the diffraction pattern and also the reconstruction of the original object from the achieved pattern may be performed if the propagation distance is longer than the critical distance:

3. Investigation of the sampling condition

3.1 Aliasing if $z<z_c$

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 in-depth, 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$z<z_c$. 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.

 

Fig. 2 Input object plane, (a) magnitude and real part (b) of the complex amplitude. If not otherwise stated, the following parameters were used for the presented simulations: Wavelength $\lambda =0.633\text{μm}$, sampling spacing in the input plane $\delta x_1=\delta y_1=0.76\text{μm}$ and in the output plane $\delta x_2=\delta y_2=0.94\text{μm}$. The pixel numbers in the$x$ and$y$ axis as well as in the input and output plane are the same $N_{1,2,x,y}=265$. The width of the computational domain in the input plane is $P_{x_1}=P_{y_1}=201\text{μm}$, whereas for the output plane it is $P_{x_2}=P_{y_2}=250\text{μm}$.

Download Full Size | PDF

The simulated distance between the object and the image plane is $z_{12}=20\text{μm}$, which according to Eqs. (7)-(9), is much smaller than the critical distance $z_c=534\text{μm}$.

The calculated amplitude and real part of the diffracted field $u_{12}$ at a distance $z_{12}$ 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 $u_{12}$. Thus, the reverse RSD (for $z=-z_{12}$) has to be applied. Here the term “reverse” instead of “inverse” transform will be used in order to avoid confusion.

 

Fig. 3 Reconstruction errors for $z<z_c$ (a) magnitude and (b) real part of the complex amplitude of the diffracted image in plane 2. (c) and (d) magnitude and real part of the reconstructed object in plane 1.$N_F=1000.$

Download Full Size | PDF

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 $u_1$ and $u_{121}$, i.e. the field after transforming from position 1 to 2 and back to 1 is $r=0.72$ and $r=0.75$ 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 $z<z_c$. The Fresnel number$N_F=P_1{P}_2/\left(4\lambda z_{12}\right)$ [18] for the forward and backward propagation between input and output plane is $1000$.

3.2 Propagation distance $z>z_c$

As derived in section 1, if the condition $z_{12}>z_c$ 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 $z_c$, 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 $z_{13}=730\text{μm}$ satisfies the sampling condition. To confirm the correctness of the field $u_{13}$ 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$r=0.97$ 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).

 

Fig. 4 (a), (b) magnitude and real part of the image complex amplitude, respectively for $z_{13}>z_c$. (c) and (d) corresponding reconstructed amplitude and real part of the reconstructed object in the input plane.$N_F=27.$

Download Full Size | PDF

Comparing the correlation coefficients for the reconstructed object in Fig. 3 and Fig. 4 shows a more than $20\%$ improvement. Therefore, at the distance $z_{13}$ 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 non-physically) 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.

3.3 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 $\delta x_1<\lambda /2$, 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.

In Figs. 5(a) and 5(b) the diffracted image for the same simulation parameters like in Fig. 3 can be seen $\left(z_{12}<z_c\right)$, 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 $r=0.997$ 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.

 

Fig. 5 (a) and (b) magnitude and real part of the complex amplitude according to the RSD at the same propagation distance like Fig. 3 but with a sampling spacing of the object and image smaller than the Abbe limit. Figure 5. (c) and (d) reconstruction of the object. The structures in the object are larger than the Abbe limit. $N_F=1000.$

Download Full Size | PDF

3.4 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.

 

Fig. 6 (a), (b) magnitude and the real part of the complex amplitude for the rescaled input object. (c) and (d) magnitude and real part of the complex amplitude for the reconstructed object for a sampling spacing below the Abbe limit $\delta x_1=\delta y_1=0.025\text{μm}\text{<}\lambda \text{/2}\text{=}\text{0}\text{.32}\text{μm}$. $N_F=10.$

Download Full Size | PDF

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 $r=0.88$ and $r=0.89$ for the magnitude and the real part of the complex amplitude respectively.

For a further reduction of the fine structures in the object the effect is increased as can be seen in Figs. 7(a)-7(d). The calculated correlation coefficients are only $r=0.62$ and $r=0.63$ Thus, a subwavelength sampling in RSD cannot retain the full object information, if the fine substructures are below the Abbe limit.

 

Fig. 7 (a) and (b) magnitude and real part of the complex amplitude of the input object. (c) and (d) magnitude and real part of the complex amplitude for the corresponding output by a sampling below the Abbe limit $\delta x_1=\delta y_1=0.013\text{μm}\ll \lambda \text{/2}\text{=}\text{0}\text{.32}\text{μm}$. $N_F=\mathrm{2.5.}$

Download Full Size | PDF

4. General solution for removing the limitation of the propagation distance

In practical applications the sampling of the object is restricted by a minimum spacing as a consequence of the limited pixel size of a given CCD camera. Here a general solution for an arbitrary distance from the object will be presented.

The RSD is a linear operator $\mathcal{R}$, which transforms an input field $u_1\left(\overrightarrow{r}\right)$ over a propagation distance $z_{12}$ to the field $u_{12}={\mathcal{R}}_{12}\left\{u_1\right\}$. Theoretically, for an unlimited aperture the information in the input plane $u_1$ is completely conserved in the diffracted image $u_{12}$. Therefore, the field $u_1$ can be reconstructed from the field $u_{12}$ by a reverse application of the RSD with $z\to -z$. The combination of the forward operator ${\mathcal{R}}_{12}$ and the reverse operator ${\mathcal{R}}_{21}$ $\left(u_1={\mathcal{R}}_{21}\left\{u_{12}\right\}\right)$ of the field is an identical operator ${\mathcal{R}}_{21}{\mathcal{R}}_{12}=\mathcal{I}$. Thus, it can be written that ${\mathcal{R}}_{21}{\mathcal{R}}_{12}\left\{u_1\right\}=\mathcal{I}\left\{u_1\right\}=u_1$ (for$z_{12}>z_{{12}_c}$).

If the sampling condition is not satisfied $z_{12}<z_{{12}_c}$, a complete reconstruction is not possible and${\mathcal{R}}_{21}{\mathcal{R}}_{12}\ne \mathcal{I}$. If a set $\Upsilon$ of all propagation distances satisfying the sampling condition is introduced, it follows that the reverse transform is not an inverse transform if $z\notin \Upsilon$. Although analytically the reverse and the inverse transforms are identical. Thus, a perfect reconstruction of all the information in the object is only possible if ${\mathcal{R}}_{21}{\mathcal{R}}_{12}=\mathcal{I}$. Therefore, an RSD operator, which satisfies ${\mathcal{R}}_{21}{\mathcal{R}}_{12}=\mathcal{I}$ for $z\notin \Upsilon$ has to be found.

As described in the last section, at the distance $z_{13}\in \Upsilon$, the reconstruction of the object is almost perfect but, outside the sampling condition $z_{12}\notin \Upsilon$, the whole object information cannot be retrieved from the field $u_{12}$. Thus, for a general solution, the following approach is proposed: in a first step the image at a longer distance which satisfies the sampling condition and identity relation $z_{13}\in \Upsilon$, is calculated with the additional property $z_{13}-z_{12}\in \Upsilon$. In a second step a new propagation distance $z_{23}=z_{13}-z_{12}$ will be calculated with ${\mathcal{R}}_{23}={\mathcal{R}}_{32}^{-1}$. The operator ${\mathcal{R}}_{32}$ transforms the field $u_{13}$ at $z_{13}$ to the field $u_{132}$ at a shorter distance $z_{12}$. Thus, the operator ${\mathcal{R}}_{132}={\mathcal{R}}_{32}{\mathcal{R}}_{13}$, which transforms the field $u_1$ to the field $u_{132}={\mathcal{R}}_{132}\left\{u_1\right\}$ at the distance $z_{12}$ is introduced. Although $z_{12}\notin \Upsilon$, it can be easily shown that ${\mathcal{R}}_{231}$ satisfies the identity relation as follows:

The reverse of the operator ${\mathcal{R}}_{132}$ is the operator ${\mathcal{R}}_{231}$. According to operator theory [19]:

Which means the reverse and inverse transforms are the same ${\mathcal{R}}_{231}={\mathcal{R}}_{132}^{-1}$. If the loss of information due to the limited aperture for practical applications is neglected, the new image $u_{132}$ contains all information from $u_1$. The operator ${\mathcal{R}}_{132}$ depends on two propagation variables $z_{12}$ and $z_{13}$. The first is the real variable, which determines the distance between the object and the image. The second is just an arbitrary parameter, which has to fulfill the condition $z_{13}\in \Upsilon$, $z_{13}-z_{12}\in \Upsilon$. Thus, the set $\Upsilon$ has an infinite number of elements, which are all valid. However, a cutting of diffracted field values due to the limited size of the computational plane 3 leads to a loss of information. Thus, for a fixed value of the pixel size and pixel number, the optimal choice for the propagation distance $z_{13}$ is the minimum allowed value. In Fig. 8 the calculated field $u_{132}$ at the distance $z_{12}\notin \Upsilon$ is compared with the field $u_{12}$ at the same distance. This field $u_{12}$ was calculated for a sampling spacing below the Abbe limit and can be used as a reference, as discussed in section 3-3.

 

Fig. 8 (a) and (b) magnitude and real part of the complex amplitude $u_{132}$. (c) and (d) magnitude and real part of the complex amplitude $u_{12}$ by sampling below the Abbe limit, used as a reference. $N_F=1000.$

Download Full Size | PDF

The correlation coefficient for the magnitude and real part of the amplitude are $r=0.97$ respectively $r=0.95$, which shows a remarkable improvement compared to $r=0.65$ and $r=0.49$ for the case of applying the conventional RSD ${\mathcal{R}}_{12}$. If we compare the image in Figs. 3(a) and 3(b) with the below Abbe sampling image in Figs. 8(c) and 8(d), we have a $32\%$ improvement in the magnitude and $46\%$ in the real part of the amplitude.

In Fig. 9 the reconstruction of the object $u_{13231}$ by the use of the operator ${\mathcal{R}}_{132}$ for the forward and ${\mathcal{R}}_{231}$ for the backward propagation is presented. The correlation coefficients for the magnitude and the real part of the complex amplitude are $r=0.97$. Again, the capability of the proposed approach can be seen.

 

Fig. 9 (a), (b) amplitude and real part of the input object (c), (d) magnitude and real part of the amplitude for the reconstructed object.

Download Full Size | PDF

5. 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 $r=0.97$. 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 $\delta x_1$ and $\delta x_2$ to prevent aliasing.