46 J. Opt. Soc. Am. A/Vol. 2, No. 1/January 1985

Fast iterative solution to exact equations for the two-dimensional phase-retrieval problem

K. Chalasinska-Macukow* and H. H. Arsenault

Laboratoire de Recherches en Optique et Laser, D6partement de Physique, Universit6 Laval, Ste-Foy, Qu6becGIK 7P4, Canada

Received December 12, 1983; accepted September 18, 1984

A modification of a previously published two-dimensional phase-retrieval algorithm based on the sampling theo-rem [Opt. Commun. 47, 380 (1983)] is introduced. With an appropriate choice of the intermediate sampling grid,a significant reduction in the number of computations is achieved, because an independent system of equationsmay be solved for each line of the Fourier transform of the object. The results of computer simulations on realobjects are given. The algorithm converges rapidly.

INTRODUCTION

The problem of recovering phase from intensity measure-ments arises in a number of areas, such as astronomy, electronmicroscopy, wave-front sensing, holographic imaging, co-herence theory, and inverse-source problems. The phase-retrieval problem may be formulated in several ways, withvarious constraints, in one or two dimensions. We concen-trate here on the two-dimensional problem because it is a mostimportant case in optics.

The first practical iterative algorithm for phase retrievalwas introduced by Gerchberg and Saxton.' Fienup modifiedthe Gerchberg-Saxton algorithm and proposed a more pow-erful version.2 '3

The use of the sampling theorem to obtain an exact solutionfor the one-dimensional phase problem was proposed by Ar-senault and Lowenthal.4 Bates et al.5 8 also proposed analgorithm based on the sampling theorem in one and two di-mensions, using a special interpolating function.

In our previous paper, 9 we proposed a method for two-dimensional phase retrieval based on the Whittaker-Shannonsampling theorem.1 0 The method applies to band-limitedobjects that are equal to zero outside of a finite support. Al-though it is true that an object and its spectrum cannot beboth limited to a finite support, any object with a finitenumber M X M of degrees of freedom is completely deter-mined by M X M independent samples in the Fourierplane."11 2 So an object can both have a finite support andbe band limited, if the bandwidth is understood to mean thesupport of the M X M independent Shannon samples that canexactly describe the object by means of the sampling theorem.A simple example of such an object in one dimension is thefunction rect(x)(1 - cos 7rx), which is completely determinedby the spectral samples at the three frequencies (-1, 0, +1).The number of degrees of freedom is equal to the space-bandwidth product'3

M X M = 4LxLyumaxvmax, (1)

where LX, Ly are the finite support of the object and um,_ and

vmax are the maximum spatial frequencies of the nonzerosamples in the Fourier spectrum of the object.

The phase-reconstruction task in this method is accom-plished by solving a set of M X M simultaneous nonlinearequations with the phase values as unknowns. These equa-tions are exact equations, and any exact solution to this setof equations is an exact solution to the phase-retrieval prob-lem. There is usually more than one such exact solution, anda major unsolved problem related to this method is how to finda solution that satisfies additional given constraints; for in-stance, if the object is known not to have any negativevalues.

The algorithm presented in Ref. 9 converges quickly andyields accurate solutions. However, the method suffers froma defect common to all methods based on a system of equa-tions: It requires an M2 X M2 matrix to solve a system of MX M nonlinear equations for an object with M X M degreesof freedom. Also, much computing time is required to solvethe problem for large two-dimensional images. This timegrows exponentially with the number of data points.1 4 Also,the accuracy of the solution decreases with the number ofunknowns. This is a limitation of all methods, because thelarge matrix representing the equations becomes ill condi-tioned.

In the present paper, we propose to circumvent this problemand to reduce the size of the matrix and the computing timeby reducing the two-dimensional phase-retrieval problem toa set of one-dimensional problems. We show that the modi-fied algorithm also converges quickly and after a few iterationsyields accurate solutions. In the absence of any additionalconstraints, the method does not yield a unique solution.

THE REDUCTION OF A TWO-DIMENSIONALPROBLEM TO A SET OF ONE-DIMENSIONALPROBLEMS

We concentrate here on the two-dimensional problem. Letthe object be a two-dimensional function f(x, y) with a finitesupport:

0740-3232/85/010046-05$02.00 © 1985 Optical Society of America

K. Chaldsinska-Macukow and H. H. Arsenault

Vol. 2, No. 1/January 1985/J. Opt. Soc. Am. A 47

f(xy) = ff(x,y) forlxl I L,to for lxi >L~,

IYI SLy

IYI > Ly(2)

We take Lx = Ly = 0.5 below.A diffraction-limited optical system yields a band-limited

Fourier transform F(u, v) of this object with the maximumfrequency Umax, Vmax in the u and v directions, respectively.To use the sampling theorem, we have to know both the finitesupport of f(x, y) and its highest frequency, which determinesthe finite support of F(u, v). From the sampling theoremapplied in the Fourier-transform plane for a complex function,we have' 0

K MIA(u, v)j exp[iso(u, v)] = E E IAk,mI

k=-K m=-M

X exp(iknm)sinc(u - k)sinc(v -m),

* 0 4

0o 0

* 0 4

0o 0

0 0

(3)

where IA(u, v)I and so(u, v) are the modulus and the phase,respectively, at any point with coordinates (u, v) and IAk,mland S~km are the modulus and the phase, respectively, at theindependent sampling points (h, m). We assume that theobject has (2K -t 1) (2M + 1) degrees of freedqm. From Eq.(3), the intensity at point (u, v) is equal to

K MIA(u, V)12 = ' E ' IAkmI

k=-K m=-M

X sinc(u - k)sinc(v - m)cos k,m]

+1K M 2.+ E E IAk,mI sinc(u - k)sinc(v - m)sin bk,m

.k=-K m=-M

(4)

This is a nonlinear equation with (2K + 1)(2M + 1) unknownphases Pkm- Only the modulus of the Fourier transform, i.e.,IA(u, v)I and IAk,ml, is known. We see that the -unknownphase values are coded irito the intensity at points with (u, v)coordinates, because Eq. (4) gives the intensity values for any(u, v) in terms of the values at the sample points (k, m).

In order to search for unknown phase values; that is, to solvethe phase-retrieval problem, we form' a system of equations. 9

Two sampling grids are used: the set of (k, m) independentsampling points (represented by filled circles in Fig. 1) anda second set of intermediate points with (u, v) coordinates,which are different from those of the independent set. Westated in Ref. 9 that the choice of the second grid containingthe intermediate sample points is important. This choicedetermines the complexity of the system of equations. Theintermediate grid composed of the in-line intermediatesamples, 5' 9 shown in Fig. 1 as squares, yields an independentsystem of equations for each row or column and reduces thetwo-dimensional problem to a set of one-dimensjonal prob-lems. Only the in-line intermediate samples reduce thetwo-dimensional problem to a set of one-dimensional prob-lems. For this grid, the coordinates of the intermediate pointsare (u, m) or (k, v), and Eq. (4) becomes

IA(u, m)12 = FkE IAk,mI sinc(u - k)cos (Pk,mllk=-K

+ | E Ak,mI sinc(u - k)sin (Pkmk=-K

(5a)

for each row, where m = -M, . . . +M is the number of the row,

03

I

* 0 4

0l I

* 0 4

LV

I 0 6

0 0

3 0

0

0

L.J W U W U

0

I 0

0

S

0 0

* 0

Fig. 1. Fourier-plane sampling. Ns represent the intermediatesample points A(u, v), and O's represent the independent samplepoints Ak,m.-

and

M 2IA(k, v)12 = IAk,mI sinc(v - m)cos (Pkm

m=-M. M 12

+ E IAk,ml sinc(v - m)sin (°km] (5b)m=-M I

for each column, where k = -K... +K is the number of thecolumn.

For the kth column we have the system of nonlinear equa-tions

Fkl(Rk) = 0

Fkj(Rh) = 0

Fkp(Rk) = 0 J

(6)

where P = (2M + 1), Rk = (SOk,-M ... (Pk,M) is a (2M + 1)-dimensional vector of unknown phases (Pkjm, and

Fkj(Rk) = | E IAkmI sinc(v - m)cos (Pk,mIm=-M I

rM 12+ | IAk ml sinc(v - m)sin spkm _ -IA(k, v)12.

Um -M I

The same procedure yields a system of nonlinear equationsfor each row. Each system of equations has (2K + 1) or (2M+ 1) unknown SOk,m, and each system solves the problem onlyfor one row or for one column of the sampled Fourier trans-form.

Each system of equations has an infinite number of solu-tions of the type ePk,m = const. (k, m). To avoid this problem,we set (po,o = 0 and solve the system of equations for the cen-tral row or for the central column (i.e., the row with m = 0 orthe column with k = 0). The solution for this row or column

- - .

- t=1 M- io

A

I

II

K. Chalasinska-Maeukow and H. H. Arsenault

48 J. Opt. Soc. Am. A/Vol. 2, No. 1/January 1985

yields the values required to set the central values (Pk,o and(pO,m for the system of equations for each alternate row andcolumn, respectively.

This procedure allows the reduction of the two-dimensionalproblem to a set of one-dimensional problems. The two-dimensional phase problem with M X M degrees of freedommay be solved by using M independent one-dimensionalsystems of equations with M - 1 unknowns each.

DIGITAL VERIFICATION OF THE ALGORITHM

To test our method, two real objects with 5 X 5 samples andwith 25 X 25 samples were used. The amplitude lAk,ml of theindependent samples in the Fourier-transform plane for theseobjects was calculated digitally using the fast Fourier trans-form. The intermediate sampling set IA(k, v)l was obtainedby interpolation in the Fourier-transform plane.

To solve the system of equations, we used the Newton-Raphson iterative method' 5 and the Gaussian eliminationmethod with pivoting.'6

Let us compare the modified algorithm presented abovewith the algorithm introduced in Ref. 9. We use the same realobject with 5 X 5 samples, as in the previous paper. To solvethe system of equations by the Newton-Raphson iterativemethod, initial phase values are required. The initial phasevalues were chosen close enough to the true solution to avoidthe problem of multiple solutions. The initial phase valueswere the same for the two algorithms.

Table 1 is a comparison of the results for the two algorithms.The calculation times for 10 iterations of each algorithm, theerror E,, of the reconstructed phase, and the error EF of theamplitude of the Fourier-transform plane were compared.The errors E., and EF were defined as

E 2= 1 E ('Pk,. - PO,(2K + 1)(2M + 1) km

E [IIA(u, v)I2 - IA 0(u, V)12]2

EF 4 = UV

Figure 3 shows the image corresponding to the initial phasevalues. The part of the image having the smallest contrastis barely visible. Areas with higher contrasts have intensityfluctuations caused by the phase errors in the Fourier-trans-form plane.

The solution of the system of equations yields the recon-structed phase values. These reconstructed phase values,combined with known moduli after Fourier transformation,yield the reconstructed image of Fig. 4. The quality of theimage of Fig. 4 is equal to that of the original of Fig. 2.

Table 1. Comparison of Two Proposed MethodsBased on the Sampling Theorem

Systems Time for Phaseof 10 Iterations Error

Equations (sec) (rad) EF

One system of 94.6 1 X 10-2 1 X 10-3

two-dimensional equationsSet of systems of 7.2 6 X 10-4 2 X 10-4

one-dimensional equations

(7a)

(7b)

E IA(u, v)14Uu

where (Pk,m is the reconstructed phase value, Spk,mO is originalphase value, IA0 (u, v)I is the original modulus, and IA(u, v)Iis the calculated modulus after each iteration. For the ithiteration, IA (u, v)I is calculated using Eq. (5b), where 'Pk,m =

(Pk,mi, i.e., the phase SPkm in the independent points (k, m) isequal to the phase (Pk,mi obtained after ith iteration. We seethat, for those three parameters, the modified method is muchbetter.

The second object used to investigate our modified methodis shown in Fig. 2.17 This object contained 25 X 25 samples.It is practically impossible to solve this two-dimensionalphase-retrieval problem using a single system of nonlinearequations, because there are 625 unknown phase values. Thereduction of the two-dimensional problem to a set of one-dimensional problems allows the solution to be found bysolving 25 systems of equations with 24 unknown phase valueseach.

These systems were also solved using the Newton-Raphsoniterative method; the initial phase values were again takenclose enough to the true solution to avoid multiple solutions.

Fig. 2. Original image containing 25 X 25 samples.

Fig. 3. Image corresponding to the initial phase values.

K. Chalasinska-Macukow and H. H. Arsenault

Vol. 2, No. 1/January 1985/J. Opt. Soc. Am. A 49

image corresponding to the reconstructed

Table 2. Comparison of the Parameters of the InitialImage and the Reconstructed Image for the Modified

Method

Phase ErrorImages (rad) EF Image Error

Initial image 0.3 0.36 211.0Reconstructed 0.07 0.001 15.0

image

A

100

-2

I3'U

C

EF

0

0

0

0

.

0

0

In Table 2, the error El, of the phase, the error EF of theFourier-transform modulus at the intermediate points, andthe error of the image (defined similarly to the phase error)are given.

The advantage of the modified method is the fast conver-gence of the algorithm to an accurate solution. To analyzethe convergence of the algorithm, the error EF in the Fourierplane was calculated for each iteration, with the results shownin Fig. 5.

CONCLUSION

We have described a method to solve the phase-retrievalproblem when only the modulus of the Fourier transform andthe support of the object are known. The method is based onthe two-dimensional sampling theorem. To solve thephase-retrieval problem by this method, more than onesampling grid is required to form a system of nonlinearequations. The particular choice of the second grid composedof intermediate sample points yields an independent systemof equations for each row or column, and the two-dimensionalproblem is reduced to a set of one-dimensional problems. Thesize of the matrix required to solve each system of equationsis thus reduced, along with the number of computations. Wehave demonstrated also that the modified algorithm permitsthe solution of the two-dimensional phase-retrieval problems,which are practically impossible to solve using a single systemof nonlinear equations (with 625 unknown phase values). Themodified algorithm converges quickly and yields an accuratesolution. In the absence of any additional constraints, thesolution is not generally unique when the initial phase valuesare not close to the solution. This is not surprising when oneconsiders that we have a system of equations representing onlypartial information on an object.

As mentioned above, the question of how to take into ac-count additional given constraints remains the main unsolvedproblem with this method. For the object of Fig. 2, for in-stance, a constraint is that the object has only positive values.Another constraint may be that the object is binary. Suchconstraints are not included in the system of equations, whichhas other solutions that do not satisfy the constraint. Thetrue solution is found only if the initial values taken for theiteration are close enough to the true values. At the presenttime, there is no known way to find such initial conditions.One possible approach would be to combine this method withmaximum-entropy techniques, which have been successfulwith some problems involving constraints.

In our method, the sensitivity to noise was analyzed for anoise level smaller than 1%, and it was found that low-levelnoise has no influence on the quality of the phase recon-struction. The influence of higher noise levels is understudy.

ACKNOWLEDGMENT

Portions of this paper were delivered at October 1983 AnnualII I I I I I a I I | Meeting of the Optical Society of America in New Orleans,

5 10 Louisiana.

NUMBER OF ITERATIONSFig. 5. Convergence of the modified algorithm. The time of oneiteration for the image with 25 X 25 samples is 26 sec.

* On leave from the Institute of Geophysics, University ofWarsaw, Pasteura 7, 02-093 Warsaw, Poland.

Fig. 4. Reconstructedphase values.

K. Chalasinska-Macukow and H. H. Arsenault

50 J. Opt. Soc. Am. A/Vol. 2, No. 1/January 1985

REFERENCES

1. R. W. Gerchberg and W. 0. Saxton, "A practical algorithm forthe determination of phase from image and diffraction planepictures," Optik 35, 237-246 (1972).

2. J. R. Fienup, "Reconstruction of an object from the modulus ofits Fourier transform," Opt. Lett. 3, 27-29 (1978).

3. J. R. Fienup, "Phase-retrieval algorithm: a comparison," Appl.Opt. 21, 2758-2769 (1982).

4. H. H. Arsenault and S. Lowenthal, "La restitution de la phase apartir de mesures d'6clairement," C. R. Acad. Sci. (Paris) 269,518-521 (1969).

5. R. H. T. Bates, "Fourier phase problems are uniquely solvablein more than one dimension. I: underlying theory," Optik 61,247-262 (1982).

6. K. L. Garden and R. H. T. Bates, "Fourier phase problems areuniquely solvable in more than one dimension. II: one-di-mensional considerations," Optik 62, 131-142 (1982).

7. W. R. Fright and R. H. T. Bates, "Fourier phase problems areuniquely solvable in more than one dimension. III: computa-tional examples for two dimensions," Optik 62, 219-230(1982).

8. R. H. T. Bates and W. R. Fright, "Composite two-dimensional

K. Chalasinska-Macukow and H. H. Arsenault

phase-restoration procedure," J. Opt. Soc. Am. 73, 358-365(1983).

9. H. H. Arsenault and K. Chalasinska-Macukow, "A solution to thephase-retrieval problem using the sampling theorem," Opt.Commun. 47, 380-386 (1983).

10. J. W. Goodman, Introduction to Fourier Optics (McGraw-Hill,New York, 1968), p. 25.

11. G. Toraldo di Francia, "Resolving power and information," J. Opt.Soc. Am. 45, 497 (1955).

12. G. Toraldo di Francia, "Degrees of freedom of an image," J. Opt.Soc. Am. 59, 358 (1969).

13. J. B. DeVelis and G. 0. Reynolds, "Communication theory," inHandbook of Optical Holography, H. J. Caulfield, ed. (Academic,New York, 1979), p. 71.

14. W. 0. Saxton, Computer Techniques for Image Processing inElectron Microscopy (Academic, New York, 1978), p. 96.

15. K. E. Atkinson, An Introduction to Numerical Analysis (Wiley,New York, 1978), p. 92.

16. M. L. James, G. M. Smith, and J. C. Wolford, Applied NumericalMethods for Digital Computation with FORTRAN and CSMP(Harper & Row, New York, 1977), p. 178.

17. K. Chalasinska-Macukow and H. H. Arsenault, "A solution to thephase-retrieval problem using the sampling theorem," J. Opt. Soc.Am. 73, 1875 (1983).