and Numerical Results - Semantic Scholar
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. MI-1, NO. 2, OCTOBER 1982
95
Restoration by the Method of Convex Projections: Part 2-Applications
Image
and
Numerical Results M.
I.
SEZAN
AND
Abstract-The image restoration theory discussed in a previous paper by Youla and Webb [ 11 is applied to a simulated image and the results compared with the well-known method known as the GerchbergPapoulis algorithm. The results show that the method of image restoration by projection onto convex sets, by providing a convenient technique for utilizing a priori information, performs significandy better than the Gerchberg-Papoulis method.
I. INTRODUCTION -N a number of important practical situations it is required to restore a signal or image from incomplete information. The problem of restoring an image from its phase, the reconstruction of a tomographic image from incomplete view data, the restoration of a sharp image from its low-pass projection, and the extrapolation or interpolation of power spectra beyond the range of the measuring instrumentation are some significant examples. For this reason the Gerchberg-Papoulis (GP) algorithms [2], [3] and their offsprings [4], [5] have generated considerable interest. A generalization of these linear iterative algorithms was furnished by Youla [6] in which he also addressed the important questions of ill-posedness or noise-sensitivity of the restoration procedure. The restoration theory discussed by Youla and Webb in [1] has a very significant advantage over the GP and related algorithms in that it enables a large number of a priori known constraints to be incorporated in the algorithm through the mechanism of projection onto a convex set. Not all constraints can be handled this way. For example, the operation of digitizing a signal is not equivalent to projection onto a convex set. Nor, for that matter is the problem of reconstructing a signal from a prescribed magnitude. However, numerous other constraints can be treated this way and we shall illustrate this with examples in this paper. Eleven signal and image constraints and their associated projection operators are furnished in [11 and many of these are of physical significance. II. SETS AND PROJECTIONS
Because we shall deal with space-limited objects, defined by their reflectance, transmittance, or absorbance over a region Q, the functions of interest in this study will be assumed to be
H. STARK
members of L2x2(Q2), the space of all functions f(x,y) of twovariables x, y, square-integrable over Q. The associated Hilbert space JC is L2X2(Q) with inner product and norm defined by, respectively,
ffg(x, y) h*(x, y) dx dy
(g, h)
llgll A [(g,g)] /2.
(1) (2)
The following sets will be used. 1) C : The set of allf 's that vanish a.e. outside a prescribed region SC Q. C1 is an extension of the set 'T introduced in [1] and is a closed linear manifold (CLM) without interior. Given an arbitrary f E AC its projection onto C I is realized by Plf=
(x, y) EE 8
.4
05
(3)
(X,Y)O.
2) C2: The set of all f's in JC whose Fourier transforms assume a prescribed value G over a closed region £ in the u-v Fourier plane. This set is essentially the "dual," i.e., the frequency domain analog, of Ca (g)-the linear variety introduced in [1]. The projection of an arbitrary f E XC onto C2 is realized by
G(u, v), (u, V)E.C (F(u, v), (UIV)0
(4)
-
where F(u, v) = -f[f(x, y)] etc. for G(u, v), and Y is the Fourier transform operator. 3) C3: The set of all nonnegative functions in AC that satisfy the energy constraint
f If(x,y) 12 dXdyE_p2
(5)
where p = \/Eh. To show convexity we write
llfl +(1-A)f2 11ip()llflp=p
(6) +(1 - P) 2 11 0. Then
f
~
~~ \
In -f12 > AIInl2 + IIfI12 - 21(fn,f)I > ilfnII2 + IIfIh2 - 211AI1 lifll =
(Ilf 1- IIfl1)2 > e2.
But this leads to a contradiction since e
=0.
llfn - fli -
0. Thus
Now consider the projection onto C 3 of an arbitrary f E K. We show in the Appendix that this projection is realized by the following rule: (0, fA E
E
(8)
E,lZ
where f1 is the real part of f, fl is the rectified portion off1, andEt is the energy inf+, i.e.,
El _
Jf(f )2 dx dy.
(9)
4) C4: The subset of all functions in XC that are nonnegative. This set, already discussed in [1] is a closed convex set with vertex 0 (the zero function). The projection onto this set is
If,0,
>0(f=fl +jf2) (10) otherwise. 5) C5: The set of allf's in KC whose amplitudes must lie in a prescribed closed interval [a, b] a > 0, b > 0, a < b. To show convexity of e S note that
pf
fl
A1f, + (I - p) f2 > pxa + (1 - ,u) a = a
(11)
and
C"
(12)
any fi and f2 E es. To demonstrate closure we note that3 is open, hence Cs is closed. The projection onto Cs is
realized by the following rule:
21 6 is the set consisting of the elements in a that are not in 6. 3The complement of a set A is denoted by A c, -
f(x,y) < a
(a,
sf = if(x,y), tb,
a b.
We show in subsequent sections how the operators Pi and more generally the operators T1 _ 1 + X1(P1 - 1) i = 1, * *, 5 can be used, either sequentially or as compositions, to restore the missing data in the image. -
III. OBJECT DESCRIPTION AND NUMERICAL CONSIDERATIONS The object consists of three nested rectangles, the longest being 24 X 32 pixels, centered on a 64 X 64 pixel field; it is shown in Fig. 1. The object is represented by a sequence, f(m,n),m = 1, - * * ,N, n = 1, - * * ,N(N=64)whichrepresents the gray levels of the mnth pixel. The gray levels are confined to 0 < f(m, n) < 1 and are 0.4, 0.8, and 1.0 in going from the largest to the smallest rectangle, respectively. The background level is zero. The two- dimensional (Cartesian) discrete Fourier transform (2-D DFT) of a sequence {f(m, n)} is given by
F(k,l) =
Afi + (1 - i)f2 0. Thus 0 .a2 6E. > V7- - 1 > 0. But in this region [CO/(I + c)] 2 is monotonUsing the standard methods of the calculus of variations, we ically increasing so that a minimum is reached at = Jk- 1. wish to minimize Hence for E' >E -
co
-
dxc-h [F-2] f[(f - h)2 +h2]
(A 1O)
where is a Lagrange multiplier. By differentiation, obtain the solution co
h=
f+
we
(Al 1)
1 +C
which, when inserted into (A7) gives
(A12)
h=f
=
0
(A13)
However, this solution must satisfy the constraint that
(A14)
h2dx<E. Thus, (A13) can only be a solution if
(ft+)2dX EL QE 1
(A 15)
Thus we have the result
h=P3f =fl
ifEE.
(A16)
What happens when E > E? Then the solution associated with co = 0 is not reachable and the function [w/(l + co)]2> 1 (Fig. 12). From (A8) and (Al 1) we obtain E+ (1
+
X)2 =E- a2.
Now write Et= kE,k > 1. Then
Summarizing all these results, we obtain
0o
fi
E
which was to be demonstrated. REFERENCES
giving
.
(A19)
ft
sto
min[1
The minimum clearly occurs at
h =P3f=
(A17)
[11 D. C. Youla and H. Webb, "Image restoration by the method of convex projections: Part 1-Theory," this issue, pp. 81-94. [2] R. W. Gerchberg, "Super-resolution through error reduction," Opt. Acta, vol. 21, pp. 709-720, 1974. [3] A. Papoulis, "A new algorithm in spectral analysis and bandlimited extrapolation," IEEE Trans. Circuits Syst., vol. CAS-22, pp. 735-742, 1975. [4] H. Stark, D. Cahana, and H. Webb, "Restoration of arbitrary fimite-energy optical objects from limited spatial and spectral information," J. Opt. Soc. Amer., vol. 71, no. 6, pp. 635-642,
1981. [5] J. A. Cadzow, "An extrapolation procedure for band-limited signal extrapolation: The extrapolation matrix," IEEE Trans. Circuits Syst., vol. CAS-25, pp. 74-78, 1978. [61 D. C. Youla, "Generalized image restoration by the method of alternating projections," IEEE Trans. Circuits Syst., vol. CAS-25, pp. 694-702, 1978. [7] H. Stark, J. W. Woods, I. Paul, and R. Hingorani, "Direct Fourier reconstruction in computer tomography," IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-29, pp. 237-245, 1981. [81 H. C. Andrews, Computer Techniques in Image Processing. New York: Academic, 1970. [9] M. I. Sezan and H. Stark, "An investigation of image restoration by projections onto convex sets," Image Processing Lab., Rensselaer Polytechnic Inst., Rep. IPL-TR-03 1. [10] L. G. Gubin, B. T. Polyak, and E. V. Raik, "The method of projections for finding the common point of convex sets," U.S.S.R. Comput. Math. and Math. Phys., vol. 7, no. 6, pp. 1-24, 1967. '[11] A. Lent and H. Tuy, "An iterative method for the extrapolation of band limited functions," J. Math. Anal. Appl., vol. 83, pp. 554-565, 1981.
95
Restoration by the Method of Convex Projections: Part 2-Applications
Image
and
Numerical Results M.
I.
SEZAN
AND
Abstract-The image restoration theory discussed in a previous paper by Youla and Webb [ 11 is applied to a simulated image and the results compared with the well-known method known as the GerchbergPapoulis algorithm. The results show that the method of image restoration by projection onto convex sets, by providing a convenient technique for utilizing a priori information, performs significandy better than the Gerchberg-Papoulis method.
I. INTRODUCTION -N a number of important practical situations it is required to restore a signal or image from incomplete information. The problem of restoring an image from its phase, the reconstruction of a tomographic image from incomplete view data, the restoration of a sharp image from its low-pass projection, and the extrapolation or interpolation of power spectra beyond the range of the measuring instrumentation are some significant examples. For this reason the Gerchberg-Papoulis (GP) algorithms [2], [3] and their offsprings [4], [5] have generated considerable interest. A generalization of these linear iterative algorithms was furnished by Youla [6] in which he also addressed the important questions of ill-posedness or noise-sensitivity of the restoration procedure. The restoration theory discussed by Youla and Webb in [1] has a very significant advantage over the GP and related algorithms in that it enables a large number of a priori known constraints to be incorporated in the algorithm through the mechanism of projection onto a convex set. Not all constraints can be handled this way. For example, the operation of digitizing a signal is not equivalent to projection onto a convex set. Nor, for that matter is the problem of reconstructing a signal from a prescribed magnitude. However, numerous other constraints can be treated this way and we shall illustrate this with examples in this paper. Eleven signal and image constraints and their associated projection operators are furnished in [11 and many of these are of physical significance. II. SETS AND PROJECTIONS
Because we shall deal with space-limited objects, defined by their reflectance, transmittance, or absorbance over a region Q, the functions of interest in this study will be assumed to be
H. STARK
members of L2x2(Q2), the space of all functions f(x,y) of twovariables x, y, square-integrable over Q. The associated Hilbert space JC is L2X2(Q) with inner product and norm defined by, respectively,
ffg(x, y) h*(x, y) dx dy
(g, h)
llgll A [(g,g)] /2.
(1) (2)
The following sets will be used. 1) C : The set of allf 's that vanish a.e. outside a prescribed region SC Q. C1 is an extension of the set 'T introduced in [1] and is a closed linear manifold (CLM) without interior. Given an arbitrary f E AC its projection onto C I is realized by Plf=
(x, y) EE 8
.4
05
(3)
(X,Y)O.
2) C2: The set of all f's in JC whose Fourier transforms assume a prescribed value G over a closed region £ in the u-v Fourier plane. This set is essentially the "dual," i.e., the frequency domain analog, of Ca (g)-the linear variety introduced in [1]. The projection of an arbitrary f E XC onto C2 is realized by
G(u, v), (u, V)E.C (F(u, v), (UIV)0
(4)
-
where F(u, v) = -f[f(x, y)] etc. for G(u, v), and Y is the Fourier transform operator. 3) C3: The set of all nonnegative functions in AC that satisfy the energy constraint
f If(x,y) 12 dXdyE_p2
(5)
where p = \/Eh. To show convexity we write
llfl +(1-A)f2 11ip()llflp=p
(6) +(1 - P) 2 11 0. Then
f
~
~~ \
In -f12 > AIInl2 + IIfI12 - 21(fn,f)I > ilfnII2 + IIfIh2 - 211AI1 lifll =
(Ilf 1- IIfl1)2 > e2.
But this leads to a contradiction since e
=0.
llfn - fli -
0. Thus
Now consider the projection onto C 3 of an arbitrary f E K. We show in the Appendix that this projection is realized by the following rule: (0, fA E
E
(8)
E,lZ
where f1 is the real part of f, fl is the rectified portion off1, andEt is the energy inf+, i.e.,
El _
Jf(f )2 dx dy.
(9)
4) C4: The subset of all functions in XC that are nonnegative. This set, already discussed in [1] is a closed convex set with vertex 0 (the zero function). The projection onto this set is
If,0,
>0(f=fl +jf2) (10) otherwise. 5) C5: The set of allf's in KC whose amplitudes must lie in a prescribed closed interval [a, b] a > 0, b > 0, a < b. To show convexity of e S note that
pf
fl
A1f, + (I - p) f2 > pxa + (1 - ,u) a = a
(11)
and
C"
(12)
any fi and f2 E es. To demonstrate closure we note that3 is open, hence Cs is closed. The projection onto Cs is
realized by the following rule:
21 6 is the set consisting of the elements in a that are not in 6. 3The complement of a set A is denoted by A c, -
f(x,y) < a
(a,
sf = if(x,y), tb,
a b.
We show in subsequent sections how the operators Pi and more generally the operators T1 _ 1 + X1(P1 - 1) i = 1, * *, 5 can be used, either sequentially or as compositions, to restore the missing data in the image. -
III. OBJECT DESCRIPTION AND NUMERICAL CONSIDERATIONS The object consists of three nested rectangles, the longest being 24 X 32 pixels, centered on a 64 X 64 pixel field; it is shown in Fig. 1. The object is represented by a sequence, f(m,n),m = 1, - * * ,N, n = 1, - * * ,N(N=64)whichrepresents the gray levels of the mnth pixel. The gray levels are confined to 0 < f(m, n) < 1 and are 0.4, 0.8, and 1.0 in going from the largest to the smallest rectangle, respectively. The background level is zero. The two- dimensional (Cartesian) discrete Fourier transform (2-D DFT) of a sequence {f(m, n)} is given by
F(k,l) =
Afi + (1 - i)f2 0. Thus 0 .a2 6E. > V7- - 1 > 0. But in this region [CO/(I + c)] 2 is monotonUsing the standard methods of the calculus of variations, we ically increasing so that a minimum is reached at = Jk- 1. wish to minimize Hence for E' >E -
co
-
dxc-h [F-2] f[(f - h)2 +h2]
(A 1O)
where is a Lagrange multiplier. By differentiation, obtain the solution co
h=
f+
we
(Al 1)
1 +C
which, when inserted into (A7) gives
(A12)
h=f
=
0
(A13)
However, this solution must satisfy the constraint that
(A14)
h2dx<E. Thus, (A13) can only be a solution if
(ft+)2dX EL QE 1
(A 15)
Thus we have the result
h=P3f =fl
ifEE.
(A16)
What happens when E > E? Then the solution associated with co = 0 is not reachable and the function [w/(l + co)]2> 1 (Fig. 12). From (A8) and (Al 1) we obtain E+ (1
+
X)2 =E- a2.
Now write Et= kE,k > 1. Then
Summarizing all these results, we obtain
0o
fi
E
which was to be demonstrated. REFERENCES
giving
.
(A19)
ft
sto
min[1
The minimum clearly occurs at
h =P3f=
(A17)
[11 D. C. Youla and H. Webb, "Image restoration by the method of convex projections: Part 1-Theory," this issue, pp. 81-94. [2] R. W. Gerchberg, "Super-resolution through error reduction," Opt. Acta, vol. 21, pp. 709-720, 1974. [3] A. Papoulis, "A new algorithm in spectral analysis and bandlimited extrapolation," IEEE Trans. Circuits Syst., vol. CAS-22, pp. 735-742, 1975. [4] H. Stark, D. Cahana, and H. Webb, "Restoration of arbitrary fimite-energy optical objects from limited spatial and spectral information," J. Opt. Soc. Amer., vol. 71, no. 6, pp. 635-642,
1981. [5] J. A. Cadzow, "An extrapolation procedure for band-limited signal extrapolation: The extrapolation matrix," IEEE Trans. Circuits Syst., vol. CAS-25, pp. 74-78, 1978. [61 D. C. Youla, "Generalized image restoration by the method of alternating projections," IEEE Trans. Circuits Syst., vol. CAS-25, pp. 694-702, 1978. [7] H. Stark, J. W. Woods, I. Paul, and R. Hingorani, "Direct Fourier reconstruction in computer tomography," IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-29, pp. 237-245, 1981. [81 H. C. Andrews, Computer Techniques in Image Processing. New York: Academic, 1970. [9] M. I. Sezan and H. Stark, "An investigation of image restoration by projections onto convex sets," Image Processing Lab., Rensselaer Polytechnic Inst., Rep. IPL-TR-03 1. [10] L. G. Gubin, B. T. Polyak, and E. V. Raik, "The method of projections for finding the common point of convex sets," U.S.S.R. Comput. Math. and Math. Phys., vol. 7, no. 6, pp. 1-24, 1967. '[11] A. Lent and H. Tuy, "An iterative method for the extrapolation of band limited functions," J. Math. Anal. Appl., vol. 83, pp. 554-565, 1981.
