Abstract:
A method of retrieving phase information from input intensity information, representative of a target image, in which a Fourier transform is performed on data and the result used in forming a phase estimate, the phase estimate being inverse Fourier transformed, thereby producing magnitude and phase replay, and wherein not only is the phase reply component but also data derived from the magnitude replay component, iteratively fed back.

Description:
RELATED APPLICATIONS 
     This application is a Continuation of International Application No. PCT/EP2007/003973, which designated the United States and was filed on May 4, 2007, published in English, which claims priority under 35 U.S.C. §119 or 365 to Great Britain, Application No. 0609365.2, filed May 11, 2006. 
     The entire teachings of the above applications are incorporated herein by reference. 
    
    
     SUMMARY OF THE INVENTION 
     The present invention relates to a method of phase retrieval, a method for real-time holographic projection and an apparatus for producing real-time holograms. 
     A number of algorithms, many based on the Gerchberg Saxton algorithm, use Fourier transforms to derive phase information from a target image. Such phase information, when implemented on a spatial light modulator (SLM) or the like, can simulate physical kinoforms so that when the SLM is illuminated by collimated laser light, a replay field corresponding generally to the target image is provided. 
     A number of other algorithms exist for providing phase information. 
     The Gerchberg Saxton algorithm and derivatives thereof are often much faster than the other “non-Fourier transform” algorithms. However, the iterative Gerchberg Saxton algorithm lacks the quality of the other algorithms, such as direct binary search algorithms, especially where relatively low numbers of iterations have been performed. 
     The Gerchberg Saxton algorithm considers the phase retrieval problem when intensity cross-sections of a light beam, I A (x,y) and I B (x,y), in the planes A and B respectively, are known and I A (x,y) and I B (x,y) are related by a single Fourier transform. With the given intensity cross-sections, an approximation to the phase distribution in the planes A and B, Φ A (x,y) and Φ B (x,y) respectively, can be found by this method. The Gerchberg-Saxton algorithm finds good solutions to this problem by following an iterative process. 
     The Gerchberg-Saxton algorithm iteratively applies spatial and spectral constraints while repeatedly transferring a data set (amplitude and phase), representative of I A (x,y) and I B (x,y), between the spatial domain and the Fourier (spectral) domain. The spatial and spectral constraints are I A (x,y) and I B (x,y) respectively. The constraints in either the spatial or spectral domain are imposed upon the amplitude of the data set and the phase information converge through a series of iterations. 
     Either or both constraints may be the phase information and, in this case, it would be the amplitude information that is desired. 
     It is also known that the Gerchberg-Saxton algorithm may begin in either the spatial domain or the Fourier domain. 
     It is desirable to provide a method of phase retrieval which can be implemented in a way that provides convergence more rapidly than the prior art. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The present invention will be better understood after reading the following description in conjunction with the drawings, in which: 
         FIG. 1  shows a diagram of the Gerchberg-Saxton algorithm; 
         FIG. 2  shows a diagram of a derivative of the Gerchberg-Saxton algorithm; 
         FIG. 3  shows a first algorithm embodying the present invention; 
         FIG. 4  shows a second algorithm embodying the present invention; 
         FIG. 5  shows a third algorithm embodying the present invention. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     The foregoing will be apparent from the following more particular description of preferred embodiments of the invention, as illustrated in the accompanying drawings in which like reference characters refer to the same parts throughout the different views. The drawings are not necessarily to scale, emphasis instead being placed upon illustrating embodiments of the present invention. 
     Referring to  FIG. 1 , the Gerchberg-Saxton algorithm establishes a discrete phase distribution at an image plane  100  and a corresponding diffraction (Fourier) plane  110  from known discrete amplitude distributions at the image and diffraction planes respectively. The image plane and diffraction plane are related through a single Fourier transform since the latter is the far-field diffraction pattern of the former. With both the amplitude and phase information, full wavefront reconstruction is achieved at both locations. 
     The method iterates a series of operations and has an input and an output. A data set having a plurality of elements, each element containing amplitude and phase information, is received at the input. After completing an iteration, the method outputs a new approximation of the received data set, and this approximation forms the basis for the input to the next iteration. It is intended that each iteration is a better approximation than the last iteration. 
     Referring to  FIG. 1 , for an n th  iteration, phase information  182  from the previous (n−1) iteration is multiplied in multiplication block  130  by the spatially corresponding image plane amplitudes  102  from the target image  100 . The resultant data set  132  is input to processing block  140  which applies a fast Fourier transform (FFT) to provide a second data set  142  in the frequency domain. Processing block  150  extracts the phase information  152  resulting from the FFT. Phase information  152  is then multiplied in processing block  160  by the stored diffraction plane amplitudes  112  (which are by definition in the frequency domain) to provide a third data set  162 ; this is an estimate of the complex diffraction pattern (it now has phase as well as amplitude). The third data set  162  is inverse Fourier transformed by processing block  170  into a fourth data set  172  in the spatial domain. The phase information  182  of data set  172  is extracted by processing block  180 . The data set with phase information  182  and amplitude information  102  provides the complex input to the second iteration. When the phase distributions converge sufficiently, the algorithm stops. 
     In the absence of phase information from the preceding iteration, the first iteration uses a random phase generator  120  to supply phase information  122  as a starting point. 
     It is found that this algorithm provides convergence on phase distributions in the spatial and Fourier domains which give a good approximation to the sampled image and its Fourier transform. 
     A known modification of the Gerchberg-Saxton algorithm will now be described with respect to  FIG. 2 . This algorithm retrieves the phase distribution at the diffraction (Fourier) plane which gives rise to full wavefront reconstruction of the corresponding image, at an image plane, when illuminated by a suitable light beam and viewed through a Fourier lens (or inverse Fourier transformed). 
     Referring to  FIG. 2 , for an n th  iteration, the input data set  202  is in the Fourier domain. It consists of amplitude information and phase information. The amplitude information is equal to the amplitude information of the Fourier transform of the target image and the phase information in the frequency domain is from the previous (n−1) iteration. This input data set is inverse Fourier transformed by processing block  220  to produce a second data set  222  in the spatial domain. The amplitude information of the second data set  222  is set to unity by processing block  230  and the phase is quantized by processing block  240  to produce a modified data set  242  having unit magnitude. The modified data set  242  represents a phase distribution that approximates to the Fourier transform of the target image and can be used to reconstruct a phase-only holographic representation of the target image. Modified data set  242  is then Fourier transformed back into the frequency domain in processing block  260  and the phase information output from the block  260  is supplied as an input to processing block  270  which, in turn, supplies the input to the next iteration. 
     For the first iteration there is no phase information from any preceding iteration, and hence the first iteration uses a random phase generator  280  to supply a starting set of phase information to processing block  270 . 
     With each iteration, the algorithm outputs phase information having a Fourier transform R[x,y] (in the replay field) which is an approximation to T[x;y] (target image). The difference between the replay field and target image gives a measure of convergence of the phase information ψ[x,y] and is assessed by an error function. 
     A first embodiment of the present invention is shown in  FIG. 3 . The figure shows a modified algorithm which retrieves the phase information ψ[x,y] of the Fourier transform of the data set which gives rise to a known amplitude information T[x,y]  362 . Amplitude information T[x,y]  362  is representative of a target image (e.g. a photograph). The phase information ψ[x,y] is used to produce a holographic representative of the target image at an image plane. 
     Since the magnitude and phase are intrinsically combined in the Fourier transform, the transformed magnitude (as well as phase) contains useful information about the accuracy of the calculated data set. Thus, embodiments of the present invention provide the algorithm with feedback on both the amplitude and the phase information. 
     The algorithm shown in  FIG. 3  can be considered as having a complex wave input (having amplitude information  301  and phase information  303 ) and a complex wave output (also having amplitude information  311  and phase information  313 ). For the purpose of this description, the amplitude and phase information are considered separately although they are intrinsically combined to form a data set. It should be remembered that both the amplitude and phase information are themselves functions of the spatial coordinates x and y and can be considered amplitude and phase distributions. 
     Referring to  FIG. 3 , processing block  350  produces a Fourier transform from a first data set having magnitude information  301  and phase information  303 . The result is a second data set, having magnitude information and phase information ψ n [x,y]  305 . The amplitude information from processing block  350  is discarded but the phase information ψ n [x,y]  305  is retained. Phase information  305  is quantized by processing block  354  and output as phase information ψ[x,y]  309 . Phase information  309  is passed to processing block  356  and given unit magnitude by processing block  352 . The third data set  307 ,  309  is applied to processing block  356  which performs an inverse Fourier transform. This produces a fourth data set R n [x,y] in the spatial domain having amplitude information |R n [x, y]|  311  and phase information &lt;R n [x, y]  313 . 
     Starting with the fourth data set, its phase information  313  forms the phase information of a fifth data set, applied as the first data set of the next iteration  303 ′. Its amplitude information R n [x, y]  311  is modified by subtraction from amplitude information T[x,y]  362  from the target image to produce an amplitude information  315  set. Scaled amplitude information  315  (scaled by α) is subtracted from target amplitude information T[x,y]  362  to produce input amplitude information η[x,y]  301  of the fifth data set for application as first data set to the next iteration. This is expressed mathematically in the following equations:
 
 R   n+1   [x,y]=F ′{exp( iψ   n   ·[u,v ])}
 
ψ n   [u,v]=               F {η*exp( i             R   n   ,[x,y ])}
 
η= T[x,y ]−α(| R   n   [x,y]|−T[x,y ])

     Where: 
     F′ is the inverse Fourier transform; 
     F is the forward Fourier transform; 
     R is the replay field; 
     T is the target image; 
                  is the angular information;
     Ψ is the quantized version of the angular information; 
     η is the new target magnitude, η≧0; and 
     α is a gain element ˜1. 
     In this embodiment, the gain element α is predetermined based on the size and rate of the incoming target image data. 
     The algorithm shown in  FIG. 3  produces phase information ψ[x,y] used to reconstruct a complex wavefront. The complex wavefront gives rise to an accurate holographic representation of target image intensity pattern T[x,y] at a replay field. 
     A second embodiment of the present invention is shown in  FIG. 4 . This embodiment differs from the first in that the resultant amplitude information  402  from processing block  350  is not discarded. The target amplitude information  362  is subtracted from amplitude information  402  to produce a new amplitude information  404 . A multiple of amplitude information  404  is subtracted from amplitude information  362  to produce the input amplitude information  406  for processing block  356 . 
     In a third embodiment, the final image reconstruction quality is improved by increasing the number of pixels in the Fourier plane. Since the normal processing method is to use a Fast Fourier Transform (FFT), the number of pixels in the Fourier domain are increased to match the number of pixels in the spatial domain, however target image size will not be increased, with the image being padded with additional data. The same gain feedback method as the first embodiment may be used. 
     A fourth embodiment of the present invention is shown in  FIG. 5 . This embodiment is the same as the first except that the phase is not fed back in full and only a portion proportion to its change over the last two iterations is fed back. 
     The phase information            R n [x,y]  313  output by processing block  356  is not fed straight back into processing block  350  as in the first embodiment. The difference between the phase information output in the current  313  and previous  504  iterations (         R n [x,y]−         R n−1 [x,y]) are calculated to give new phase information  502 . A multiple, β, of phase information  502  is subtracted from the phase information R n−1 [x,y]  504  of the previous iteration to give new input phase information  506  which provides the phase input for processing block  350 . This may be expressed mathematically in the following equations:
 
 R   n+1   [x,y]=F′{iψ   n   [u,v ])}
 
ψ n   [u,v]=             F {η exp( I θ)}
 
η= T[x,y ]−α(| R   n   [x,y]|−T[x,y ])
 
Θ=           R   n−1   [x,y ]+β(&lt; R   n   [x,y]−             R   n−1   [x,y ])

     Where: 
     F′ is the inverse Fourier transform; 
     F if the forward Fourier transform; 
     R is the replay field; 
     T is the target image; 
                  is the angular information;
     Ψ is the quantized version of the angular information; 
     η is the new target magnitude; η≧0; 
     θ is the new phase angle to match the new target magnitude; 
     α is a gain element ˜1; and 
     β is ratio of phase acceptance ˜1. 
     In this way, the algorithm will use the amplitude and phase information to predict the future values of phase. This can significantly reduce the number of iterations required. The gain values α and β are chosen to provide optimized speed and quality performance when the algorithm is used to sequentially process a sequence of image frames at conventional video rates. 
     Embodiments of the present invention may be used to dynamically change the phase pattern on a SLM in response to the output of the algorithm. The algorithms have been optimized to dynamically output phase patterns which produce a sufficient quality holographic representation of a received intensity pattern as perceived by a human viewer. The received intensity patterns may be supplied by a convention video camera. In essence, such a system would comprise: a computer-controlled SLM; video-capture apparatus, such as a video camera and frame-grabber, to provide the source intensity patterns (target images); the phase-retrieval algorithm as described in the above embodiments; and a suitably-chosen light source for illuminating the SLM and reconstructing the holographic image. 
     The skilled person will appreciate that the algorithm is not sensitive to the source of the input target images, only the rate at which images are received. The algorithm is optimized to output phase patterns on the SLM leading to dynamic holographic image of “acceptable” quality, time-correlated with the input intensity patterns. Embodiments of the present invention make the optimum trade-off between the quality of holographic images and the speed at which an “acceptable” quality holographic image is produced for series of input intensity patterns arriving at conventional video rates, for example 50 or 60 Hz). 
     Embodiments of the present invention are suitable for producing real-time phase patterns on a SLM for purposes other than real-time holography. For example, the SLM may be dynamically modified using calculations from the algorithm, using any of the described embodiments, to redirect incoming EM waves in a predetermined direction. This may be useful in applications where it is advantageous to redirect arbitrary RF waves in specific directions such as towards particular antenna. 
     The gain factors α and β may be fixed values or dynamically changing with time in a predetermined manner or in response to the rate of convergence, for example. 
     Some embodiments of the present invention may sufficiently converge on a phase information in fewer iterations than other embodiments but the time taken to converge may increase owing to the greater computational time required for each iteration. 
     Embodiments of the present invention are particularly suitable for phase retrieval of MPEGs since this video format only records changes in the image from frame to frame and the feedback factor |R n [x,y]|−T[x,y] necessitates this calculation. Thus, computation time may be saved in the phase retrieval if this calculation has already been made. 
     Embodiments of the present invention find application in small projectors and head-up displays. 
     Embodiments may provide phase information at high rates—for example more than one image per frame rate, e.g. to allow for grey scale reproduction, or color reproduction. 
     The invention has now been described with reference to a number of embodiments thereof. The invention is not however to be restricted to features of the described embodiments. 
     While this invention has been particularly shown and described with references to preferred embodiments thereof, it will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the scope of the invention encompassed by the appended claims.