Abstract

The phase-retrieval problem, fundamental in applied physics and engineering, addresses the question of how to determine the phase of a complex-valued function from modulus data and additional a priori information. Recently we identified two important methods for phase retrieval, namely, Fienup’s basic input–output and hybrid input–output (HIO) algorithms, with classical convex projection methods and suggested that further connections between convex optimization and phase retrieval should be explored. Following up on this work, we introduce a new projection-based method, termed the hybrid projection–reflection (HPR) algorithm, for solving phase-retrieval problems featuring nonnegativity constraints in the object domain. Motivated by properties of the HPR algorithm for convex constraints, we recommend an error measure studied by Fienup more than 20 years ago. This error measure, which has received little attention in the literature, lends itself to an easily implementable stopping criterion. In numerical experiments we found the HPR algorithm to be a competitive alternative to the HIO algorithm and the stopping criterion to be reliable and robust.

© 2003 Optical Society of America

Full Article  |  PDF Article

References

  • View by:
  • |
  • |
  • |

  1. R. W. Gerchberg, W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Stuttgart) 35, 237–246 (1972).
  2. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982).
    [CrossRef] [PubMed]
  3. A. Levi, H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984).
    [CrossRef]
  4. A. Levi, H. Stark, “Restoration from phase and magnitude by generalized projections,” in Image Recovery: Theory and Application, H. Stark, ed. (Academic, Orlando, Fla., 1987), pp. 277–320.
  5. P. L. Combettes, “The foundations of set theoretic estimation,” Proc. IEEE 81, 182–208 (1993).
    [CrossRef]
  6. D. C. Youla, H. Webb, “Image restoration by the method of convex projections: Part I—theory,” IEEE Trans. Med. Imaging MI-1, 81–94 (1982).
    [CrossRef]
  7. J. P. Boyle, R. L. Dykstra, “A method for finding projections onto the intersection of convex sets in Hilbert spaces,” Lect. Notes Stat. 37, 28–47 (1986).
  8. P.-L. Lions, B. Mercier, “Splitting algorithms for the sum of two nonlinear operators,” SIAM J. Numer. Anal. 16, 964–979 (1979).
    [CrossRef]
  9. H. H. Bauschke, P. L. Combettes, D. R. Luke, “Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization,” J. Opt. Soc. Am. A 19, 1334–1345 (2002).
    [CrossRef]
  10. R. P. Millane, “Phase retrieval in crystallography and optics,” J. Opt. Soc. Am. A 7, 394–411 (1990).
    [CrossRef]
  11. In general, the object-domain space need not be restricted to real-valued signals. For a review of a phase-retrieval application in which the iterates are complex valued, see Ref. 12.
  12. D. R. Luke, J. V. Burke, R. Lyon, “Optical wavefront reconstruction: theory and numerical methods,” SIAM Rev. 44, 169–224 (2002).
    [CrossRef]
  13. H. Stark, ed., Image Recovery: Theory and Application (Academic, Orlando, Fla.1987).
  14. Typically, 0.5≤β≤1.
  15. It was pointed out in Remark 4.1 of Ref. 9that a more literal reformulation of Eq. (13) leads to (31)(∀t∈X)xn+1(t)=(PM(xn))(t)if t∈D  or (PM(xn))(t)=0xn(t)-β(PM(xn))(t)otherwise.Under the assumption that the zero crossings of PM(xn)outside of Dare negligible, Eq. (31) reduces to Eq. (14). In digital computing, this assumption is justified by the fact that the probability of obtaining zero numbers is virtually zero.
  16. If m=0,then x=0is the unique solution of the corresponding phase-retrieval problem.
  17. For the BIO algorithm, this was already pointed out in Remark 5.6 of Ref. 9.
  18. We stress that monitoring the sequences (PM(xn))n∈Nand (PS+PM(xn))n∈Nis well motivated from the convex consistent setting. Replace the nonconvex set Mand its corresponding projector PMwith the convex set Band the corre-sponding projector PB.If S+∩B≠∅,then, using the results in Ref. 9, one can prove that ES+(xn)→0with equality precisely when PB(xn)∈S+∩B.However, if S+∩B=∅,which is likely a better approximation of the geometry of the phase-retrieval problem, then minimizing ES+(xn)corresponds to finding a displacement vector for S+and Bin the sense of Ref. 19. If the problem is feasible but ES+(xn)is positive, then the algorithm has stagnated, and the value of ES+(xn)is an indicator of the quality of the stagnation point.
  19. H. H. Bauschke, J. M. Borwein, “On the convergence of von Neumann’s alternating projection algorithm for two sets,” Set-Valued Anal. 1, 185–212 (1993).
  20. J. C. Dainty, J. R. Fienup, “Phase retrieval and image reconstruction for astronomy,” in Image Recovery: Theory and Application, H. Stark, ed. (Academic, Orlando, Fla., 1987), pp. 231–275.
  21. T. R. Crimmins, J. R. Fienup, B. J. Thelen, “Improved bound on object support from autocorrelation support and application to phase retrival,” J. Opt. Soc. Am. A 7, 3–13 (1990).
    [CrossRef]
  22. J. R. Fienup, C. C. Wackerman, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3, 1897–1907 (1986).
    [CrossRef]

2002

1993

H. H. Bauschke, J. M. Borwein, “On the convergence of von Neumann’s alternating projection algorithm for two sets,” Set-Valued Anal. 1, 185–212 (1993).

P. L. Combettes, “The foundations of set theoretic estimation,” Proc. IEEE 81, 182–208 (1993).
[CrossRef]

1990

1986

J. R. Fienup, C. C. Wackerman, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3, 1897–1907 (1986).
[CrossRef]

J. P. Boyle, R. L. Dykstra, “A method for finding projections onto the intersection of convex sets in Hilbert spaces,” Lect. Notes Stat. 37, 28–47 (1986).

1984

1982

D. C. Youla, H. Webb, “Image restoration by the method of convex projections: Part I—theory,” IEEE Trans. Med. Imaging MI-1, 81–94 (1982).
[CrossRef]

J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982).
[CrossRef] [PubMed]

1979

P.-L. Lions, B. Mercier, “Splitting algorithms for the sum of two nonlinear operators,” SIAM J. Numer. Anal. 16, 964–979 (1979).
[CrossRef]

1972

R. W. Gerchberg, W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Stuttgart) 35, 237–246 (1972).

Bauschke, H. H.

H. H. Bauschke, P. L. Combettes, D. R. Luke, “Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization,” J. Opt. Soc. Am. A 19, 1334–1345 (2002).
[CrossRef]

H. H. Bauschke, J. M. Borwein, “On the convergence of von Neumann’s alternating projection algorithm for two sets,” Set-Valued Anal. 1, 185–212 (1993).

Borwein, J. M.

H. H. Bauschke, J. M. Borwein, “On the convergence of von Neumann’s alternating projection algorithm for two sets,” Set-Valued Anal. 1, 185–212 (1993).

Boyle, J. P.

J. P. Boyle, R. L. Dykstra, “A method for finding projections onto the intersection of convex sets in Hilbert spaces,” Lect. Notes Stat. 37, 28–47 (1986).

Burke, J. V.

D. R. Luke, J. V. Burke, R. Lyon, “Optical wavefront reconstruction: theory and numerical methods,” SIAM Rev. 44, 169–224 (2002).
[CrossRef]

Combettes, P. L.

Crimmins, T. R.

Dainty, J. C.

J. C. Dainty, J. R. Fienup, “Phase retrieval and image reconstruction for astronomy,” in Image Recovery: Theory and Application, H. Stark, ed. (Academic, Orlando, Fla., 1987), pp. 231–275.

Dykstra, R. L.

J. P. Boyle, R. L. Dykstra, “A method for finding projections onto the intersection of convex sets in Hilbert spaces,” Lect. Notes Stat. 37, 28–47 (1986).

Fienup, J. R.

Gerchberg, R. W.

R. W. Gerchberg, W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Stuttgart) 35, 237–246 (1972).

Levi, A.

A. Levi, H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984).
[CrossRef]

A. Levi, H. Stark, “Restoration from phase and magnitude by generalized projections,” in Image Recovery: Theory and Application, H. Stark, ed. (Academic, Orlando, Fla., 1987), pp. 277–320.

Lions, P.-L.

P.-L. Lions, B. Mercier, “Splitting algorithms for the sum of two nonlinear operators,” SIAM J. Numer. Anal. 16, 964–979 (1979).
[CrossRef]

Luke, D. R.

Lyon, R.

D. R. Luke, J. V. Burke, R. Lyon, “Optical wavefront reconstruction: theory and numerical methods,” SIAM Rev. 44, 169–224 (2002).
[CrossRef]

Mercier, B.

P.-L. Lions, B. Mercier, “Splitting algorithms for the sum of two nonlinear operators,” SIAM J. Numer. Anal. 16, 964–979 (1979).
[CrossRef]

Millane, R. P.

Saxton, W. O.

R. W. Gerchberg, W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Stuttgart) 35, 237–246 (1972).

Stark, H.

A. Levi, H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984).
[CrossRef]

A. Levi, H. Stark, “Restoration from phase and magnitude by generalized projections,” in Image Recovery: Theory and Application, H. Stark, ed. (Academic, Orlando, Fla., 1987), pp. 277–320.

Thelen, B. J.

Wackerman, C. C.

Webb, H.

D. C. Youla, H. Webb, “Image restoration by the method of convex projections: Part I—theory,” IEEE Trans. Med. Imaging MI-1, 81–94 (1982).
[CrossRef]

Youla, D. C.

D. C. Youla, H. Webb, “Image restoration by the method of convex projections: Part I—theory,” IEEE Trans. Med. Imaging MI-1, 81–94 (1982).
[CrossRef]

Appl. Opt.

IEEE Trans. Med. Imaging

D. C. Youla, H. Webb, “Image restoration by the method of convex projections: Part I—theory,” IEEE Trans. Med. Imaging MI-1, 81–94 (1982).
[CrossRef]

J. Opt. Soc. Am. A

Lect. Notes Stat.

J. P. Boyle, R. L. Dykstra, “A method for finding projections onto the intersection of convex sets in Hilbert spaces,” Lect. Notes Stat. 37, 28–47 (1986).

Optik (Stuttgart)

R. W. Gerchberg, W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Stuttgart) 35, 237–246 (1972).

Proc. IEEE

P. L. Combettes, “The foundations of set theoretic estimation,” Proc. IEEE 81, 182–208 (1993).
[CrossRef]

Set-Valued Anal.

H. H. Bauschke, J. M. Borwein, “On the convergence of von Neumann’s alternating projection algorithm for two sets,” Set-Valued Anal. 1, 185–212 (1993).

SIAM J. Numer. Anal.

P.-L. Lions, B. Mercier, “Splitting algorithms for the sum of two nonlinear operators,” SIAM J. Numer. Anal. 16, 964–979 (1979).
[CrossRef]

SIAM Rev.

D. R. Luke, J. V. Burke, R. Lyon, “Optical wavefront reconstruction: theory and numerical methods,” SIAM Rev. 44, 169–224 (2002).
[CrossRef]

Other

H. Stark, ed., Image Recovery: Theory and Application (Academic, Orlando, Fla.1987).

Typically, 0.5≤β≤1.

It was pointed out in Remark 4.1 of Ref. 9that a more literal reformulation of Eq. (13) leads to (31)(∀t∈X)xn+1(t)=(PM(xn))(t)if t∈D  or (PM(xn))(t)=0xn(t)-β(PM(xn))(t)otherwise.Under the assumption that the zero crossings of PM(xn)outside of Dare negligible, Eq. (31) reduces to Eq. (14). In digital computing, this assumption is justified by the fact that the probability of obtaining zero numbers is virtually zero.

If m=0,then x=0is the unique solution of the corresponding phase-retrieval problem.

For the BIO algorithm, this was already pointed out in Remark 5.6 of Ref. 9.

We stress that monitoring the sequences (PM(xn))n∈Nand (PS+PM(xn))n∈Nis well motivated from the convex consistent setting. Replace the nonconvex set Mand its corresponding projector PMwith the convex set Band the corre-sponding projector PB.If S+∩B≠∅,then, using the results in Ref. 9, one can prove that ES+(xn)→0with equality precisely when PB(xn)∈S+∩B.However, if S+∩B=∅,which is likely a better approximation of the geometry of the phase-retrieval problem, then minimizing ES+(xn)corresponds to finding a displacement vector for S+and Bin the sense of Ref. 19. If the problem is feasible but ES+(xn)is positive, then the algorithm has stagnated, and the value of ES+(xn)is an indicator of the quality of the stagnation point.

J. C. Dainty, J. R. Fienup, “Phase retrieval and image reconstruction for astronomy,” in Image Recovery: Theory and Application, H. Stark, ed. (Academic, Orlando, Fla., 1987), pp. 231–275.

In general, the object-domain space need not be restricted to real-valued signals. For a review of a phase-retrieval application in which the iterates are complex valued, see Ref. 12.

A. Levi, H. Stark, “Restoration from phase and magnitude by generalized projections,” in Image Recovery: Theory and Application, H. Stark, ed. (Academic, Orlando, Fla., 1987), pp. 277–320.

Cited By

OSA participates in CrossRef's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (13)

Fig. 1
Fig. 1

Original images and corresponding data used for the comparison of the HIO and HPR algorithms. (a) A 128×128 pixel image of a satellite used in Ref. 20. The center of (d) is a 38×38 pixel section of the standard Lena image, zero padded to a 128×128 matrix. (b) and (e) The Fourier magnitude data m corresponding to images (a) and (d), respectively. The same object-domain support constraint of size 64×64 pixels, shown in (c), is used for each example. An initial image x0 is shown in (f).

Fig. 2
Fig. 2

Asymptotic behavior of the HIO and HPR algorithms (β=0.75 and 1.0) over 100 random trials for the original image shown in Fig. 1(a). Shown here is the mean value of the performance measure ES+ defined by Eq. (30) at each iteration.

Fig. 3
Fig. 3

Asymptotic behavior of the HIO and HPR algorithms (β=0.75 and 1.0) over 100 random trials for the original image shown in Fig. 1(d). Shown here is the mean value of ES+ defined by Eq. (30) at each iteration.

Fig. 4
Fig. 4

Typical images recovered with 600 iterations of the HIO and HPR algorithms for the data corresponding to Fig. 1(a) and a random initial guess similar to that shown in Fig. 1(f). <./>(a) and (b) Reconstructed with the HIO algorithm for β=0.75 and β=1.0, respectively. (c) and (d) Images reconstructed with the HPR algorithm for β=0.75 and β=1.0, respectively. The rotation of all recovered images from the original image, with the exception of (d), reflects the nonuniqueness in recovery from phase.

Fig. 5
Fig. 5

Worst images recovered by the HIO and HPR algorithms for the data corresponding to Fig. 1(a) and a random initial guess similar to that shown in Fig. 1(f). (a) and (b) Images reconstructed with the HIO algorithm for β=0.75 and β=1.0, respectively. (c) and (d) Images reconstructed with the HPR algorithm for β=0.75 and β=1.0, respectively. The rotation of all recovered images from the original image, with the exception of (b), reflects the nonuniqueness in recovery from phase.

Fig. 6
Fig. 6

Typical images recovered with 600 iterations of the HIO and HPR algorithms for the data corresponding to Fig. 1(d) and a random initial guess similar to that shown in Fig. 1(f). (a) and (b) Images reconstructed with the HIO algorithm for β=0.75 and β=1.0, respectively. (c) and (d) Images reconstructed with the HPR algorithm for β=0.75 and β=1.0, respectively.

Fig. 7
Fig. 7

By use of the stagnated HIO (β=1.0) iterate shown in Fig. 6(b) as the initial point, the HPR (β=1.0) algorithm finds a stagnation point with a much better picture quality.

Fig. 8
Fig. 8

Behavior of the HIO and HPR algorithms (β=0.75 and 1.0) with noise (SNR=34 dB) over 100 random trials for the original image shown in Fig. 1(a). Shown here is the mean value of ES+ at each iteration.

Fig. 9
Fig. 9

Behavior of the HIO and HPR algorithms (β=0.75 and 1.0) with noise (SNR=34 dB) over 100 random trials for the original image shown in Fig. 1(d). Shown here is the mean value of ES+ at each iteration.

Fig. 10
Fig. 10

Typical images recovered by the HIO (β=1) and HPR (β=1) algorithms with noise (SNR=34 dB) for the data corresponding to Fig. 1(d) and the normalized zero-phase initial guess shown in Fig. 1(c). (a), (b), and (c) Images reconstructed with 100, 200, and 400 iterations of the HIO algorithm, respectively. (d), (e), and (f) Images reconstructed with 100, 200, and 400 iterations of the HPR algorithm, respectively.

Fig. 11
Fig. 11

Typical images recovered by the HIO (β=0.75) and HPR (β=0.75) algorithms with noise (SNR=34 dB) for the data corresponding to Fig. 1(d) and the normalized zero-phase initial guess shown in Fig. 1(c). (a), (b), and (c) Images reconstructed with 100, 200, and 400 iterations of the HIO algorithm, respectively. (d), (e), and (f) Images reconstructed with 100, 200, and 400 iterations of the HPR algorithm, respectively.

Fig. 12
Fig. 12

Typical images recovered by coupling the HIO (β=1), HPR (β=1), and ER algorithms for the noisy data corresponding to Fig. 1(d) (SNR=34 dB) and the normalized zero-phase initial guess shown in Fig. 1(c). (a) Images reconstructed with 200 iterations of HIO followed by 200 iterations of ER. (b) Images reconstructed with 200 iterations of HIO followed by 200 iterations of HPR. (c) Images reconstructed with 200 iterations of HPR followed by 200 iterations of ER. (d) For reference, the reconstruction obtained with 200 iterations of HPR.

Fig. 13
Fig. 13

Typical images recovered by coupling the HIO (β=0.75), HPR (β=0.75), and ER algorithms for the noisy data corresponding to Fig. 1(d) (SNR=34 dB) and the normalized zero-phase initial guess shown in Fig. 1(c). (a) Images reconstructed with 200 iterations of HIO followed by 200 iterations of ER. (b) Images reconstructed with 200 iterations of HIO followed by 200 iterations of HPR. (c) Images reconstructed with 200 iterations of HPR followed by 200 iterations of ER. (d) For reference, the reconstruction obtained with 200 iterations of HPR.

Equations (38)

Equations on this page are rendered with MathJax. Learn more.

xn+1=T(xn),
findxSM,
M={yL||yˆ|=m},
PM(x)=F-1(y^0),
y^0(ω)=m(ω) xˆ(ω)|xˆ(ω)|ifxˆ(ω)0m(ω)otherwise.
{tX|x*(t)0}D.
S={yL|y1D=0}.
PS(x)=x1D.
S+={yL|y1D=0andy0}.
(tX)
(PS+(x))(t)=max{0, x(t)}iftD0otherwise.
findxS+M.
γn={tX|PM(xn) violates the object-domain constraint at t}.
(tX)
xn+1(t)=(PM(xn))(t)iftγnxn(t)-β(PM(xn))(t)iftγn.
(tX)xn+1(t)=(PM(xn))(t)iftDxn(t)-β(PM(xn))(t)otherwise.
xn+1=12(RS(RM+(β-1)PM)+I+(1-β)PM)(xn).
xn+1=1DPM(xn)+1D(xn-βPMxn)=1DPM(xn)+(1-1D)(xn-βPMxn)=1D((1+β)PM-I)(xn)+xn-βPM(xn)=PS((1+β)PM-I)(xn)+xn-βPM(xn)=12((2PS-I)((1+β)PM-I)+I+(1-β)PM)(xn)=12(RS(RM+(β-1)PM)+I+(1-β)PM)(xn),
xn+1=12(RSRM+I)(xn).
(tX)
xn+1(t)=(PM(xn))(t)iftDand(PM(xn))(t)0xn(t)-β(PM(xn))(t)otherwise.
xn+1=12(RS+RM+I)(xn).
xn+1=12(RS+(RM+(β-1)PM)+I+(1-β)PM)(xn).
xn+1=12(RS+(RM+(β-1)PM)+I+(1-β)PM)(xn)=12((2PS+-I)((1+β)PM-I)+I+(1-β)PM)(xn)=PS+((1+β)PM-I)(xn)+xn-βPM(xn)=(1D((1+β)PM-I)(xn))++xn-βPM(xn)=1DPM(xn)+1D(xn-βPMxn)-(1D((1+β)PM-I)(xn))-.
xn+1(t)=(PM(xn))(t)-(((1+β)PM-I)(xn))-(t)=(PM(xn))(t)if(1+β)(PM(xn))(t)xn(t)xn(t)-β(PM(xn))(t)otherwise.
xn+1(t)=xn(t)-β(PM(xn))(t).
(tX)xn+1(t)=(PM(xn))(t)iftDand(RM(xn))(t)(1-β)(PM(xn))(t)xn(t)-β(PM(xn))(t)otherwise.
(tX)
xn+1(t)=(PM(xn))(t)iftDand(RM(xn))(t)0xn(t)-(PM(xn))(t)otherwise.
 
(PM(xn))(t)0,
(RM(xn))(t)0.
Erms(xn)=|x^n|-mm,
(PM(xn))nN
(PS+(PM(xn)))nN.
ES+(xn)=PS+(PM(xn))-PM(xn)2PM(xn)2=PS+(PM(xn))-PM(xn)2m2.
(tX)
xn+1(t)=(PM(xn))(t)iftDor(PM(xn))(t)=0xn(t)-β(PM(xn))(t)otherwise.

Metrics