814 Chapter 18.Integral Equations and Inverse Theory ·known bounds(i.e.,uL(r)≤t(z)≤uw(x)for specified functions uL and uv). (In this last case,the bounds might be related to an initial estimate and its error bars, e.g.,uo()+(),where y is of order I or 2.)Notice that these,and similar, constraints can be either in the image space,or in the Fourier transform space,or(in fact)in the space of any linear transformation of u. If C;is a convex set,then Pi is called a nonexpansive projection operator onto that set if(i)Pi leaves unchanged any u already in Ci,and (ii)Pi maps any u outside Ci to the closest element of Ci.in the sense that 81 |P:i-t创≤lia-for all in C. (18.5.24) While this definition sounds complicated,examples are very simple:A nonexpansive B projection onto the set of positive u's is "set all negative components of u equal to zero."A nonexpansive projection onto the set of (z)'s bounded by uL()< ICAL ()<u()is "set all values less than the lower bound equal to that bound,and set all values greater than the upper bound equal to that bound."A nonexpansive RECIPES projection onto functions with compact support is "zero the values outside of the region of support.” 9 The usefulness of these definitions is the following remarkable theorem:Let C be the intersection of m convex sets C1,C2,....Cm.Then the iteration )=(PiP2...Pm)() (18.5.25 巴6 will converge to C from all starting points,as k-oo.Also,if C is empty (there is no intersection),then the iteration will have no limit point.Application of this theorem is called the method of projections onto comvex sets or sometimes POCS [7]. A generalization of the POCS theorem is that the Pi's can be replaced by a set of Ti's, T≡1+0：(P-1)0<6:<2 (18.5.26) A well-chosen set of Bi's can accelerate the convergence to the intersection set C. 会 Numerica 10621 Some inverse problems can be completely solved by iteration(18.5.25)alone! For example,a problem that occurs in both astronomical imaging and X-ray diffraction work is to recover an image given only the modulus of its Fourier transform (equivalent to its power spectrum or autocorrelation)and not the phase. 台第 Here three convex sets can be utilized:the set of all images whose Fourier transform has the specified modulus to within specified error bounds:the set of all positive images;and the set of all images with zero intensity outside of some specified region. In this case the POCS iteration(18.5.25)cycles among these three,imposing each constraint in turn;FFTs are used to get in and out of Fourier space each time the Fourier constraint is imposed. The specific application of POCS to constraints alternately in the spatial and Fourier domains is also known as the Gerchberg-Saxton algorithm [81.While this algorithm is non-expansive,and is frequently convergent in practice,it has not been proved to converge in all cases [9].In the phase-retrieval problem mentioned above, the algorithm often "gets stuck"on a plateau for many iterations before making sudden,dramatic improvements.As many as 104 to 105 iterations are sometimes814 Chapter 18. Integral Equations and Inverse Theory Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). • known bounds (i.e., uL(x) ≤ u(x) ≤ uU (x) for specified functions uL and uU ). (In this last case, the bounds might be related to an initial estimate and its error bars, e.g., u0(x) ± γσ(x), where γ is of order 1 or 2.) Notice that these, and similar, constraints can be either in the image space, or in the Fourier transform space, or (in fact) in the space of any linear transformation of u. If Ci is a convex set, then Pi is called a nonexpansive projection operator onto that set if (i) Pi leaves unchanged any u already in Ci, and (ii) Pi maps any u outside Ci to the closest element of Ci, in the sense that |Piu − u|≤|ua − u| for all ua in Ci (18.5.24) While this definition sounds complicated, examples are very simple: A nonexpansive projection onto the set of positive u’s is “set all negative components of u equal to zero.” A nonexpansive projection onto the set of u(x)’s bounded by u L(x) ≤ u(x) ≤ uU (x) is “set all values less than the lower bound equal to that bound, and set all values greater than the upper bound equal to that bound.” A nonexpansive projection onto functions with compact support is “zero the values outside of the region of support.” The usefulness of these definitions is the following remarkable theorem: Let C be the intersection of m convex sets C1, C2,...,Cm. Then the iteration u(k+1) = (P1P2 ···Pm)u(k) (18.5.25) will converge to C from all starting points, as k → ∞. Also, if C is empty (there is no intersection), then the iteration will have no limit point. Application of this theorem is called the method of projections onto convex sets or sometimes POCS [7]. A generalization of the POCS theorem is that the Pi’s can be replaced by a set of Ti’s, Ti ≡ 1 + βi(Pi − 1) 0 < βi < 2 (18.5.26) A well-chosen set of βi’s can accelerate the convergence to the intersection set C. Some inverse problems can be completely solved by iteration (18.5.25) alone! For example, a problem that occurs in both astronomical imaging and X-ray diffraction work is to recover an image given only the modulus of its Fourier transform (equivalent to its power spectrum or autocorrelation) and not the phase. Here three convex sets can be utilized: the set of all images whose Fourier transform has the specified modulus to within specified error bounds; the set of all positive images; and the set of all images with zero intensity outside of some specified region. In this case the POCS iteration (18.5.25) cycles among these three, imposing each constraint in turn; FFTs are used to get in and out of Fourier space each time the Fourier constraint is imposed. The specific application of POCS to constraints alternately in the spatial and Fourier domains is also known as the Gerchberg-Saxton algorithm [8]. While this algorithm is non-expansive, and is frequently convergent in practice, it has not been proved to converge in all cases [9]. In the phase-retrieval problem mentioned above, the algorithm often “gets stuck” on a plateau for many iterations before making sudden, dramatic improvements. As many as 10 4 to 105 iterations are sometimes