Abstract

The method of projections onto convex sets can be used to solve many problems in image restoration, e.g., restoration from phase, spectral extrapolation, and signal recovery in computer-aided tomography. However, image-restoration problems involving nonconvex constraints cannot be handled by the method of projection onto convex sets in a fashion that ensures convergence. The restoration-from-magnitude (RFM) problem is such a case. To handle the RFM as well as other nonconvex constraints, we describe an algorithm known as generalized projections and discuss its properties. When sets are nonconvex, it is possible for the algorithm to exhibit pathological behavior that is never manifest in convex projections. We introduce an error criterion called the summed-distance error (SDE) and show under what circumstances the SDE is a monotonically decreasing function of the number of iterations. Near-optimum performance of the algorithm is achieved by relaxation parameters. Comparisons with other RFM methods are furnished for synthetic imagery.

© 1984 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 35, 237–246 (1972).
  2. R. Gerchberg, “Super resolution through error energy reduction,” Opt. Acta 21, 709–720 (1974).
    [CrossRef]
  3. A. Papoulis, “A new algorithm in spectral analysis and band-limited extrapolation,” IEEE Trans. Circuits Syst. CAS-22, 735–742 (1975).
    [CrossRef]
  4. D. C. Youla, “Generalized image restoration by the method of alternating orthogonal projections,” IEEE Trans. Circuits Syst. CAS-25, 694–702 (1979).
  5. H. Stark, D. Cahana, H. Webb, “Restoration of arbitrary finite-energy optical objects from limited spatial and spectral information,” J. Opt. Soc. Am. 71, 635–642 (1981).
    [CrossRef]
  6. H. Stark, S. Cruze, G. Habetler, “Restoration of optical objects subject to nonnegative spatial or spectral constraints,” J. Opt. Soc. Am. 72, 993–1000 (1982).
    [CrossRef]
  7. L. M. Bregman, “Finding the common point of convex sets by the method of successive projections,” Dokl. Akad. Nauk SSSR 162, 487–490 (1965).
  8. L. G. Gubin, B. T. Polyak, E. V. Raik, “The method of projections for finding the common point of convex sets,” USSR Computational Math. Math. Phys. 7(6), 1–24 (1967).
    [CrossRef]
  9. A. Lent, H. Tuy, “An iterative method for the extrapolation of band-limited functions,” J. Math. Anal. Appl. 83, 554–565 (1981).
    [CrossRef]
  10. D. C. Youla, H. Webb, “Image restoration by the method of projections onto convex sets. Part I,” IEEE Trans. Med. Imaging TMI-1, 81–94 (1982).
    [CrossRef]
  11. We work in a Hilbert space ℋ with the norm of g ∈ ℋ denoted by ||g|| and the inner product of g ∈ ℋ and f ∈ ℋ denoted by (g, f). A sequence (fx) is said to converge weakly to f if (fK, g) → (f, g) for all g ∈ ℋ and to converge strongly to f if ||f− fK|| → 0.
  12. Youla showed weak convergence in the general case, but in some cases it is also strong convergence, as in the case when we deal with finite discrete sequences. (See Ref. 10).
  13. M. I. Sezan, H. Stark, “Image restoration by the method of projections onto convex sets. Part II,” IEEE Trans. Med. Imaging TMI-1, 95–101 (1982).
    [CrossRef]
  14. A. Levi, H. Stark, “Signal restoration from phase by projection onto convex sets,” J. Opt. Soc. Am. 73, 810–822 (1983).
    [CrossRef]
  15. M. H. Hayes, “The reconstruction of a multidimensional sequence from the phase of its fourier transform,” IEEE Trans. Acoust. Speech Signal Process. ASSP-30, 140–154 (1982).
    [CrossRef]
  16. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982).
    [CrossRef] [PubMed]
  17. A Fourier-transform pair is denoted by f(x) ↔ F(w).
  18. It should be understood that, in general, λ1and λ2depend on the iteration number n and may be denoted by λ1,n and λ2,n, respectively. However, in order to simplify notation throughout the paper we omit the subscript n. The optimal values for λi, i= 1, 2 are denoted by λi,optin the per-step optimization case and by λi* in the per-cycle optimization case.

1983

1982

H. Stark, S. Cruze, G. Habetler, “Restoration of optical objects subject to nonnegative spatial or spectral constraints,” J. Opt. Soc. Am. 72, 993–1000 (1982).
[CrossRef]

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

M. I. Sezan, H. Stark, “Image restoration by the method of projections onto convex sets. Part II,” IEEE Trans. Med. Imaging TMI-1, 95–101 (1982).
[CrossRef]

M. H. Hayes, “The reconstruction of a multidimensional sequence from the phase of its fourier transform,” IEEE Trans. Acoust. Speech Signal Process. ASSP-30, 140–154 (1982).
[CrossRef]

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

1981

H. Stark, D. Cahana, H. Webb, “Restoration of arbitrary finite-energy optical objects from limited spatial and spectral information,” J. Opt. Soc. Am. 71, 635–642 (1981).
[CrossRef]

A. Lent, H. Tuy, “An iterative method for the extrapolation of band-limited functions,” J. Math. Anal. Appl. 83, 554–565 (1981).
[CrossRef]

1979

D. C. Youla, “Generalized image restoration by the method of alternating orthogonal projections,” IEEE Trans. Circuits Syst. CAS-25, 694–702 (1979).

1975

A. Papoulis, “A new algorithm in spectral analysis and band-limited extrapolation,” IEEE Trans. Circuits Syst. CAS-22, 735–742 (1975).
[CrossRef]

1974

R. Gerchberg, “Super resolution through error energy reduction,” Opt. Acta 21, 709–720 (1974).
[CrossRef]

1972

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

1967

L. G. Gubin, B. T. Polyak, E. V. Raik, “The method of projections for finding the common point of convex sets,” USSR Computational Math. Math. Phys. 7(6), 1–24 (1967).
[CrossRef]

1965

L. M. Bregman, “Finding the common point of convex sets by the method of successive projections,” Dokl. Akad. Nauk SSSR 162, 487–490 (1965).

Bregman, L. M.

L. M. Bregman, “Finding the common point of convex sets by the method of successive projections,” Dokl. Akad. Nauk SSSR 162, 487–490 (1965).

Cahana, D.

Cruze, S.

Fienup, J. R.

Gerchberg, R.

R. Gerchberg, “Super resolution through error energy reduction,” Opt. Acta 21, 709–720 (1974).
[CrossRef]

Gerchberg, R. W.

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

Gubin, L. G.

L. G. Gubin, B. T. Polyak, E. V. Raik, “The method of projections for finding the common point of convex sets,” USSR Computational Math. Math. Phys. 7(6), 1–24 (1967).
[CrossRef]

Habetler, G.

Hayes, M. H.

M. H. Hayes, “The reconstruction of a multidimensional sequence from the phase of its fourier transform,” IEEE Trans. Acoust. Speech Signal Process. ASSP-30, 140–154 (1982).
[CrossRef]

Lent, A.

A. Lent, H. Tuy, “An iterative method for the extrapolation of band-limited functions,” J. Math. Anal. Appl. 83, 554–565 (1981).
[CrossRef]

Levi, A.

Papoulis, A.

A. Papoulis, “A new algorithm in spectral analysis and band-limited extrapolation,” IEEE Trans. Circuits Syst. CAS-22, 735–742 (1975).
[CrossRef]

Polyak, B. T.

L. G. Gubin, B. T. Polyak, E. V. Raik, “The method of projections for finding the common point of convex sets,” USSR Computational Math. Math. Phys. 7(6), 1–24 (1967).
[CrossRef]

Raik, E. V.

L. G. Gubin, B. T. Polyak, E. V. Raik, “The method of projections for finding the common point of convex sets,” USSR Computational Math. Math. Phys. 7(6), 1–24 (1967).
[CrossRef]

Saxton, W. O.

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

Sezan, M. I.

M. I. Sezan, H. Stark, “Image restoration by the method of projections onto convex sets. Part II,” IEEE Trans. Med. Imaging TMI-1, 95–101 (1982).
[CrossRef]

Stark, H.

Tuy, H.

A. Lent, H. Tuy, “An iterative method for the extrapolation of band-limited functions,” J. Math. Anal. Appl. 83, 554–565 (1981).
[CrossRef]

Webb, H.

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

H. Stark, D. Cahana, H. Webb, “Restoration of arbitrary finite-energy optical objects from limited spatial and spectral information,” J. Opt. Soc. Am. 71, 635–642 (1981).
[CrossRef]

Youla, D. C.

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

D. C. Youla, “Generalized image restoration by the method of alternating orthogonal projections,” IEEE Trans. Circuits Syst. CAS-25, 694–702 (1979).

Appl. Opt.

Dokl. Akad. Nauk SSSR

L. M. Bregman, “Finding the common point of convex sets by the method of successive projections,” Dokl. Akad. Nauk SSSR 162, 487–490 (1965).

IEEE Trans. Acoust. Speech Signal Process.

M. H. Hayes, “The reconstruction of a multidimensional sequence from the phase of its fourier transform,” IEEE Trans. Acoust. Speech Signal Process. ASSP-30, 140–154 (1982).
[CrossRef]

IEEE Trans. Circuits Syst.

A. Papoulis, “A new algorithm in spectral analysis and band-limited extrapolation,” IEEE Trans. Circuits Syst. CAS-22, 735–742 (1975).
[CrossRef]

D. C. Youla, “Generalized image restoration by the method of alternating orthogonal projections,” IEEE Trans. Circuits Syst. CAS-25, 694–702 (1979).

IEEE Trans. Med. Imaging

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

M. I. Sezan, H. Stark, “Image restoration by the method of projections onto convex sets. Part II,” IEEE Trans. Med. Imaging TMI-1, 95–101 (1982).
[CrossRef]

J. Math. Anal. Appl.

A. Lent, H. Tuy, “An iterative method for the extrapolation of band-limited functions,” J. Math. Anal. Appl. 83, 554–565 (1981).
[CrossRef]

J. Opt. Soc. Am.

Opt. Acta

R. Gerchberg, “Super resolution through error energy reduction,” Opt. Acta 21, 709–720 (1974).
[CrossRef]

Optik

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

USSR Computational Math. Math. Phys.

L. G. Gubin, B. T. Polyak, E. V. Raik, “The method of projections for finding the common point of convex sets,” USSR Computational Math. Math. Phys. 7(6), 1–24 (1967).
[CrossRef]

Other

We work in a Hilbert space ℋ with the norm of g ∈ ℋ denoted by ||g|| and the inner product of g ∈ ℋ and f ∈ ℋ denoted by (g, f). A sequence (fx) is said to converge weakly to f if (fK, g) → (f, g) for all g ∈ ℋ and to converge strongly to f if ||f− fK|| → 0.

Youla showed weak convergence in the general case, but in some cases it is also strong convergence, as in the case when we deal with finite discrete sequences. (See Ref. 10).

A Fourier-transform pair is denoted by f(x) ↔ F(w).

It should be understood that, in general, λ1and λ2depend on the iteration number n and may be denoted by λ1,n and λ2,n, respectively. However, in order to simplify notation throughout the paper we omit the subscript n. The optimal values for λi, i= 1, 2 are denoted by λi,optin the per-step optimization case and by λi* in the per-cycle optimization case.

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 (7)

Fig. 1
Fig. 1

A demonstration that the algorithm fn+1 = P1P2P3fn may not possess the set-distance reduction property.

Fig. 2
Fig. 2

Illustration of a trap and a tunnel for an algorithm of the form fn+1 = P1P2fn (a) Starting at the point f0, the sequence {fn} converges to a trap point T, whereas the true solution must belong to C 1 C 2. (b) Starting at the point f0, the algorithm enters into a long tunnel toward the solution at the point S.

Fig. 3
Fig. 3

Original image 1.

Fig. 4
Fig. 4

Original image 2.

Fig. 5
Fig. 5

Restoration of image 1 after 30 iterations with initial point f0 = 0. (a) Restoration by fn+1 = P1P2fn. (b) Restoration by fn+1 = P1T2fn. (c) Restoration by fn+1 = P1LP2fn. (d) Restoration by fn+1 = P1LT2fn.

Fig. 6
Fig. 6

Restoration of image 2 after 100 iterations with f0 = 0. (a) Restoration by fn+1 = P1P2fn (b) Restoration by fn+1 = P1T2fn. (c) Restoration by fn+1 = P1LP2fn. (d) Restoration by fn+1 = P1LT2fn.

Fig. 7
Fig. 7

Restoration of image 2 after 40 iterations with initial point f01(x, y) = 0.72 for 20 × 20 pixels in the center of the 30 × 30 pixels of the support S p of the function f. f01(x, y) = 0.36 for all other points of S p. (a) Restoration by fn+1 = P1LP2fn. (b) Restoration by fn+1 = P1LT2fn.

Tables (3)

Tables Icon

Table 1 Restoration of Image 1a

Tables Icon

Table 2 Restoration of Image 2a

Tables Icon

Table 3 Normalized SDE, J(fn), and λ2* as a Function of the Iteration Number n for the Restoration of Image 2a

Equations (72)

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

g - h = min y - h ,             over all y C i
f n + 1 = P 1 P 2 P m f n ,             f 0 arbitrary
f n + 1 = T 1 T 2 T m f n ,             f 0 arbitrary ,
T i = 1 + λ i ( P i - 1 ) ,             i = 1 , 2 , , m
f n + 1 = P 1 P 2 f n ,
e n = f n - f
f n + 1 = T 1 T 2 f n ,             f 0 arbitrary
d ( h , C ) = Min y - h ,             over all y C .
d ( h , C ) = P h - h .
J ( f n ) P 1 f n - f n + P 2 f n - f n ,
J ( T 2 f n ) = P 1 T 2 f n - T 2 f n + P 2 T 2 f n - T 2 f n .
J ( f n + 1 ) J ( T 2 f n ) J ( f n )
0 λ i A i 2 + A i A i 2 + A i - ½ ( A i + B i ) a i             i = 1 , 2 ,
A 1 P 1 T 2 f n - T 2 f n P 2 T 2 f n - T 2 f n ,
A 2 P 2 f n - f n P 1 f n - f n ,
B 1 ( P 2 T 2 f n - T 2 f n , P 1 T 2 f n - T 2 f n ) P 2 T 2 f n - T 2 f n 2 ,
B 2 ( P 1 f n - f n , P 2 f n - f n ) P 1 f n - f n 2 .
f n + 1 = P 1 P 2 f n ,             f 0 arbitrary
f n + 1 = P 1 P 2 P 3 f n ,             f 0 arbitrary .
C 1 = { g ( x ) :             g ( x ) = 0 , x S 1 c and g ( x ) = P ( x ) , x S 2 } ,
P 1 f = { 0 x S 1 c f ( x ) x S 1 S 2 c P ( x ) x S 2 , f ( x ) > 0 - P ( x ) x S 2 , f ( x ) 0
C 2 = { g ( x ) G ( ω ) :             ϕ ( ω = θ ( ω )             for ω S 1
G ( ω ) = M ( ω )             for ω S 2 } .
P 2 f ( x ) H ( ω ) = { F ( ω ) ω ( S 1 S 2 ) c M ( ω ) exp [ j θ ( ω ) ] ω S 1 S 2 exp [ j θ ( ω ) ] F ( ω ) cos [ θ ( ω ) - ϕ ( ω ) ] ω S 1 S 2 c , cos [ θ ( ω ) - ϕ ( ω ) [ 0 0 ω S 1 S 2 c , cos [ θ ( ω ) - ϕ ( ω ) ] < 0 M ( ω ) exp [ j ϕ ( ω ) ] ω S 2 S 1 c .
C 1 = { g ( x ) :             g ( x ) = 0             for x > a } ,
C 2 = { g ( x ) G ( ω ) :             G ( ω ) = M ( ω )             for all ω } .
P 1 g ( x ) = { g ( x ) x < a 0 x a
P 2 g ( x ) M ( ω ) exp [ j ϕ ( ω ) ] ,
f n + 1 = P 2 T 1 f n ,             f 0 arbitrary .
λ i , opt 1 ,             i = 1 , 2.
λ 1 * 1.
f n + 1 = T 2 T 1 f n ,             f 0 arbitrary .
f n + 1 = P 1 T 2 f n ,             f 0 arbitrary ,
J ( f n + 1 ) = P 2 P 1 T 2 f n - P 1 T 2 f n + P 1 P 1 T 2 f n - P 1 T 2 f n = P 2 P 1 T 2 f n - P 1 T 2 f n .
J 2 ( f n + 1 ) = - [ M ( ω ) - F n + 1 ( ω ) ] 2 d ω 2 π ,
F n + 1 ( ω ) = F ( P 1 { f n ( x ) + λ 2 [ P 2 f n ( x ) - f n ( x ) ] } )
F n + 1 ( ω ) = ( 1 - λ 2 ) F n ( ω ) + λ 2 F [ P 1 P 2 f n ( x ) ] ,
f n + 1 = T 1 T 2 T m f n = f n .
f n + 1 = P 1 P 2 f n ,             f 0 arbitrary .
P 1 P 2 f n = f n .
P 2 f n ( x ) rect ( x / 2 a ) = f n ( x ) ,
- M ( α ) exp [ j ϕ n ( α ) ] 2 sin [ ( ω - α ) a ] ω - α d α = A ( ω ) exp [ j ϕ n ( ω ) ] .
C 1 L = { g ( x ) :             g ( x ) = 0             for x > a             and             0 g ( x ) 1             for x a } ,
f = [ 1 2 π - M 2 ( ω ) d ω ] 1 / 2 .
J ( T 2 f n ) P 1 f n - T 2 f n + P 2 f n - T 2 f n .
J ( T 2 f n ) [ P 1 f n - f n 2 - 2 λ 2 ( P 1 f n - f n , P 2 f n - f n ) + λ 2 2 P 2 f n - f n 2 ] 1 / 2 + 1 - λ 2 P 2 f n - f n .
J ( T 2 f n ) - J ( f n ) [ ( 1 - 2 λ 2 B 2 + λ 2 2 A 2 2 ) 1 / 2 - 1 ] P 1 f n - f n + ( 1 - λ 2 - 1 ) P 2 f n - f n ,
( 1 - 2 λ 2 B 2 + λ 2 2 A 2 2 ) 1 / 2 A 2 + 1 - 1 - λ 2 A 2 .
- 1 A 2 λ 2 2 + 1 A 2 .
λ 2 ( A 2 2 - B 2 ) ( A 2 2 + A 2 ) ( 1 - 1 - λ 2 ) ,
2 A 2 > A 2 + B 2 0.
0 λ 2 A 2 2 + A 2 A 2 2 + A 2 - ½ ( A 2 + B 2 ) .
A 2 2 + A 2 A 2 2 + A 2 - ½ ( A 2 + B 2 ) A 2 2 + A 2 A 2 2 + 1 + 1 A 2 < 2 + 1 A 2 .
J ( T 2 f n ) = ( λ 2 + 1 - λ 2 ) P 2 f n - f n ,
λ 2 + 1 - λ 2 - 1 0
0 λ 2 1.
0 λ 1 A 1 2 + A 1 A 1 2 + A 1 - ½ ( A 1 + B 1 ) ,
J ( T i g ) J ( P i g )             for all λ i 1 , i = 1 , 2 ,
T i g - g + T i g - P i g = ( λ i 1 - λ i ) P i g - g .
T i g - g + T i g - P i g = P i g - g ,             0 λ i 1.
T i g - y < T i g - P i g .
P i g - g > T i g - y + T i g - g y - g .
T i g - y T i g - P i g             for all y C i
P i T i g = P i g             for 0 λ i 1.
J ( T i g ) = P 1 T i g - T i g + P 2 T i g - T i g P 1 T i g - P 2 T i g
J ( T 1 g ) P 1 g - P 2 T 1 g             0 λ 1 1 P 1 g - P 2 P 1 g = J ( P 1 g ) ,
J ( T 1 g ) P 1 T 1 g - P 2 T 1 g = P 1 g + λ 1 P 1 ( P 1 g - g ) - P 2 T 1 g .
J ( T 1 g ) P 1 g - P 2 P 1 g = J ( P 1 g )             for all λ 1 .
J 2 ( f n ) = P 2 f n - f n 2 = 1 2 π - [ M ( ω ) - A ( ω ) ] 2 d ω ,
e n 2 = f - f n 2 = 1 2 π - M ( ω ) exp [ j ϕ ( ω ) ] - A ( ω ) exp [ j ϕ n ( ω ) ] 2 d ω .
e n 2 = 1 2 π - [ M ( ω ) - A ( ω ) ] 2 d ω + 4 2 π - M ( ω ) A ( ω ) sin 2 [ ϕ ( ω ) - ϕ n ( ω ) 2 ] d ω .
e n 2 J 2 ( f n ) .

Metrics