Reconstruction and synthesis applications of an iterative algorithm
Reconstruction Reconstruction and and synthesis synthesis applications applications of an iterative algorithm R. Fienup J. R.
Abstract. This paper paper reviews reviews the the GerchbergGerchberg-Saxton Saxton algorithm and and variavaria-
Environmental Environmental Research Research Institute Institute of Michigan P.O. P.O. Box 8618 Ann Arbor, Arbor, Michigan 48107 48107 Ann
tions thereof that have have been been used used to solve a number number of difficult difficult reconreconstruction and and synthesis synthesis problems problems in in optics and and related related fields. It can can be be used information (including (including both both used on on any any problem problem in in which which only partial information measurements the wavefront wavefront or or signal signal isis available available in in measurements and and constraints) of the one information isis available available inin another another domain one domain domain and and other other partial partial information domain (usually (usually the the Fourier Fourier domain). domain). The Thealgorithm algorithm combines combines the the information in both arrive at of the the wavefront wavefront or both domains domains to arrive at the complete description of or signal. are reviewed, reviewed, including synthesis of Fourier signal. Various Various applications are Fourier transform pairs having having desirable desirable properties properties as as reconstruction transform pairs as well as reconstruction problems. and the the convergence convergence properties problems. Variations of the algorithm and properties of of the are discussed. the algorithm are
1. INTRODUCTION INTRODUCTION 1. There very difficult There exist exist many many problems problemsthat that are are very to solve difficult to solve in astronomy, x x-ray crystallography, electron -ray crystallography, electron microscopy, microscopy, specspectroscopy, wavefront sensing, holography, particle scattering, wavefront sensing, superresolution, synthesis, filter filter design, design, superresolution, radar radar signal signal and and antenna synthesis, and other disciplines disciplines that and other that share share an an important feature. feature. These These are are problems that involve the problems that involve the reconstruction or synthesis of aa reconstruction or synthesis of wavefront signal, etc.) when partial information wavefront (or (or an an object object or a signal, or constraints constraints exists exists in in each each of of two two different different domains. domains. The second domain is usually Fourier transform domain is transform domain. domain. This This paper usually the the Fourier paper describes combining all the available available information in describesaa method method of combining in the complete description, description, thereby thereby solvsolvthe two two domains domains to to arrive at a complete ing ing the problems. The problems fall into two general categories: categories: (1) (1) reconstruct reconstruct the entire (an image, image, wavefront, wavefront, signal, entire information information about about a function (an signal, etc.) is available etc.) when when only only partial partial information is available in in each each of two two domains; and (2) (2) synthesize mains; and (Fourier) transform transform pair pair having having synthesize aa (Fourier) desirable domains. A reconstruction desirable properties propertiesinin both both domains. reconstruction problem problem arises arises when when only only partial partial information information isis measured measured in in one one domain, and and in in the other domain either partial information information isis measured measured or or certain certain constraints constraints are are known known aa priori. priori. The information available in any one reconstruct the any one domain domain isis insufficient insufficient to to reconstruct the function function or its its A synthesis synthesis problem typically arises when one wants the transform. A to have have certain certain desirable desirable properties properties (such (such as transform of a function to as uniform spectrum, spectrum, low low sidelobes, sidelobes, etc.) while the function function itself itself uniform etc.) while satisfy certain must satisfy certain constraints or have certain desirable properties. Because arbitrary constraints can conBecause arbitrary sets sets of of properties properties and and constraints can be contradictory, there exist a transform pair that tradictory, there may may not exist that isis completely completely desirable and Nevertheless, one seeks a desirable and satisfies satisfies all all the the constraints. Nevertheless, transform pair comes as close close as as possible possible to transform pair that that comes to having having the the desirable properties desirable properties and and satisfying satisfying the the constraints constraints in in both domains. Both the the reconstruction reconstruction and the the synthesis synthesis problems problems can be be exexBoth pressed word "constraints" ''constraints" isis pressed asas follows, follows,ifif the the meaning meaning of of the the word broadened broadened to include include any any kind kind of measured measured data, desirable proper-
ties, or a priori conditions: conditions: ties, Given aa set constraints placed placed on Given set of constraints on a function and another set of constraints constraints placed placed on its transform, find find aa transform set transform pair (i.e., a function and and its its transform) transform) that that satisfies satisfies both sets of constraints. Once aa solution solution is problem, the the question question often often Once is found found to such a problem, remains: is the solution solution unique? unique? For Forsynthesis synthesis problems, problems, the the is the uniqueness one isis satisfied satisfied with with any any solusoluuniqueness isis usually usually unimportant unimportant -one tion satisfies all tion that satisfies all the the constraints; constraints; often often a more more important probproblem lem isis whether whether there there exists exists any any solution solution that satisfies satisfies what may be arbitrary and conflicting conflicting constraints. For reconstruction reconstruction problems, problems, the uniqueness uniqueness properties properties of the solution are of central importance. importance. If many different functions functions satisfying satisfying the the constraints constraintscould couldgive give rise rise to the same same measured measured data, data, then a solution that is found could not be guaranteed guaranteed to be be the the correct correct solution. solution. The The question question of ofuniqueuniqueness studied for for each each problem. problem. Fortunately, ness must must be be studied Fortunately, as as will will be be described some important reconstruction reconstruction problems described later, later, for for some problems the solution solution usually is unique. An effective effective approach solving the the large large class class of of problems problems An approach to solving described above iterative algorithms algorithms related described above isis the the use use of iterative related to to the Gerchberg-Saxton algorithms involve involve the the iterative GerchbergSaxtonalgorithm. algorithm.'1 The algorithms iterative back and and forth forth between between the the two two domains, domains, with with the transformation back the repetitively in each each domain. known constraints applied repetitively is presented in Sec. The basic algorithm is Sec. 2. 2. A number of different applications having types of of constraints constraints are are described, described, applications having different different types examples are shown in Sec. 3. In In Sec. Sec. 44 the the convergence convergence propand examples erties of the algorithm are discussed, discussed, and improved versions versions of erties of the the algorithm are reviewed. A summary and comments comments are are inalgorithm are reviewed. A brief summary Sec. 5. 5. cluded in Sec.
THE BASIC BASIC ITERATIVE ITERATIVE ALGORITHM ALGORITHM 2. THE account of ofthe theiterative iterativealgorithm algorithmwas was its its use The first published account use by by Gerchberg and Saxton' Saxton 1 to tosolve solve the Gerchberg theelectron electronmicroscopy microscopy problem. problem. this problem problem both both the themodulus modulus (magnitude) (magnitude) of For this ofaacomplexcomplexSPIE Signal Processing (1981) 147 SPIE Vol. Vol.373 373 Transformations Transformations in Optical Signal (1981 )// 147
Downloaded From: http://opticalengineering.spiedigitallibrary.org/ on 11/25/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx
are transform are Fourier transform its Fourier modulus of its valued image valued image and and the the modulus measured, measured, and the goal is is to to reconstruct the phase in both domains. was inApparently unknown to Gerchberg Gerchberg and Saxton, the method was inLesem2 to solve vented vented somewhat somewhat earlier earlierby by Hirsch, Hirsch, Jordan, Jordan, and Lesem2 solve aa synthesis problem problem for for computer -generatedholograms hologramsthat that has has a computer-generated synthesis more later in more described later will be described (This will constraints. (This similar similar set set of constraints. in problem in similar problem detail.) detail.) The method was again reinvented for aa similar 3 The computer computer holography holographyby by Gallagher Gallagherand and Liu. Liu.3 The fact fact that that the simplicity and efalgorithm algorithm was was invented invented repeatedly repeatedly testifies testifies to its simplicity fectiveness. Gerchberg-Saxton 2.1. 2.1. GerchbergSaxton algorithm is described in In what immediately immediately follows, follows, the iterative algorithm is reconstruction microscopy reconstruction electron microscopy the electron terms terms of its application to the microscopy phase electron microscopy the electron excellent treatment of the problem. problem. An excellent problem and its its solution solution by by this this and and other methods can be be found in principles to a large Ref. 4. 4. Later it is shown how to apply the same principles large Ref. class class of problems. image plane is wave function electron wave Suppose that that the the electron function in an image Suppose (2-D) two-dimensional described by the two -dimensional (2 -D) complex-valued complex -valued function function f(x) f(x) = f(x)
(1)
ei>G(x)
far-field in aa far function in wave function the wave Fourier transform, the Its Fourier Its -field diffraction diffraction plane, is given by
F(u) = = | ei_ where gk gk(x) 0, and gk _,+_ j(x) 1(x) 0. In (x) ), is is the the In spectroscopy, Fourier transform transform of the the complex complex degree degree of of temporal temporal coherence, coherence, Fourier 7(7), which 1X7)1 most easily easily measured. measured. In xx-ray Y(r), of which 11(r) is most -ray crystallography, the nonnegative electron electron density density function, crystallography, function, Q(x, @(x, y, y, z), which which isis periodic, z), periodic, isis the the Fourier Fourier transform transform of the structure facfacFhkl , of of which which Fhkl | is is measured measured by diffractometer. The tor Fhkh by a diffractometer. The will be astronomy problem will be described described in in more more detail detail later.
3.2.1. Uniqueness Uniqueness of of solutions solutions 3.2.1. one-dimensional For the one -dimensionalproblem, problem, use use of of the iterative algorithm (or any other other method) method) to reconstruct reconstruct the function function from from its any its Fourier Fourier since the solution solution in modulus is of limited interest since in the the general general case case isis usually not unique.26,27 unique. 26'27 The solution for for the theone one-The uniqueness uniqueness of the solution dimensional problem dimensional problem can can be be analyzed analyzed using using the the theory theory of analytic functions, additional solutions solutions can functions, from from which which one one finds finds that that additional can be be generated by of the the Fourier Fourier transform transform analytically generated by "flipping zeros" of analytically extended over plane. 26'27 The extended over the complex plane.26,27 Theadditional additional"solutions" "solutions" have the same support the original original function, function, but are not not have the same support as as the but are guaranteed to nonnegative; therefore therefore one one could could reduce reduce the guaranteed to be nonnegative; the degree of degree of ambiguity ambiguity by bygenerating generating all allpossible possible"solutions" "solutions" and and then then nonnegative ones.28 ones. 28 keeping only the nonnegative one-dimensional For certain special types of one -dimensional functions, there is is a high probability For aa function high probability that that the solution is unique. For function having having separated intervals separated by an interval two separated intervals of of support, support, being separated which the function function isis zero, zero, the thesolution solutionusually usuallyisisunique,29,3° over which unique, 29' 30 only if the two intervals of of support supportare aresufficiently sufficientlyseparated.31 but only separated. 31 special type function for for which which the the solution solutionisis usually Another special type of function usually unique is consisting of a summation summation of of aanumber numberofofdelta unique is one consisting delta-functions randomly space; for such such functions, functions, one one functions randomly distributed distributed in in space; does not need need the they can can be does the iterative iterative method method -they be reconstucted reconstucted by by aa noniterative method method involving involving the product of the product simple noniterative of three three
7
(b) (b) Fig. (a) and (b) (b) having modulus. Fig. 6. 6. Functions (a) having the same Fourier modulus.
the autocorrelation autocorrelation function.32 function. 32 translates of the In the event that multiple solutions do exist, exist, it would not appear the algorithm algorithm would would be biased toward one over another, another, and that the to converge converge to different different solutions, solutions, one would expect the algorithm to depending on the initial input to to the the algorithm. algorithm. For Forexample, example, Fig. Fig. 66 shows two same Fourier Fourier modulus. modulus. In In aa comcomshows two functions having the same puter experiment experiment using using the iterative iterative reconstruction reconstruction algorithm algorithm on functions' Fourier Fourier modulus, modulus, itit converged converged to one one of of the the solusoluthe functions' tions in about half half of of the the trials trials and andconverged converged to to the the other other solution solution of the the trials, trials, depending depending on on the the random random number number sesein the other half of quences used initial input to the the algorithm. algorithm. quences used as as the initial For the problem in two or more more dimensions, dimensions, it appears appears that the the is usually usually unique. Considering Considering sampled sampled functions functions defined defined solution is on a rectangular grid of points, points, Bruck Bruck and and Sodin33 Sodin33 showed showed that that the the existence of additional solutions solutions is is equivalent equivalent to to the the factorability factorability of of existence polynomial representation representation of the Fourier Fourier transform. transform. Since a polynomial of the Since aa of one one variable variable of of degree degree M M can always be polynomial of be factored into prime factors, factors,there thereare are2M 2M-1 solutions ininthe theone one-dimensional M prime -1 solutions -dimensional case. Once Once again, again, only only some some of the the "solutions" "solutions" may may be be non non-case. negative. On the other other hand, hand,polynomials polynomialsof oftwo twoor ormore morevariables variables negative. coefficients are are only only rarely factorable; consequentconsequenthaving arbitrary coefficients the two two-dimensional Attempts have ly, the -dimensional problem is usually unique. Attempts have also been made to extend this concept to continuous, as opposed to functions. 34 Although itit is is always always possible to make up exdiscrete, functions.34 amples in dimensions that are are not not unique,35 unique, 35 it appears appears to be amples in two dimensions be
152 SPIE Vol Vol. 373 Transformations Transformations in Optical Optical Signal Signal Processing (1981) (1981) 152 //SPIE
Downloaded From: http://opticalengineering.spiedigitallibrary.org/ on 11/25/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx
real world, the real two-dimensional for two true that for -dimensional functions drawn from the two-the two of the uniqueness of general uniqueness The general the solution is usually unique. The results dimensional case case isis indicated indicated by by experimental experimental reconstruction results Fourier the Fourier noise in the algorithm. 36 Furthermore, iterative algorithm.36 using using the the iterative Furthermore, noise noise to the adding noise the effect of adding had the modulus data modulus data has has had conalgorithm to concausing the reconstructed function rather rather than causing the algorithm reconstructed function solution. 37 different solution.37 verge to a radically different verge
Astronomical reconstruction 3.2,2, Astronomical 3.2.2. nonnegative two-dimensional reconstructing a atwo ofreconstructing problem of The problem The -dimensional nonnegative
(a)
(b)
(c)
(d)
(e)
(f)
(9) (g)
(h)
in arises in function from from the the modulus of its transform arises Fourier transform its Fourier modulus of function atresolution atthe resolution turbulence, the atmospheric turbulence, astronomy. Due astronomy. Due to atmospheric tainable from from large large optical optical telescopes telescopeson on earth earth isis only only about about one one tainable
second of arc, many times worse worse than the diffraction limit limit imposed a five-meter Fora five by -meter aperture. For telescope aperture. thetelescope of the diameter of the diameter by the diffraction-limited the diffraction telescope aperture, aperture, the -limited resolution resolution would would be be telescope atmospheric Despite atmospheric finer. Despite times finer. fifty times about about 0.02 0.02 seconds seconds of of arc arc -fifty turbulence, itit isis possible possible to to measure measure the the modulus modulus of of the Fourier turbulence, the limit of the diffraction limit the diffraction space object transform transform of of aa space object out out to the techniques. 38"41 telescope telescope using using interferometric interferometric techniques.38 -41The Theautocorrelation autocorrelation allowing the object object can can be be computed computed from the Fourier modulus, allowing of the the diameter diameter of of the the object object to be determined. the unless the However, unless determined. However, the previously not was previously measured, it was Fourier transform transform phase phase isis also also measured, not Fourier possible to determine special some special for some except for itself, except object itself, the object determine the possible to cases. Previous Previous attempts attempts to to solve solve this this problem problem had had not not proven proven to to be be cases. objects. two-dimensional practical for complicated -dimensional objects. complicated two interferometer The an object from interferometer object from reconstructing an problem of reconstructing The problem Fourier-The Fourier method. 42' 36 The iterativemethod.42,36 theiterative bythe solved by be solved can be data data can Fourier modulus equal the Fourier Fourier modulus domain constraint isis that that the Fourier domain constraint function-domain the function modulus measured measured by by an an interferometer, and the -domain modulus Figure 77 nonnegative. Figure be nonnegative. function be object function constraint that the the object constraint isis that computer-synthesized shows aa computer7(a) shows shows an an example. example. Fig. 7(a) synthesized object shows sun-like used for the experiment experiment -aa sun -like disk disk having having "solar "solar flares" and Fourier transform its Fourier of its modulus of Themodulus bright and and dark dark "sunspots." "sunspots." The random of random square of shows aa square 7(c) shows Figure 7(c) 7(b). Figure Fig, 7(b). shown in is is shown in Fig. numbers used used as as the the initial initial input input for the iterative algorithm. Figures numbers 7(d), 230, 20, 230, results after 20, reconstruction results 7(f) show the reconstruction 7(e), and 7(f) 7(d), 7(e), input initial input 7(g) shows the initial Figure 7(g) respectively. Figure 600 iterations, respectively. and 600 215 and 215 after 22 and results after reconstruction results the reconstruction second trial, and the for for a second iterations are shown in Comparing respectively. Comparing 7(i), respectively. and 7(i), 7(h) and Figs, 7(h) in Figs. sees that one sees Figs. 7(i) with with the the original object in Fig. 7(a), one 7(f) and 7(i) Figs, 7(f) object original object the original images match the reconstructed images the reconstructed for both trials, the very closely. very closely. Note Note that that inverted inverted solutions such as Fig. 7(f) are permitted for this problem since since the the modulus modulus of of the the Fourier Fourier transform transform f(-x) of f( -x) equals equalsthe the modulus modulusof of the the Fourier Fourier transform transform of of f(x) f(x) for for real-valued real -valuedf(x). f(x).Other Other successful successful reconstruction reconstruction experiments experiments have have noise presof noise types of the types have the simulated to have performed on data simulated been been performed 39 and speckle interferometry, stellar speckle in stellar ent in interferometry,39 anditit appears appears that that under diffraction-realistic levels levelsof ofphoton photon noise noise for for fairly fairly bright objects, diffraction realistic 37 Initial limited limited images images can can be be reconstructed. reconstructed.37 Initial expements expements have have also telescopes. 43 fromtelescopes.43 datafrom on data been carried out on
3.2.3. synthesis and synthesis reconstruction and Pupil reconstruction 3.2. 3. Pupil two-reconstruct aatwo to reconstruct want to may want one may which one Another case in in which Another case modulus is dimensional nonnegative dimensional nonnegativefunction function from from its its Fourier Fourier modulus is in in optical diffraction-limited In a adiffraction determination. In function determination. pupil function pupil -limited optical
system, the pointspread function function isis the the squared Fourier modulus modulus point-spread system, the transfer optical transfer Equivalently, the optical pupil function. Equivalently, system's pupil of the system's function. 44 Given the pupil function.'" function function isis the the autocorrelation autocorrelation of of the the pupil one image plane, one pointspreadfunction functionatat aa given given location location in in an image point-spread could use the iterative algorithm to retrieve the corresponding pupil the equivalent to mathematically equivalent function, inin a way to the that isismathematically way that function, astronomy problem. astronomy problem. Turning Turning this this problem problem around, around, one could use synthesize (design) the iterative iterative algorithm algorithm to to synthesize (design)aa pupil pupil function function that the possibly while possibly function while point-spread desired pointgiven, desired would spread function yield a given, would yield well. satisfying satisfying other other desirable desirable constraints as well. aperture over part of an aperture extent—measurement Finite extent 3.3. 3.3. Finite -measurement over
In a number problems, there there isis aa function function of of reconstruction problems, number of reconstruction In
(i)
Fourier Fig. 7.7. Reconstruction of a function from from its Fourier nonnegative function a nonnegative Reconstruction of Fig. Test object; (b) (b) modulus modulus of its Fourier (c) initransform; (c) Fourier transform; modulus, modulus. (a) Test -(f) reconstruction results results (d)-(f) test);(d) (firsttest); object(first the object of the estimate of tial estimate - number (g) initial initial estimate of of 600; (g) (f) 600; 230, (f) (e) 230, 20,(e) (d)20, iterations:(d) numberofoffiterations: the object (second -(i)reconstruction reconstructionresults results-number number of of (h)-(i) (second test); (h) the iterations: (h) 215. (i) 215. 2, (i) (h) 2, iterations:
reconstruct the wishes to reconstruct one wishes and one known known finite finite extent extent (or support) and Fourier in. the Fourier function with with resolution resolution appropriate appropriate to an aperture in function domain more were measurements were which measurements overwhich oneover theone than the complete than more complete actually taken. In larger simply larger aperture isis simply desired aperture the desired cases, the some cases, In some measurements were taken, and so one which measurements over which aperture over the aperture than the wishes to extrapolate extrapolate the the function's function's Fourier Fourier transform, transform, i.e., to obwishes to tain superresolution has made one has cases, one other cases, In other superresolution of the function. In one case one which case in which aperture, in filled aperture, measurements over measurements over aa partially filled and function, and wishes to interpolate interpolate the the Fourier Fourier transform transform of the function, wishes to dofunction dothe function response in the impulse response thereby improved impulse an improved thereby obtain an main.
superresolution or superresolution Extrapolation or 3.3.1. 3.3.1. Extrapolation The error-reduction error -reduction algorithm algorithm was was first first applied applied to to the the extrapolation extrapolation (or superresolution) problem writbeen writhas been Much has Gerchberg. 45 Much by Gerchberg.45 problem by error-reduction theerror specifically the algorithm,specifically ten about -reduction iterative algorithm, the iterative about the ways of algorithm, various ways including various problem, including this problem, relates to this algorithm, as it relates proofs of 2) and Sec. 2) understanding (see the the end end of Sec. and proofs algorithm (see understanding the algorithm convergence.1o4243 °46-48 this particularproblem, problem, the the nature nature of particular this ForFor convergence. 10' 12 ' 13 '46"48 algorithm by implement the algorithm possible to implement the constraints makes it possible by aa taking on the order seconds 10"9 seconds of10-9 order of processor49'50 taking feedback optical processor49,5o Smith case. Marks and Smith two-dimensional thetwo for the per iteration even -dimensional case. even for describe these these matters in detail volume. elsewhere in this volume. detail elsewhere Interpolation 3.3.2. Interpolation 3.3.2. In tomographic systems, many many projections projections of of the the object object imaging systems, tomographic imaging In slice yielding information about aa slice projection yielding each projection measured, each are measured, through the Fourier transform measurements When measurements object. When the object. of the transform of over only a limited cone of angles effective aperture the effective angles are made, the SP /EVol. Vol373 373 Transformations Transformations In in Optical S/gna/ Signal Processing 153 1981) )// 153 Processing f(1981 SPIE
Downloaded From: http://opticalengineering.spiedigitallibrary.org/ on 11/25/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx
in Fourier domain domain has gaps, and and the the impulse impulse response response of in the Fourier of the the system system isis highly highlyirregular. irregular. In In applying applying the iterative algorithm to this this problem, 51 ' 52 the problem,51,52 thefunction-domain function -domainconstraint constraintisis the the finite finite extent extent and and nonnegativity nonnegativity of of the the object, object, and and the the Fourier Fourier domain domain constraint constraint isis that that the the Fourier Fourier transform transform equal equal the the measured measured Fourier Fourier transform transform over over the the measurement measurement aperture. problem similar similar to tomography problem A problem to the tomography in radio problem arises arises in astronomy. sky brightness brightness map map isisaatwo two-dimensional astronomy. The The radio sky -dimensional real, real, nonnegative nonnegativefunction function which whichisisthe the Fourier Fourier transform transform of of the the complex complex visibility visibilityfunction. function. The The visibility visibility function function isis measured measured by by radio interferornetry, interferometry, and in the case case of -baseline inoflong long-baseline ter ferometry, the visibility visibility function is measured only over a limited terferometry, "tracks"ininthe theFourier Fourierdomain, domain,resulting resultinginina apartially partially-filled set of "tracks" -filled effective aperture. err or-reduction has been been used used to to effective aperture. The error -reduction algorithm has obtain in effect, effect, interpolating interpolating the obtain improved improved maps maps by, by, in the visibility visibility function fill in the area between between the tracks.53 function to fill tracks. 53 For this problem, that itit be be nonnegative nonnegative and the constraints on the brightness map are that field of view. view. In the visibility visibility plane, be zero zero outside outside the the known known field plane, the the constraint is that the the complex complex visiblity visiblity function equal equal the the measured measured value within value within the the area area of of the the tracks. 3.4. Modulus quantized values 3.4. Modulusquantized values As spectrum, shaping, in in comcomAs mentioned mentioned earlier earlier in in connection with spectrum puter holography one may wish wish to to encode encode the the Fourier transform of an image as computer-generated an image as a computer -generatedhologram, hologram,but but some some types types of of computer-generated computer -generatedholograms hologramscan can encode encode only only certain, certain quantized quantized complex kinoform example example discussed discussed earlier earlierisis aa special complex values. values. The kinoform special type quantization. A A more more general general example example is is the the Lohmann Lohrnann type of quantization. hologram, 54 54for forwhich whichthe themodulus modulusand and phase phase of a complex complex sample sample are determined respectively, of of an determined by the area and relative position, respectively, aperture within a sampling sampling cell. cell. The The number number of of allowable allowable quantized quantized values resolution elements, elements, of the values isis determined determined by by the the number number of resolution recording device recording deviceused usedto to fabricate fabricate the the hologram, hologram, used used to form form, one one cell. synthesis problem, function-domain cell. For this synthesis problem, the functiondomain constraint is is that the modulus of the function function equal equal the the desired desired image image modulus and the Fourier Fourier-domain and the -domain constraint constraint isis that that the the complex complex Fourier Fourier coefficients prescribed set set of coefficientsfall fall on on aa prescribed of quantized quantized values. values. ExExperiments synthesizing such transform periments have have shown shown that synthesizing such aa Fourier transform pair is possible using using the iterative algorithm.55,7 algorithm,. 55 ' 7 For example, Fig. 8(a) image produced produced by Lohrnann 8(a) shows shows aa simulation simulation of of an an image by a Lohmann hologram having modulus and four four phase hologram having only only four four modulus phase quantization quantization levels levels when, when the the image image was wasrandom, random phase coded. Figure Figure 8(b) 8(b) shows the image after 13 1,3 iterations, iterations, aa considerable considerable improvement. improvement. This This the image problem general class class of problems problem isis one one of a more general the problems regarding regarding the transmission transmission of of coded coded data. 3.5. Finite Finite extent extent -phase phase 3.5. Finally, the iterative algorithm Finally, the iterative algorithm has has been been used used to reconstruct reconstruct the modulus of a band band-limited from its its phase.56.57 phase. 56 ' 57 Or, -limited signal from Or, looking at itit in way, given given that a function has has finite finite extent extent and and given in another way, given the reconstruct the the phase phase of of its its Fourier Fourier transform, transform, reconstruct the modulus its modulus of its Fourier shown that Fourier transform. transform. For this application, it has been shown that for a wide is unique.56 unique. 56 This application wide class class of of conditions conditions the solution is will Sec. 4. will be be discussed discussed further further in Sec.
4. ALGORITHM CONVERGENCE AND 4. ALGORITHM ACCELERATED ACCELERATED ALGORITHMS ALGORITHMS As mentioned mentioned in Sec. 2, the the basic basic iterative iterative algorithm algorithm depicted As in Sec. depicted in in Fig. 1, referred has been been shown Fig. 1, referred to as the error-reduction error -reduction algorithm, has shown to converge section, the the convergence convergence converge for for some some applications. applications. In this section, is proven is proven for for all all applications. applications. In In addition, modified algorithms that often error-reduction often converge converge much much faster faster than than the error -reduction algorithm algorithm are are discussed. Convergence of the the errorerror-reduction algorithm 4.1. Convergence reduction algorithm For the error error-reduction mean-squared -reduction algorithm, algorithm, the mean -squared error error can be defined in It is is aa normalized normalized version version defined in general general by by Eq. Eq. (7) (7) or or Eq. (8). It the integral integral over the square of of the the amount amount by bywhich which the of the thecorncorn-
Fig. Computer-simulated Fig. 8. 8. Computer -simulated images images from from hologram hologram with with four magnitude magnitude and and four phase phase quantized quantized levels. levels, (a) (a) Object Object random random phased phased coded; coded; (b) (b) after after 13 13 iterations iterations of the the iterative iterative method. method.
puted function function (or (or the the computed computed Fourier Fourier transform) transform) violates violates the the puted appropriatedomain. domain.When Whenthe constraints in the appropriate -squared erthemean mean-squared eris zero, zero, then then aa Fourier Fourier transform transform pair pair has has been ror is been found found that that satisfies all satisfies all the the constraints constraints in in both domains. Consider steps in in the error-reduction the error -reduction algorithm algorithm Consider again again the the steps described k th iteration starts starts with with an an estimate estimate gk(x) described in in Sec. Sec. 2. 2. The kth g^(x) that satisfies satisfies the function-domain the function -domain constraints. constraints. For any coordinate, the complex complex values values that g(x) g(x) can have have that that satisfy satisfy the the function function-x, the domain space. For domain constraints constraints form form some some set set of of points points in in phasor space. For example, modulus must equal equal f(x) f(x)1, then the set |, then example, ifif the the modulus set of of such such points circle of radius | f(x) f(x)1| in phasor space; space; if if the the function function points isis a circle must must be be nonnegative, nonnegative, then then the set of such points is line on on is the half line the nonnegative real axis. The The function function estimate estimate gk(x) gk (x) is is Fourier Fourier the nonnegative real axis. transformed, yielding G k (u). The transformed, yielding Gk(u). The next next step step in the algorithm algorithm is is to form Gk(u) G^(u) by changing changing G^(u) Gk(u) by by the the smallest smallest possible possibleamount amount that that Fourier-domain constraints. Gk(u) G£(u)isis then then ininallows it to satisfy the Fourierdomain constraints. verse Fourier transformed, yielding yielding gk(x) g£(x) in the function verse function domain. domain. In the final final step, step, gk g^ +1(x) + jM is is formed In the by the formed by by changing changing gk(x) g£(x) by the smallest amount that that allows allows itit to tosatisfy satisfythe smallest -domain conthefunction function-domain constraints. Now Now consider consider the unnormalized unnorrnalized squared squared error, error, given given by by straints. (7) and (8). the numerators in Eqs. (7) (8). In the Fourier domain, domain, the ununnormalized squared k th iteration iteration is is normalized squared error at the kth co 2 _ eFk -
G k (u)IGk(u) - Gk(u) 1 2 du du _00
154 SPIE Vol. Vol. 373 373 Transformations Transformations in in Optical Optical Signal Signal Processing Processing (1981) 154 //SPIE
Downloaded From: http://opticalengineering.spiedigitallibrary.org/ on 11/25/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx
(14) (14)
•/
oc
v^/
dx,, |gk(x)-g£(x)| = I I gk(x) - gk(x) 122dx -00
Parseval's from Parseval's results from equation results this equation second line where the the second line in in this where theorem. The The unnormalized unnormalized squared squared error error in in the function domain theorem. given by is given at the kth iteration is w 2
= eOk = e0k
fI gk
dx . + i(x)-g£(x)| + l(x) - gk(x) I 22dx.
INPUT 9 g
(15)
-03
function-domain the function satisfy the -domain i(x) by definition satisfy + i(x) gk + Both gk(x) gk(x) and gk j(x) is gk++ 1(x) coordinatex,x,gk given coordinate any given at any constraints. Also at is the point in function-domain satisfying the function phasor -domain constraints constraints that isis space satisfying phasor space closest to gk(x). gk(x). Therefore, Therefore, for all values of x, closest to +1W II gk gk + 1(x) --
g£ gk(x) -- gk(x) gk WI I
Abstract. This paper paper reviews reviews the the GerchbergGerchberg-Saxton Saxton algorithm and and variavaria-
Environmental Environmental Research Research Institute Institute of Michigan P.O. P.O. Box 8618 Ann Arbor, Arbor, Michigan 48107 48107 Ann
tions thereof that have have been been used used to solve a number number of difficult difficult reconreconstruction and and synthesis synthesis problems problems in in optics and and related related fields. It can can be be used information (including (including both both used on on any any problem problem in in which which only partial information measurements the wavefront wavefront or or signal signal isis available available in in measurements and and constraints) of the one information isis available available inin another another domain one domain domain and and other other partial partial information domain (usually (usually the the Fourier Fourier domain). domain). The Thealgorithm algorithm combines combines the the information in both arrive at of the the wavefront wavefront or both domains domains to arrive at the complete description of or signal. are reviewed, reviewed, including synthesis of Fourier signal. Various Various applications are Fourier transform pairs having having desirable desirable properties properties as as reconstruction transform pairs as well as reconstruction problems. and the the convergence convergence properties problems. Variations of the algorithm and properties of of the are discussed. the algorithm are
1. INTRODUCTION INTRODUCTION 1. There very difficult There exist exist many many problems problemsthat that are are very to solve difficult to solve in astronomy, x x-ray crystallography, electron -ray crystallography, electron microscopy, microscopy, specspectroscopy, wavefront sensing, holography, particle scattering, wavefront sensing, superresolution, synthesis, filter filter design, design, superresolution, radar radar signal signal and and antenna synthesis, and other disciplines disciplines that and other that share share an an important feature. feature. These These are are problems that involve the problems that involve the reconstruction or synthesis of aa reconstruction or synthesis of wavefront signal, etc.) when partial information wavefront (or (or an an object object or a signal, or constraints constraints exists exists in in each each of of two two different different domains. domains. The second domain is usually Fourier transform domain is transform domain. domain. This This paper usually the the Fourier paper describes combining all the available available information in describesaa method method of combining in the complete description, description, thereby thereby solvsolvthe two two domains domains to to arrive at a complete ing ing the problems. The problems fall into two general categories: categories: (1) (1) reconstruct reconstruct the entire (an image, image, wavefront, wavefront, signal, entire information information about about a function (an signal, etc.) is available etc.) when when only only partial partial information is available in in each each of two two domains; and (2) (2) synthesize mains; and (Fourier) transform transform pair pair having having synthesize aa (Fourier) desirable domains. A reconstruction desirable properties propertiesinin both both domains. reconstruction problem problem arises arises when when only only partial partial information information isis measured measured in in one one domain, and and in in the other domain either partial information information isis measured measured or or certain certain constraints constraints are are known known aa priori. priori. The information available in any one reconstruct the any one domain domain isis insufficient insufficient to to reconstruct the function function or its its A synthesis synthesis problem typically arises when one wants the transform. A to have have certain certain desirable desirable properties properties (such (such as transform of a function to as uniform spectrum, spectrum, low low sidelobes, sidelobes, etc.) while the function function itself itself uniform etc.) while satisfy certain must satisfy certain constraints or have certain desirable properties. Because arbitrary constraints can conBecause arbitrary sets sets of of properties properties and and constraints can be contradictory, there exist a transform pair that tradictory, there may may not exist that isis completely completely desirable and Nevertheless, one seeks a desirable and satisfies satisfies all all the the constraints. Nevertheless, transform pair comes as close close as as possible possible to transform pair that that comes to having having the the desirable properties desirable properties and and satisfying satisfying the the constraints constraints in in both domains. Both the the reconstruction reconstruction and the the synthesis synthesis problems problems can be be exexBoth pressed word "constraints" ''constraints" isis pressed asas follows, follows,ifif the the meaning meaning of of the the word broadened broadened to include include any any kind kind of measured measured data, desirable proper-
ties, or a priori conditions: conditions: ties, Given aa set constraints placed placed on Given set of constraints on a function and another set of constraints constraints placed placed on its transform, find find aa transform set transform pair (i.e., a function and and its its transform) transform) that that satisfies satisfies both sets of constraints. Once aa solution solution is problem, the the question question often often Once is found found to such a problem, remains: is the solution solution unique? unique? For Forsynthesis synthesis problems, problems, the the is the uniqueness one isis satisfied satisfied with with any any solusoluuniqueness isis usually usually unimportant unimportant -one tion satisfies all tion that satisfies all the the constraints; constraints; often often a more more important probproblem lem isis whether whether there there exists exists any any solution solution that satisfies satisfies what may be arbitrary and conflicting conflicting constraints. For reconstruction reconstruction problems, problems, the uniqueness uniqueness properties properties of the solution are of central importance. importance. If many different functions functions satisfying satisfying the the constraints constraintscould couldgive give rise rise to the same same measured measured data, data, then a solution that is found could not be guaranteed guaranteed to be be the the correct correct solution. solution. The The question question of ofuniqueuniqueness studied for for each each problem. problem. Fortunately, ness must must be be studied Fortunately, as as will will be be described some important reconstruction reconstruction problems described later, later, for for some problems the solution solution usually is unique. An effective effective approach solving the the large large class class of of problems problems An approach to solving described above iterative algorithms algorithms related described above isis the the use use of iterative related to to the Gerchberg-Saxton algorithms involve involve the the iterative GerchbergSaxtonalgorithm. algorithm.'1 The algorithms iterative back and and forth forth between between the the two two domains, domains, with with the transformation back the repetitively in each each domain. known constraints applied repetitively is presented in Sec. The basic algorithm is Sec. 2. 2. A number of different applications having types of of constraints constraints are are described, described, applications having different different types examples are shown in Sec. 3. In In Sec. Sec. 44 the the convergence convergence propand examples erties of the algorithm are discussed, discussed, and improved versions versions of erties of the the algorithm are reviewed. A summary and comments comments are are inalgorithm are reviewed. A brief summary Sec. 5. 5. cluded in Sec.
THE BASIC BASIC ITERATIVE ITERATIVE ALGORITHM ALGORITHM 2. THE account of ofthe theiterative iterativealgorithm algorithmwas was its its use The first published account use by by Gerchberg and Saxton' Saxton 1 to tosolve solve the Gerchberg theelectron electronmicroscopy microscopy problem. problem. this problem problem both both the themodulus modulus (magnitude) (magnitude) of For this ofaacomplexcomplexSPIE Signal Processing (1981) 147 SPIE Vol. Vol.373 373 Transformations Transformations in Optical Signal (1981 )// 147
Downloaded From: http://opticalengineering.spiedigitallibrary.org/ on 11/25/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx
are transform are Fourier transform its Fourier modulus of its valued image valued image and and the the modulus measured, measured, and the goal is is to to reconstruct the phase in both domains. was inApparently unknown to Gerchberg Gerchberg and Saxton, the method was inLesem2 to solve vented vented somewhat somewhat earlier earlierby by Hirsch, Hirsch, Jordan, Jordan, and Lesem2 solve aa synthesis problem problem for for computer -generatedholograms hologramsthat that has has a computer-generated synthesis more later in more described later will be described (This will constraints. (This similar similar set set of constraints. in problem in similar problem detail.) detail.) The method was again reinvented for aa similar 3 The computer computer holography holographyby by Gallagher Gallagherand and Liu. Liu.3 The fact fact that that the simplicity and efalgorithm algorithm was was invented invented repeatedly repeatedly testifies testifies to its simplicity fectiveness. Gerchberg-Saxton 2.1. 2.1. GerchbergSaxton algorithm is described in In what immediately immediately follows, follows, the iterative algorithm is reconstruction microscopy reconstruction electron microscopy the electron terms terms of its application to the microscopy phase electron microscopy the electron excellent treatment of the problem. problem. An excellent problem and its its solution solution by by this this and and other methods can be be found in principles to a large Ref. 4. 4. Later it is shown how to apply the same principles large Ref. class class of problems. image plane is wave function electron wave Suppose that that the the electron function in an image Suppose (2-D) two-dimensional described by the two -dimensional (2 -D) complex-valued complex -valued function function f(x) f(x) = f(x)
(1)
ei>G(x)
far-field in aa far function in wave function the wave Fourier transform, the Its Fourier Its -field diffraction diffraction plane, is given by
F(u) = = | ei_ where gk gk(x) 0, and gk _,+_ j(x) 1(x) 0. In (x) ), is is the the In spectroscopy, Fourier transform transform of the the complex complex degree degree of of temporal temporal coherence, coherence, Fourier 7(7), which 1X7)1 most easily easily measured. measured. In xx-ray Y(r), of which 11(r) is most -ray crystallography, the nonnegative electron electron density density function, crystallography, function, Q(x, @(x, y, y, z), which which isis periodic, z), periodic, isis the the Fourier Fourier transform transform of the structure facfacFhkl , of of which which Fhkl | is is measured measured by diffractometer. The tor Fhkh by a diffractometer. The will be astronomy problem will be described described in in more more detail detail later.
3.2.1. Uniqueness Uniqueness of of solutions solutions 3.2.1. one-dimensional For the one -dimensionalproblem, problem, use use of of the iterative algorithm (or any other other method) method) to reconstruct reconstruct the function function from from its any its Fourier Fourier since the solution solution in modulus is of limited interest since in the the general general case case isis usually not unique.26,27 unique. 26'27 The solution for for the theone one-The uniqueness uniqueness of the solution dimensional problem dimensional problem can can be be analyzed analyzed using using the the theory theory of analytic functions, additional solutions solutions can functions, from from which which one one finds finds that that additional can be be generated by of the the Fourier Fourier transform transform analytically generated by "flipping zeros" of analytically extended over plane. 26'27 The extended over the complex plane.26,27 Theadditional additional"solutions" "solutions" have the same support the original original function, function, but are not not have the same support as as the but are guaranteed to nonnegative; therefore therefore one one could could reduce reduce the guaranteed to be nonnegative; the degree of degree of ambiguity ambiguity by bygenerating generating all allpossible possible"solutions" "solutions" and and then then nonnegative ones.28 ones. 28 keeping only the nonnegative one-dimensional For certain special types of one -dimensional functions, there is is a high probability For aa function high probability that that the solution is unique. For function having having separated intervals separated by an interval two separated intervals of of support, support, being separated which the function function isis zero, zero, the thesolution solutionusually usuallyisisunique,29,3° over which unique, 29' 30 only if the two intervals of of support supportare aresufficiently sufficientlyseparated.31 but only separated. 31 special type function for for which which the the solution solutionisis usually Another special type of function usually unique is consisting of a summation summation of of aanumber numberofofdelta unique is one consisting delta-functions randomly space; for such such functions, functions, one one functions randomly distributed distributed in in space; does not need need the they can can be does the iterative iterative method method -they be reconstucted reconstucted by by aa noniterative method method involving involving the product of the product simple noniterative of three three
7
(b) (b) Fig. (a) and (b) (b) having modulus. Fig. 6. 6. Functions (a) having the same Fourier modulus.
the autocorrelation autocorrelation function.32 function. 32 translates of the In the event that multiple solutions do exist, exist, it would not appear the algorithm algorithm would would be biased toward one over another, another, and that the to converge converge to different different solutions, solutions, one would expect the algorithm to depending on the initial input to to the the algorithm. algorithm. For Forexample, example, Fig. Fig. 66 shows two same Fourier Fourier modulus. modulus. In In aa comcomshows two functions having the same puter experiment experiment using using the iterative iterative reconstruction reconstruction algorithm algorithm on functions' Fourier Fourier modulus, modulus, itit converged converged to one one of of the the solusoluthe functions' tions in about half half of of the the trials trials and andconverged converged to to the the other other solution solution of the the trials, trials, depending depending on on the the random random number number sesein the other half of quences used initial input to the the algorithm. algorithm. quences used as as the initial For the problem in two or more more dimensions, dimensions, it appears appears that the the is usually usually unique. Considering Considering sampled sampled functions functions defined defined solution is on a rectangular grid of points, points, Bruck Bruck and and Sodin33 Sodin33 showed showed that that the the existence of additional solutions solutions is is equivalent equivalent to to the the factorability factorability of of existence polynomial representation representation of the Fourier Fourier transform. transform. Since a polynomial of the Since aa of one one variable variable of of degree degree M M can always be polynomial of be factored into prime factors, factors,there thereare are2M 2M-1 solutions ininthe theone one-dimensional M prime -1 solutions -dimensional case. Once Once again, again, only only some some of the the "solutions" "solutions" may may be be non non-case. negative. On the other other hand, hand,polynomials polynomialsof oftwo twoor ormore morevariables variables negative. coefficients are are only only rarely factorable; consequentconsequenthaving arbitrary coefficients the two two-dimensional Attempts have ly, the -dimensional problem is usually unique. Attempts have also been made to extend this concept to continuous, as opposed to functions. 34 Although itit is is always always possible to make up exdiscrete, functions.34 amples in dimensions that are are not not unique,35 unique, 35 it appears appears to be amples in two dimensions be
152 SPIE Vol Vol. 373 Transformations Transformations in Optical Optical Signal Signal Processing (1981) (1981) 152 //SPIE
Downloaded From: http://opticalengineering.spiedigitallibrary.org/ on 11/25/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx
real world, the real two-dimensional for two true that for -dimensional functions drawn from the two-the two of the uniqueness of general uniqueness The general the solution is usually unique. The results dimensional case case isis indicated indicated by by experimental experimental reconstruction results Fourier the Fourier noise in the algorithm. 36 Furthermore, iterative algorithm.36 using using the the iterative Furthermore, noise noise to the adding noise the effect of adding had the modulus data modulus data has has had conalgorithm to concausing the reconstructed function rather rather than causing the algorithm reconstructed function solution. 37 different solution.37 verge to a radically different verge
Astronomical reconstruction 3.2,2, Astronomical 3.2.2. nonnegative two-dimensional reconstructing a atwo ofreconstructing problem of The problem The -dimensional nonnegative
(a)
(b)
(c)
(d)
(e)
(f)
(9) (g)
(h)
in arises in function from from the the modulus of its transform arises Fourier transform its Fourier modulus of function atresolution atthe resolution turbulence, the atmospheric turbulence, astronomy. Due astronomy. Due to atmospheric tainable from from large large optical optical telescopes telescopeson on earth earth isis only only about about one one tainable
second of arc, many times worse worse than the diffraction limit limit imposed a five-meter Fora five by -meter aperture. For telescope aperture. thetelescope of the diameter of the diameter by the diffraction-limited the diffraction telescope aperture, aperture, the -limited resolution resolution would would be be telescope atmospheric Despite atmospheric finer. Despite times finer. fifty times about about 0.02 0.02 seconds seconds of of arc arc -fifty turbulence, itit isis possible possible to to measure measure the the modulus modulus of of the Fourier turbulence, the limit of the diffraction limit the diffraction space object transform transform of of aa space object out out to the techniques. 38"41 telescope telescope using using interferometric interferometric techniques.38 -41The Theautocorrelation autocorrelation allowing the object object can can be be computed computed from the Fourier modulus, allowing of the the diameter diameter of of the the object object to be determined. the unless the However, unless determined. However, the previously not was previously measured, it was Fourier transform transform phase phase isis also also measured, not Fourier possible to determine special some special for some except for itself, except object itself, the object determine the possible to cases. Previous Previous attempts attempts to to solve solve this this problem problem had had not not proven proven to to be be cases. objects. two-dimensional practical for complicated -dimensional objects. complicated two interferometer The an object from interferometer object from reconstructing an problem of reconstructing The problem Fourier-The Fourier method. 42' 36 The iterativemethod.42,36 theiterative bythe solved by be solved can be data data can Fourier modulus equal the Fourier Fourier modulus domain constraint isis that that the Fourier domain constraint function-domain the function modulus measured measured by by an an interferometer, and the -domain modulus Figure 77 nonnegative. Figure be nonnegative. function be object function constraint that the the object constraint isis that computer-synthesized shows aa computer7(a) shows shows an an example. example. Fig. 7(a) synthesized object shows sun-like used for the experiment experiment -aa sun -like disk disk having having "solar "solar flares" and Fourier transform its Fourier of its modulus of Themodulus bright and and dark dark "sunspots." "sunspots." The random of random square of shows aa square 7(c) shows Figure 7(c) 7(b). Figure Fig, 7(b). shown in is is shown in Fig. numbers used used as as the the initial initial input input for the iterative algorithm. Figures numbers 7(d), 230, 20, 230, results after 20, reconstruction results 7(f) show the reconstruction 7(e), and 7(f) 7(d), 7(e), input initial input 7(g) shows the initial Figure 7(g) respectively. Figure 600 iterations, respectively. and 600 215 and 215 after 22 and results after reconstruction results the reconstruction second trial, and the for for a second iterations are shown in Comparing respectively. Comparing 7(i), respectively. and 7(i), 7(h) and Figs, 7(h) in Figs. sees that one sees Figs. 7(i) with with the the original object in Fig. 7(a), one 7(f) and 7(i) Figs, 7(f) object original object the original images match the reconstructed images the reconstructed for both trials, the very closely. very closely. Note Note that that inverted inverted solutions such as Fig. 7(f) are permitted for this problem since since the the modulus modulus of of the the Fourier Fourier transform transform f(-x) of f( -x) equals equalsthe the modulus modulusof of the the Fourier Fourier transform transform of of f(x) f(x) for for real-valued real -valuedf(x). f(x).Other Other successful successful reconstruction reconstruction experiments experiments have have noise presof noise types of the types have the simulated to have performed on data simulated been been performed 39 and speckle interferometry, stellar speckle in stellar ent in interferometry,39 anditit appears appears that that under diffraction-realistic levels levelsof ofphoton photon noise noise for for fairly fairly bright objects, diffraction realistic 37 Initial limited limited images images can can be be reconstructed. reconstructed.37 Initial expements expements have have also telescopes. 43 fromtelescopes.43 datafrom on data been carried out on
3.2.3. synthesis and synthesis reconstruction and Pupil reconstruction 3.2. 3. Pupil two-reconstruct aatwo to reconstruct want to may want one may which one Another case in in which Another case modulus is dimensional nonnegative dimensional nonnegativefunction function from from its its Fourier Fourier modulus is in in optical diffraction-limited In a adiffraction determination. In function determination. pupil function pupil -limited optical
system, the pointspread function function isis the the squared Fourier modulus modulus point-spread system, the transfer optical transfer Equivalently, the optical pupil function. Equivalently, system's pupil of the system's function. 44 Given the pupil function.'" function function isis the the autocorrelation autocorrelation of of the the pupil one image plane, one pointspreadfunction functionatat aa given given location location in in an image point-spread could use the iterative algorithm to retrieve the corresponding pupil the equivalent to mathematically equivalent function, inin a way to the that isismathematically way that function, astronomy problem. astronomy problem. Turning Turning this this problem problem around, around, one could use synthesize (design) the iterative iterative algorithm algorithm to to synthesize (design)aa pupil pupil function function that the possibly while possibly function while point-spread desired pointgiven, desired would spread function yield a given, would yield well. satisfying satisfying other other desirable desirable constraints as well. aperture over part of an aperture extent—measurement Finite extent 3.3. 3.3. Finite -measurement over
In a number problems, there there isis aa function function of of reconstruction problems, number of reconstruction In
(i)
Fourier Fig. 7.7. Reconstruction of a function from from its Fourier nonnegative function a nonnegative Reconstruction of Fig. Test object; (b) (b) modulus modulus of its Fourier (c) initransform; (c) Fourier transform; modulus, modulus. (a) Test -(f) reconstruction results results (d)-(f) test);(d) (firsttest); object(first the object of the estimate of tial estimate - number (g) initial initial estimate of of 600; (g) (f) 600; 230, (f) (e) 230, 20,(e) (d)20, iterations:(d) numberofoffiterations: the object (second -(i)reconstruction reconstructionresults results-number number of of (h)-(i) (second test); (h) the iterations: (h) 215. (i) 215. 2, (i) (h) 2, iterations:
reconstruct the wishes to reconstruct one wishes and one known known finite finite extent extent (or support) and Fourier in. the Fourier function with with resolution resolution appropriate appropriate to an aperture in function domain more were measurements were which measurements overwhich oneover theone than the complete than more complete actually taken. In larger simply larger aperture isis simply desired aperture the desired cases, the some cases, In some measurements were taken, and so one which measurements over which aperture over the aperture than the wishes to extrapolate extrapolate the the function's function's Fourier Fourier transform, transform, i.e., to obwishes to tain superresolution has made one has cases, one other cases, In other superresolution of the function. In one case one which case in which aperture, in filled aperture, measurements over measurements over aa partially filled and function, and wishes to interpolate interpolate the the Fourier Fourier transform transform of the function, wishes to dofunction dothe function response in the impulse response thereby improved impulse an improved thereby obtain an main.
superresolution or superresolution Extrapolation or 3.3.1. 3.3.1. Extrapolation The error-reduction error -reduction algorithm algorithm was was first first applied applied to to the the extrapolation extrapolation (or superresolution) problem writbeen writhas been Much has Gerchberg. 45 Much by Gerchberg.45 problem by error-reduction theerror specifically the algorithm,specifically ten about -reduction iterative algorithm, the iterative about the ways of algorithm, various ways including various problem, including this problem, relates to this algorithm, as it relates proofs of 2) and Sec. 2) understanding (see the the end end of Sec. and proofs algorithm (see understanding the algorithm convergence.1o4243 °46-48 this particularproblem, problem, the the nature nature of particular this ForFor convergence. 10' 12 ' 13 '46"48 algorithm by implement the algorithm possible to implement the constraints makes it possible by aa taking on the order seconds 10"9 seconds of10-9 order of processor49'50 taking feedback optical processor49,5o Smith case. Marks and Smith two-dimensional thetwo for the per iteration even -dimensional case. even for describe these these matters in detail volume. elsewhere in this volume. detail elsewhere Interpolation 3.3.2. Interpolation 3.3.2. In tomographic systems, many many projections projections of of the the object object imaging systems, tomographic imaging In slice yielding information about aa slice projection yielding each projection measured, each are measured, through the Fourier transform measurements When measurements object. When the object. of the transform of over only a limited cone of angles effective aperture the effective angles are made, the SP /EVol. Vol373 373 Transformations Transformations In in Optical S/gna/ Signal Processing 153 1981) )// 153 Processing f(1981 SPIE
Downloaded From: http://opticalengineering.spiedigitallibrary.org/ on 11/25/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx
in Fourier domain domain has gaps, and and the the impulse impulse response response of in the Fourier of the the system system isis highly highlyirregular. irregular. In In applying applying the iterative algorithm to this this problem, 51 ' 52 the problem,51,52 thefunction-domain function -domainconstraint constraintisis the the finite finite extent extent and and nonnegativity nonnegativity of of the the object, object, and and the the Fourier Fourier domain domain constraint constraint isis that that the the Fourier Fourier transform transform equal equal the the measured measured Fourier Fourier transform transform over over the the measurement measurement aperture. problem similar similar to tomography problem A problem to the tomography in radio problem arises arises in astronomy. sky brightness brightness map map isisaatwo two-dimensional astronomy. The The radio sky -dimensional real, real, nonnegative nonnegativefunction function which whichisisthe the Fourier Fourier transform transform of of the the complex complex visibility visibilityfunction. function. The The visibility visibility function function isis measured measured by by radio interferornetry, interferometry, and in the case case of -baseline inoflong long-baseline ter ferometry, the visibility visibility function is measured only over a limited terferometry, "tracks"ininthe theFourier Fourierdomain, domain,resulting resultinginina apartially partially-filled set of "tracks" -filled effective aperture. err or-reduction has been been used used to to effective aperture. The error -reduction algorithm has obtain in effect, effect, interpolating interpolating the obtain improved improved maps maps by, by, in the visibility visibility function fill in the area between between the tracks.53 function to fill tracks. 53 For this problem, that itit be be nonnegative nonnegative and the constraints on the brightness map are that field of view. view. In the visibility visibility plane, be zero zero outside outside the the known known field plane, the the constraint is that the the complex complex visiblity visiblity function equal equal the the measured measured value within value within the the area area of of the the tracks. 3.4. Modulus quantized values 3.4. Modulusquantized values As spectrum, shaping, in in comcomAs mentioned mentioned earlier earlier in in connection with spectrum puter holography one may wish wish to to encode encode the the Fourier transform of an image as computer-generated an image as a computer -generatedhologram, hologram,but but some some types types of of computer-generated computer -generatedholograms hologramscan can encode encode only only certain, certain quantized quantized complex kinoform example example discussed discussed earlier earlierisis aa special complex values. values. The kinoform special type quantization. A A more more general general example example is is the the Lohmann Lohrnann type of quantization. hologram, 54 54for forwhich whichthe themodulus modulusand and phase phase of a complex complex sample sample are determined respectively, of of an determined by the area and relative position, respectively, aperture within a sampling sampling cell. cell. The The number number of of allowable allowable quantized quantized values resolution elements, elements, of the values isis determined determined by by the the number number of resolution recording device recording deviceused usedto to fabricate fabricate the the hologram, hologram, used used to form form, one one cell. synthesis problem, function-domain cell. For this synthesis problem, the functiondomain constraint is is that the modulus of the function function equal equal the the desired desired image image modulus and the Fourier Fourier-domain and the -domain constraint constraint isis that that the the complex complex Fourier Fourier coefficients prescribed set set of coefficientsfall fall on on aa prescribed of quantized quantized values. values. ExExperiments synthesizing such transform periments have have shown shown that synthesizing such aa Fourier transform pair is possible using using the iterative algorithm.55,7 algorithm,. 55 ' 7 For example, Fig. 8(a) image produced produced by Lohrnann 8(a) shows shows aa simulation simulation of of an an image by a Lohmann hologram having modulus and four four phase hologram having only only four four modulus phase quantization quantization levels levels when, when the the image image was wasrandom, random phase coded. Figure Figure 8(b) 8(b) shows the image after 13 1,3 iterations, iterations, aa considerable considerable improvement. improvement. This This the image problem general class class of problems problem isis one one of a more general the problems regarding regarding the transmission transmission of of coded coded data. 3.5. Finite Finite extent extent -phase phase 3.5. Finally, the iterative algorithm Finally, the iterative algorithm has has been been used used to reconstruct reconstruct the modulus of a band band-limited from its its phase.56.57 phase. 56 ' 57 Or, -limited signal from Or, looking at itit in way, given given that a function has has finite finite extent extent and and given in another way, given the reconstruct the the phase phase of of its its Fourier Fourier transform, transform, reconstruct the modulus its modulus of its Fourier shown that Fourier transform. transform. For this application, it has been shown that for a wide is unique.56 unique. 56 This application wide class class of of conditions conditions the solution is will Sec. 4. will be be discussed discussed further further in Sec.
4. ALGORITHM CONVERGENCE AND 4. ALGORITHM ACCELERATED ACCELERATED ALGORITHMS ALGORITHMS As mentioned mentioned in Sec. 2, the the basic basic iterative iterative algorithm algorithm depicted As in Sec. depicted in in Fig. 1, referred has been been shown Fig. 1, referred to as the error-reduction error -reduction algorithm, has shown to converge section, the the convergence convergence converge for for some some applications. applications. In this section, is proven is proven for for all all applications. applications. In In addition, modified algorithms that often error-reduction often converge converge much much faster faster than than the error -reduction algorithm algorithm are are discussed. Convergence of the the errorerror-reduction algorithm 4.1. Convergence reduction algorithm For the error error-reduction mean-squared -reduction algorithm, algorithm, the mean -squared error error can be defined in It is is aa normalized normalized version version defined in general general by by Eq. Eq. (7) (7) or or Eq. (8). It the integral integral over the square of of the the amount amount by bywhich which the of the thecorncorn-
Fig. Computer-simulated Fig. 8. 8. Computer -simulated images images from from hologram hologram with with four magnitude magnitude and and four phase phase quantized quantized levels. levels, (a) (a) Object Object random random phased phased coded; coded; (b) (b) after after 13 13 iterations iterations of the the iterative iterative method. method.
puted function function (or (or the the computed computed Fourier Fourier transform) transform) violates violates the the puted appropriatedomain. domain.When Whenthe constraints in the appropriate -squared erthemean mean-squared eris zero, zero, then then aa Fourier Fourier transform transform pair pair has has been ror is been found found that that satisfies all satisfies all the the constraints constraints in in both domains. Consider steps in in the error-reduction the error -reduction algorithm algorithm Consider again again the the steps described k th iteration starts starts with with an an estimate estimate gk(x) described in in Sec. Sec. 2. 2. The kth g^(x) that satisfies satisfies the function-domain the function -domain constraints. constraints. For any coordinate, the complex complex values values that g(x) g(x) can have have that that satisfy satisfy the the function function-x, the domain space. For domain constraints constraints form form some some set set of of points points in in phasor space. For example, modulus must equal equal f(x) f(x)1, then the set |, then example, ifif the the modulus set of of such such points circle of radius | f(x) f(x)1| in phasor space; space; if if the the function function points isis a circle must must be be nonnegative, nonnegative, then then the set of such points is line on on is the half line the nonnegative real axis. The The function function estimate estimate gk(x) gk (x) is is Fourier Fourier the nonnegative real axis. transformed, yielding G k (u). The transformed, yielding Gk(u). The next next step step in the algorithm algorithm is is to form Gk(u) G^(u) by changing changing G^(u) Gk(u) by by the the smallest smallest possible possibleamount amount that that Fourier-domain constraints. Gk(u) G£(u)isis then then ininallows it to satisfy the Fourierdomain constraints. verse Fourier transformed, yielding yielding gk(x) g£(x) in the function verse function domain. domain. In the final final step, step, gk g^ +1(x) + jM is is formed In the by the formed by by changing changing gk(x) g£(x) by the smallest amount that that allows allows itit to tosatisfy satisfythe smallest -domain conthefunction function-domain constraints. Now Now consider consider the unnormalized unnorrnalized squared squared error, error, given given by by straints. (7) and (8). the numerators in Eqs. (7) (8). In the Fourier domain, domain, the ununnormalized squared k th iteration iteration is is normalized squared error at the kth co 2 _ eFk -
G k (u)IGk(u) - Gk(u) 1 2 du du _00
154 SPIE Vol. Vol. 373 373 Transformations Transformations in in Optical Optical Signal Signal Processing Processing (1981) 154 //SPIE
Downloaded From: http://opticalengineering.spiedigitallibrary.org/ on 11/25/2015 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx
(14) (14)
•/
oc
v^/
dx,, |gk(x)-g£(x)| = I I gk(x) - gk(x) 122dx -00
Parseval's from Parseval's results from equation results this equation second line where the the second line in in this where theorem. The The unnormalized unnormalized squared squared error error in in the function domain theorem. given by is given at the kth iteration is w 2
= eOk = e0k
fI gk
dx . + i(x)-g£(x)| + l(x) - gk(x) I 22dx.
INPUT 9 g
(15)
-03
function-domain the function satisfy the -domain i(x) by definition satisfy + i(x) gk + Both gk(x) gk(x) and gk j(x) is gk++ 1(x) coordinatex,x,gk given coordinate any given at any constraints. Also at is the point in function-domain satisfying the function phasor -domain constraints constraints that isis space satisfying phasor space closest to gk(x). gk(x). Therefore, Therefore, for all values of x, closest to +1W II gk gk + 1(x) --
g£ gk(x) -- gk(x) gk WI I
