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

238 APODIZATION, SUPERRESOLUTION, AND RECOVERY OF MISSING INFORMATION

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

^u ¼ Hþv

ð14:15-4Þ

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

Hþ ¼ ðHtHÞ 1Ht

ð14:15-5Þ

Hþ is of size N2 N1. Hþ satisfies

HþH ¼ I

ð14:15-6Þ

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

HHþ ¼ I

ð14:15-7Þ

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

Hþ ¼ HtðHHtÞ 1

ð14:15-8Þ

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

H ¼	1	cmfmt	ð14:16-1Þ
	lm2

m¼1

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

COMPUTATION OF Hþ BY SINGULAR VALUE DECOMPOSITION (SVD)

239

The SVD representation of the pseudoinverse Hþ is given by

X

r

Hþ ¼	1	fmcmt	ð14:16-2Þ
	lm2

m¼1

Using Eq. (14.16-2) in Eq. (14.15-4), the pseudoinverse solution can be written as

r

^u ¼		1		ð14:16-3Þ
	lm2fmcmt v
m¼1
X
which is the same as
	r
^u ¼	m¼1	omfm		ð14:16-4Þ
	X
where
om ¼		cmt v		ð14:16-5Þ
		lm1=2

Equation (14.16-4) shows that ^u is in the vector space spanned by the eigenvectors of HtH.

The use of Eqs. (14.16-4) and (14.16-5) is relatively simple for small problems. However, for large scale problems, it becomes very difficult to do so. For N1 ¼ N2 ¼ 256, H is of size 65;536 65;536. Then, it becomes prohibitive.

If the degradation point spread function (PSF) is separable so that v ¼ H1uH2,

the generalized inverse is also separable, and
	^u		¼		HþvHþ				ð	14:16-6	Þ
					1	2
EXAMPLE 14.7 Determine the least squares solution of
2	2	1	3		x1	¼	2	2	3
	1	2		2				1
4			5		x		4		5
	1	3						3
by determining the pseudoinverse.
Solution: H and A ¼ HtH are given by
	1		2				6		7
	4				5
H ¼ 2 2 1					3A ¼ 7 14

1 3

240 APODIZATION, SUPERRESOLUTION, AND RECOVERY OF MISSING INFORMATION

singular, Hþ is given by		l1 ¼	10	þ			l2	¼	10
The eigenvalues of A are				p65 and							p65. Since A is non-
Hþ ¼ A 1Ht ¼ 35		7 6			2	1	3	¼ 35			5	8 11
1		14	7		1	2	1		1		0	21	7

The least-squares solution is given by

" #

0:6

^u ¼ Hþy ¼

0:628

EXAMPLE 14.8 Suppose that the solution is constrained in some sense in the signal domain and the constraint is applied as a matrix acting on ^u. For example, bounding the image to be nonzero in a certain region would be such an operator. Determine the least squares solution in this case.

Solution: The error measure to be minimized can be written as

E ¼ jv HTuj2

where T is the constraint operator in matrix form. In the case of boundedness, T is a diagonal matrix whose diagonal elements are 1s and 0s. We have

E ¼ ðv HTuÞtðv HTuÞ

¼ vtv 2utTtHtv þ utTtHtHTu

The gradient of E is given by

g ¼ 2TtHtv þ 2TtHtHTu

Setting g equal to 0 yields

^u ¼ ½TtHtHT& 1TtHtv

provided that the matrix inverse exists.

14.17THE STEEPEST DESCENT ALGORITHM

When the goal is to obtain ^u rather than Hþ, the iterative gradient method can be used. This method is based on the fact that the steepest descent of E ¼ jv Huj2 is in the direction of the negative gradient of E with respect to u.

THE STEEPEST DESCENT ALGORITHM

241

Suppose that uk is the vector u at the kth iteration. The gradient of E with respect to uk is given by

gk ¼ 2Htðv HukÞ	ð14:17-1Þ
The iterative gradient algorithm is given by
ukþ1 ¼ uk akgk with u0 ¼ 0	ð14:17-2Þ

where ak is a scalar. Using Eq. (14.17-2) in Eq. (14.17-1), gk can be written as

gk ¼ gk 1 ak 1HtHgk 1			ð14:17-3Þ
It can be shown that uk converges to ^u if
2
0 ak			ð14:17-4Þ
	lmax

where lmax is the maximum eigenvalue of A ¼ HtH.

If ak is chosen constant, the method is called the one-step gradient algorithm. The optimum value of a for fastest convergence is given by

2
aopt ¼ lmax þ lmin	ð14:17-5Þ

where lmin is the minimum eigenvalue of A.

If A is highly ill-conditioned, the condition number of A equal to lmax=lmin is large, and convergence may be very slow.

If ak is optimized at each iteration, the method is called the method of steepest descent. The optimal value of ak is given by

gkt gk
ak ¼ gkt Agk	ð14:17-6Þ

When lmax=lmin is large, convergence may still be slow when ak is used.

EXAMPLE 14.9 Determine the approximate least squares solution of Example 14.8 by (a) the one-step gradient algorithm in 12 iterations, (b) the steepest descent algorithm in three iterations.

Solution: (a) In the one-step gradient algorithm, sopt is found as

2

aopt ¼ l1 þ l2 ¼ 0:1

242 APODIZATION, SUPERRESOLUTION, AND RECOVERY OF MISSING INFORMATION

In the first iteration, we have
	g0 ¼		8		0:8
			13 ;	u1 ¼ 1:3
and the error E1 ¼ jv	Hu1j2 equals 9.46.
After 12 such iterations, we get
	u11 ¼	0:685		u12 ¼	0:555
		1:253 ;			0:581
with the error E12 ¼ jv	Hu12j2 ¼ 1:101.

(b) With the steepest descent algorithm, ak is optimized at each iteration according to Eq. (14.17-6). The result after three iterations is

u3 ¼ 0:5592 0:629

where the error E3 ¼ jv Hu3j2 ¼ 1:0285.

14.18THE CONJUGATE GRADIENT METHOD

In the conjugate gradient method, the scalar ak as well as the correction vector is optimized at each iteration. This allows faster convergence. For this purpose, Eq. (14.17-2) is replaced by

ukþ1 ¼ uk þ akrk	ð14:18-1Þ
where the correction vector rk replaces the gradient vector gk.
The correction vector rk has the following properties:
gkt þ1rk ¼ 0	ð14:18-2Þ
½Hrk&t½Hrm& ¼ rkt Arm ¼ 0 k ¼6 m; 0 k; m < N	ð14:18-3Þ

where N is the rank of A.

Due to these properties, the algorithm converges in N iterations. The iterative equations in addition to Eq. (14.18-1) are as follows:

rk ¼ gk þ bkrk 1 r0 ¼ g0					ð14:18-4Þ
	gt	Ark 1
bk ¼	k				ð14:18-5Þ
	rt	1	Ark 1
	k

THE CONJUGATE GRADIENT METHOD				243
ak ¼	gkt rk			ð14:18-6Þ
	rkt Ark
gk ¼ Htv þ Auk ¼ gk 1 þ ak 1Ark 1				ð14:18-7Þ
The final solution ^u can also be written as
		X
		N	1	ð14:18-9Þ
		^u ¼	amrm

m¼0

When N is a large number, the algorithm can be stopped earlier than N iterations. If the rank of H is M < N, it is satisfactory to run the algorithm with M iterations.

EXAMPLE 14.10 Determine the conjugate gradient method for the least squares problem of Example 14.8.

Solution: The least-squares solution satisfies

TtHtHTu ¼ TtHtv

This is in the same form as before with the replacement of H by H0 ¼ HT. The only iterative equation which needs to be modified is the equation for the gradient, which becomes

gk ¼ gk 1 ak 1TtHtHrk 1

The initial image u0 is chosen as satisfying the constraint represented by T:

Tu0 ¼ u0

g0 rm and r0 are given by

g0 ¼ r0 ¼ TtHtðv Hu0Þ

15

Diffractive Optics I

15.1INTRODUCTION

Diffractive optics involves creation of holograms by computation in a digital computer, followed by actual implementation by a recording system such as a laser film writer, scanning electron beam writer, high resolution printer, usually followed by photoreduction, and the like. Diffractive optics is very versatile since any type of wave can be considered for computation within the computer. Digital holograms created with such technology are more commonly called diffractive optical element (DOE). Another name for DOE is computer-generated hologram (CGH). Because of these equivalent terminologies, the words hologram, DOE, and CGH will be used equivalently in this chapter.

In communication engineering, it is well known that a complex signal must be coded in some fashion before it can be transmitted over a single channel. This can be done by modulating the amplitude of a carrier by the modulus of the complex signal and/or by modulating the phase of the carrier by the phase of the complex signal. The real part of this modulated carrier then can be transmitted, and through proper demodulation procedures, the complex signal can be recovered at the receiver. Both holography and diffractive optics are based on such modulation and demodulation principles.

DOEs have found applications in many areas such as wave shaping, laser machining, 3-D displays [Yatagai, 1975], optical pattern recognition, optical interconnects, security devices [Dittman, et al., 2001], optical communications/networking, and spectroscopy. Some important advantages of DOEs are as follows:

A DOE can perform more than one function, for example, have multiple focal points, corresponding to multiple lenses on a single element. It can also be designed for use with multiple wavelengths.

DOEs are generally much lighter and occupy less volume than refractive and reflective optical elements.

With mass manufacturing, they can be manufactured less expensively for a given task.

Diffraction, Fourier Optics and Imaging, by Okan K. Ersoy

Copyright # 2007 John Wiley & Sons, Inc.

244

INTRODUCTION

245

Since DOEs are based on diffraction, they are highly dispersive, namely, their properties depend on wavelength. Hence they are more commonly used with monochromatic light. However, the dispersive properties can also be used to advantage. For example, light at one wavelength can be focused at one point, and light at another wavelength can be focused at another point in space. In devices such as optical phased arrays used in optical communications and networking, dispersion is used to separate different wavelengths at different focal points. In addition, refractive and diffractive optical elements can be combined in such a way that wavelength dispersion can be removed, and goals such as removal of spherical aberration can be achieved.

There are several steps in the creation of a DOE. These are sampling of waves, computation of wave propagation, usually with the FFT, coding of complex wave information on the hologram, for example, as a real and positive image on the hologram plane, and recording of the resulting DOE.

In 3-D wave coding, phase is usually much more important to be represented correctly than amplitude. Hence, phase is given much more attention. In applications in which the output image intensity is used, output amplitude variations can be significantly reduced by multiplying the output image with a random phase factor (diffuser). This does not change the intensity of the image [Burckhardt, 1970].

This chapter consists of ten sections. The first few methods discussed utilizes the Fourier transformation property of a lens and generate a DOE design for a Fourier transform hologram. The first such method is the Lohmann method discussed in Section 15.2. In this method, amplitude modulation is achieved by controlling the sizes of the apertures, and phase modulation is achieved by controlling the positions of the apertures generated. Section 15.3 discusses the approximations involved in implementing DOEs based on the Lohmann method. These approximations are compensated for in a method called the Lohmann-ODIFIIT method in Section 16.8. Section 15.4 simplifies the Lohmann method by choosing constant aperture size, implying constant amplitude. In order to compensate for the error generated, iterative optimization discussed in Chapter 14 is utilized. Section 15.5 describes computer experiments with the Lohmann methods. Another Fourier method based on hard-clipping of the hologram transmittance function resulting in a binary hologram is discussed in Section 15.6.

The methods discussed so far assume that the reconstructed image is on an output plane. The method discussed in Section 15.7 is capable of generating output image points at arbitrary locations in 3-D. It is a simple method to use. It was the method used in demonstrating for the first time how to make DOEs with a scanning electron microscope system. Section 15.8 extends the method of Section 15.7 and makes it capable of generating many object points in 3-D space.

Essentially all DOE methods involve nonlinear coding of amplitude and phase. As a result, they generate harmonic images. When the FFT is used, the implicit assumption is that the object is periodic. This also causes the images to repeat in a periodic fashion on the reconstruction plane. The one-image-only holography method discussed in Section 15.9 generates a single image by using semi-irregular sampling and a spherical reference wave whose origin is located near the hologram

246		DIFFRACTIVE OPTICS I
DOE plane	Lens	Output plane

Figure 15.1. The Fourier transform system used in DOE design.

plane. This method is later utilized in Chapter 19 for the design of phasars for dense wavelength division multiplexing in order to achieve many more wavelength channels in optical communications and networking than what is possible with the present phasar methods.

The final section goes back in time and introduces a classical DOE that functions as a flat lens. It is called a binary Fresnel zone plate.

15.2LOHMANN METHOD

The Lohmann method as well as a number of other methods for DOEs are developed for the Fourier transform arrangement shown in Figure 15.1.

In the Lohmann method (also called the detour phase method), the hologram plane is divided into smaller rectangles each containing an aperture [Lohmann, 1970]. An example of a Lohmann cell is shown in Figure 15.2. The size of the aperture is used to control amplitude, and its position is changed to adjust phase. This results in a binary transmission pattern.

Figure 15.2. The (n,m)th cell of a Lohmann hologram.

APPROXIMATIONS IN THE LOHMANN METHOD

247

Let hðx; yÞ be the output amplitude from the hologram Hðnx; nyÞ. It should be proportional to the desired image f ðx; yÞ. The binary transmission function of the hologram can be written as

H

;

nyÞ ¼

X X

rect

nx ðn þ PnmÞdn

rect

ny mdn

ð

15:2-1

Þ

ðnx

n m

cdn

Wnmdn

When a tilted plane wave, expð2pix0nxÞ is incident on the binary hologram, the complex amplitude after the hologram is Hðnx; nyÞ expð2pix0nxÞ. The complex

amplitude at the image plane is the Fourier transform of Hðnx; nyÞ expð2pix0nxÞ:

ðð

Hðnx; nyÞe2pi½ðxþx0Þnxþyny&dnxdny

X X

¼ cðdnÞ2sinc½cdnðx þ x0Þ& WnmsincðyWnmdnÞ ð15:2-2Þ

n m

expf2pi½dnððx þ x0Þðn þ PnmÞ þ ymÞ&g

The parameters Wnm, Pnm, and the two constants x0,c are chosen such that the complex amplitude in the image plane matches the desired image f ðx; yÞ.

Equation (15.2.2) can be compared to the desired image by writing f ðx; yÞ in the form

f ðx; yÞ ¼ ðð Fðnx; nyÞe2piðxnxþynyÞdnxdny ¼	n m	Fðndn; mdnÞe2pi½dnðxnþymÞ&
	X X	ð15:2-3Þ

The two sinc functions and the factor exp½2piðxPnmdnÞ& in Eq. (15.2-2) can be assumed to be close to unity. The validity of this assumption will be discussed in Section 15.3. Equating the Fourier coefficients in Eqs. (15.2-2) and (15.2-3) yields

cðdnÞ2Wnm expf2pi½x0dnðn þ PnmÞ&g / Fðndn; mdnÞ;							fnm=2p		; ð15:2-4Þ
F	ndn; mdn	Þ /	c	dn	Þ	2Anm exp 2pi		Þ&
ð			ð			½ ð
Wnm Anm; Pnm þ n fnm=2px0dn:
By choosing x0dn equal to an integer M, we obtain
	Pnm fnm=2pM:								ð15:2-5Þ

Equations (15.2-4) and (15.2-5) show that the height W and the position P of the aperture in each cell can be chosen proportional to the amplitude A and the phase f of the complex amplitude F at the cell, respectively.

15.3APPROXIMATIONS IN THE LOHMANN METHOD

In the previous section, the three approximations (a) sinc½cdnðx þ x0Þ& const, (b) sincðyWnmdnÞ 1, and (c) exp½2piðxPnmdxÞ& 1 were made for simplicity. The