Computational algorithm from the Huygens-Fresnel’s diffraction integral for two-dimensional holographic reconstruction

Abstract

While most common holographic methods of digital reconstruction are based on the convolution theory, for the ease in the mathematical approach, here we present an algorithm by a discretization of the Huygens-Fresnel integral from a Taylor series expansion to produce a bidimensional Fourier transform. Compared to the digital convolution method, the algorithm presented here is more concise and generates a reduction in processing time, since the Fourier transform appears only once in the discretization. Another advantage is associated with the production of results in the frequency domain, allowing the optical information to be obtained directly.

Keywords
Digital holography; Fresnel method; Digital reconstruction; Fresnel transform

1. Introduction

Holography was proposed in 1948 by the hungarian physicist Dennis Gabor in the work: A new microscope principle [¹]. His initial idea was to improve the quality of the image in an electron microscope. Gabor achieved his goal when he discovered that a second reference wave, interfering with the object wave, produced a significant improvement in the overall quality of the image.

In the early 1960s, holography became of great interest due to the invention of the laser, theorized by Arthur Schawlow and built by Theodore H. Maiman [²]. The coherent laser light, coupled with the adjustment of a holographic technique called the off-axis recording technique developed by Emmett Leith and Juris Upatnieks in 1962 [³], provided an increase in the quality of the holographic images and motivated, even more, the research of the time. With this, the Holography found relevant importance in the scientific universe taking Dennis Gabor to the Nobel Prize of physics in 1971.

Digital holography (HD) was introduced in the early 1970s by Yaroslavski et al. [⁴]. Currently, some techniques use digital cameras to record holograms and computational methods for intensity and phase reconstructions. There are several methods to digitally reconstruct holograms, the main two are known as the convolution method and Fresnel method [⁵, ⁶, ⁷].

The wieldy use of HD throughout the years has resulted in great advances in both techniques and methods inside the scientific and technological universes, as shown by the works of Wang (2018), Narayanamurthy (2017), Da Silva (2017), e Zhang (2012) [⁸, ⁹, ¹⁰, ¹¹]. Currently, digital holography is a fast and efficient tool that can be very useful for the analysis of mechanical effects on materials [¹²].

This paper shows the construction of an algorithm to perform digital holographic reconstruction, based on a discretization of the Fresnel approximation for the Huygens-Fresnel diffraction integral. When compared to discretization by convolution, this discretization is considerably simpler, since the Fourier transform, which appears only once, implies a reduction in processing time. Furthermore, the result obtained is given in the frequency domain, facilitating works that use this domain directly, as opposed to the convolution method that produces results in the spatial domain.

2. Theoretical Analysis

2.1. Two-dimensional Fresnel transform

The method presented in this work is based on the Fresnel approximation to the Huygens-Fresnel integral of the function (1) [¹³, ¹⁴], and so it is called the Fresnel method.

(1) $${\overrightarrow{\psi }}_0(\overrightarrow{\mathrm{r}})=\frac{1}{\mathrm{i}\lambda }{\iint }_A\frac{{\overrightarrow{\psi }}_{\mathrm{dif}}({\overrightarrow{\mathrm{r}}}^{\prime })}{{\mathrm{r}}_{\mathrm{P}}}{\mathrm{e}}^{-\mathrm{i}\overrightarrow{\mathrm{k}}\cdot {\overrightarrow{\mathrm{r}}}_{\mathrm{P}}}\cos \theta \mathrm{ds}$$

Expression (1) represents the calculation of the diffracted image field, ${\overrightarrow{\psi }}_0(\overrightarrow{\mathrm{r}})$, in HD. While ${\overrightarrow{\psi }}_{\mathrm{dif}}({\overrightarrow{\mathrm{r}}}^{\prime })$ represents the diffraction between the incident wave with wavelength $\lambda$ and the hologram with area A. The term $\cos \theta =\mathrm{z}/\mathrm{r}$ is associated with the geometry of the reconstruction, as shown in Figure1.

Figure 1
Fresnel’s holographic reconstruction geometry.

The goal of this work is to reconstruct a digital image of a sample object by holography. The diffraction will be the result of a reference wave, ${\overrightarrow{\psi }}_{\mathrm{ref}}({\overrightarrow{\mathrm{r}}}^{\prime })$, which will focus on the hologram, represented by its intensity ${\mathrm{I}}_{\mathrm{hol}}({\overrightarrow{\mathrm{r}}}^{\prime })$. The final product of this method will be a discrete expression, which is the basis of the digital holographic reconstruction algorithm. With these considerations, the initial expression of Huygens-Fresnel is

(2) $${\overrightarrow{\psi }}_{\mathrm{rec}}(\overrightarrow{\mathrm{r}})=\frac{1}{\mathrm{i}\lambda }\iint {\mathrm{I}}_{\mathrm{hol}}({\overrightarrow{\mathrm{r}}}^{\prime }){\overrightarrow{\psi }}_{\mathrm{ref}}({\overrightarrow{\mathrm{r}}}^{\prime })\frac{{\mathrm{e}}^{-\mathrm{i}\overrightarrow{\mathrm{k}}\cdot {\overrightarrow{\mathrm{r}}}_{\mathrm{P}}}}{{\mathrm{r}}_{\mathrm{P}}}\cos \theta \mathrm{ds}$$

where ${\overrightarrow{\psi }}_{\mathrm{rec}}(\overrightarrow{\mathrm{r}})$ is the complex amplitude of the reconstructed wave. Using the holographic reconstruction geometry outlined in Fig.1, where xy represents the plane of the object, $\xi \eta$ the plane of the hologram, uv the reconstruction plane and z the direction of propagation, function (2) can be rewritten as

(3) $${\overrightarrow{\psi }}_{\mathrm{rec}}(u,v)={\displaystyle \frac{1}{\mathrm{i}\lambda }}{\displaystyle \iint {\mathrm{I}}_{\mathrm{hol}}(\xi ,\eta ){\overrightarrow{\varphi }}_{\mathrm{ref}}(\xi ,\eta )\frac{{\mathrm{e}}^{-\mathrm{i}\overrightarrow{\mathrm{k}}\cdot \overrightarrow{\mathrm{r}}(\xi ,\eta ,u,v)}}{r(\xi ,\eta ,u,v)}}$$ $$\mathrm{}\times \cos \theta \mathrm{d}\xi \mathrm{d}\eta$$

If $\mathrm{r}(\xi ,\eta ,u,v)={[{(\mathrm{u}-\xi )}^2+{(\mathrm{v}-\eta )}^2+{\mathrm{z}}^2]}^{1/2}$ is the distance between a pixel of the hologram plane and a pixel of the reconstruction plane, and z = r cos $\theta$ is the distance between the two planes with $\theta$ very small, then, after using approximations by Taylor series, the function (^2.1) can be rewritten as follows [³, ⁴],

(4) $${\overrightarrow{\psi }}_{\mathrm{rec}}(u,v)={\displaystyle \frac{{\mathrm{e}}^{-\mathrm{ik}\mathrm{z}}}{\mathrm{i}\lambda \mathrm{z}}}{\mathrm{e}}^{-\frac{\mathrm{ik}}{2\mathrm{z}}({\mathrm{u}}^2+{\mathrm{v}}^2)}{\displaystyle \iint \mathrm{f}(\xi ,\eta ){\mathrm{e}}^{\mathrm{i}({\mathrm{k}}_{\xi }\xi +{\mathrm{k}}_{\eta }\eta )}d\xi d\eta }$$

where were defined: ${\mathrm{k}}_{\xi }=-\frac{\mathrm{k}}{\mathrm{z}}\mathrm{u}$, ${\mathrm{k}}_{\eta }=-\frac{\mathrm{k}}{\mathrm{z}}\mathrm{v}$ and,

$$\mathrm{f}(\xi ,\eta )={\mathrm{I}}_{\mathrm{hol}}(\xi ,\eta ){\overrightarrow{\psi }}_{\mathrm{ref}}(\xi ,\eta ){\mathrm{e}}^{-\frac{\mathrm{ik}}{2\mathrm{z}}({\xi }^2+{\eta }^2)}{\mathrm{e}}^{-\frac{\mathrm{ik}}{2\mathrm{z}}({\mathrm{u}}^2+{\mathrm{v}}^2)}.$$

The term ${\mathrm{e}}^{-\frac{\mathrm{ik}}{2\mathrm{z}}({\xi }^2+{\eta }^2)}$ is called the chirp function, and represents the quadratic approximation of the spherical wave generated by the encounter of the reference wave with the hologram. The Fourier transform of the function $\mathrm{f}(\xi ,\eta )=\frac{1}{2\pi }\iint ℱ({\mathrm{k}}_{\xi },{\mathrm{k}}_{\eta }){\mathrm{e}}^{-\mathrm{i}({\mathrm{k}}_{\xi }\xi +{\mathrm{k}}_{\eta }\eta )}{\mathrm{dk}}_{\xi }{\mathrm{dk}}_{\eta }$ can be defined as

(5) $$ℱ(k_{\xi },k_{\eta })={\displaystyle \iint \mathrm{f}(\xi ,\eta ){\mathrm{e}}^{\mathrm{i}({\mathrm{k}}_{\xi }\xi +{\mathrm{k}}_{\eta }\eta )}d\xi d\eta }$$

Thus, function (4) can be rewritten as

(6) $${\overrightarrow{\psi }}_{\mathrm{rec}}(u,v)=\frac{{\mathrm{e}}^{-\mathrm{ik}\mathrm{z}}}{\mathrm{i}\lambda \mathrm{z}}{\mathrm{e}}^{-\frac{\mathrm{ik}}{2\mathrm{z}}({\mathrm{u}}^2+{\mathrm{v}}^2)}ℱ(k_{\xi },k_{\eta })$$

with the field of the reconstructed image expressed in the universe of frequencies. Function (6) is called the Two-Dimensional Fresnel Transform.

2.2. Discretization of the Fresnel transform

To discretize the Two-Dimensional Fresnel transform, expression (6), $\mathrm{\Delta }\mathrm{u}=\frac{\lambda \mathrm{z}}{\mathrm{N}\mathrm{\Delta }\xi }$ and $\mathrm{\Delta }\mathrm{v}=\frac{\lambda \mathrm{z}}{\mathrm{M}\mathrm{\Delta }\eta }$ were defined, based on the scheme of Figure2.

Figure 2
Discretization scheme of the two-dimensional fresnel transform.

But as, $\mathrm{u}=\mathrm{n}\mathrm{\Delta }\mathrm{u}=\mathrm{n}\frac{\lambda \mathrm{z}}{\mathrm{N}\mathrm{\Delta }\xi }$ and $\mathrm{v}=\mathrm{m}\mathrm{\Delta }\mathrm{v}=\mathrm{m}\frac{\lambda \mathrm{z}}{\mathrm{M}\mathrm{\Delta }\eta }$. Therefore,

(7) $${\psi }_{\mathrm{nm}}=\frac{{\mathrm{e}}^{-\mathrm{ik}\mathrm{z}}}{\mathrm{i}\lambda \mathrm{z}}{\mathrm{e}}^{-\frac{\mathrm{ik}}{2\mathrm{z}}{\lambda }^2\left[{\left(\frac{\mathrm{n}}{\mathrm{N}\mathrm{\Delta }\xi }\right)}^2+{\left(\frac{\mathrm{m}}{\mathrm{M}\mathrm{\Delta }\eta }\right)}^2\right]}ℱ(\mathrm{\Delta }{k}_{\xi },\mathrm{\Delta }{k}_{\eta })$$

where $\mathrm{\Delta }{\mathrm{k}}_{\xi }=-\frac{\mathrm{k}}{\mathrm{z}}\mathrm{\Delta }\mathrm{u}=-\frac{\mathrm{k}}{\mathrm{z}}\frac{\lambda \mathrm{z}}{\mathrm{N}\mathrm{\Delta }\xi }=-\frac{2\pi }{\mathrm{N}\mathrm{\Delta }\xi }$ and, $\mathrm{\Delta }{\mathrm{k}}_{\eta }=-\frac{\mathrm{k}}{\mathrm{z}}\mathrm{\Delta }\mathrm{u}=-\frac{\mathrm{k}}{\mathrm{z}}\frac{\lambda \mathrm{z}}{\mathrm{M}\mathrm{\Delta }\eta }=-\frac{2\pi }{\mathrm{M}\mathrm{\Delta }\eta }$, $\mathrm{n}ϵ[1,\mathrm{N}]$ and $\mathrm{m}ϵ[1,\mathrm{M}]$.

Thus, the intensity and phase of the reconstructed field can be determined, respectively, as follows,

(8a) $${\mathrm{I}}_{\mathrm{rec}}={⟨\psi _2^{\mathrm{nm}}⟩}_{\tau }={\displaystyle \frac{1}{2}}\Re e\{\psi _{*}^{\mathrm{nm}}{\psi }_{\mathrm{nm}}\}$$

(8b) $${\mathrm{\Phi }}_{\mathrm{rec}}=\mathrm{arc}\mathrm{tg}\left[{\displaystyle \frac{\Im \mathrm{m}({\psi }_{\mathrm{nm}})}{\Re e({\psi }_{\mathrm{nm}})}}\right]$$

As the recordings of holograms are carried out on the sensor of a digital camera, it must be far from the interference plane so that there is an angular difference ${\theta }_{\max }$, between the reference wave and that of the object.

This angular difference is important, as it allows the identification of diffraction orders due to the simulation of reference waves with holograms, during digital reconstruction [¹⁵]. By the sampling theorem, at least two sensor pixels are required for each period of spatial variation in the hologram intensity distribution [¹⁶, ¹⁷]. Thus,

(9) $${\theta }_{\max }<2\mathrm{arcsen}\left(\frac{\lambda }{4\mathrm{\Delta }\xi }\right).$$

3. Computational Program

Here we describe the general steps for the creation of a computer program, based on the Fresnel method, capable of reconstructing the hologram picture. The determinations of intensity and phase maps were performed utilizing the mathematical algorithms shown in expressions (7), (8a), and (8b).

The program input is composed of two digital images: (i) a hologram, obtained by any experimental technique of Holography, and (ii) an image of the same reference wave used during the production of the hologram.

3.1. Structure:

PART 1 – Initial data capture and treatment:

• (1.1)

  Read input images, hologram and reference, and assign images to two 8-bit matrix variables: one for the hologram and one for the reference;

• (1.2)

  Convert the arrays of 8-bit images to floating-point arrays and calculate the mean intensities of the converted matrices;

• (1.3)

  Calculate the intensities for each cell of the matrices using the following expression:

  (10) $$\mathit{\text{New Int.}}=\frac{\mathit{\text{Old intensity}}\cdot \mathit{\text{Intensity factor}}}{\mathit{\text{Mean intensity}}}$$

  The intensity factor is the value calculated from the measured power that arrives on the digital camera sensor;

• (1.4)

  Store the height and width of the images in variables. The images must have the same area.

PART 2 – Select region of interest of the hologram:

• (2.1)

  Filter and centralize hologram and reference image;

• (2.2)

  Calculate the logarithm of the absolute value of the hologram transform;

• (2.3)

  Select the region for reconstruction, in the hologram picture;

• (2.4)

  Round the values of the selection regions for construction of new regions, and store the new images for the final transform.

PART 3 – Construct a matrix of the image plane from the propagation matrix of the hologram plane:

• (3.1)

  Center the selected images from the dimensions of the hologram;

• (3.2)

  Include only the amplitudes of the selected images, in the complex format, in the central regions of the matrices.

PART 4 – Calculate the Chirp function and carry out the reconstruction:

• (4.1)

  Calculate the coefficient of the Chirp function, by approximation of the Taylor series in first order.

• (4.2)

  Calculate the reconstruction matrices from the hologram and the reference pictures;

• (4.3)

  Calculate the intensities matrix from the reconstruction matrices;

• (4.4)

  Calculate the phase matrix from the reconstruction matrices.

PART 5 – Display the reconstruction:

• (5.1)

  Show the reference and hologram images.

• (5.2)

  Show the phase and intensity maps.

4. Application

The program was used to reconstruct a hologram of a half-wave plate placed in one of the arms of a Mach-Zehnder interferometer. With the addition of two polarizers, one for the object arm and another for the reference arm of the interferometer, it was possible to control the intensity of the images and, consequently, improve the quality of the reconstruction.

Figure3 shows the experimental configuration used to register the images.

Figure 3
Mach-Zehnder interferometer. Where L is the laser; M1, M2, M3 and M4 are plane mirrors; F is a spatial filter; C is a plane-concave lens; I is an iris; B1 and B2 are beam splitters; P1 and P2 are linear polarizers; A is a polarizer analyzer; S is the sample (half-wave plate); D is the CMOS digital camera.

The best reconstruction was obtained when the intensities of both the object and reference arms were approximated to the same value.

The results of intensity and phase reconstructions, with the developed program, are presented in Figure4.

Figure 4
A – hologram part of a half-wave plate; B – Fresnel transform of the hologram; C – intensity map; D – demodulated phase map.

During data processing, a larger region was selected in the Fourier space, frame B, to improve the definition of the resulting image. In frame D, the outline represented by the darker curved line divides the regions of the half-wave plate, on the left, and the free space, on the right.

5. Conclusion

The program showed efficiency in the reconstruction of the images, both intensity and phase, by the Fresnel method, observed in the images of Figure4. Good data obtained from property calibrated experimental configurations are important for good holographic reconstructions. A suitable selection of the region of the phase space, Fresnel space, showed improvement in the definition of the reconstructed images.

It is important to state that, one can improve the observations of stress and deformation of materials by making the holographic reconstruction algorithm include a Third dimension in the reconstructed image, thus increasing the range of values to be analyzed.

References

• ¹
  D.A. Gabor, Nature 161, 4098 (1948).
• ²
  B.E.A. Saleh and M.C. Teich, Fundamentals of Photonics (John Wiley & Sons Inc., Ottawa, 2007), 2ª ed.
• ³
  E.N. Leith and J. Upatnieks, JOSA 52, 1123 (1962).
• ⁴
  M.A. Kronrod, N.S. Merzlyakov and L.P. Yaroslavski, Sov. Phys.-Tech. Phys. 17, 329 (2018).
• ⁵
  M.S.Z. Taietti, Métodos recursivos para o cálculo da integral de convolução. Masters Dissertation, Universidade Federal do Rio Grande do Sul, Porto Alegre (2002).
• ⁶
  U. Schnars and W.P. Jüptner, Measurement science and technology 13, R85 (2002).
• ⁷
  A.V. Fantin, Holografia digital complexa utilizando um interferômetro shearing . Doctoral Thesis, Universidade Federal de Santa Catarina, Florianópolis (2003).
• ⁸
  F. Wang, Y. Zhang, H. Wang, W. Xu, Y. Zhang and C. Li, Journal of Optics 47, 547 (2018).
• ⁹
  C.S. Narayanamurthy, G. Pedrini, and W. Osten, Applied Optics 56, F213 (2017).
• ¹⁰
  S.L. Silva, Estudo de tensões em amostras fotoelásticas com Holografia Digital (Novas Edições Acadêmicas, São Paulo, 2017), 1ª ed.
• ¹¹
  Y. Zhang, M. Huang, H. Liang and F. Lao, Optics and Lasers in Engineering 50, 619 (2012).
• ¹²
  S.L. Da Silva, F.M. Prado, D.J. Toffoli and N.U. Wetter, in: Proceedings of Photonics West – Practical Holography XXXIV: Displays, Materials, and Applications, edited by H. I. Bjelkhagen (SPIE, San Francisco, 2020), p. 113060.
• ¹³
  E. Hecht, Óptica (Fundação Calouste Gulbenkian, Lisboa, 2002), 2ª ed.
• ¹⁴
  M.A. Kronrod, N.S. Merzlyakov and L.P. Yaroslavsii, Sov. Phys. Tech. Phys. 17, 333 (1972).
• ¹⁵
  I. Naydenova, Advanced Holography: Metrology and Imaging (InTech, Rijeka, 2011).
• ¹⁶
  E. Carcolé, J. Campos and S. Bosch. Applied Optics 33, 162 (1994).
• ¹⁷
  M. Jacquot, P. Sandoz and G. Tribillon, Optics communications 190, 87 (2001).

Publication Dates

• Publication in this collection
  06 Dec 2021
• Date of issue
  2022

History

• Received
  24 May 2021
• Reviewed
  20 Oct 2021
• Accepted
  05 Nov 2021