Fresnel diffraction in the case of an inclined image plane

Peter Modregger; Daniel Lübbert; Peter Schäfer; Rolf Köhler; Timm Weitkamp; Michael Hanke; Tilo Baumbach

doi:10.1364/OE.16.005141

1. Introduction

The well-known Fresnel integral relates a known complex wave defined in the object plane (the input wave field) to the observable complex wave (the output wave field) defined in the image plane after free-space propagation [1]. In this particular case, the corresponding imaging system is said to be shift-invariant if the object and image plane are parallel to each other. This means that if the object is shifted, the same image is received, except for a corresponding shift in the image plane [2]. In this case, the Fresnel integral can be quickly and easily calculated by Fourier transform. This advantageous property was essential for the development of phase sensitive imaging techniques such as x-ray in-line holography [3], phase retrieval with visible-light [4] or electron microscopy [5].

However, if the image plane is inclined with respect to the incident beam, the effective propagation distance will vary over the image plane. Consequently, the imaging system is not shift-invariant and the Fresnel integral can no longer be calculated as a Fourier transform in straight-forward way. In [6] it is stated that the use of a Fourier-based approach to the theoretical formalism requires the imaging system to be linear and shift-invariant. In this letter we show that the imaging system corresponding to Fresnel diffraction with an inclined image plane can still be treated by a Fourier-based approach. This is accomplished by replacing the original coordinate system and thus showing the shift-invariance of the imaging system.

Several imaging techniques face the problem of varying propagation distances and calculating the corresponding Fresnel diffraction integral. One example is in-line holography applied for mirror characterization. The problem is explicitly mentioned in [7] and worked around by using discrete wavelet transform instead of Fourier transform.

Closely related are energy dispersive speckle experiments used for reconstruction of the height distribution of a sample surface [8]. A pinhole located upstream of the sample is used to ensure coherent illumination. The problem is to calculate the Fresnel diffraction pattern of the pinhole at the position of the sample. In this special case the so-called Lommel formalism [9] can be used, implying basically the summation of approximately one hundred Bessel functions.

Another example is x-ray imaging with two-fold asymmetric Bragg reflection, the so called “Bragg Magnifier”. In [10] numerical simulations of the imaging process were performed by applying the Fourier transform for a set of perpendicular planes at fixed propagation distances and then interpolating to the desired positions in the image plane.

However, all mentioned approaches of calculating the Fresnel diffraction integral in case of an inclined image plane have either limited fields of application and/or are numerically slow compared to calculating a simple Fourier transform.

Another approach to the problem was published in [11] and numerically verified in [12]. Based on Rayleigh-Sommerfeld diffraction, the authors used a plane-wave decomposition of the input wave field in order to calculate the wave field on an inclined image plane. Similar to the treatment presented in this article, the plane-wave decomposition takes advantage of the computational speed of Fourier transform. However, the two approaches differ in the way how the inclined image plane is treated theoretically. Thus, this article can be viewed as an alternative treatment of the problem in the case of Fresnel diffraction.

In the subsequent sections we will first demonstrate the influence of varying propagation distances to the observable intensity for the example of imaging with the Bragg Magnifier. Then we will lay out the theoretical formalism of Fresnel diffraction in the case of an inclined image plane, showing that Fourier transform can still be used to calculate the observable intensity. The viability of the resulting algorithm as well as the reduction of computation time by several orders of magnitude is shown by a numerical example.

Naturally, the objective of phase-sensitive imaging is to reconstruct the phase distribution of the wave directly after transmission through the sample (i.e. the phase map). The advantage of greatly reduced computation time now allows us to achieve this objective in the case of imaging with the Bragg Magnifier. By adapting an iterative phase retrieval algorithm, previously developed for electron microscopy, the phase map becomes accessible while taking effects of Fresnel diffraction into account. As a first approach the discussion will be limited to one spatial dimension. The convergence and applicability of the iterative algorithm is shown by a numerical simulation under quasi-realistic imaging conditions.

2. Experimental example

Fig. 1. (Color online) (a) Experimental image of a copper mesh obtained with the Bragg Magnifier. (b) Intensity of the first and last period along the line indicated in (a). The different intensity distribution (see arrows) of identical objects features due to varying propagation distances within the image is clearly visible. The mean propagation distances between sample and the surface of the analyzer crystal are about 5mm for the first and 39mm for the last period, respectively.

Download Full Size | PDF

In this section, the Bragg Magnifier is used to experimentally demonstrate the occurrence of varying propagation distances in the image plane. The experimental setup consists of two asymmetrically cut analyzer crystals with perpendicular diffraction planes. The analyzer crystals act both as magnifying devices (each for one direction) and as phase-to-intensity contrast converters. The latter can be explained as follows. The sample introduces refraction into the incident beam, so that beam deviations from the main beam direction occur locally. During reflection of the beam at the analyzer crystal these local beam deviations are filtered according to their angular position on the reflection curve. As the refraction angle is directly proportional to the derivative of the phase, the signal obtained is then (neglecting effects of Fresnel diffraction) proportional to the derivative of the phase itself (see [13] for more details). Since free-space propagation after reflection is practically negligible [10], the surface of the analyzer crystal can be regarded as an inclined image plane.

Figure 1(a) shows an experimental image of a copper mesh as acquired with the Bragg Magnifier. The 10µm thick copper mesh is normally used for calibration and has a periodicity of 63.5µm in both directions. A photon energy of 8.048 keV corresponding to a wavelength of 1.54 Åwas used, yielding 40.4-fold magnification in vertical and 41.6-fold magnification in horizontal direction. A highly efficient CCD camera (Bruker AXS Smart Apex 2) with a pixel size of 15µm at 4096² pixels was used for image acquisition. The measurement was performed at the beamline ID19 at the European Synchrotron Radiation Facility in Grenoble (France). Figure 1(b) shows the intensity distribution of the first and last period along the horizontal line in Fig. 1(a). The intensity was averaged over seven horizontal lines to increase the signal-to-noise ratio.

In general, both Fresnel diffraction and Bragg reflection contribute to the visible intensity distribution of the Bragg Magnifier. On the one hand, a complex interaction between these two effects might be expected, but in [10] it was shown that Fresnel fringes are attenuated during reflection. This effect turns out to be asymmetric, meaning that Fresnel fringes on left hand side of object edges are more strongly attenuated than on the right hand side if the angular position of the main beam direction is chosen to be on the left slope of the reflection curve (clearly visible in the left part of Fig. 1(b)).

On the other hand, the contribution of Fresnel diffraction increases with increasing propagation distance. Therefore, as the periods of the copper mesh consist of identical units, the difference of the intensity distribution in the left and right part of figure 1(b) can be explained by different Fresnel propagation distances.

The left side of Fig. 1(a) corresponds to smaller, the right side to larger propagation distances. The averaged sample-to-analyzer distance was close to 5mm for the first period shown in Fig. 1(b) and about 39mm for the last period.

3. Theory

Fig. 2. Coordinate system for the theory.

Download Full Size | PDF

The Fresnel diffraction integral in Fourier space relates the Fourier transform D̂_in(q) of the input wave field D_in(x) to the observable output wave field D_out (x) after free-space propagation over the distance z+z₀. The Fourier transform variable q corresponds to the spatial variable x. For a one-dimensional object (the extension to two object dimensions will be discussed below) it can be written as [1]

D_{out} (x) \propto \int dq {\hat{D}}_{in} (q) \exp (- i \frac{z + z_{0}}{2 K} q^{2}) \exp (iqx)

where the modulus of the wavevector K=2π/λ is given by the wavelength λ and z₀ is an arbitrary distance between object and image plane (the mean propagation distance would be a convenient choice; see Fig. 2). As long as the image plane is perpendicular to the main beam direction (i.e. z is constant) Eq. (1) constitutes a simple Fourier transform.

But if the image plane is inclined (i.e.α≠π/2), the propagation distance becomes dependent on the position in the image (i.e. z=z(x)). As stated in the introduction, D_out can be calculated by performing Fourier transforms for many planes of constant propagation distances (typically several hundred slices) and interpolating the results to the coordinates of the inclined image plane s. However, this is obviously a time-consuming approach.

Therefore, we now introduce a more convenient approach that will require only one Fourier transform even for an inclined image plane. First, the integral (1) is rewritten to express the output wave field with respect to the position on the image plane. This yields

D_{out} (s) \propto \int dq {\hat{D}}_{in} (q) \exp (- i \frac{z_{0}}{2 K} q^{2}) \exp (is (q \sin α - \frac{\cos α}{2 K} q^{2}))

where the identities x=s sinα and z=scosα (see Fig. 2) have been used. The idea of the following substitution is to reduce the last factor to exp(is f), thus converting the integral (1) to a simple Fourier transform. Obviously, this can be achieved by choosing

f = q \sin α - \frac{\cos α}{2 K} q^{2} .

This implies

q = q (f) = K \tan α \pm \sqrt{K^{2} \tan^{2} α - \frac{2 K}{\cos α} f}

where the occurring ambiguity (q=K tanα divides the graph of Eq. (4) into two branches) is resolved as follows. The Fourier component q corresponds to a plane wave propagating with an angular deviation of ≈q/K with respect to the main beam direction. But the validity of the Fresnel diffraction integral is limited to small angular deviations from the main beam direction [1]. This implies the condition q≪K which can only be fulfilled by choosing the negative sign in Eq. (4).

For reasons of convenience, we will use the so-called propagator P̂(q) to relate input and output wave field in Fourier space. While in the usual case of shift-invariant imaging techniques the relation is given by D̂_out(q)=P̂(q)D̂_in(q), it will turn out that in the present case the relation needs to be modified. Therefore, we use the propagator in a generalized sense.

Using $dq = \frac{df}{(\sin α - \frac{\cos α}{K} q)}$ the substitution leads to

D_{out} (s) \propto \int d f \hat{P} (f) {\hat{D}}_{in} (q (f)) e^{ifs}

where the generalized propagator for Fresnel diffraction in the case of an inclined image plane is introduced as

\hat{P} (f) = \frac{1}{\sin α - \frac{\cos α}{K} q (f)} \exp (- i \frac{z_{0}}{2 K} q {(f)}^{2}) .

Please note that all occurring q’s in equations (5) and (6) are functions of f according to Eq. (4) with negative sign.

Equations (4)–(6) can now be used to numerically calculate the Fresnel diffraction pattern on an inclined image plane with one Fourier transform. In the actual computation the Fourier components f form a regular grid so that sample values of the input wave field in Fourier space D̂_in(q(f)) have to be interpolated onto that grid. This is a necessary, potentially time-consuming step during the numerical calculations. However, compared to the calculation with many Fourier transforms plus interpolation as described above, a time saving in the order of the desired number of points along the image plane can still be expected. In order to numerically prove the viability of the proposed algorithm and to determine the time saving, Fig. (3) shows a numerical example (details see caption).

Fig. 3. (Color online) Simulated Fresnel diffraction in case of an inclined image plane (α=0.5°). The figure compares the intensity distribution of a double rectangular input wave field (5µm width; 10µm apart) as computed by many Fourier transforms for planes of constant propagation distance (here 512) plus interpolation and by the proposed algorithm according to equations (4)–(6). The effect of varying propagation distance is clearly visible. While the intensity distributions show negligible difference, the computation time of the proposed algorithm is reduced by a factor of about 200.

Download Full Size | PDF

The extension of the proposed algorithm to two object dimensions (x and y direction) is simple if the image plane is only rotated around the y-axis and the orientations of the x,y coordinate axes in the object plane are chosen in a convenient way. By setting the y-direction to be perpendicular to the normal of the image plane, the propagation distance in y-direction is constant in the image plane. Hence, the y-direction can be treated as usual, yielding

D_{out} (s, y) \propto \int d f d p {\hat{D}}_{in} (q (f), p) \hat{P} (f) \exp (i f s + i p y - i \frac{z_{0}}{2 K} p^{2}) .

The Fourier transform variable p corresponds to the spatial variable y.

4. Iterative phase retrieval algorithm for the Bragg Magnifier

The standard post-detection analysis of images obtained with analyzer-based imaging techniques is the diffraction enhanced imaging algorithm (DEI) first introduced by Chapman et al. [13]. Using DEI it is possible to separate the contribution of absorption and refraction contrast to images obtained with a one analyzer crystal setup. Recently, we proposed an extension to imaging with the Bragg Magnifier and its two analyzer crystals with perpendicular diffraction planes called 2D-DEI [14]. This is necessary to account for a non-linearity of the orthogonal refraction directions contributing to the observable intensity.

However, the underlying theory of DEI is based on a geometrical optics approach, thus neglecting the coherent influences of free-space propagation and reflection. These influences act typically on a submicrometer scale. In order to avoid artifacts on this scale an alternative procedure based on an iterative phase retrieval algorithm introduced by Coene et al. for electron microscopy [5] is proposed in the following.

In this first investigation, we restrict the discussion to one-dimensional objects and perfectly collimated, monochromatic incident plane waves (i.e. the case of perfect coherence). The possible inclusion of finite coherence and the extension to the two-dimensional case will be discussed in the conclusions.

Including the reflection curve (see [10]) in the propagator in Eq. (6) yields the 1D-propagator for the Bragg Magnifier

{\hat{P}}_{B M} (f) = \frac{1}{\sin α - \frac{\cos α}{K} q (f)} \hat{R} (- \frac{q (f)}{K} + Δ θ) \exp (- i \frac{z_{0}}{2 K} q {(f)}^{2})

where Δθ defines the angular position of the main beam direction on the analyzer reflection curve (working point). The observable intensity of a known input wave field is then given by

I (s) = {∣ \int d f {\hat{P}}_{B M} (q (f)) \hat{D} (q (f)) e^{i f s} ∣}^{2} .

A challenge of phase-sensitive imaging is the loss of phase information during detection. As both absorption and phase information contribute to the visible intensities one needs at least two independent experimental images to separate these contributions. Mathematically speaking, independent means that the propagators corresponding to the experimental images have to differ. This can be achieved by varying a free experimental parameter (e.g. propagation distance in case of in-line holography [3] or focal series in electron microscopy [5]). In the case of the Bragg Magnifier Δθ is the natural choice of the free experimental parameter realizing the independence of the corresponding experimental images.

A general approach of reconstructing the complex input wave field (corresponding to absorption and phase information about the sample) is to minimize the quadratic difference between the experimental and calculated images. Let I^exp_m be the mth experimentally obtained intensity of a total of M pictures(e.g. M=2 for two pictures recorded on opposite slope positions of the reflection curve). Further, let I_m(s) be the calculated intensity according to Eq. (9). Then the quadratic difference is given by

S = \sum_{m} \int ds {(I_{m}^{\exp} (s) - I_{m} (s))}^{2} .

The proposed phase retrieval algorithm starts with an arbitrary input wave field in Fourier space D̂ ₀(f). The j+1th iteration is performed by calculating (a direct adaptation of Eq. (6) in [5])

{\hat{D}}_{j + 1} (f) = {\hat{D}}_{j} (f) + \frac{γ}{M} \sum_{m = 1}^{M} {\hat{P}}_{m}^{*} (f) F {(I_{m}^{\exp} (s) - I_{m j} (s)) F^{- 1} {{\hat{P}}_{m} (f) {\hat{D}}_{j} (f)}}

with γ a convergence parameter, F{…} (F ^-1{…}) indicating the (inverse) Fourier transform of its argument and I_mj(s) the intensity (calculated according to Eq. 9) corresponding to the current wave field D̂_j.

Obviously, the algorithm has a fixed point at the desired solution (I^exp_m(s)=I_mj(s)) while the function of the iteration can be understood by reading Eq. (11) from inside to outside. First, the current wave field D̂_j is propagated and transformed to real space (F ^-1{D̂_j(f)P̂_m(f)}), yielding the current input wave field in real space. This wave field is now weighted with the difference between experiment and calculation (I^exp_m(s)-I_mj(s)) and the result is then back-transformed to Fourier space. The multiplication with the conjugated propagator P̂*_m(f) corresponds to a back-propagation of the modified wave field. Finally, the new wave field D̂ _j+1 is obtained by the weighted sum over all modified wave fields.

The result of this algorithm is the wave field as a function f. So, before obtaining the wave field in real space it has to be rewritten in terms of the original Fourier argument q. In this first attempt, this is done by simply calculating q=f/sinα, which is a linearized version of Eq. (3) avoiding the need of interpolation.

A demonstration of the feasibility

Fig. 4. (Color online) (a) The Shepp-Logan phantom, which is used as a realistic model sample. It is illuminated in the direction indicated by the arrows. (b) Numerically simulated intensities after free-space propagation (mean distance 30mm) and asymmetric reflection according to eq. (9). The two intensity distributions corresponding to opposite slope position (top left; bottom right) of the reflection curve, are used as input for the iteration algorithm (see Fig. 5).

Download Full Size | PDF

Except for the restriction to perfect coherence and to one spatial dimension, we have chosen as realistic imaging conditions as possible for the following forward and backward simulations.

A Shepp-Logan phantom [15], usually used to check algorithms for tomographic reconstruction, was chosen to provide a realistic model sample (see Fig. 4a). The total sample width was set to 1mm to allow for significantly varying propagation distances in the intensity distribution. The refractive index of the sample was adapted to simulate a typical biological object and summation along the direction of the arrows indicated in Fig. 4(a) provided the input wave field D_in(x). In order to give a smooth appearance and to minimize the discontinuity at the edges, the input wave field was convoluted with a small rectangular filter of 2.5µm width.

The forward simulation of the observable intensities (Fig. 4b) was performed by numerically calculating the integral (9). The mean propagation distance was 30mm (similar to the actual experimental setup). The reflection curve was calculated for a Si-224σ reflection at 8.048 keV and 40-fold magnification corresponding to the first analyzer crystal of the Bragg Magnifier.

The iterative phase retrieval algorithm was implemented according to eq. (11) using MAT-LAB (Mathworks). Figure 5(a) shows the comparison between the original and reconstructed phase distribution of the Shepp-Logan phantom. We have used 2048 pixels and 50,000 iterations steps were performed taking 14 minutes on a standard 2 GHz PC. From the very good agreement between original and reconstructed phase we conclude on the convergence and feasibility of the proposed algorithm. Furthermore, Fig. 5(b) shows a comparison between the input intensity on the right slope (same as bottom curve in Fig. 4(b)) and the output intensity corresponding to the final input wave field that was delivered by the iterative phase reconstruction algorithm. The excellent agreement underlines the consistency of the reconstructed data with the input intensity.

Fig. 5. (Color online) (a) The comparison of the original and reconstructed phase distribution of the Shepp-Logan phantom (see Fig. 4a) shows a very good agreement. (b) The original intensity distribution on the right slope (same as bottom line in Fig. 4(b)) and the intensity distribution corresponding to the input wave field after reconstruction calculated according to eq. (9). The comparison shows the consistency of the reconstructed data with the input intensity.

Download Full Size | PDF

5. Conclusion

We have extended the theoretical formalism of Fresnel diffraction to the case of an inclined image plane for one- and two-dimensional objects. Thereby, we have shown that the Fresnel integral can still be calculated by simple Fourier transform. The resulting numerical algorithm speeds up computation times by typically several orders of magnitude compared to calculating many Fourier transforms plus interpolation. This could prove beneficial for forward simulations as well as image analysis procedures for experiments facing the problem of an inclined image plane.

An example for this kind of experiment is x-ray imaging with the Bragg Magnifier. Due to the time-saving of calculating the Fresnel integral it was possible to adapt an iterative phase retrieval algorithm to this imaging technique. Numerical forward and backward simulations demonstrated the convergence and usability of the algorithm in the case of one-dimensional objects and coherent illumination.

The next step will be to extend the theoretical formalism of the image formation process and consequently the iterative phase retrieval algorithm to two-dimensional objects and partially coherent illumination. Since the second spatial direction can be treated identically to the first [10], the extension to two-dimensional objects is straightforward. Partial coherence on the other hand can be included in a way similar to the one presented in [16]. Using both extensions will lead to a full description of the image formation process and thus open the possibility of quantitative tomography on a sub-micrometer scale with the Bragg Magnifier.

We would like to acknowledge the assistance of Jane Richter and Rainer Schurbert in experiment preparation as well as Jürgen Härtwig and Lukas Helfen during the experiments.

References and links

1. M. Born and E. Wolf, Principles of Optics (7th edn; 1999), Cambridge University Press, Cambridge.

2. C.E. Metz and K. Doi, “Transfer function analysis of radiographic imaging systems,” Phys. Med. Biol. 24, 1079–1106 (1979). [CrossRef] [PubMed]

3. P. Cloetens, W. Ludwig, J. Baruchel, D. van Dyck, J. van Landuyt, J. P. Guigay, and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron raditaion x-rays,” Appl. Phys. Lett. 75, 2912–2914 (1999). [CrossRef]

4. V. Yu. Ivanov, V. P. Sivokon, and M. A. Vorontsov, “Phase retrieval from a set of intensity measurements: theory and experiment,” J. Opt. Soc. Am. A 9,1515–1524 (1992). [CrossRef]

5. W. Coene, G. Janssen, M. Op de Beeck, and D. Van Dyck, “Phase retrieval through focus variation for ultra-resolution in field-emission transmission electron microscopy,” Phys. Rev. Lett. 69, 3743–3746 (1992). [CrossRef] [PubMed]

6. I. A. Cunningham and R. Shaw, “Signal-to-noise optimization of medical imaging systems,” J. Opt. Soc. Am. A 16621–632 (1999). [CrossRef]

7. A. Souvorov, M. Yabashi, K. Tamasaku, T. Ishikawa, Y. Mori, K. Yamauchi, K. Yamamura, and A. Saito, “Deterministic retrieval of surface waviness by means of topography with coherent X-rays,” J. Synchrotron Rad. 9, 223–228 (2002). [CrossRef]

8. T. Panzner, G. Dleber, T. Sant, W. Leitenberger, and U. Pietsch, “Coherence experiments at the white-beam beamline of BESSY II,” Thin Sol. Films 5155563–5567 (2007). [CrossRef]

9. K. D. Mielenz, “Algorithms for Fresnel diffraction at rectangular and circular apertures,” J. Res. Natl. Inst. Stand. Technol. 103497–509 (1998). [CrossRef]

10. P. Modregger, D. Lübbert, P. Schäfer, and R. Köhler, “Magnified phase imaging using asymmetric Bragg reflection: experiment and theory,” Phys. Rev. B 74054107-1–054107-10 (2006). [CrossRef]

11. N. Delen and B. Hooker, “Free-space beam propagation between arbitrarily oriented planes based on full diffraction theory: a fast Fourier transform approach,” J. Opt. Soc. Am. A 15, 857–867 (1998). [CrossRef]

12. N. Delen and B. Hooker, “Verification and comparison of a fast Fourier transform-based full diffraction method for tilted and offset planes,” Appl. Opt. 40, 3525–3531 (2001). [CrossRef]

13. D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pisano, N. Gmur, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced x-ray imaging,” Phys. Med. Biol. 422015–2025 (1997). [CrossRef] [PubMed]

14. P. Modregger, D. Lübbert, P. Schäfer, and R. Köhler, “Two dimensional diffraction enhanced imaging algorithm,” Appl. Phys. Lett. 90, 193501-1–193501-3 (2007). [CrossRef]

15. A.C. Kak and M. Slaney “Principles of Computerized Tomographic Imaging,” IEEE Press. (1988).

16. Ya. I. Nesterets, T. E. Gureyev, and S. W. Wilkins, “Polychromaticity in the combined propagation-based/analyser-based phase-contrast imaging,” J. Appl. Phys. D: Appl. Phys. 384259–4271 (2005). [CrossRef]

Fresnel diffraction in the case of an inclined image plane

Abstract

1. Introduction

2. Experimental example

3. Theory

4. Iterative phase retrieval algorithm for the Bragg Magnifier

A demonstration of the feasibility

5. Conclusion

References and links

Cited By

Figures (5)

Equations (11)

Optics Express