Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Computer-generated holography in the intermediate domain

Open Access Open Access

Abstract

Iterative Fourier transform algorithms are widely used for hologram generation for phase-modulating spatial light modulators. In this paper, we introduce a new technique called the “intermediate domain,” which decomposes the Fourier transforms used into multiple subtransforms, the combination of which can offer major performance benefits over traditional approaches. To demonstrate this, we introduce ID-GS, an implementation of the intermediate domain technique for possibly the best known hologram generation algorithm, Gerchberg–Saxton. We discuss the performance of this across a wide range of configurations with a focus on computational performance.

Published by Optica Publishing Group under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

1. INTRODUCTION

Phase retrieval problems are commonly encountered in a wide range of areas, including x-ray crystallography [13], telescopy [46], and microscopy [79]. In these contexts the experimentalist is able to measure the intensity of an incident wave but is also interested in its phase. This problem becomes tractable if the intensity of the incident wave is measured in two different planes, on the proviso that the transformation undergone by the wave when traveling from one plane to the next is well known. For example, the intensity of an x-ray beam can be measured, left to propagate to the far field (equivalent to taking the two-dimensional Fourier transform of the field), and its intensity measured again. Knowledge of the field intensity before and after propagation allows the problem to be sufficiently constrained so as to allow the phase of the field to be inferred, although the solution may not be unique.

In a modern holographic projector, a beam of incident coherent light is modulated by a hologram displayed on a spatial light modulator (SLM) to form the diffraction field, which propagates through the optical system to form the replay field. Many commonly used SLMs can only modulate the phase of the incident field, and so the amplitude of the diffraction field must equal that of the illuminating beam. Propagation through the optical system is described by a linear transform ${\cal F}$ relating the hologram field $H$ to the replay field $R$ [Eq. (1)]. The desired target field $T$ is defined in the replay field plane, and its counterpart $G$ is defined in the hologram plane. $G$ and $T$ are also related by the linear transform ${\cal F}$ [Eq. (1)]. We denote the amplitude of the beam illuminating the SLM as $I$ (see Fig. 1):

 figure: Fig. 1.

Fig. 1. Computer-generated holography coordinate systems.

Download Full Size | PDF

$$G\underset{{{\mathcal{F}}^{-1}}}{\overset{\mathcal{F}}{\mathop{\rightleftarrows }}}\,T,\quad H\underset{{{\mathcal{F}}^{-1}}}{\overset{\mathcal{F}}{\mathop{\rightleftarrows }}}\,R.$$

The aim of computer-generated holography (CGH) is then to determine what hologram displayed on the SLM yields a replay field that resembles the target field as closely as possible. In this manner, the problem can be cast in the guise of a phase retrieval problem. CGH holds promise for a variety of applications, including virtual reality [10], augmented reality [11], lithography [12], and optical tweezing [13]. Although holography is the specific application considered in this paper, the results presented are also more widely applicable to other fields where phase retrieval problems are encountered.

In this manner, the problem is cast in the guise of a phase retrieval problem. Although holography is the specific application considered in this paper, the results presented are more widely applicable to fields where phase retrieval problems are encountered.

One of the earliest algorithms addressing the phase retrieval problem was introduced by Gerchberg and Saxton in 1972 [14]. The proposed algorithm has since become a staple for tackling phase retrieval problems and will be referred to here as the GS algorithm. GS operates by alternately projecting between the hologram and replay field planes as shown in Fig. 2 (note the use of the $\angle$ operator here, which returns a unit phasor of given angle, i.e, $\angle X \equiv X/| X |$). After each projection, the amplitude of the hologram or the replay field is constrained, but the phase profile is conserved. It has been shown that, in this manner, the algorithm converges on a solution where $H$ and $R$ both adhere to their known magnitude constraints [15].

 figure: Fig. 2.

Fig. 2. Gerchberg–Saxton algorithm.

Download Full Size | PDF

This paper introduces, in Section 2, a groundbreaking new approach to generating holograms based on a novel concept known as the intermediate domain (ID). To aid comprehension, only the necessary concepts are introduced, and the underlying mathematical justification is deferred until Section 5. We then introduce, in Section 3, an example of using the ID for a new iterative algorithm that offers significant performance improvements over industry staples such as GS [14]. The performance envelope of this algorithm is presented in Section 4. The remaining sections deal with more in-depth aspects of the intermediate domain and its applications. Section 5 presents the mathematical justification for the assumptions previously made in Section 2 before Section 6 presents alternative formulations and applications. Finally, conclusions are drawn in Section 8.

2. INTERMEDIATE DOMAIN

In this section, we introduce the concept of the ID.

A. 2D Transform Decomposition

A 2D linear transform ${\cal F}$ is used to determine the replay field from the diffraction field. Depending on the optical configuration used, a wide array of linear transforms may describe the propagation of the wavefront, such as Fourier, Fresnel, and Fraunhofer transforms as well as less common transforms such as Gyrator [16] and Wavelet transforms [17]. We decompose ${\cal F}$ into partial transforms ${{\cal T}_1}$ and ${{\cal T}_2}$:

$${\cal F}\{H\} \equiv {{\cal T}_2}\{{{\cal T}_1}\{H\} \} ,\quad {\rm{and}}\quad {{\cal F}^{- 1}}\{R\} \equiv {\cal T}_2^{- 1}\{{\cal T}_1^{- 1}\{R\} \} .$$

By decomposing the 2D linear transformation into two separate partial transformations, we can introduce intermediate domains $C$ and $D$ as illustrated in Fig. 3:

$$G\underset{\mathcal{T}_{1}^{-1}}{\overset{{{\mathcal{T}}_{1}}}{\mathop{\rightleftarrows }}}\,C\underset{\mathcal{T}_{2}^{-1}}{\overset{{{\mathcal{T}}_{2}}}{\mathop{\rightleftarrows }}}\,T,\quad H\underset{\mathcal{T}_{1}^{-1}}{\overset{{{\mathcal{T}}_{1}}}{\mathop{\rightleftarrows }}}\,D\underset{\mathcal{T}_{2}^{-1}}{\overset{{{\mathcal{T}}_{2}}}{\mathop{\rightleftarrows }}}\,R.$$
 figure: Fig. 3.

Fig. 3. Relationship of the “intermediate domain” to the diffraction and replay fields, assuming a 2D discrete Fourier transform (DFT) decomposed into a column-wise DFT (${{\cal T}_1}$) and a row-wise DFT (${{\cal T}_2}$).

Download Full Size | PDF

Previously, the authors have shown that Parseval’s theorem entails that the replay field can be represented in the hologram plane and vice versa. This allows an error metric such as the phase-sensitive mean-squared error (${{\rm MSE}_{{\rm PS}}}$) to be calculated in any plane [18]. In this case, this allows both the diffraction field and the replay field to be represented in the intermediate domain without loss of generality.

B. Example Case

A 2D linear transform can be decomposed in an infinite number of ways, and the choice of ${{\cal T}_1}$ and ${{\cal T}_2}$ is arbitrary For simplicity, the bulk of this work will use the example of a far-field hologram where the sampled replay field is given by the discrete Fourier transform (DFT) of the hologram. ${{\cal T}_1}$ is chosen to be the DFT of image columns and ${{\cal T}_2}$ to be the DFT of image rows. The order of row and column Fourier transforms does not affect the combined result:

$$\mathop {{\rm{FFT}}}\limits_{2D} \equiv \mathop {{\rm{FFT}}}\limits_{{\rm columns}} \circ \mathop {{\rm{FFT}}}\limits_{{\rm rows}} \equiv \mathop {{\rm{FFT}}}\limits_{{\rm rows}}{\circ}\mathop {{\rm{FFT}}}\limits_{{\rm columns}} ,$$
where “$\circ$” represents function chaining. This particular decomposition is chosen as it closely represents the way that the DFT of 2D images is calculated in practice. This is not the only subdivision of the 2D DFT available, and we return to this in more detail in Section 6.B.

In the next section we introduce an algorithm that uses this ID concept to significantly improve algorithmic performance for a wide range of phase hologram cases.

Tables Icon

Algorithm 1. Intermediate Domain Gerchberg-Saxton

 figure: Fig. 4.

Fig. 4. ID-GS. The “$\bullet$” and “$\diamondsuit$” symbols represent combination and decision, respectively.

Download Full Size | PDF

3. ID-GS ALGORITHM

Defining the error metric in the intermediate domain means that, instead of updating the entire ${N^2}$ pixels of $H$ or $R$ in a single step followed by an $O({N^2}\log (N))$ fast Fourier transform (FFT), instead we can update $N$ pixels in either a row or column followed by an $O(N\log (N))$ FFT.

This allows us to define an intermediate domain Gerchberg–Saxton (ID-GS) algorithm shown in Algorithm 1 and Fig. 4. By iterating through all the rows and columns in a random permutation, we end up with an algorithm of identical algebraic/computational performance to GS. We show below that this improves convergence and reliability over other iterative Fourier transform algorithms (IFTAs), particularly in cases where they are known to struggle.

4. PERFORMANCE/RESULTS

ID-GS does not offer a significant advantage over GS in the standard case of a continuous-phase SLM, uniform illumination, and full target region of interest. Improvements of 5% to 10% are observed for continuous-phase SLM, uniform illumination, full target fields.

More significant improvements are observed in more constrained scenarios. It is noted that an SLM is a digital device and is only capable of discrete phase modulation levels. This is accounted for in both the GS and ID-GS algorithms by quantizing the hologram pixels to the nearest allowed value. Results for this are shown in Fig. 5 (left), where ID significantly outperforms GS after only a few iterations on a 16-level phase SLM (the Matlab code to produce this and other convergence graphs can be found in Demo_Graph.m found in Supplement 1). The effect on visual quality is shown in Fig. 6 (the Matlab code to produce this and other image comparisons can be found in Demo_Image.m found in Supplement 1).

 figure: Fig. 5.

Fig. 5. Convergence of ID-GS against GS for 16 level phase SLM without a region of interest (left) and a 60% fill factor region of interest (right). Error bars show 2 standard deviations. Values are taken as being the normalized mean of 20 runs per image of the six test images shown in Fig. 8 with independent random phase profiles. The converged images are shown in Figs. 6 and 7.

Download Full Size | PDF

 figure: Fig. 6.

Fig. 6. Example of the Mandrill test image, Fig. 8(a), calculated in Fig. 5 after 300 iterations of GS (left) against ID algorithm (right) for a 16-level phase SLM with no region of interest and uniform illumination.

Download Full Size | PDF

Next, a scenario is considered where a region of interest is incorporated into the replay field. In this case, even greater improvements are observed as shown in Fig. 5 (right). The effect on visual quality is shown in Fig. 7.

 figure: Fig. 7.

Fig. 7. Example of the Mandrill test image, Fig. 8(a), calculated in Fig. 5 after 300 iterations of GS (left) against ID algorithm (right) for a 16-level phase SLM with a 60% fill factor region of interest and uniform illumination.

Download Full Size | PDF

 figure: Fig. 8.

Fig. 8. Standard test images from the USC-SIPI image database [19]. From left to right: Mandrill, Peppers, Man, Camera Man, Aerial, Landscape.

Download Full Size | PDF

All of these runs are taken as the normalized mean of 20 runs per image in the set of test images in Fig. 8. The phase-insensitive mean-squared error (${{\rm MSE}_{{\rm PI}}}$) and $\rm SSIM$ error metrics used are given by Eqs. (5) and (6) [20], respectively. In both cases the sum is performed over the region of interest. For the structural similarity index (SSIM) metric, ${\mu _T}$ and ${\mu _R}$ are the window means, ${\sigma _T}$ and ${\sigma _R}$ are the window variances, and ${\sigma _{{\rm TR}}}$ is the covariance between the two windows. ${c_1}$ and ${c_2}$ are functions of pixel dynamic range $L$, such that ${c_1} = (0.01L{)^2}$ and ${c_2} = (0.03L{)^2}$:

$${{\rm{MSE}}_{{\rm{PI}}}}(T,R) = {\left\| {\left| T \right| - \left| R \right|} \right\|^2},$$
and SSIM error metric given by [21]
$${\rm{SSIM}}(T,R) = \underbrace {\frac{{\left({2{\mu _T}{\mu _R} + {c_1}} \right)}}{{\left({\mu _T^2 +\mu_R^2 + {c_1}} \right)}}}_{{S_1}}\underbrace {\frac{{\left({2{\sigma _{\textit{TR}}} + {c_2}} \right)}}{{\left({\sigma _T^2 + \sigma _R^2 + {c_2}} \right)}}}_{{S_2}}.$$

In both cases, the sum is performed over the region of interest. For the SSIM metric, ${\mu _T}$ and ${\mu _R}$ are the window means, ${\sigma _T}$ and ${\sigma _R}$ are the window variances, and ${\sigma _{{\rm TR}}}$ is the covariance between the two windows. ${c_1}$ and ${c_2}$ are functions of pixel dynamic range $L$, such that ${c_1} = (0.01L{)^2}$ and ${c_2} = (0.03L{)^2}$.

5. JUSTIFICATION

Convergence of GS is proven in [15], where it is demonstrated that the ${{\rm MSE}_{{\rm PS}}}$ between $R$ and $T$ is reduced at each iteration. Within the same iteration, the error between $G$ and $I$ is also reduced. By a similar argument, it can be shown that a single iteration of ID-GS reduces the ${{\rm MSE}_{{\rm PS}}}$ between either $R$ and $T$ if a row $q$ is chosen, or between $G$ and the $I$ if a column $p$ is chosen. We provide a proof of this in Appendix A.

As shown by our results, this additional flexibility favors ID-GS. The algorithm’s branching at the point where a row or a column is selected adds more potential for refinement. Heuristically, it makes sense to select more (replay-side) rows than (hologram-side) columns during the course of the algorithm to, one, reduce the overall amount of nonlinear phase quantization steps, in cases where SLM quantisation is used, and two, place greater emphasis on the side of the optimization where there are more complex amplitude profiles, the target plane.

6. FURTHER CONSIDERATIONS

A. Parallelization

Perhaps the most apparent advantage of ID-GS is the ease of parallelization. Rows can be updated independently of each other; columns likewise. Changes in one type only effect future modifications of the other type. This allows for multi-core systems to operate on multiple rows or columns in parallel provided the data is synced before moving onto the other type. This offers significant benefits in multi-core systems where the problem can be broken down more straightforwardly than with GS.

B. Alternative Decompositions

In Section 2.B we defined the intermediate domain as being the column-wise DFT of the diffraction field and row-wise IDFT of the replay field. Obviously, problem symmetry allows definition of the ID as being the row-wise DFT of the diffraction field and column-wise IDFT of the replay field. The question arises, therefore, of what other decompositions are available.

 figure: Fig. 9.

Fig. 9. Efficiency savings in the ID using a region of interest.

Download Full Size | PDF

Given two invertible, orthogonal transforms ${{\cal T}_1}$ and ${{\cal T}_2}$, the requirements are that we can write

$${\cal F}\{H\} \equiv {{\cal T}_2}\{{{\cal T}_1}\{H\} \} ,\quad {\rm{and}}\quad {{\cal F}^{- 1}}\{R\} \equiv {\cal T}_2^{- 1}\{{\cal T}_1^{- 1}\{R\} \}$$
and that both ${{\cal T}_1}$ and ${{\cal T}_2}$ follow Parseval’s theorem, i.e.,
$$\sum a\overline b = \sum {{\cal T}_1}\{a\} \overline {{{\cal T}_1}\{b\}} ,\quad {\rm{and}}\quad \sum a\overline b = \sum {{\cal T}_2}\{a\} \overline {{{\cal T}_2}\{b\}} .$$

The most obvious solution to these is the row-/column-wise transforms used earlier, but this could be split in many other ways. 2D FFT algorithms exist for block-/stride-wise approaches that also offer $O({N_r}{N_c}\log ({N_r}{N_c}))$ performance.

C. All Intermediate Domains Are Created Equal

Provided the conditions of Section 6.B are kept, the choice of intermediate domain is arbitrary. Phase-sensitive MSE is the same, independent of the choice of ID ${{\cal T}_1}$ and ${{\cal T}_2}$.

This result is important, as it means that the ID can be placed anywhere in the transform, even at the ends. As suggested by Section 5, choosing an ID with approximately similar size transforms in each direction resulted in the best computational performance.

D. Multiple Intermediate Domains Are Possible

Another observation is that it is possible to use multiple IDs. For example, we could redefine Eq. (2) as

$${\cal F}\{H\} \equiv {{\cal T}_3}\{{{\cal T}_2}\{{{\cal T}_1}\{H\} \} \} ,\quad {\rm{and}}\quad {{\cal F}^{- 1}}\{R\} \equiv {\cal T}_3^{- 1}\{{\cal T}_2^{- 1}\{{\cal T}_1^{- 1}\{R\} \} \} .$$

This could be used to keep all transforms to smaller than $\sqrt {{N_r}{N_c}}$. The authors would be interested to explore whether this could be used to offer benefits sufficient to justify the extra overhead.

E. Fresnel Holograms, Aberration Correction, and 3D

Not only can we define alternative decomposition strategies as in Section 6.B, we can also choose transforms other than Fourier transforms provided they meet the requirements. Perhaps the most common alternative transform is the Fresnel transform, which adds a phase quadratic to both sides of the relationship

$${R_{u,v}} = \mathop {\cal F}\limits_{{\rm{Fresnel}}} \{{H_{x,y}}\} = {\Phi _1}\mathop {\cal F}\limits_{{\rm{Fourier}}} \{{H_{x,y}}{\Phi _2}\} ,$$
where ${\Phi _1} = \exp \frac{{i\pi}}{{\lambda z}}({x^2} + {y^2})$ and ${\Phi _2} = \exp \frac{{i\pi}}{{\lambda z}}({u^2} + {v^2})$.

We can also include linear modifiers to transforms such as Zernike polynomials [22] and Seidel aberrations [23]. Finally, while this paper discussed 2D holograms, its concepts are equally applicable to any iterative technique, including those for 3D holography.

F. Effect of Region of Interest on Performance

ID-GS is algebraically identical to GS in terms of number of arithmetic operations when there is no region of interest. Adding a region of interest, however, offers differential benefits for ID-GS over GS. The reason for this is that, in GS, the whole field is always required to be transformed. In ID-GS, however, only the rows of $R$ that contain constrained pixels require transforming. This is schematically shown in Fig. 9. For example, in Fig. 5 (right) and Fig. 7, 22% of target rows did not ever require transforming, with an overall time saving of an additional 11% over GS. This can be improved upon as well. For example, using a block/stride decomposition as discussed in Section 6.B saved over 40% of transforms between $D$ and $R$ with a total runtime saving of over 20%.

 figure: Fig. 10.

Fig. 10. Region of interest and supersampling regimes: (a) no supersampling, no region of interest; (b) no supersampling, region of interest; and (c) supersampling, no region of interest.

Download Full Size | PDF

 figure: Fig. 11.

Fig. 11. Efficiency savings in the ID with supersampling.

Download Full Size | PDF

G. Effect of Supersampling on Performance

Similar to the region of interest case, CGH is often carried out with supersampling where a larger target image is used than the available SLM pixels. This concept is shown in Fig. 10(c) in comparison to the basic case in Fig. 10(a) and the region of interest case Fig. 10(b).

 figure: Fig. 12.

Fig. 12. Convergence of ID-GS with a 1:1 ratio between columns and rows and a 3:1 ratio between columns and rows. Values shown are taken as the normalized MSE/SSIM. Error bars show 2 standard deviations. Values are taken as being the mean of 20 runs per image of the six test images shown in Fig. 8 with independent random phase profiles.

Download Full Size | PDF

Unlike the region of interest case, where the area outside the region of interest is allowed to assume any value, in the supersampling case it must be held to zero. This can be done for many columns in the intermediate domain without reference back to the SLM plane. Figure 11 shows this in action.

H. Hologram/Replay Preference

So far we have assumed that rows and columns are equally weighted in the algorithm. Intuitively, however, the greater complexity of the amplitude profiles in the replay fields suggests that a preferential treatment could be given to the row transforms on the replay side over the column transforms on the hologram side. Figure 12 shows the effect of taking a weighting chance of 3:1 between rows and columns when compared to a 1:1 weighting chance. This gives better convergence in the longer term but worse behavior in the short term. The selection priority is given by

$$\begin{split}{\rm{row}} &= {\rm randi}\,\,[1,{N_r}]\,\,{\rm{where}}\,\,{\rm rand}\,[ 0 ,1] \\&\lt \frac{{w{N_r}}}{{w{N_r} + {N_c}}}{\rm{otherwisec}}\,\,{\rm{olumn}} = {\rm randi}\,[1,{N_c}],\end{split}$$
where $w$ represents the weighting value for preferring rows over columns.

Even greater improvements were found using more complex probability profiles, which favored replay-side rows at the start of the algorithm and hologram-side columns toward the end. This could nearly double the performance differences in Fig. 12 but required additional overhead and loss of generality. The authors believe that the primary factor here is the variation in complexity between target replay and hologram amplitudes, but we were unable to fully develop a theory predicting this behavior. We will pursue this further in our research.

I. Randomization Strategy

With the exception of the previous section, this paper has assumed that each random selection of row or column is independently chosen from a uniform distribution of the available options. We found a number of significant factors that further improved the algorithm at the expense of additional overhead.

The first is to define a permutation of the ${N_r} + {N_c}$ rows and columns for each of the $N$ top-level iterations and to iterate through that so that each row and column is updated exactly $N$ times. This gave better performance than the described approach of randomly selecting from all available options $N({N_r} + {N_c})$ times with no guarantee of how many times each row or column will be updated. This comes at the expense of additional overhead. In some cases, we found that this could be further improved by using the same permutation order in each of the $N$ top-level iterations, though we failed to determine the exact cause for this.

The independence of row updates from each other and column updates from each other means that the randomization works best when rows and columns are alternated in the process. This comes at the expense of requiring a more advanced randomization algorithm.

Finally, we found evidence that weighting the rows and columns according to importance of local error improved convergence but found no evidence that this was sufficient to merit the additional error measurement calculation step it required.

J. ID Location Preference

One less obvious corollary of the work here is that while a centrally located ID offers the best performance, the decrease in performance away from the center is symmetric. ID converges equally well for an image that is $1024 \times 8$ as it does for an image that is $8 \times 1024$. This is in contrast to the performance improvements in weighting in favor of replay-side rows over hologram-side columns as discussed in Section 6.H and is due to the fundamental equality of various IDs as discussed in Section 6.C.

7. FURTHER WORK

In this work we have introduced the ID concept and one algorithmic implementation, ID-GS, that uses this concept. Operation in the ID offers lots of exciting opportunities, as it allows the hologram generation problem to be split into smaller, more manageable chunks. In this paper we have only discussed applying the ID technique to the GS algorithm as probably the best known hologram generation algorithm in the field. A wide array of alternative algorithms have been developed since 1972 when GS was first published, however, and the authors propose a more thorough exploration of how ID may be applied to newer algorithms including Yang-Gu [24,25], FIDOC [26,27], and Fienup [15,28].

8. CONCLUSION

This paper has introduced a new concept in computer-generated holography, that of the intermediate domain (ID). This offers a number of exciting opportunities for faster generation algorithms with improved convergence. We demonstrated the ID concept in action by presenting a variation of the well known Gerchberg–Saxton (GS) algorithm known as ID-GS.

ID-GS is comparable to GS in the standard use case of a continuous-phase SLM, uniform illumination, and full target region of interest. This changes for more realistic use cases including with discrete levelled SLMs, non-uniform illumination, regions of interest, supersampling, and non-Fourier operations where ID can offer convergent errors significantly lower than GS and other IFTA algorithms.

APPENDIX A

The forward and backward branches of the ID-GS algorithm can be analyzed symmetrically. For any given point in the algorithm, let $D$ be the intermediate domain field at the start, and let $D^\prime $ be the field obtained by performing a single iteration of ID-GS. The relevant steps of Algorithm 1 can then be written in parallel as follows:

$$\begin{split}\begin{array}{*{20}{l}}H &\leftarrow {\cal T}_1^{- 1}\{D\}\qquad R \leftarrow {{\cal T}_2}\{D\} , \\ H^\prime &\leftarrow |I|\angle H\qquad\,\, R^\prime \leftarrow |T|\angle R, \\ D^\prime &\leftarrow {{\cal T}_1}\{H^\prime \}\qquad D^\prime \leftarrow {\cal T}_2^{- 1}\{R^\prime \} .\end{array}\end{split}$$

We can then write the following steps:

$$\begin{split}\begin{array}{*{20}{l}}{\rm{MSE}} &= \left\| {|{\cal T}_1^{- 1}\{D^\prime \} | - |I|} \right\|\qquad\quad{\rm{MSE}} = \left\| {|{{\cal T}_2}\{D^\prime \} | - |T|} \right\| \\\\ &\le \left\| {{\cal T}_1^{- 1}\{D^\prime \} - |I|\angle {\cal T}_1^{- 1}\{D\}} \right\|\qquad \le \left\| {{{\cal T}_2}\{D^\prime \} - |T|\angle {{\cal T}_2}\{D\}} \right\|\end{array}\end{split}$$
$$\begin{split}\quad\,\,\,\begin{array}{*{20}{l}}&\le \left\| {{\cal T}_1^{- 1}\{D\} - |I|\angle {\cal T}_1^{- 1}\{D\}} \right\| \qquad\, \le \left\| {{{\cal T}_2}\{D\} - |T|\angle {{\cal T}_2}\{D\}} \right\|\end{array}\end{split}$$
$$\begin{split}\,\,\,\,\begin{array}{*{20}{l}}&= \left\| {|{\cal T}_1^{- 1}\{D\} | - |I|} \right\|\qquad \qquad= \left\| {|{{\cal T}_2}\{D\} | - |T|} \right\|.\end{array}\end{split}$$

Step (A2) follows from application of the triangle inequality, while (A3) follows from the definition of $D^\prime $. The above analysis shows that a single iteration of the ID-GS algorithm improves (or leaves unchanged) error performance in the hologram plane (left side above) or in the target plane (right side above). Concretely, the intermediate domain algorithm can be thought of as relaxing the requirement of a conventional GS algorithm of simultaneous error reduction in both the input and output planes. Instead, ID-GS reduces and counterbalances the error in each plane sequentially. The analysis supports the heuristic intuition that the ID approach has more flexibility than the typical GS algorithm.

Further use of the triangle inequality provides an intuitive interpretation of the relationship of hologram plane, ID, and target plane in the context of the algorithm:

$$\begin{split}\left\| {|{\cal T}\{H^\prime \} | - |T|} \right\| &\le \left\| {{\cal T}\{H^\prime \} - |T|\angle {{\cal T}_2}\{D\}} \right\| \\ &\le \left\| {{\cal T}\{H^\prime \} - {{\cal T}_2}\{D\}} \right\| + \left\| {|{{\cal T}_2}\{D\} | - |T|} \right\|.\end{split}$$

Of the two terms on the right-hand side of the inequality in Eq. (A5), iterations on the hologram side of the ID reduce the first, and iterations on the target side of the ID reduce the second. Hence iterations on either side of the ID can be thought of as constraining this composite upper bound on the distance between target and image.

Funding

Engineering and Physical Sciences Research Council (EP/T008369/1, EP/V055003/1); Biotechnology and Biological Sciences Research Council (BB/M011194/1).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Supplemental document

See Supplement 1 for supporting content.

REFERENCES

1. R. P. Millane, “Phase retrieval in crystallography and optics,” J. Opt. Soc. Am. A 7, 394–411 (1990). [CrossRef]  

2. R. W. Harrison, “Phase problem in crystallography,” J. Opt. Soc. Am. A 10, 1046–1055 (1993). [CrossRef]  

3. J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of X-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature 400, 342–344 (1999). [CrossRef]  

4. C. Fienup and J. Dainty, “Phase retrieval and image reconstruction for astronomy,” in Image Recovery: Theory and Application (1987), Vol. 231, p. 275.

5. J. R. Fienup, J. C. Marron, T. J. Schulz, and J. H. Seldin, “Hubble space telescope characterized by using phase-retrieval algorithms,” Appl. Opt. 32, 1747–1767 (1993). [CrossRef]  

6. J. E. Krist and C. J. Burrows, “Phase-retrieval analysis of pre-and post-repair hubble space telescope images,” Appl. Opt. 34, 4951–4964 (1995). [CrossRef]  

7. W. Coene, G. Janssen, M. O. 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]  

8. J. Zuo, I. Vartanyants, M. Gao, R. Zhang, and L. Nagahara, “Atomic resolution imaging of a carbon nanotube from diffraction intensities,” Science 300, 1419–1421 (2003). [CrossRef]  

9. Y. Zhang, G. Pedrini, W. Osten, and H. J. Tiziani, “Phase retrieval microscopy for quantitative phase-contrast imaging,” Optik 115, 94–96 (2004). [CrossRef]  

10. A. Maimone, A. Georgiou, and J. S. Kollin, “Holographic near-eye displays for virtual and augmented reality,” ACM Trans. Graph. 36, 85 (2017). [CrossRef]  

11. Z. He, X. Sui, G. Jin, and L. Cao, “Progress in virtual reality and augmented reality based on holographic display,” Appl. Opt. 58, A74–A81 (2019). [CrossRef]  

12. R. M. Von Bunau, H. Fukuda, and T. Terasawa, “Phase retrieval from defocused images and its applications in lithography,” Jpn. J. Appl. Phys. 36, 7494–7498 (1997). [CrossRef]  

13. J. Liesener, M. Reicherter, T. Haist, and H. J. Tiziani, “Multi-functional optical tweezers using computer-generated holograms,” Opt. Commun. 185, 77–82 (2000). [CrossRef]  

14. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35, 237–246 (1972).

15. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef]  

16. J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Gyrator transform: properties and applications,” Opt. Express 15, 2190–2203 (2007). [CrossRef]  

17. M. Antonini, M. Barlaud, P. Mathieu, and I. Daubechies, “Image coding using wavelet transform,” IEEE Trans. Image Process. 1, 205–220 (1992). [CrossRef]  

18. P. J. Christopher, R. Mouthaan, M. E. Guendy, and T. D. Wilkinson, “Linear-time algorithm for phase-sensitive holography,” Opt. Eng. 59, 085104 (2020). [CrossRef]  

19. USC-SIPI Database, 2020, https://sipi.usc.edu/database/.

20. A. Hore and D. Ziou, “Image quality metrics: PSNR vs. SSIM,” in 20th International Conference on Pattern Recognition (IEEE, 2010), pp. 2366–2369.

21. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process. 13, 600–612 (2004). [CrossRef]  

22. M. Born and E. Wolf, Principles of Optics, 6th ed. (Pergamon, 1980).

23. D. Claus, J. Watson, and J. Rodenburg, “Analysis and interpretation of the Seidel aberration coefficients in digital holography,” Appl. Opt. 50, H220–H229 (2011). [CrossRef]  

24. L. Wang, B. Dong, and G. Yang, “Phase retrieval from two intensity measurements in an optical system involving nonunitary transformation,” Appl. Opt. 29, 3422–3427 (1990). [CrossRef]  

25. G. Yang, B. Dong, B. Gu, J. Zhuang, and O. K. Ersoy, “Gerchberg–Saxton and Yang–Gu algorithms for phase retrieval in a nonunitary transform system: a comparison,” Appl. Opt. 33, 209–218 (1994). [CrossRef]  

26. H. Akahori, “Spectrum leveling by an iterative algorithm with a dummy area for synthesizing the kinoform,” Appl. Opt. 25, 802–811 (1986). [CrossRef]  

27. A. Georgiou, J. Christmas, N. Collings, J. Moore, and W. Crossland, “Aspects of hologram calculation for video frames,” J. Opt. A Pure Appl. Opt. 10, 035302 (2008). [CrossRef]  

28. H. H. Bauschke, P. L. Combettes, and D. R. Luke, “Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization,” J. Opt. Soc. Am. A 19, 1334–1345 (2002). [CrossRef]  

Supplementary Material (1)

NameDescription
Supplement 1       Supplemental document.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (12)

Fig. 1.
Fig. 1. Computer-generated holography coordinate systems.
Fig. 2.
Fig. 2. Gerchberg–Saxton algorithm.
Fig. 3.
Fig. 3. Relationship of the “intermediate domain” to the diffraction and replay fields, assuming a 2D discrete Fourier transform (DFT) decomposed into a column-wise DFT (${{\cal T}_1}$) and a row-wise DFT (${{\cal T}_2}$).
Fig. 4.
Fig. 4. ID-GS. The “$\bullet$” and “$\diamondsuit$” symbols represent combination and decision, respectively.
Fig. 5.
Fig. 5. Convergence of ID-GS against GS for 16 level phase SLM without a region of interest (left) and a 60% fill factor region of interest (right). Error bars show 2 standard deviations. Values are taken as being the normalized mean of 20 runs per image of the six test images shown in Fig. 8 with independent random phase profiles. The converged images are shown in Figs. 6 and 7.
Fig. 6.
Fig. 6. Example of the Mandrill test image, Fig. 8(a), calculated in Fig. 5 after 300 iterations of GS (left) against ID algorithm (right) for a 16-level phase SLM with no region of interest and uniform illumination.
Fig. 7.
Fig. 7. Example of the Mandrill test image, Fig. 8(a), calculated in Fig. 5 after 300 iterations of GS (left) against ID algorithm (right) for a 16-level phase SLM with a 60% fill factor region of interest and uniform illumination.
Fig. 8.
Fig. 8. Standard test images from the USC-SIPI image database [19]. From left to right: Mandrill, Peppers, Man, Camera Man, Aerial, Landscape.
Fig. 9.
Fig. 9. Efficiency savings in the ID using a region of interest.
Fig. 10.
Fig. 10. Region of interest and supersampling regimes: (a) no supersampling, no region of interest; (b) no supersampling, region of interest; and (c) supersampling, no region of interest.
Fig. 11.
Fig. 11. Efficiency savings in the ID with supersampling.
Fig. 12.
Fig. 12. Convergence of ID-GS with a 1:1 ratio between columns and rows and a 3:1 ratio between columns and rows. Values shown are taken as the normalized MSE/SSIM. Error bars show 2 standard deviations. Values are taken as being the mean of 20 runs per image of the six test images shown in Fig. 8 with independent random phase profiles.

Tables (1)

Tables Icon

Algorithm 1. Intermediate Domain Gerchberg-Saxton

Equations (16)

Equations on this page are rendered with MathJax. Learn more.

GFF1T,HFF1R.
F{H}T2{T1{H}},andF1{R}T21{T11{R}}.
GT1T11CT2T21T,HT1T11DT2T21R.
FFT2DFFTcolumnsFFTrowsFFTrowsFFTcolumns,
MSEPI(T,R)=|T||R|2,
SSIM(T,R)=(2μTμR+c1)(μT2+μR2+c1)S1(2σTR+c2)(σT2+σR2+c2)S2.
F{H}T2{T1{H}},andF1{R}T21{T11{R}}
ab¯=T1{a}T1{b}¯,andab¯=T2{a}T2{b}¯.
F{H}T3{T2{T1{H}}},andF1{R}T31{T21{T11{R}}}.
Ru,v=FFresnel{Hx,y}=Φ1FFourier{Hx,yΦ2},
row=randi[1,Nr]whererand[0,1]<wNrwNr+Ncotherwisecolumn=randi[1,Nc],
HT11{D}RT2{D},H|I|HR|T|R,DT1{H}DT21{R}.
MSE=|T11{D}||I|MSE=|T2{D}||T|T11{D}|I|T11{D}T2{D}|T|T2{D}
T11{D}|I|T11{D}T2{D}|T|T2{D}
=|T11{D}||I|=|T2{D}||T|.
|T{H}||T|T{H}|T|T2{D}T{H}T2{D}+|T2{D}||T|.
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.