Optimization of Phase-Only Computer-Generated Holograms Based on the Gradient Descent Method

Abstract

The Gerchberg-Saxton (GS) algorithm is a Fourier iterative algorithm that can effectively optimize phase-only computer-generated holograms (CGHs). This study proposes a new optimization technique for phase-only CGHs based on the gradient descent method. The proposed technique evaluates the intensity distributions of reconstructed images to directly obtain the phase distributions of the CGHs, whereas the GS algorithm equivalently evaluates the amplitude distributions of reconstructed images and extracts phase distributions from complex-amplitude distributions of the holograms using a constant-amplitude constraint. The proposed technique can reduce the errors in the reconstructed images with fewer iterations than the GS algorithm.

1. Introduction

Holography [1,2] is an ideal technique for producing three-dimensional (3D) images based on wave front reconstruction. An object wave and a reference wave are optically superimposed to generate an interference pattern, which is then recorded on a photosensitive film to produce a hologram. Computer-generated holograms (CGHs) [3,4,5,6] are produced as the interference pattern which is calculated by a computer and recorded on a film. The CGH technique enables the generation of 3D images of virtual objects stored on a computer. When CGHs are displayed on spatial light modulators (SLMs), moving 3D images can be produced.

The holography technique transforms the amplitude and phase distributions of an object wave into the intensity distribution of an interference pattern. When a reconstruction wave illuminates the interference pattern, the amplitude and phase distributions of the object wave are reconstructed. Phase-only holograms, which modulate the phase distribution but not the amplitude distribution of the reconstruction wave, were developed to obtain bright reconstructed images. Optically recorded holograms can be bleached to produce phase-only holograms [7,8,9]. Because the amplitude information is lost, the reconstructed images are degraded. When phase-only holograms are produced using the CGH technique, the phase distributions can be optimized using a computer to reduce the degradation of the reconstructed images. The interference patterns do not need to be calculated because the shape of the light wavefront can be directly controlled by the phase distributions. The complex-conjugate image does not appear because the interference pattern is not recorded. Because of their high light efficiency, phase-only CGHs have been used for a wide range of applications, such as optical tweezers [10], optical interconnections [11], and optogenetic techniques [12].

Because neither the amplitude nor the phase distributions of the reconstructed object wave can be fully controlled using only phase distributions of the holograms, phase-only CGHs are generally designed to control the intensity distributions of reconstructed images. Because the phase distributions of phase-only holograms cannot be directly calculated from the intensity distributions of target images, a number of optimization techniques have been developed. Iterative algorithms, such as simulated annealing [13], genetic algorithms [14,15,16], and direct binary searches [17,18,19], have been developed to optimize phase-only CGHs. As the CGH resolution increases, the number of iterations required for these algorithms increases rapidly. Although genetic algorithms need few iterations, the large quantity of calculations is required at each iteration step. The Gerchberg-Saxton (GS) algorithm [20], which iteratively performs Fourier transforms and inverse Fourier transforms, can more rapidly optimize phase-only CGHs than the above algorithms. In particular, the number of iterations required for the GS algorithm does not increase substantially as the CGH resolution increases. Several variants of the GS algorithm have also been developed. The mixed-region amplitude freedom algorithm [21] reduced the optimization time; however, it suffered from a low diffraction efficiency. The offset mixed-region amplitude freedom algorithm [22] provided highly accurate reconstructed images, but the black region could not be controlled efficiently. Wang et al. [23] proposed yet another revised GS algorithm, which yielded high-quality reconstructed images in the preset area, nevertheless, it left the remaining regions filled with speckles. A modified GS algorithm (MGSA) [24] has also been proposed, in which the Fourier transform is replaced by the Fresnel transform in order to generate reconstructed images at any depth position. The MGSA has been used to generate multiple reconstructed images at different depth positions to create 3D images with a wide depth range [7,25].

Fienup [26] proved that the GS algorithm is mathematically equivalent to the gradient descent method, which reduces the mean square error (MSE) defined for the amplitude distribution of a reconstructed image. The gradient descent method has been widely adopted in many scientific fields, such as machine learning [27]. Although the GS algorithm reduces the error of the reconstructed image during the iteration process, the rate of error reduction decreases after several iterations.

In this paper, we adopt the gradient descent method to generate phase-only CGHs. The proposed technique calculates the gradient of the MSE defined for the intensity distributions of reconstructed images. Therefore, the proposed technique directly improves the intensity distributions of the reconstructed images. Moreover, we also modify the technique to increase the optimization speed. The proposed technique directly provides the phase distributions of phase-only CGHs, whereas the GS algorithm provides the complex-amplitude distributions of CGHs, which are then approximated for phase-only distributions.

2. Theory

2.1. Gradient Descent Method for Optimizing Phase-Only CGHs

In the gradient descent method, an iterative optimization algorithm is applied to identify a local minimum of a given function. To find a local minimum, one takes steps toward the opposite direction of the gradient of the function with a certain step length. Thus, there are two essential components in the gradient descent method. First, the direction of the iteration must be determined, and second, the step length of the iteration should be determined. The formulated description of the gradient descent method is given by

$${\theta }_{t+1}={\theta }_t-{\gamma }_t\nabla F\left({\theta }_t\right)$$

(1)

where t is the iteration number, ${\gamma }_t$ is the step length coefficient for iteration, and $\nabla F\left({\theta }_t\right)$ is the gradient of the function $F\left(\theta \right)$.

In this study, we optimize Fourier-transform-type phase-only CGHs, in which the Fourier transform of the CGH phase distributions provides the reconstructed images. The phase distribution of the CGH is denoted by ${\theta }_{pq}$, where p and q represent the coordinates of the hologram pixels and P and Q are the numbers of rows and columns of the hologram. The amplitude and phase distributions of the reconstructed image are denoted by $a_{mn}$ and ${\varphi }_{mn}$, respectively, where m and n represent the coordinates of the pixels in the reconstructed image, and M and N are the numbers of rows and columns of the reconstructed image. The relationships between the phase-only CGH and the reconstructed image are given by:

$$a_{mn}{e}^{j{\varphi }_{mn}}=\mathcal{F}\left\{e^{j{\theta }_{pq}}\right\}=\sum _p^P\sum _q^Q{e}^{j{\theta }_{pq}}{e}^{-j2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)}$$

(2)

$${\mathcal{F}}^{-1}\left\{a_{mn}{e}^{j{\varphi }_{mn}}\right\}=\frac{1}{MN}\sum _m^M\sum _n^N{a}_{mn}{e}^{j{\varphi }_{mn}}{e}^{j2\pi \left(\frac{mp}{M}+\frac{nq}{N}\right)}$$

(3)

The intensity distribution of the reconstructed image, denoted by $I_{mn}$, is given by:

$$I_{mn}=a_{mn}^2$$

(4)

The GS algorithm is illustrated in Figure 1. Two constraints are applied in each iteration; one constraint applies to the reconstructed image, and the second constraint applies to the hologram. As mentioned in Section 1, the GS algorithm is mathematically equivalent to the gradient descent method, which minimizes the MSE defined for the amplitude distribution of the reconstructed image [24]:

$$E_{GS}=\frac{1}{MN}\sum _m^M\sum _n^N{({\widehat{a}}_{mn}-a_{mn})}^2$$

(5)

where ${\widehat{a}}_{mn}$ is the amplitude distribution of the target image. The gradient of the MSE is calculated with respect to the complex-amplitude of the phase-only CGH, which is $g_{pq}=e^{j{\theta }_{pq}}$. Thus, the gradient is given by $\partial E_{GS}/\partial g_{pq}$. A constant value is used for the step length coefficient ${\gamma }_t$. The phase distribution of the CGH ${\theta }_{pq}$ is extracted from the complex-amplitude distribution $g_{pq}$.

In this study, the MSE is defined for the intensity distribution of the reconstructed image and the gradient of the MSE is calculated with respect to the phase distribution of the CGHs. Because the intensity distribution is optimized directly, the proposed technique might control the intensity distribution of the reconstructed image more precisely than the GS algorithm. The phase distribution of the CGHs is directly provided. When the intensity distribution of the target image is denoted by ${\widehat{I}}_{mn}$, the MSE is defined as:

$$E=\frac{1}{MN}\sum _m^M\sum _n^N{({\widehat{I}}_{mn}-I_{mn})}^2$$

(6)

The gradient of the MSE with respect to the phase distribution, $\partial E/\partial {\theta }_{pq}$, is calculated by the chain rule.

$$\frac{\partial E}{\partial {\theta }_{pq}}=\sum _m^M\sum _n^N\frac{\partial E}{\partial I_{mn}}\frac{\partial I_{mn}}{\partial {\theta }_{pq}}$$

(7)

$\partial E/\partial I_{mn}$ is calculated as follows:

$$\frac{\partial E}{\partial I_{mn}}=-\frac{2}{MN}({\widehat{I}}_{mn}-I_{mn})$$

(8)

The following process shows the calculation of $\partial I_{mn}/\partial {\theta }_{pq}$:

$$I_{mn}={\left\{\sum _p^P\sum _q^{Q}cos\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\right\}}^2+{\left\{\sum _p^P\sum _q^{Q}sin\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\right\}}^2$$

(9)

$$\begin{array}{cc}\hfill \frac{\partial I_{mn}}{\partial {\theta }_{pq}}& =\frac{\partial {\left\{{\displaystyle \sum _p^P}{\displaystyle \sum _q^Q}cos\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\right\}}^2}{\partial {\theta }_{pq}}+\frac{\partial {\left\{{\displaystyle \sum _p^P}{\displaystyle \sum _q^Q}sin\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\right\}}^2}{\partial {\theta }_{pq}}\hfill \\ & =-2\left\{\sum _p^P\sum _q^{Q}cos\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\right\}sin\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\hfill \\ & +2\left\{\sum _p^P\sum _q^{Q}sin\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\right\}cos\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\hfill \\ & =-2\mathrm{Re}\left\{a_{mn}{e}^{j{\varphi }_{mn}}\right\}sin\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\hfill \\ & +2\mathrm{Im}\left\{a_{mn}{e}^{j{\varphi }_{mn}}\right\}cos\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\hfill \\ & =-2a_{mn}cos{\varphi }_{mn}sin\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]+2a_{mn}sin{\varphi }_{mn}cos\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\hfill \\ & =2a_{mn}sin\left[{\varphi }_{mn}-{\theta }_{pq}+2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\hfill \end{array}$$

(10)

Then, the gradient of the MSE is given as follows:

$$\frac{\partial E}{\partial {\theta }_{pq}}=\frac{4}{MN}\sum _m^M\sum _n^N\left({\widehat{I}}_{mn}-a_{mn}^2\right)a_{mn}sin\left[{\varphi }_{mn}-{\theta }_{pq}+2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]$$

(11)

The above calculation can be described using the Fourier transform.

$$\frac{\partial E}{\partial {\theta }_{pq}}=\frac{4}{MN}\mathrm{Im}\left\{e^{-j{\theta }_{pq}}\sum _m^M\sum _n^N\left({\widehat{I}}_{mn}-a_{mn}^2\right)a_{mn}{e}^{j\left[{\varphi }_{mn}+2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]}\right\}$$

(12)

Substitute the Fourier transform (3) into (12)

$$\frac{\partial E}{\partial {\theta }_{pq}}=-4\mathrm{Im}\left\{e^{-j{\theta }_{pq}}{\mathcal{F}}^{-1}\left[\left({\widehat{I}}_{mn}-a_{mn}^2\right)a_{mn}{e}^{j{\varphi }_{mn}}\right]\right\}$$

(13)

where ${\mathcal{F}}^{-1}$ represents the inverse Fourier transform. When the fast Fourier transform algorithm is used, the calculation time can be reduced.

According to the gradient descent method, the phase distribution is provided at each iteration as follows:

$${\theta }_{pq}^{{}^{\prime }}={\theta }_{pq}+\mathsf{\Delta }{\theta }_{pq}={\theta }_{pq}-\gamma \frac{\partial E}{\partial {\theta }_{pq}}$$

(14)

2.2. Determination of Step Length Coefficient

Next, the determination of the step length coefficient is explained. In this study, the coefficient is determined using the estimated MSE of the next iteration.

The MSE of the next iteration is estimated as below:

$$E^{{}^{\prime }}=\sum _m^M\sum _n^N\frac{{\left[{\widehat{I}}_{mn}-\left(I_{mn}+\mathsf{\Delta }{I}_{mn}\right)\right]}^2}{MN}$$

(15)

where $\mathsf{\Delta }{I}_{mn}$ is the estimated change of the intensity distribution for the phase change $\mathsf{\Delta }{\theta }_{pq}$. The step length coefficient $\gamma$ is determined so that $E^{{}^{\prime }}$ becomes minimum, i.e., $\frac{\partial E^{{}^{\prime }}}{\partial \gamma }$.

$$E^{{}^{\prime }}=\frac{1}{MN}\sum _m^M\sum _n^N{\left[\left({\widehat{I}}_{mn}-a_{mn}^2\right)-\mathsf{\Delta }{I}_{mn}\right]}^2$$

(16)

$$\frac{\partial E^{{}^{\prime }}}{\partial \gamma }=\frac{2}{MN}\sum _m^M\sum _n^N\left[\left({\widehat{I}}_{mn}-a_{mn}^2\right)-\mathsf{\Delta }{I}_{mn}\right]\left(-\frac{\partial \mathsf{\Delta }{I}_{mn}}{\partial \gamma }\right)=0$$

(17)

We can adopt the first order approximation to calculate $\mathsf{\Delta }{I}_{mn}$

$$\mathsf{\Delta }{I}_{mn}\approx \sum _p^P\sum _q^Q\frac{\partial I_{mn}}{\partial {\theta }_{pq}}\mathsf{\Delta }{\theta }_{pq}$$

(18)

Substrate (10) into (18)

$$\begin{array}{cc}\hfill \mathsf{\Delta }{I}_{mn}& \approx \sum _p^P\sum _q^{Q}2a_{mn}\mathsf{\Delta }{\theta }_{pq}sin\left[{\varphi }_{mn}-{\theta }_{pq}+2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]\hfill \\ & \approx 2\mathrm{Im}\left\{a_{mn}{e}^{j{\varphi }_{mn}}\sum _p^P\sum _q^Q\mathsf{\Delta }{\theta }_{pq}{e}^{-j\left[{\theta }_{pq}-2\pi \left(\frac{mp}{P}+\frac{nq}{Q}\right)\right]}\right\}\hfill \end{array}$$

(19)

Substrate (2) into (19)

$$\mathsf{\Delta }{I}_{mn}\approx 2\mathrm{Im}\left\{a_{mn}{e}^{j{\varphi }_{mn}}{\left[\mathcal{F}\left(\mathsf{\Delta }{\theta }_{pq}{e}^{j{\theta }_{pq}}\right)\right]}^{*}\right\}=2\gamma \mathrm{Im}\left\{a_{mn}{e}^{j{\varphi }_{mn}}{\left[\mathcal{F}\left(-\frac{\partial E}{\partial {\theta }_{pq}}{e}^{j{\theta }_{pq}}\right)\right]}^{*}\right\}$$

(20)

Substrate (20) into (17).

$$\frac{2}{MN}{\sum }_m^M{\sum }_n^N\left[\left({\widehat{I}}_{mn}-a_{mn}^2\right)-2\gamma \mathrm{Im}\left\{a_{mn}{e}^{j{\varphi }_{mn}}{\left[\mathcal{F}\left(-\frac{\partial E}{\partial {\theta }_{pq}}{e}^{j{\theta }_{pq}}\right)\right]}^{*}\right\}\right]\left[-2\mathrm{Im}\left\{a_{mn}{e}^{j{\varphi }_{mn}}{\left[\mathcal{F}\left(-\frac{\partial E}{\partial {\theta }_{pq}}{e}^{j{\theta }_{pq}}\right)\right]}^{*}\right\}\right]=0$$

(21)

Finally, $\gamma$ is determined as follows:

$$\gamma =\frac{{\displaystyle \sum _m^M}{\displaystyle \sum _n^N}({\widehat{I}}_{mn}-a_{mn}^2)\mathrm{Im}\left\{a_{mn}{e}^{j{\varphi }_{mn}}{\left[\mathcal{F}\left(-\frac{\partial E}{\partial {\theta }_{pq}}{e}^{j{\theta }_{pq}}\right)\right]}^{*}\right\}}{2{\displaystyle \sum _m^M}{\displaystyle \sum _n^N}{\left[\mathrm{Im}\left\{a_{mn}{e}^{j{\varphi }_{mn}}{\left[\mathcal{F}\left(-\frac{\partial E}{\partial {\theta }_{pq}}{e}^{j{\theta }_{pq}}\right)\right]}^{*}\right\}\right]}^2}$$

(22)

The proposed technique calculates the step length coefficient for each iteration by following the above equation.

The proposed algorithm can be extended for applications in the Fresnel regime for the optimization of the phase-only holograms. The MSE function used in this study can then be replaced by other error MSE functions.

3. Proposed Calculation Process

Figure 2 shows the whole calculation process of the proposed technique. A random phase distribution is used as an initial phase distribution of the CGH. The phase-only CGH is Fourier transformed to obtain the reconstruct the image, which is multiplied by the difference of the intensity distributions of the reconstructed image and the target image. Next, the inverse Fourier transform is performed and the phase distribution of the CGH of the previous iteration is multiplied. Then, the imaginary part is extracted and –4 is multiplied to obtain the gradient. The step length coefficient is calculated to obtain the increment of the phase distribution. Finally, the phase distribution of the CGH is revised. This process is performed iteratively.

4. Verification of Proposed Technique

This section presents a verification of the proposed technique. Figure 3 shows three target images used for the verification, including one grayscale image (“mandrill”) and two binary images (“hikari” and “circle”).

The improvements in the reconstructed images obtained via the optimization process are shown in Figure 4. Here, the reconstructed images for “hikari” obtained by the proposed technique after 1, 2, 5, 10, 50, and 100 iterations are shown. For comparison, the reconstructed images obtained by the GS algorithm are shown in Figure 5. These results demonstrate that the proposed technique required few iterations than the GS algorithm for the phase optimization task.

Figure 6 shows the reconstructed images for the three target images using the proposed technique and the GS algorithm after 100 iterations. The changes in the MSE defined for the intensity distribution of the reconstructed image during the optimization processes are shown in Figure 7. These results show that lower MSE values were obtained with the proposed algorithm than with the GS algorithm. The SNRs of the reconstructed images obtained after 100 iterations were about 30 dB, when the proposed technique was adopted.

For 100 iterations, 3D intensity profiles of the reconstructed images for the “circle” target image obtained by the proposed algorithm and the GS algorithm are shown in Figure 8a,b, respectively. As shown by the results, the proposed technique more precisely controlled the intensity distribution of the reconstructed image than the GS algorithm. The intensity variation corresponding to the white region is much smaller for the proposed algorithm than for the GS algorithm, while that corresponding to the black region is larger for the proposed algorithm than for the GS algorithm. The GS algorithm reduces evenly the MSE based on the amplitude distribution of the reconstructed image, while the proposed algorithm reduces the MSE defined for the intensity distribution. Therefore, the GS algorithm can uniformly control the amplitude distribution, whereas, the proposed technique can evenly influence the intensity distribution. Because the GS algorithm cannot impact the intensity distribution uniformly, it can be used for the more precise control of low intensities but not for high intensities. Therefore, the GS algorithm could control the black region more strictly, while the proposed technique could manage the white region more accurately.

The effectiveness of the proposed approach for determining the step length coefficient was verified. Figure 9a–c show changes in the MSE during the optimization processes when constant step length coefficients were used ($\gamma$ = 0.20, 0.25, and 0.30) and when an optimized step length coefficient was used. When constant step length coefficients were used, the decrease in the MSE was dependent on the test image. When an optimized step length coefficient was used, the MSE was minimized for all test images. Figure 9d–f show the changes in the coefficient $\gamma$ corresponding to the results in Figure 9a–c.

5. Experiments

The proposed technique was experimentally verified utilizing the experimental system shown in Figure 10. A PLUTO-2 phase-only SLM (HOLOEYE Photonics AG) was utilized. The resolution was 1920 × 1080, the pixel pitch was 8 $\mathsf{\mu }$m, and the active area size was $15.36\times 8.64$ mm. A He-Ne laser with a wavelength of 633 nm was used as the light source. The focal length of the collimator lens was 150 mm, and that of the Fourier transform lens was 100 mm. The reconstructed images were captured using the digital camera (Canon EOS kiss X9 camera, Canon Inc.) The resolution was 6000 × 4000 and the sensor size was $22.3\times 14.9{\mathrm{mm}}^2$.

Figure 11a–c show the phase distributions used for the experiments, which correspond to the reconstructed images shown in Figure 6a–c and Figure 11d–f are the phase distributions correspond to Figure 6d–f. The experimentally obtained reconstructed images by proposed algorithm are shown in Figure 12a–c, and the reconstructed images by GS algorithm are shown in Figure 12d–f.

From Figure 12, the test images were reconstructed. Figure 12a–c show the reconstructed images obtained using the phase-only holograms calculated by the proposed method, while Figure 12d–f illustrate those obtained using the holograms calculated by the GS algorithm. From these results, the intensity variations in the white regions were smaller for the proposed algorithm than the GS algorithm. However, the difference was not so obvious because of some noises caused by the optical system and speckles in the reconstructed images. The reconstructed images contained the speckles because the phase distributions of the reconstructed images were not controlled. Although the images “hikari” and “circle” were reconstructed clearly, the image “mandrill” was not as clear as them. Because the “hikari” and “circle” were binary images, the image contrast was not affected so much by the speckles. On the contrary, because the “mandrill” was a gray-scale image, the image contrast was heavily affected.

6. Conclusions

We proposed a new technique for calculating phase-only CGHs using the gradient descent method, which minimizes the MSE defined for the intensity distributions of reconstructed images. We confirmed that the proposed technique optimizes the CGH phase distribution more rapidly than the GS algorithm and provides higher-precision reconstructed images. The proposed technique was also experimentally verified. The proposed technique could control the intensity of the reconstructed images evenly, because the MSE was directly defined for intensity distribution. The GS algorithm could control lower intensities more precisely than the proposed technique because it focused on controlling the amplitude distribution. For the binary images, the GS algorithm could manage the black region more accurately, but the proposed technique could control the white region more precisely. The proposed algorithm based on the gradient descent method can be utilized for the optimization of the Fresnel-type phase-only holograms.