Файл: Ersoy O.K. Diffraction, Fourier optics, and imaging (Wiley, 2006)(ISBN 0471238163)(427s) PEo .pdf

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

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

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,

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

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

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

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

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

where j j2 denotes the L2 norm.

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

238 APODIZATION, SUPERRESOLUTION, AND RECOVERY OF MISSING INFORMATION

If Ht½H of size N2 N2 is nonsingular, the least-square solution is given by

where Hþ is called the generalized inverse (pseudo inverse) of H, and is given by

Hþ is of size N2 N1. Hþ satisfies

However, it is not true that HHþ equals I. A necessary condition for HtH to be nonsingular is that N1 N2 and the rank r of H is N2.

If N1 < N2, and the rank of H is N1, Hþ is defined such that

Hþ satisfying Eq. (14.15-7) is not unique. Uniqueness is achieved by constraining the solution given by Eq. (14.15-4) to have minimum norm. In other words, among all possible solutions, the one with minimum j^xj2 is chosen. Then, the pseudo inverse is given by

In conclusion, Hþ always exists uniquely as discussed above, and the least squares solution is Hþv.

Hþ has the following additional properties:

1.HHþ ¼ ðHHþÞt In other words, HHþ is symmetric.

2.HþH ¼ ðHþHÞt In other words, HþH is symmetric.

3.HHþH ¼ H

4.HþHHt ¼ Ht

14.16 COMPUTATION OF Hþ BY SINGULAR VALUE DECOMPOSITION (SVD)

The SVD representation of the matrix H of size N1 N2 and rank r can be written as

X

r

m¼1

where cm and fm are the eigenvectors of HHt and HtH, respectively, with the singular values lm, 1 m r.