The phase distribution of coherent light is difficult to measure directly, because phase information is lost in the measured intensity. Phase retrieval is used to determine the phase distribution on one plane by the intensity distribution measured from several other planes, and thereby is widely used both for wavefront sensing and wavefront modulation. For wavefront sensing, as in x-ray crystallography, electron microscopy,^{1, 2} and astronomy,³ the intensity distributions on several output planes are directly measured, while the coherent light distribution on the input plane is “sensed” by a retrieval algorithm. For wavefront modulation, as in optical tweezers⁴ and 2-D video frame projection,^{5, 6} the coherent light on the input plane is directly “modulated” to generate a special image pattern on the output plane.

A phase-only spatial light modulator (POSLM) is widely used for wavefront modulation. Since only the phase distribution is modulated, leaving the amplitude intact, no energy is lost and thus very high diffraction efficiency is achieved. Recently, liquid crystal on silicon (LCOS) with $2\pi$ range of phase modulation was developed to be used as POSLM.⁷ By modulating the “pure” phase on LCOS, the amplitude on the output plane with arbitrary distribution was realized, even at high video frame rates.⁶

To generate a special image pattern on the output plane, the required phase distribution on the LCOS modulator is deduced by either iteration algorithms^{8, 9} or direct algorithms.¹⁰ These algorithms are based on scalar diffraction theory,¹¹ in which the coherent light propagates from one plane to another, and therefore they are efficient for deducing the light distribution in a 2-D plane. Modulation of coherent light in 3-D space will extend the application field of phase retrieval, for example, 3-D holographic displays¹² and 3-D manipulation of biomedical particles.⁴ 3-D light distribution can be realized by simultaneously modulating coherent light in several parallel 2-D planes. Two kinds of multiplane iteration methods have been introduced to modulate 3-D light distribution. One is based on consecutive propagation between multiplanes;^{13, 14} the other is based on the addition of complex distribution directly propagated from multiplanes.^{15, 16} However, to converge a sharp image, the distance between planes cannot be very small in the first category method, and the reconstructed image in the second category method sometimes looks dirty due to the addition of complex distribution propagated from several planes. The purpose of this study is to develop a multiplane iteration method that can modulate the intensity distribution on several closely spaced planes with sharp images.

As mentioned before, a 3-D light distribution can be discretized into several 2-D planes. The light distribution on each plane is efficiently calculated by fast Fourier transform (FFT) based on the scalar diffraction theory. Our iteration method is based on Gerchberg-Saxton’s iteration method.^{8, 9} To retrieve the pure phase distribution on a POSLM, forward and backward propagation is calculated from plane to plane, and the amplitude distribution on each output plane is modified according to the amplitude constraint, which is the image pattern that will be modulated on each plane.

Figure 1 shows the proposed iteration method to modulate images on two output planes by a POSLM. The choice of two output planes is arbitrary and does not affect the generality of the method and results. The method consists of four steps: 1. the light starts at the POSLM and propagates to the first plane $\left(P_1\right)$ ; 2. the amplitude distribution on plane $P_1$ is modified according to the amplitude constraint and then propagates to plane $P_2$ ; 3. the amplitude distribution on the final plane $P_2$ is inversely propagated to plane $P_1$ , after applying the amplitude constraint; 4. the light distribution on plane $P_1$ is inversely propagated to the POSLM, and finally the amplitude distribution on POSLM is modified to a uniform light distribution. During each step, the phase distribution on each plane is calculated by scalar diffraction from the previous plane, combined with the modified amplitude distribution, and is then propagated to the next plane.

Fig. 1

Scheme of proposed multiplane iteration algorithm.

The angular spectrum method,¹¹ based on the scalar diffraction theory, is used to propagate light between output planes, because it is suitable for light propagation over a short distance, and thus high resolution 3-D light distributions can be realized by discretizing in small steps. Moreover, the sampling space is invariant between two propagation planes. The distance between the POSLM and the first plane $\left(P_1\right)$ is not limited to very short distances. For propagation over a long distance, single FFT Fresnel diffraction is used. It should be noted that the sampling space during Fresnel diffraction is changed, and thus forward and inverse propagation should take place between the same two planes. Therefore, during inverse propagation, the light is first propagated back to plane $P_1$ in step 3 and then finally propagated back to the POSLM plane in step 4.

The numerical simulation is based on the scheme depicted in Fig. 1. The wavelength used in the simulation is $\lambda =5\times {10}^{-7}\mathrm{m}$ . The distance between the POSLM and the first plane $\left(P_1\right)$ is $L_1=0.5\mathrm{m}$ . The distances between the two output planes is $L_2=1\times {10}^{-2}\mathrm{m}$ . The sampling space on the output plane is $dx=dy=2\times {10}^{-5}\mathrm{m}$ , and the sampling numbers in horizontal and vertical directions are both $N=512$ . The sampling space on the POSLM plane is determined by Fresnel diffraction.

Figures 2 and 2 show the amplitude constraints that are used in the simulation on two output planes, i.e., planes $P_1$ and $P_2$ , respectively. They are generated by segmenting the original image in Fig. 2. The resolution of these three images is $256\times 256$ , and is padded to $512\times 512$ by zeros during the calculations to avoid circular convolution. Figure 2 shows the constraint area on plane $P_2$ . When light propagates to this plane, the amplitude distribution inside the constraint area (the white area) is substituted by the image content, and the amplitude distribution in the dark area is intact during propagation. The constraint area on plane $P_1$ is an inverse one that is the dark area in Fig. 2.

Fig. 2

(a) The original image used to generate constraint. (b) and (c) are amplitude constraints on planes $P_1$ and $P_2$ , respectively. (d) The constraint area on plane $P_2$ .

Figures 3 and 3 show the simulated amplitude distribution on planes $P_1$ and $P_2$ after 512 iterations based on our proposed method. Inside the constraint area on each plane, the amplitude distribution converges to a sharp image. At the same area on the adjacent planes, this sharp image becomes a smeared image due to diffraction. For comparison, Figs. 3 and 3 show the simulated amplitude distribution based on the method in Ref. 16. The edges look blurred due to the light distribution being close to the edge on one plane propagated into the other side on other plane. The convergence of two algorithms presented in Fig. 3 is also different. Our proposed algorithm shows a fast convergence due to the do-not-care area outside the constraint.⁶ The root mean square error (rms) of the calculated amplitude distribution on plane $P_2$ is lower than 10 after 4 times of iteration. In contrast, the rms for the method in Ref. 16 is higher than 35 after 512 times iteration on the same plane due to the blurred edges.

Fig. 3

The simulated amplitude distributions on planes $P_1$ and $P_2$ , respectively. (a) and (b) are the results by our method. (c) and (d) are the results by the complex distribution addition method.

As mentioned before, the output planes are placed close to each other in our simulations; this enables the generation of high resolution 3-D light distribution. However, they also display smeared images on the other planes. Fortunately, this is not an issue in some applications. For instance, when the 3-D light distribution is discretized into several 2-D planes and placed close to human eyes,¹² people will clearly perceive a sharp image on one plane, and the images on other planes look blurred due to accommodation. Another example is a multiplane projection: when we use a POSLM to project an image on an arbitrarily shaped surface, e.g., a sphere, this 3-D surface can be discretized into several 2-D planes. On each plane a sharp image is projected on one area and stops further propagation.

In conclusion, we develop a new iteration method to project a 3-D image on several output planes by one POSLM. The numerical simulation implies that this method is very suitable for modulation of 3-D images with very fine structure. The method described can be used to generate 3-D holographic displays, and also has potential application in optical tweezers.