Diffraction patterns from opaque planar objects simulated with Maggi-Rubinowicz method and angular spectrum theory

The Maggi–Rubinowicz method (MRM) is a useful tool to compute diffraction patterns from opaque planar objects. We adapted the MRM to planar rectangles. In the first part of this study, differences between diffraction patterns, both the intensity and the phase distributions, from a tilted rectangle and from the square having the same orthogonal projection on the observation plane, have been analyzed. In the second part, we compared results obtained with the MRM to those obtained with angular spectrum theory (AST) coupled to fast Fourier transform (FFT). The main novelty of this work is the fact that MRM is particularly well suited for evaluating anti-aliasing procedures applied to AST-FFT calculations. © 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
The phenomenon known as diffraction has been studied for almost four centuries. Most of these studies are based on the Huygens-Fresnel principle which states that every point of a wavefront may be considered as a source of secondary spherical wavelets, and the wavefront at any later instant may be regarded as the envelope of these wavelets [1]. This principle led to the prediction of the existence of a bright spot at the center of the shadow of an opaque disc: the "Poisson's spot" [2].
The rigorous theory of diffraction is based on Maxwell's equations; the effect of radiation incident on a perfectly conducting body may be interpreted in terms of the induced surface currents (see Chap. 11 in [1]). Such techniques can be considered as emission of waves by a specified fixed surface as a whole. At the same time, the diffraction pattern produced by an object with large dimensions and at large distances compared to the wavelength can be computed by the scalar diffraction theory resulting from Green's theorem applied to the Helmholtz wave equation [2].
The Helmholtz-Kirchhoff diffraction integral, the Fresnel-Kirchhoff formula, the Rayleigh-Sommerfeld formulation, and the angular spectrum theory (AST) can be considered as different formulations of the general scalar diffraction theory [2,3]. The Maggi-Rubinowicz representation of the Helmholtz-Kirchhoff diffraction integral, hereafter the Maggi-Rubinowicz method (MRM), along with the generalization of [4] has the same domain of validity. The key peculiarity of the MRM consists in the fact that it is based on a boundary diffraction approach. The MRM can be adapted to different planar objects. However, the method is quite time consuming and the analytical parameterization of the contour has to be developed for each object.
Fast Fourier transform (FFT) can be used for numerical implementation of the AST. FFT-based algorithms are easy to implement and effective. That is why the AST is extensively employed in different domains, including simulations of diffraction patterns (e.g., see [5,6]). According to our computations, the AST-FFT modeling is several orders of magnitude faster than MRM-based calculations (all other things being the same). At the same time, the AST-FFT has shortcomings inherent in calculations employing 2D-FFT such as artifacts due to aliasing. The first objective of this paper is to present simulation results obtained with the MRM for opaque planar objects (including tilted rectangles). A second objective is to compare diffraction patterns computed with the MRM and the AST-FFT in order to evaluate performance of a low-pass filter employed along with the AST-FFT calculations. The main novelty of this work is the fact that MRM is particularly well suited for evaluating anti-aliasing procedures applied to AST-FFT calculations. Throughout this work, it is assumed that the dimensions of opaque planar objects and the distance from an object to the corresponding diffraction pattern are much larger than the wavelength. Section 2 is devoted to analytical formulations of the MRM and simulation results. Diffraction patterns obtained with the AST-FFT are compared to MRM results in Section 3. Figure 1(a) schematically illustrates the geometric setup for the diffraction study of a monochromatic plane wave U(Q) of wavelength λ interacting with an opaque planar object. The contour Γ represents the boundary of the object. For simplicity sake, it is assumed that the wave is of unit amplitude. Within the accuracy of the Kirchhoff diffraction theory, the diffracted field U d (Q) behind the object may be expressed as a superposition of the boundary wave U b (Q) with the "geometrical wave" U g (Q):

Maggi-Rubinowicz method
where and S is the distance PQ, R is the distance OP, k = 2π λ is the wavenumber, and ì k, ì p, ì s, and ì l are the unit vectors shown in Fig. 1(a). Equations (1)-(3) are the Maggi-Rubinowicz representation of the Helmholtz-Kirchhoff diffraction integral [4]. Note that ì s has the opposite direction in [4]. Hence, the minus sign before ì s in the denominator of Eq. (1). In the numerator, the sign before ì s is unchanged because Eq. (3) is written for an opaque object (not for an aperture) and the Babinet's principle is taken into account. If the plane wave is perpendicular to the opaque planar object, then ì k · ì p = 0 and ì s ∧ ì k · ì l = ì k ∧ ì l · ì s = ì n · ì s with ì n the inward unit vector normal to the contour Γ. Moreover, if the considered diffraction pattern at the distance Z is parallel to the opaque planar object, thus perpendicular to the plane wave, then ì s · ì k = Z S . In this case, Eq. (3) can be simplified: For any known shape, it is possible, using Eqs. (2) and (3) or Eqs. (2) and (4), to compute the corresponding diffraction patterns at any distance. Throughout this work, the z-axis is normal to the observation plane. In the following subsections we apply this theory to circular and rectangular planar objects.

Opaque disc
If the opaque planar object is a disc with radius R ( Fig. 1(b)), then ì n · ì s = R−r cos α S and integrating along the contour (dl = Rdα), Eq. (4) becomes The angle α in Eqs. (5) and (6) can be defined to start from the direction O Q rather than from the direction Ox of the object plane due to the fact that the integrand of Eq. (5) is a periodic function with the period 2π. Equations (5) and (6) taken together lead to another form of Eq. (5) by straightforward algebraic operations Generally, S depends on α, so do all the terms of the integrand in Eq. (7). If point Q belongs to the axis OO , r = 0 and the on-axis amplitude behind the disc is expressed as follows: (8) is identical to Eq. (9) of the work by Lucke [7] where the on-axis amplitude was obtained using the Fresnel-Kirchhoff diffraction integral. Figure 2(a) shows the diffraction pattern from an opaque disc for the Fresnel number The diffraction pattern is presented in light intensity percentage of the incident light at the planar object. Figure axes ξ x and ξ y are x and y coordinates normalized by R. It is seen that the Poisson's spot, i.e., a bright spot in the center of the disc shadow is well reproduced by the MRM. If we consider the outer limit of a 50 % intensity threshold (orange dashed line) as the edge of the disc shadow, we notice that the disc size is 30 % larger than the real disc object for this specific case when N F = 1.

Opaque rectangle
If the opaque planar object is a rectangle with the width 2a and the length 2b ( Fig. 1(c)), then integrating along the contour Γ of Eq. (4) yields: with where y P = −a in the first integral, x P = b in the second integral, y P = a in the third integral, and x P = −b in the fourth integral. We simulated a huge set of diffraction patterns from opaque rectangles of different sizes in a large range of distances Z for the wavelength of 0.783 µm. The patterns were computed for the screen of the size of 1000 × 1000 µm at the observation plane. The shorter and the longer sides of the rectangles were within the ranges 2a > 30λ and 2a < 2b < 400 µm.
In general, a diffraction pattern primarily depends on the Fresnel number N F = b 2 /(Zλ). However, it also depends on the aspect ratio (length-to-width ratio). This is why the values of N F should be considered as approximate estimates.
For N F about 10 or higher, the patterns are quite close to the shape of the rectangle. For 2 < N F < 10, the patterns are very complex and a number of small Poisson's spots can be observed. For N F < 2, the patterns are less complex and more smoothed. Figure 2(b) shows an example of the diffraction pattern for N F = 1, i.e. the Fresnel diffraction regime. There are two large Poisson's spots. Basically, the shape is an oval rather than a rectangle. For N F < 0.3, the patterns resemble those computed using the Fraunhofer approximation. It should be mentioned that patterns in the focal plane of a lens, i.e. for Z → ∞, are not considered in this work. The Fraunhofer diffraction regime was only considered for large but limited distances Z, i.e. N F values were not smaller than 0.03. For example, if the Fraunhofer approximation is used (see [2], Chap. 4.4.1) to simulate diffraction from an opaque square of 100 × 100 µm at the distance of 2.5 cm and the Babinet's principle is applied to the amplitude, the obtained intensity pattern has much in common with Fig. 3(a)

Tilted opaque rectangle
Diffraction calculations provide a basis for computer-generated holograms (e.g., see [8] and references therein). Fast algorithms are required for such calculations. That is why, the overwhelming majority of algorithms employ the AST-FFT modeling. Simulations of diffraction patterns from tilted planar objects are not an exception (e.g., see [5]). As mentioned above, the MRM is time consuming. However, it has the following advantages. It is a straightforward matter to develop an analytical parameterization of a tilted planar rectangle. Calculations can be performed on the equidistant sampling grid.
In this section, using the MRM, we treat the single-axis rotation around the y-axis (Fig. 1(c)). The diffraction patterns and the phase distributions from the tilted opaque planar rectangle of The results for these two case studies can be seen in the video files (Visualization 1, Visualization 2, Visualization 3, and Visualization 4). When the z-axis is normal to the rectangle, θ = 0°. The orthogonal projection of a 60°tilted rectangle on the observation plane is equal to a square of 100 × 100 µm.
It is instructive to compare diffraction properties of the tilted 100 × 200 µm rectangle when θ = 60°and the corresponding square (100 × 100 µm) when the z-axis is normal to the square. show the difference ∆ = data rectangle − data square . For Z = 2.5 cm, the differences between the intensities ∆ I (Fig. 3(b)) and between the phase distributions ∆ P (Fig. 3(d)) are quite small. For Z = 0.5 cm, the differences ∆ I (Fig. 4(b)) and ∆ P (Fig. 4(d)) are well pronounced (different color scales compared to Figs. 3(b) and 3(d)). In our opinion, the differences ∆ I (Fig. 4(b)) between the intensity patterns are sufficient to distinguish the tilted rectangle from the square if the patterns will be used as holograms for reconstructing 3D object information.
The important property of diffraction patterns computed with the MRM is their high quality. There are no artifacts such as aliasing. Figures 2-4 were obtained without any smoothing. If some applications need a higher resolution, this can be easily achieved by grid refinement. To complete this Sect. 2.3, we would like to point out that analytical parameterizations of laterally shifted or tilted laterally-shifted rectangles are straightforward as well. Diffraction patterns of the same quality for such planar objects can be obtained without extra effort.

AST fundamental relationships
Other than the Maggi-Rubinowicz method, the angular spectrum theory (AST) can be used to compute the diffraction pattern from an illuminated opaque object. We recall the milestones of the AST following the textbooks [2,3]. Suppose that a monochromatic plane wave is propagating in the positive Z direction and the complex field across the plane Z = 0 is represented by U i (x, y, 0). A diffracting structure is introduced in the plane Z = 0. The amplitude transmittance function t A (x, y) is the ratio of the transmitted field amplitude U t (x, y, 0) to the incident field amplitude U i (x, y, 0) at each position (x, y) in the Z = 0 plane, with: In this work, the amplitude transmittance function t A (x, y) is defined as follows: In the Fourier domain, the angular spectrum of U t (x, y, 0) is: The angular spectrum of the wave field U(x, y, z) at Z = z is given by a solution of the differential equation which represents the Helmholtz equation in the Fourier domain [3]. The solution can be written as: where under the condition of homogeneous waves (k 2 > 4π 2 ( f 2 x + f 2 y )), which is the case in this study. Applying the inverse Fourier transform gives the resulting wave field at the distance Z = z: Note that despite the apparent differences of their approaches, "the angular spectrum approach and the first Rayleigh-Sommerfeld solution yield identical predictions of diffracted fields" (see [2], p. 61). A theoretical demonstration of the interrelationship between the AST and the general Rayleigh-Sommerfeld diffraction formula can be found in the work by Harvey [9].

Comparison of AST and MRM simulations
Fast Fourier transform (FFT) algorithms [10] allow very fast computations of discrete Fourier transform. In this study, we implement AST-FFT to compute diffraction patterns produced from planar opaque objects. The resulting wave field can be computed by two approaches. The first one is direct: the amplitude transmittance function t A (x, y) corresponds to an opaque planar object and Eqs. (11)-(15) can be employed. The second approach is based on Babinet's principle: the resulting wave field U(x, y, z) is computed for the corresponding aperture and the complex field U i (x, y, z) is subtracted (amplitude and phase), not only the intensity (amplitude squared) (e.g., see [1], Section 8.3). Both approaches give the same results for the spatial distributions of the amplitude and of the phase when the AST-FFT modeling is used. Figure 5(a) shows the diffraction pattern from the opaque planar rectangle 2a × 2b = 200 × 400 µm at the distance of 2.5 cm (N F = 2.04) computed with the AST-FFT. Detailed analysis of the pattern (e.g., see the magnified image in the inset) reveals some artifacts, i.e. aliasing noise, appearing as very close horizontal and vertical lines. When looking at the intensity profiles along the y-axis (Fig. 5(b)) and x-axis (Fig. 5(c)), these artifacts correspond to sharp peaks on the intensity profiles (red curves) compared to the smooth profiles from Maggi-Rubinowicz method simulations (blue curves). It is well known that aliasing noise is common to algorithms using the FFT (see e.g, [8,11,12] and references therein). A number of anti-aliasing procedures can be found in the literature. The effectiveness of anti-aliasing procedures has to be evaluated by FFT-free means. In our opinion, the MRM is particularly well suited for such an evaluation.
It is seen that the aliasing noise belongs to the high-frequency domain in the case of Figs. 5(a)-5(c). A cut-off low-pass filter, which is a very simple anti-aliasing procedure, has been tested on the AST-FFT patterns. The cut-off frequency ( f c = 0.06 µm −1 ) was chosen using the root mean square deviation (RMSD) between the MRM and the AST-FFT patterns. When the low-pass filter was applied, the AST-FFT provided the diffraction pattern shown in Fig. 5(d). The pattern now appears smooth (see the magnified image in the inset). Nevertheless, there are smaller differences remaining between the MRM and the filtered AST-FFT data (see the profiles in Figs. 5(e) and 5(f)). Still, it can be seen that the AST-FFT profiles are not excessively smoothed.

Conclusions
The simulations performed in this study demonstrate that the Maggi-Rubinowicz method is a flexible tool to compute diffraction patterns from opaque planar objects. The MRM has the same domain of validity as the angular spectrum theory and the first Rayleigh-Sommerfeld solution and can provide valuable data within a large range of Fresnel-number values. An important feature of the MRM is the high quality of simulated diffraction patterns (both intensity and phase 2D distributions). The patterns are artifact-free and the resolution can be increased just by grid refinement. The MRM has two major shortcomings. First, the method is time consuming (according to our estimations, it is about three orders of magnitude slower than the AST). And second, the analytical parameterization of the contour has to be developed for each object. Direct computing of the Rayleigh-Sommerfeld integral should be free of issues inherent in FFT calculations. Compared to such an approach, the MRM has the advantage that a surface integral is replaced by a line integral. Generally, that leads to more efficient calculations.
We adapted the MRM to planar rectangles. It follows from our simulations that the intensity diffraction patterns are quite close to the shape of the rectangle for the Fresnel number N F = 10 or higher. For 2 < N F < 10, the patterns are very complex and a number of small Poisson's spots can be observed. For N F < 2, the patterns are less complex and more smoothed and their shape is an oval rather than a rectangle. For N F < 0.3, the patterns resemble those computed using the Fraunhofer approximation.
We compared diffraction patterns from a tilted 100 × 200 µm opaque rectangle when tilt angle θ = 60°and from a 100 × 100 µm opaque square when tilt angle θ = 0°, which is equal to the projection of the rectangle on the observation plane. The patterns (both intensity and phase 2D distributions) look very similar. Still, the detailed analysis revealed some differences, which should be sufficient to distinguish between the tilted rectangle and the non-tilted square if the patterns might be used as holograms for reconstructing 3D object information.
We simulated the diffraction pattern from the opaque planar rectangle (N F = 2.04) using AST-FFT and compared it with MRM results. This comparison demonstrates that the Maggi-Rubinowicz method is particularly well suited for evaluating anti-aliasing procedures applied to AST-FFT calculations. This fact is the major novelty of this work.