Diffraction algorithm suitable for both near and far field with shifted destination window and oblique illumination

Cheng-Shan Guo; Yi-Yan Xie; Bei Sha

doi:10.1364/OL.39.002338

With the development of computer, image sensing, and modulating technology, numerical simulation of free-space diffraction in computer is more and more required in many research fields for optical information processing, such as lensless diffractive imaging [1,2], reconstruction of digital holograms [3,4], design of computer-generated holograms or diffractive optical elements [5], and so on. Many algorithms for simulation of free-space diffraction have been reported up to now [6–22], in which more commonly interested and adopted are the algorithms based on the angular spectrum (AS) theory of diffraction and fast Fourier transform (FFT) [12–22]. According to the AS theory, the diffraction field of a beam after propagating through a distance of $z$ in free-space can be simply expressed as [23]
随着计算机、图像传感和调制技术的发展，在光学信息处理的许多研究领域，如无透镜衍射成像[1, 2]、数字图像重建等，越来越需要计算机对自由空间衍射进行数值模拟。全息图 [3, 4]、计算机生成全息图或衍射光学元件的设计 [5] 等等。迄今为止，已经报道了许多模拟自由空间衍射的算法[6-22]，其中更普遍感兴趣和采用的是基于衍射角谱（AS）理论和快速傅里叶变换（FFT）的算法[ 12-22]。根据AS理论，光束在自由空间中传播 $z$ 距离后的衍射场可以简单地表示为[23]

u_{o} (x, y) = F^{- 1} {F {u_{i} (x, y)} H (ξ, η)},

in which,

F {}

and

F^{- 1} {}

are forward and inverse Fourier transform, respectively,

u_{i} (x, y)

is the complex amplitude distribution of the object beam at the input plane.

H (ξ, η)

is a function in spatial frequency domain, which has the following form:
其中，

F {}

和

F^{- 1} {}

分别为正傅里叶变换和逆傅里叶变换，

u_{i} (x, y)

为输入平面处物光束的复振幅分布。

H (ξ, η)

是空间频域的函数，其形式如下：

H (ξ, η) = \exp {j k z \sqrt{1 - λ^{2} (ξ^{2} + η^{2})}},

where

k = 2 π / λ

,

λ

is the wavelength. Equation (1) reveals the linearity of the free-space diffraction transform from the input field

u_{i} (x, y)

to the output field

u_{o} (x, y)

. So we often named the function

H (ξ, η)

as the transfer function (TF) of the free-space diffraction transform. In principle, Eq. (1) is suitable for calculating diffraction of both near and far field, because it is derived directly from the scalar Helmholtz equation without approximation. In the numerical simulations, however, the conventional FFT-based AS algorithms are suitable only for near-field regions because of the aliasing error induced from sampling the TF. Another limitation is that the sampling rate on the source plane is the same as that on the destination plane due to the property of FFT. For overcoming the limitations, some ways have been reported in recent years [14–22]. For example, in [14] Matsushima and Shimobaba present a band-limited angular spectrum (BLAS) method for resolving the sampling problem by limiting the bandwidth and truncating unnecessary high-frequency signals in the input source field. However, when the diffraction distance is relatively large, the accuracy of BLAS would also rapidly decrease. The method recently proposed by Xiao et al. [20] is much more effective, in which a large enough window size is chosen to make sure that the spectral sampling interval is small enough and all the sampling in the frequency domain is effective. But a drawback of the method is the increased computational complexity. And the methods are not suitable for both near and far field with shifted destination window and large oblique source illumination at the same time.
其中

k = 2 π / λ

、

λ

是波长。方程(1)揭示了从输入场

u_{i} (x, y)

到输出场

u_{o} (x, y)

的自由空间衍射变换的线性。因此我们常将函数

H (ξ, η)

命名为自由空间衍射变换的传递函数（TF）。原则上，等式。式（1）适用于计算近场和远场的衍射，因为它是直接从标量亥姆霍兹方程推导出来的，无需近似。然而，在数值模拟中，传统的基于 FFT 的 AS 算法仅适用于近场区域，因为对 TF 进行采样会产生混叠误差。另一个限制是，由于 FFT 的特性，源平面上的采样率与目标平面上的采样率相同。为了克服这些限制，近年来报道了一些方法[14-22]。例如，在[14]中，Matsushima和Shimobaba提出了一种限带角谱（BLAS）方法，通过限制带宽并截断输入源场中不必要的高频信号来解决采样问题。然而，当衍射距离较大时，BLAS的精度也会迅速下降。肖等人最近提出的方法。 [20]更有效，其中选择足够大的窗口大小以确保频谱采样间隔足够小并且频域中的所有采样都是有效的。但该方法的缺点是计算复杂度增加。并且该方法不适用于同时具有移动目标窗口和大斜光源照明的近场和远场。

In this Letter, we propose a method for simulation of both near- and far-field diffractions with shifted destination window and large oblique source illumination at the same time. The method is based on a so-called generalized transfer function (GTF) of free-space diffraction and the corresponding generalized point-spread function (GPSF) derived from the angular spectrum theory under the condition with a shifted destination window and a large oblique illumination. One of the advantages of the method is that it has the same low computing cost by using the FFT twice both for near-field and far-field diffraction. In addition, a large oblique illumination source can be introduced, even if the sampling interval at the input plane does not satisfy the sampling condition for such an illumination. This feature may be useful in the situations involved in an incident beam with large oblique angle, such as in super-resolved digital holography.

To explain our method, let us begin by considering the typical diffraction geometry as shown in Fig. 1. An object with a transmittance function of $t_{i} (x_{i}, y_{i})$ is illuminated by an oblique plane wave with a tilt angle of $θ$ . The center of the input sampling area is placed at the origin of the coordinate system, while the region of interest in the output plane is located apart from the origin. The diffracted filed of the object wave after propagating a distance $z$ through free space to the output plane can be calculated by

u_{o} (x, y) = F^{- 1} {F {t (x_{i}, y_{i}) \exp (j k x_{i} \sin θ)} H (ξ, η)} .

Fig. 1. Diffraction geometry with shifted destination window and an oblique illumination.

Download Full Size | PDF

In Eq. (3), a default setting is that the coordinates $(x, y)$ on the output plane are the same as the coordinates $(x_{i}, y_{i})$ on the input plane except a translation $z$ along the optical axis. In conventional FFT-based numerical methods for this situation, the input sampling window should be extended by increasing some zero pads so that it includes the region of interest on the output plane at the same time, which would obviously lead to huge computational complexity, especially in cases where the region of interest is far from the optical axis. In addition, the oblique incident beam with a large tilt angle would further aggravate the complexity because of the sampling problem if it is sampled together with the illuminated object function.

To overcome these limitations, we first redefine the origin of the coordinates $(x, y)$ by shifting a displacement on the output plane. For simplicity, suppose the displacement is $x_{0}$ along $x$ direction, as shown in Fig. 1. Then Eq. (3) should be changed into

u_{o} (x, y) = F^{- 1} {F {t (x_{i} + x_{0}, y_{i}) \exp (j k (x_{i} + x_{0}) \sin θ)} H (ξ, η)} = F^{- 1} {T (ξ - \frac{\sin θ}{λ}, η) \exp (j 2 π x_{0} ξ) H (ξ, η)},

where

T (ξ, η)

is the Fourier transform of the object function

t_{i} (x_{i}, y_{i})

. Equation (4) can be further simplified after a coordinates transformation of the spatial frequency as
其中

T (ξ, η)

是对象函数

t_{i} (x_{i}, y_{i})

的傅里叶变换。对空间频率进行坐标变换后，式(4)可以进一步简化为

u_{o} (x, y) = \exp {j k (x + x_{0}) \sin θ} F^{- 1} {T (ξ, η) H_{g} (ξ, η)},

in which 其中

H_{g} (ξ, η) = \exp {j k [λ x_{0} ξ + z \sqrt{1 - {(λ ξ + \sin θ)}^{2} - λ^{2} η^{2}}]} .

Equations (5) and (6) provide us with an effective algorithm for calculating the diffraction with both oblique illumination and shifted destination window without notable increase in computational complexity. For the sake of distinction, we name the function $H_{g} (ξ, η)$ the GTF and the algorithm based on Eqs. (5) and (6) the GTF sampling algorithm.
方程（5）和（6）为我们提供了一种有效的算法，用于计算倾斜照明和移动目标窗口的衍射，而不会显着增加计算复杂性。为了区分，我们将函数 $H_{g} (ξ, η)$ 命名为GTF，并将基于方程的算法命名为 $H_{g} (ξ, η)$ 。 (5)和(6)GTF采样算法。

Nevertheless, the GTF algorithm is also only suitable for near-field diffraction because of the sampling problems caused in sampling the GTF. The suitable scope of the method can be found out by using Nyquist condition to the GTF. Because the local “spatial frequency” of the GTF is equal to
然而，由于GTF采样时带来的采样问题，GTF算法也只适用于近场衍射。通过对GTF使用奈奎斯特条件可以找出该方法的适用范围。因为 GTF 的局部“空间频率”等于

\frac{\partial φ (ξ, η)}{2 π \partial ξ} = x_{0} - \frac{z (λ ξ + \sin θ)}{\sqrt{1 - {(λ ξ + \sin θ)}^{2} - λ^{2} η^{2}}},

to avoid the sampling problems, the sampling interval

δ u

in spatial frequency domain should satisfy the following condition (here we restrict ourselves to one-dimension situation just for simplicity):
为了避免采样问题，空间频域的采样间隔

δ u

应满足以下条件（这里为了简单起见，我们将自己限制为一维情况）：

\frac{1}{δ u} \geq 2 {| [x_{0} - \frac{z (λ ξ + \sin θ)}{\sqrt{1 - {(λ ξ + \sin θ)}^{2}}}] |}_{\max} .

Supposing the sampling number and the sampling interval of the input window is respectively $M$ and $δ x$ , and the FFT technique is adopted in numerical simulation, the sampling number in the spatial frequency domain will also be equal to $M$ and the sampling interval must satisfy $δ u = 1 / (M δ x)$ . Because in most situations we have $z (λ ξ_{\max} + \sin θ) / \sqrt{1 - {(λ ξ_{\max} + \sin θ)}^{2}} \geq x_{0}$ , Eq. (8) can be finally rewritten as
假设输入窗口的采样数和采样间隔分别为 $M$ 和 $δ x$ ，数值模拟中采用FFT技术，则空间频域的采样数也将等于 $M$ 并且采样间隔必须满足 $δ u = 1 / (M δ x)$ 。因为在大多数情况下我们有 $z (λ ξ_{\max} + \sin θ) / \sqrt{1 - {(λ ξ_{\max} + \sin θ)}^{2}} \geq x_{0}$ ，等式： (8) 最终可改写为

z \leq Z_{c},

where 在哪里

Z_{c} = \frac{M {(δ x)}^{2}}{λ} \frac{(1 + \frac{2 x_{0}}{M δ x})}{(1 + \frac{2 δ x \sin θ}{λ})} \sqrt{1 - \frac{λ^{2}}{4 {(δ x)}^{2}} {(1 + \frac{2 δ x \sin θ}{λ})}^{2}},

by which the allowed diffraction distance of the FFT-based GTF sampling method can be quantitatively determined. It can be seen that the allowed distance should be less than or equal to a special distance

Z_{c}

. Here we name

Z_{c}

determined by Eq. (10) the criteria distance of free-space diffraction. For example, if the simulation parameters are taken as

M = 512

,

δ x = 5 μm

,

λ = 0.5 μm

,

θ = 30 d eg

, and

x_{0} = 9.0 mm

, the criteria distance is

Z_{c} = 15.61 mm

, that is to say, the above FFT-based GTF sampling method is only suitable for a short diffraction distance less than 15.61 mm in this case.
由此可以定量确定基于FFT的GTF采样方法的允许衍射距离。可见，允许的距离应该小于或等于特殊距离

Z_{c}

。这里我们命名

Z_{c}

，由等式确定。 (10)自由空间衍射的标准距离。例如，如果模拟参数为

M = 512

、

δ x = 5 μm

、

λ = 0.5 μm

、

θ = 30 d eg

和

x_{0} = 9.0 mm

，标准距离为

Z_{c} = 15.61 mm

，也就是说，本例中上述基于FFT的GTF采样方法仅适用于小于15.61mm的短衍射距离。

Although Xiao et al. [20] have demonstrated an improved TF algorithm also suitable for a long diffraction distance, it is difficult to be expanded to the situations with both an oblique illumination and a shifted destination window. Here we resolve the problem based on the convolution representation of Eq. (5). According to the convolution theorem of Fourier transform, Eq. (5) can be expressed as
尽管肖等人。 [20]已经证明了一种改进的TF算法也适用于长衍射距离，但很难扩展到同时具有倾斜照明和移动目标窗口的情况。这里我们基于式（1）的卷积表示来解决这个问题。（5）。根据傅里叶变换的卷积定理，方程： (5) 可以表示为

u_{o} (x, y) = \frac{z}{j λ} \exp {j k x \sin θ} {t (x_{i}, y_{i}) * h_{g} (x_{i}, y_{i})},

where the sign “

*

” represents the convolution operator and
其中符号“

*

”代表卷积算子，

h_{g} (x_{i}, y_{i}) = \exp (j k x \sin θ) F^{- 1} {H_{g} (ξ, η)} \approx \frac{\exp {j k [- x_{i} \sin θ + \sqrt{z^{2} + {(x_{i} + x_{0})}^{2} + y_{i}^{2}}]}}{[z^{2} + {(x_{i} + x_{0})}^{2} + y_{i}^{2}]} .

Although Eqs. (11) and (12) are totally equivalent to Eqs. (5) and (6) in mathematics, they have a completely different sampling problem and so a different suitable scope, because the sampling to the function $h_{g} (x_{i}, y_{i})$ is carried out in real space, but not in spatial frequency domain. Accordingly, we name function $h_{g} (x_{i}, y_{i})$ presented by Eq. (12) the GPSF and the algorithm based on Eqs. (11) and (12) the GPSF sampling algorithm.

Similar analysis to the sampling problem based on the Nyquist sampling theorem indicates that the suitable scope of the GPSF sampling algorithm could be determined by the following formula:

\frac{1}{δ x} \geq 2 {| \frac{1}{λ} [- \sin θ + \frac{(x + x_{0})}{\sqrt{{(x + x_{0})}^{2} + z^{2}}}] |}_{\max} .

It can be further derived from Eq. (13) that the allowed diffraction distance of the FFT-based GPSF sampling algorithm must satisfy (for

(x + x_{0}) / \sqrt{{(x + x_{0})}^{2} + z^{2}} > \sin θ

)
可以从等式进一步推导出来。 (13) 基于FFT的GPSF采样算法的允许衍射距离必须满足（对于

(x + x_{0}) / \sqrt{{(x + x_{0})}^{2} + z^{2}} > \sin θ

）

z \geq Z_{c} .

Equation (14) reveals that the GSPF sampling algorithm is only suitable for the situations with a long diffraction distance, and interestingly, the allowed distance is just larger than or equal to the criteria distance

Z_{c}

also determined by Eq. (10).
由式(14)可知，GSPF采样算法仅适用于衍射距离较长的情况，有趣的是，允许的距离恰好大于或等于式(14)确定的标准距离

Z_{c}

。（10）。

Obviously, the optimal solution for our purpose is to combine both the GTF and the GPSF sampling algorithms into a hybrid algorithm based on the criteria distance $Z_{c}$ derived above, which could be called the GTF–GPSF algorithm. Figure 2 shows the flow chart of the GTF–GPSF algorithm for programming.

Fig. 2. Flow chart of the GTF-GSPF algorithm for programming.

Download Full Size | PDF

To verify the feasibility of the GTF–GPSF algorithm, one-dimensional diffraction by a rectangular aperture is first simulated. The simulation parameters are taken as follows: the wavelength of the source $λ = 500 nm$ , the sampling interval $δ x = 10 λ$ , the sampling number $M = 512$ , the width of the rectangular $a = 256 δ x$ , and the tilt angle of the incident beam $θ = 30 d eg$ . The shift of the destination window is set to be $x_{0} = z t g θ$ (where $z$ is the diffraction distance), to ensure that the main diffraction spot does not shift out of the destination window for any diffraction distance $z$ because of the oblique illumination beam. The simulation accuracy obtained respectively by GTF and GSPF sampling algorithms is shown in Fig. 3 as a function of the diffraction distance. The accuracy is evaluated by the SNR of the diffraction field got by the evaluated methods in comparison with a reference field provided by rigorous numerical integration [20,22]. From Fig. 3 it can be seen that, as the theoretical prediction, the GTF sampling method is valid for a short diffraction distance, while the simulation results using GSPF sampling are more accurate for the long diffraction distance. The criteria distance $Z_{c}$ provide us a reliable criterion for choice between these two sampling algorithms.

Fig. 3. Accuracy of the GTF and GSPF.

Download Full Size | PDF

We further simulate the diffraction in two dimensions. The simulation parameters are: the wavelength $λ = 500 nm$ , the sampling intervals $δ x = δ y = 10 λ$ , the sampling number $M = N = 512$ , and the tilt angle of the incident beam $θ = 30 d eg$ . The object on the input plane is an annular aperture with the inner diameter of 0.32 mm and the outer diameter of 0.64 mm. The shift of the destination window is also taken as $x_{0} = z t g θ$ . Figure 4 shows some examples of the simulation. Figures 4(a), 4(b), and 4(c) are the diffraction intensities simulated by the GTF sampling algorithm when the diffraction distance is taken as 3, 15, and 150 mm, respectively. It can be seen that some “grid noise” appears on the diffraction pattern with the long diffraction distance of 150 mm. Figures 4(d), 4(e), and 4(f) are the simulated diffraction intensities by the GSPF sampling algorithm using the same distances. In contrast, the severe error emerges in the diffraction when the short distance of 3 mm is taken. Figures 4(g), 4(h), and 4(i) give the results calculated by a program based on the GTF-GSPF sampling algorithm as shown in Fig. 2, which are all consistent with the results obtained by rigorous numerical integration.

Fig. 4. Diffraction intensities simulated by different sampling algorithms when the diffraction distance is taken as 3, 15, and 150 mm, respectively.

Download Full Size | PDF

In summary, based on the angular spectrum theory, we derived the GTF [Eq. (6)] and the GSPF [Eq. (12)] suitable for free-space diffraction simulation when both a large oblique illumination and a shifted destination window are taken into account. The theoretical analysis and the simulation results have proved that the GTF sampling is valid for the diffraction distance less than or equal to the criteria distance $Z_{c}$ , while the GSPF sampling is only suitable for the distance larger than or equal to $Z_{c}$ . This criteria distance $Z_{c}$ can be quantitatively calculated by an analytical formula as shown in Eq. (10). If a work needs rapidly simulating a free-space diffraction with a wide range of diffraction distance and with both large oblique illumination and shifted destination window at the same time, the hybrid GTF-GSPF algorithm as charted in Fig. 2 could be a new choice with some advantages in low computing cost and a wide range of adaptability.

The work is supported in part by the NSFC under grant 11074152 and the Research Foundation for the Doctoral Program of Higher Education of China under grant 20113704110002.

References

1. G. J. Williams, H. M. Quiney, B. B. Dhal, C. Q. Tran, K. A. Nugent, A. G. Peele, D. Paterson, and M. D. de Jonge, Phys. Rev. Lett. 97, 025506 (2006). [CrossRef]

2. T. Shimobaba, M. Makowski, T. Kakue, M. Oikawa, N. Okada, Y. Endo, R. Hirayama, and T. Ito, Opt. Express 21, 25285 (2013). [CrossRef]

3. W. Bishara, T.-W. Su, A. F. Coskun, and A. Ozcan, Opt. Express 18, 11181 (2010). [CrossRef]

4. D. W. E. Noom, K. S. E. Eikema, and S. Witte, Opt. Lett. 39, 193 (2014). [CrossRef]

5. N. Okada, T. Shimobaba, Y. Ichihashi, R. Oi, K. Yamamoto, M. Oikawa, T. Kakue, N. Masuda, and T. Ito, Opt. Express 21, 9192 (2013). [CrossRef]

6. L. Yaroslavsky, Proc. SPIE 6252, 625216 (2006). [CrossRef]

7. R. P. Muffoletto, J. M. Tyler, and J. E. Tohline, Opt. Express 15, 5631 (2007). [CrossRef]

8. J. F. Restrepo and J. G. Sucerquia, Appl. Opt. 49, 6430 (2010). [CrossRef]

9. J. Lin, X.-C. Yuan, S. S. Kou, C. J. R. Sheppard, O. G. Rodríguez-Herrera, and J. C. Dainty, Opt. Lett. 36, 1341 (2011). [CrossRef]

10. G. Bora Esmer, Appl. Opt. 52, A18 (2013). [CrossRef]

11. P. Lobaz, Opt. Express 21, 2795 (2013). [CrossRef]

12. S. De Nicola, A. Finizio, G. Pierattini, P. Ferraro, and D. Alfieri, Opt. Express 13, 9935 (2005). [CrossRef]

13. K. Falaggis, T. Kozacki, and M. Kujawinska, Appl. Opt. 52, 3288 (2013). [CrossRef]

14. K. Matsushima and T. Shimobaba, Opt. Express 17, 19662 (2009). [CrossRef]

15. K. Matsushima, Opt. Express 18, 18453 (2010). [CrossRef]

16. S. Odate, C. Koike, H. Toba, T. Koike, A. Sugaya, K. Sugisaki, K. Otaki, and K. Uchikawa, Opt. Express 19, 14268 (2011). [CrossRef]

17. P. Lobaz, Opt. Express 19, 32 (2011). [CrossRef]

18. T. Shimobaba, K. Matsushima, T. Kakue, N. Masuda, and T. Ito, Opt. Lett. 37, 4128 (2012). [CrossRef]

19. T. Kozacki, K. Falaggis, and M. Kujawinska, Appl. Opt. 51, 7080 (2012). [CrossRef]

20. Y. Xiao, X. Tang, Y. Qin, H. Peng, and W. Wang, Opt. Lett. 37, 4943 (2012). [CrossRef]

21. K. Yamamoto, Y. Ichihashi, T. Senoh, R. Oi, and T. Kurita, Opt. Express 20, 12949 (2012). [CrossRef]

22. T. Shimobaba, T. Kakue, M. Oikawa, N. Okada, Y. Endo, R. Hirayama, and T. Ito, Opt. Lett. 38, 5130 (2013). [CrossRef]

23. J. W. Goodman, Introduction to Fourier Optics, 3rd ed. (Roberts & Company, 2005).

Diffraction algorithm suitable for both near and far field with shifted destination window and oblique illumination
适用于具有移动目标窗口和倾斜照明的近场和远场的衍射算法

Abstract 抽象的

References

Cited By

Figures (4)

Equations (14)

Optics Letters