Ersoy O.K. Diffraction, Fourier optics, and imaging (Wiley, 2006)(ISBN 0471238163)(427s) PEo

228 APODIZATION, SUPERRESOLUTION, AND RECOVERY OF MISSING INFORMATION

Thus, ui is also the fixed point of Ti. It can also be shown that Ti is contractive if 0 < li < 2. lis are called the relaxation parameters. Note that li ¼ 1 makes Ti equal to Pi. Thus, Ti is a more general operator than Pi. Proper choice of li usually speeds up convergence.

Let T be defined by

T ¼ TmTm 1 . . . T1					ð14:8-4Þ
The following can be shown to be true [Stark]:
1. The iterations defined by
u0 arbitrary
ukþ1 ¼ Tuk
converge weakly to a fixed point in S0. Weak convergence means klim ðuk; wÞ ¼
				!1
ðu ; wÞ, 8w 2 S0 where ða; bÞ is the inner product of the vectors a and b. Under
k		j	un	u	j ¼
certain conditions, the convergence is strong, meaning lim						0.
	!1

2.li can be changed at each iteration. Let lik denote the value of li at iteration k. For every choice of lik, provided that 0 < lik < 2, the iterations defined by

u0 arbitrary

ukþ1 ¼ Tuk

converge weakly to a fixed point u .

Quite often, the POCS algorithm is used when M ¼ 2. This is especially true when one projection is in the signal domain and the other one in the transform domain. This case is shown in Figure 14.7.

uk		Transform	UK		Transform
					domain
					constraints

uK+1	Space		u ′	Inverse
	domain		k
				transform
	constraints

Figure 14.7. The POCS algorithm in the transform and signal domains.

OTHER POCS ALGORITHMS

229

There are two possible ways to estimate the optimal values of l1k and l2k for fast convergence. In the first case, l1k is first optimized, and T1 is applied. This is followed by the optimization of l2k and the application of T2. In general, it can be shown that for a linear projector lik ¼ 1, and lik 0 otherwise. In the second case, l1k and l2k can be optimized simultaneously so that T ¼ T2T1 results in minimum error [Stark].

14.9GERCHBERG–PAPOULIS (GP) ALGORITHM

The original GP algorithm is a special case of the POCS algorithm, and uses projectors P1 and P2 in the signal and the Fourier domains, respectively [Gerchberg]. Hence

ukþ1 ¼ P2P1uk

ð14:9-1Þ

where P1 and P2 are defined as follows:

Let A be a region in R2, and S1 be the set of all functions in S that vanish almost everywhere (except, possibly, over a point set of measure zero) outside A. P1 is

defined by
P1u ¼ (	u	ðx; yÞ 2 A	ð14:9-2Þ
	0	otherwise

Let S2 be the set of all functions whose Fourier transforms assume a prescribed function Vðfx; fyÞ over a closed region B in the Fourier plane. P2 is defined by

Uðfx; fy			Þ ð	fx; fy			2 B	ð14:9-3Þ
P2u $ ( V	ð	fx; fy			fx; fyÞ		=B
			Þ ð			Þ2

It can be shown that both S1 and S2 are convex sets.

14.10OTHER POCS ALGORITHMS

In the GP algorithm, the two projection operators P1 and P2 are with respect to the convex sets S1 and S2 discussed above, respectively. Now, we introduce three more projection operators utilizing convex sets.

Consider the set S3 introduced in Example 14.3. Assume u is complex. It can be written as

u ¼ uR þ juI

ð14:10-1Þ

230 APODIZATION, SUPERRESOLUTION, AND RECOVERY OF MISSING INFORMATION

where uR; uI denote the real and complex parts of u, respectively. xþR is also defined as

				u	uR > 0			ð14:10-2Þ
		uRþ		¼ 0R	otherwise
The projection operator P3 is defined by
P3u			0		uR < 0		E		14:10-3
	¼	8 uRþ			uR > 0; ERþ			ð		Þ
		>
		> p
		:		E=ERþuRþ	uR > 0; ERþ > E
		<
where ERþ is given by
		ERþ ¼ ðð½uRþðx; yÞ&2dxdy						ð14:10-4Þ

and E is the estimated or known energy of the desired image. ERþ is bounded. Consider the set S4 defined in Example 14.4. The projection operator P4 is

defined by

>

a

u x y

< a

P4u

¼

ð

x; y

Þ

ð

Þ

b

ð

14:10-5

Þ

8 u

að

; uÞx; y

:

ð

Þ

>

x; y

b

< b

u

Finally, consider the subset S5 of all real nonnegative functions in S. S5 can be

shown to be convex. The projection operator P5 is defined by
P	5	u	¼	uR	uR 0	ð	14	:	10-6	Þ
				0	otherwise

where u ¼ uR þ juI as before.

Some examples of the applications of the theory presented above is given in the following sections.

14.11RESTORATION FROM PHASE

For the sake of simplicity, the 1-D case will be discussed below. The signal uðxÞ will be assumed to have a compact support ½ a; a&. The known phase of the FT of uðxÞ is fðf Þ. The amplitude of the signal is to be estimated. It can be shown that the necessary condition for uðxÞ to be completely specified by fðf Þ within a scale factor

RESTORATION FROM PHASE	231
is that there is no point of symmetry x0 such that
uðx0 þ xÞ ¼ uðx0 xÞ	ð14:11-1Þ

The two projections of interest can be written as follows:

		u x	uðxÞj a
P1u ¼ ( ð0 Þ			jotherwise				ð14:11-2Þ
		8 Aðf Þ cos½fðf Þ cðf Þ&ejfðf Þ			p
				jfðf Þ	cðf Þj <
P2u						2	ð	14:11-3	Þ
		:	0	otherwise
	$ <
where A(f) is the FT amplitude, and cðf Þ				is the FT phase. It is observed that
Aðf Þ cosðfðf Þ		cðf ÞÞ is the projection of Uðf Þ in the direction given by fðf Þ. In a

number of papers, the relative directions given by fðf Þ and cðf Þ are not considered, and P2 is expressed as

P2u $ Aðf Þejfðf Þ

ð14:11-4Þ

It is possible that fðf Þ is not known for all f. Suppose that fðf Þ is known for f in a set S. Then P2 can be expressed as

	8 Aðf Þ cos½fðf Þ cðf Þ&ejfðf Þ		jfðf Þ cðf Þj		<	p		and f 2 S
P2u						2				14:11-5
	>			2
	<
	$ >	0	otherwise; f		S				ð		Þ
	>
	:	Uðf Þ	f2=S
	>

where Uðf Þ is the FT of u(t). It is observed that P1 is a linear operator. Then, the optimal l1k is 1. A lower bound for l2k can be derived as [Stark ]

	l2k lL ¼	1	þ	jP2uk P1P2ukj2	ð	14:11-6	Þ
				jP1P2uk ukj2

lL can be used as an estimate of l2k. However, lL may be greater than 2. In later iterations, it is necessary to restrict l2k to be less than 2 to guarantee convergence. For example, minðlL; 1:99Þ can be chosen.

Since l1k equals 1, one cycle of iterations can be written as

ukþ1	¼ P1½1 þ l2kðP2	1Þ&uk	ð14:11-7Þ
	¼ ð1 l2kÞuk þ l2kP1P2uk
where P1 uk is replaced by uk	since uk ¼ P1T2uk	1 and therefore uk 2 D1.

232 APODIZATION, SUPERRESOLUTION, AND RECOVERY OF MISSING INFORMATION

14.12 RECONSTRUCTION FROM A DISCRETIZED PHASE FUNCTION BY USING THE DFT

The discrete signal of length N will be denoted by u½n&. The phase of the DFT of u½n& as well as the support of u½n& are assumed to be known. The support of u½n& will be defined as follows:

u½n& ¼ 0 outside 0 n < M and u½0& ¼6 0

It can be shown that N 2M for the reconstruction algorithm to be satisfactory [Hayes et al., 1980]. The DFT of u½n& is given by

U½k& ¼ A½k&ejj½k&

ð14:12-1Þ

j½k& is assumed to be known.

The procedure of the algorithm is as follows:

1. Initialize A½k& with an initial guess, A1½k&. Then, the initial U1½k& is given by

U1½k& ¼ A1½k&ejj½k&

ð14:12-2Þ

2.Compute the inverse DFT of U1½k&, yielding the first estimate u1½n&. Let i equal 1.

3.Define a new sequence vi½n& as

vi½n& ¼	u n	&	0		n < M	ð14:12-3Þ
	i0½		M		n < N

4. Compute the DFT Vi½k& of vi½n&. Vi½k& can be written as

Vi½k& ¼ Ai½k&ejci½k&

ð14:12-4Þ

The rest of the procedure can be written in two ways with respect to two different algorithms as follows, respectively.

Algorithm 1 [Levi and Stark]

5. Project Vi½k& as follows:

(

Ai k

cos

j

k

ci k

&Þ

ejf½k&

if j½k& ci k

&j

<

p

U

k

2

ð

14:12-5

Þ

iþ1

½ & ¼

½ &

ð

½

0&

½

otherwise ½

This operation corresponds to P2. This equation can be modified so that T2 is used by

Uiþ1½k& ¼ ½1 þ l2iðP2 1Þ&½Vi½k&&

ð14:12-6Þ

l2i is computed by using the relation (14.11-6).

RECONSTRUCTION FROM A DISCRETIZED PHASE FUNCTION BY USING THE DFT 233

6.Generate a new estimate of uiþ1½n& by computing the inverse DFT of

Uiþ1½k&.

7.Repeat steps 2–6 until convergence.

Algorithm 2 [Hayes et al.]

The only different step is the following: 8. Project Vi½k& as follows:

Uiþ1½k& ¼ Ai½k&ejf½k&

ð14:11-7Þ

All the other steps are the same as in Algorithm 1.

EXAMPLE 14.5 Let the error Ei at the ði þ 1Þth iteration be defined as

X
N 1	u½n&	uiþ1	½n& 2
Eiþ1 ¼ n¼0

where u½n& is the desired signal to be reconstructed. Prove that Ei is nonincreasing as i increases.

Solution: By Parseval’s theorem, the following is true:

		X
Ei ¼	1 N 1		U½k&	Ui½k&	2
	N k¼0

Since Z(k) and Zi(k) have the same phase, the last equation can be written as

		X
Eiþ1 ¼	1	N 1			2
	N k¼0 U½k&			Uiþ1½k&

We also define

				X
Ei0	¼		1 N 1		U½k&	Vi½k&	2
		N k¼0

The triangle equation for vector differences can be used to note the following:

Eiþ1 Ei0

with equality if Uiþ1½k& ¼ Vi½k&. The inequality above is due to the fact that Uiþ1½k& has the same phase as U½k& whereas Vi½k& has some other phase.

234 APODIZATION, SUPERRESOLUTION, AND RECOVERY OF MISSING INFORMATION

Again using the Parseval’s theorem on this inequality results in

X
N	1
Ei	ju½n& vi½n&j2

n¼0

Since vi½n& is before the projection to create uiþ1½n&, the above equation shows that the projection reduces the error.

14.13GENERALIZED PROJECTIONS

There are many signal recovery problems involving nonconvex constraints. For example, digitizing a signal and constraining a signal or its Fourier transform to have a specified magnitude are nonconvex projections.

When Si is a nonconvex set, there may be more than one point satisfying the definition of a projector Pi. This discrepancy can be removed by specifying additional constraint(s). Another problem with nonconvex sets is that a proof for the existence of a projection has not been found.

In the discussion below, the number of projections M will be limited to 2. In this case, the following error measure is useful:

EðukÞ ¼ jP1uk ukj þ jP2uk ukj

ð14:13-1Þ

The following theorem determines the convergence properties of generalized projections with m ¼ 2 [Levi–Stark]:

Theorem: The recursion given by

ukþ1 ¼ T1T2uk;				u0 arbitrary				ð14:13-2Þ
has the property
			Eðukþ1Þ EðukÞ					ð14:13-3Þ
for every l1 and l2 which satisfy
0	li		Ai2 þ Ai					ð	14:13-4	Þ
			2	1		ðAi
			Ai þ Ai				þ BiÞ
					2
	A		jP1T2uk		T2ukj			ð	14:13-5	Þ
		1	¼ jP2T2uk		T2ukj
	A		jP2uk	ukj				ð	14:13-6	Þ
		2	¼ jP1uk	ukj

RESTORATION FROM MAGNITUDE							235
B1	¼	ðP2T2uk	T2uk; P1T2uk			T2ukÞ	ð14:13-7Þ
			2
		jP2T2uk		T2ukj
B2	¼	ðP1uk	uk; P2uk	ukÞ			ð14:13-8Þ
			2
		jP1uk ukj

It can be shown that the upper limit in Eq. (14.13-4) includes 1. Hence, the theorem applies with P1; P2 instead of T1; T2.

The projections P1 and P2 can involve multiple constraints. Quite often, one set of constraints is in the signal domain, and the other set is in the transform domain.

EXAMPLE 14.6 A projector P maps a function uðxÞ such that the resulting function gðxÞ satisfies

gðxÞ ¼

where sðxÞ for all x. Specify P.

Solution:

P Satisfies

sðxÞ x 2 D D0 0 x 2= D0

j

Pu

uj ¼

v jv

uj

inf

Since gðxÞ ¼ PuðxÞ, we have

jPu uj ¼ jg uj ¼

ð

jgðxÞ

uðxÞj2dx þ

ð

c

jgðxÞ

uðxÞj2dx

x2=D0

ð

x2D0\D

ð

þ

jgðxÞ uðxÞj2dx þ

jgðxÞ uðxÞj2dx

x2D;sðxÞ>0

x2D;sðxÞ<0

jg uj is minimum if gðxÞ is chosen as

¼

>

0

x=D0

x2

>

P1u

8 u

x

D0

Dc

> s

ðx Þ

x 2

D; x\t > 0

>

ð Þ

2

ð Þ

<

>

:sðxÞ x 2 D; xðtÞ 0

14.14RESTORATION FROM MAGNITUDE>

This problem is also called phase retrieval. It occurs in a number of fields such as astronomy, optics, and x-ray crystallography, where only the intensity equal to the

236 APODIZATION, SUPERRESOLUTION, AND RECOVERY OF MISSING INFORMATION

square of the amplitude is available at each measurement point, and the phase is to be determined.

Phase retrieval is a much more difficult problem than restoration from phase. In the 2-D case, all the functions uðx; yÞ, uðx; yÞ, uðx x0; y y0Þ, and uð x; yÞ have the same Fourier transform magnitude. Apart from these four cases, it can be shown that a timeor space-limited sequence can be uniquely determined from its Fourier magnitude if the z-transform of the sequence is irreducible (cannot be written in terms of first-degree polynomials in z or z 1) [Hayes et al.]. This result is not interesting in the 1-D case since all polynomials can be factored in terms of first-degree polynomials according to the fundamental theorem of algebra. However, in the 2-D case, most polynomials are irreducible. Consequently, restoration of images from Fourier magnitude information is feasible.

Let the known Fourier transform magnitude be jUðf Þj. The two projections can be written as

P1u ¼ (	u	u < a
	0	otherwisej j	ð14:14-5Þ
P2u $ jUðf Þjejfðf Þ			ð14:14-6Þ

where fðf Þ is the currently available Fourier phase. We observe that P1 is convex whereas P2 is not. P1 may also involve other constraints.

The general restoration algorithm can be written as

ukþ1 þ T1T2uk; u0 arbitrary	ð14:14-7Þ
where Ti ¼ 1 þ liðPi 1Þ, and T1 can be chosen equal to P1	since P1 is linear.

When both li are chosen equal to 1, the algorithm is known as Gerchberg–Saxton algorithm as discussed in Section 14.9.

When uk ¼ P1T2uk 1, uk belongs to the space D1 whose projector is P1. Hence,

P1uk ¼ uk, and EðukÞ of Eq. (14.13-1) reduces to
EðukÞ ¼ jP2ukþ1 ukþ1j	ð14:14-8Þ
Equation (14.13-2) also reduces to
ukþ1 ¼ ð1 l2kÞuk þ l2kP1P2uk	ð14:14-9Þ

Computing the Fourier transform of this equation and denoting FTs with the operator Fð Þ yields

Ukþ1 ¼ ð1 l2kÞUk þ l2kFðP1P2ukÞ

ð14:14-10Þ

IMAGE RECOVERY BY LEAST SQUARES AND THE GENERALIZED INVERSE

237

Due to Parseval’s relation, Eq. (14.14-8) can be written as

	1
Eðukþ1Þ ¼	ð	½Uðf Þ Ukþ1&2df	ð14:14-11Þ
	1

l2k can be estimated by a simple search in the range 0 < l2k < a (where a is typically a number such as 3) so that Eðukþ1Þ is minimized.

14.14.1Traps and Tunnels

When at least one of the projections is nonconvex, the convergence may occur to a point which is not a fixed point of every individual Ti. Such a point um is called a trap. xt fails to satisfy one or more of the a priori constraints but satisfies

um ¼ T1T2um

ð14:14-12Þ

Traps cannot occur with convex sets. On the other hand, tunnels may occur with box convex and nonconvex sets. A tunnel occurs when uk changes very little from iteration to iteration.

A trap can be detected by observing no change in EðukÞ > 0 from iteration to iteration. If P1 is linear and P1T2uk ¼ uk, it can be shown that the correct solution u lies in a hyperplane orthogonal to P2uk uk [Stark].

14.15 IMAGE RECOVERY BY LEAST SQUARES AND THE GENERALIZED INVERSE

The techniques discussed to this point involved iterative optimization. Another powerful approach for image recovery involves the method of least squares. For this approach, the measured image is modeled as

v ¼ Hu þ n

ð14:15-1Þ

where H is a matrix of size N1 N2, u is the desired image, and n is the noise image. In Eq. (14.15-1), the images are row or column-ordered.

The unconstrained least squares estimate of u is obtained by minimizing

E ¼ jv Huj2 ¼ ðv HuÞt ðv HxÞ

ð14:15-2Þ

where j j2 denotes the L2 norm.

Computing the partial derivative of E with respect to u and setting it equal to

zero results in
Htv ¼ HtHu	ð14:15-3Þ