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

318				COMPUTERIZED IMAGING TECHNIQUES I: SYNTHETIC APERTURE RADAR

Figure 17.5. SAR imaging geometry [Courtesy of Soumekh, 1999].

SAR imaging is usually discussed in one of two modes. In the squint mode, the center of the target area is at ðxc; ycÞ where yc is nonzero. This is what is shown in Figure 17.5. In the spotlight mode, yc equals zero such that the target region appears perpendicular to the direction of flight.

Assume that the target region consists of stationary point reflectors with reflectivity n at coordinates ðxn; ynÞ, n ¼ 1; 2; 3; . . . . A pulse signal pð Þ is used to illuminate the target area. A radar receiver at ð0; yÞ receives the echo signal sðt; yÞ reflected back from the targets as

sðt; yÞ ¼	X	ð17:8-1Þ
	npðt tnÞ

n

A SIMPLIFIED THEORY OF SAR IMAGING				319
where tn is the round-trip delay from the radar to the nth target, given by
	tn ¼	2rn		ð17:8-2Þ
		c
and rn is the distance to the nth target, namely,
rn ¼ ½xn2 þ ðyn yÞ2&21				ð17:8-3Þ
The desired image to be reconstructed can be represented as
uðx; yÞ ¼	X			ð17:8-4Þ
	ndðx xn; y ynÞ

n

So the problem is how to obtain an approximation of uðx; yÞ from the measured signal sðt; yÞ. The 2-D Fourier transform of uðx; yÞ is given by

X	ð17:8-5Þ
Uðfx; fyÞ ¼ n exp½ j2pðfxxn þ fyynÞ&

n

The 1-D Fourier transform of sðt; yÞ with respect to the time variable t is given by

X	X	ð17:8-6Þ
Sðf ; yÞ ¼ Pðf Þ ¼ n exp½ j2pftn& ¼ Pðf Þ	n exp½ j2krn&
n	n

where k ¼ 2pf =c is the wave number, and Pðf Þ is the Fourier transform of pðtÞ. The 2-D Fourier transform of sðt; yÞ with respect to both t and y can now be

computed as the 1-D Fourier transform of Sðf ; yÞ with respect to y, and is given by

X							jðkxxn þ kyynÞ&		ð17:8-7Þ
Sðf ; fyÞ ¼ Pðf Þ	n	¼ n exp½
where
ky ¼ 2pfy					2 21				ð17:8-8Þ
kx ¼ ½4k			2
					ky & ¼ 2pfx
Note that Uðfx; fyÞ can be obtained from Sðf ; fyÞ by
U	f	; f			Sðf ; fyÞ				ð	17:8-9	Þ
	ð x			yÞ ¼ P		f	Þ
					ð

This operation is called source deconvolution where the source is Pðf Þ. However, using Eq. (17.8-9) usually yields erroneous results since it is an ill-conditioned

320	COMPUTERIZED IMAGING TECHNIQUES I: SYNTHETIC APERTURE RADAR

operation, especially in the presence of noise. Instead matched filtering is used in the form

U0ð fx; fyÞ ¼ Sðf ; fy				ÞP ðf Þ							ð17:8-10Þ
¼ j	P f	Þj	2	X	½	j2p	fxx	þ	fyy	Þ&
				n exp
	ð					ð

n

In the space-time domain, it can be shown that this operation corresponds to the

convolution of sðt; yÞ with p ð tÞ:
u0ðt; yÞ ¼ sðt; yÞ p ð tÞ
¼	X	nh t	2cn		ð	17	:	8-11	Þ
	n
			r
where hð Þ is the point spread function given by
hðtÞ ¼ F 1½jPðf Þj2&					ð17:8-12Þ

Equation (17.8-10) as compared to Eq. (17.8-4) shows that the spectrum of the reconstructed signal is different from the spectrum of the desired signal by the factor of jPðf Þj2. Since this is typically a slowly varying amplitude function of f with zero phase, reconstruction by matched filtering usually gives much better results than source deconvolution in the presence of noise.

17.9 IMAGE RECONSTRUCTION WITH FRESNEL APPROXIMATION

It is possible to simplify the reconstruction algorithm by using the Fresnel approximation. Suppose that the center of the target area is at ðxc; 0Þ. The Taylor series expansion of rn around ðxc; 0Þ with xn ¼ xc þ xn yields

rn ¼ xc þ xn þ	ðyn yÞ2	þ	ð17:9-1Þ
	2xc

Fresnel approximation corresponds to keeping the terms shown. Sðf ; yÞ of Eq. (17.8-6) can now be written as

Sðf ; yÞ ’ Pðf Þ expð j2kxcÞ	X	n exp"	j2kxn þ	ðyn yÞ2	#	ð17:9-2Þ
	n			xc
Consider the 1-D Fourier transform of the ideal target function						given by
Eq. (17.8-4) with respect to the x-coordinate:
X				ynÞ		ð17:9-3Þ
Uðkx; yÞ ¼	n expð jkxxnÞdðy

n

IMAGE RECONSTRUCTION WITH FRESNEL APPROXIMATION

321

where kx ¼2w2pfx and dðx xn; y ynÞ ¼ dðx xnÞdðy

ynÞ are

used.

Letting

k

2k

that S

f

is

x ¼

¼ c , comparison of Eqs. (17.9-2) and (17.9-3) shows

ð

; yÞ

2

approximately equal to the convolution of Pðf ÞUð2k; yÞ with expð jky

=xcÞ in the y-

coordinate:

Sðf ; yÞ ’ Pðf ÞUð2k; yÞ exp

jky2

ð17:9-4Þ

xc

In order to recover uðx; yÞ, matched filtering is performed with respect to both Pðf Þ and the chirp signal expð jky2=xcÞ. The reconstruction signal spectrum can be written as

U0ðfx; yÞ ¼ P ðf ÞSðf ; yÞ exp j	ky2		ð17:9-5Þ
	xc

In general, the radar bandwidth is much smaller than its carrier frequency kc so that jk kcj kc is usually true. This is referred to as narrow bandwidth. In addition, the y-values are typically much smaller than the target range. This is

referred	2to as narrow beamwidth. Then, expðjky2=xcÞ can							be approximated by
expðjkcy =xcÞ; and Eq. (17.9-5) becomes
	U0ðfx; yÞ ¼ P ðf ÞSðf ; yÞ exp					jkcy2		ð17:9-6Þ
						xc
where
	fx ¼	k	¼	w				ð17:9-7Þ
		p		pc

Note that there are two filtering operations above. The first one is with the filter transfer function P ðf Þ, which corresponds to p ð tÞ in the time domain. The second one is with the filter impulse function expðjkcy2=xcÞ.

We recall that the correspondence between the time variable t and the range variable x with the stated approximations is given by

x ¼	ct	ð17:9-8Þ
	2

In summary, when narrow bandwidth and narrow beamwidth assumptions are valid, image reconstruction is carried out as follows:

1.Filter with the impulse response p ð tÞ in the time direction.

2.Filter with the impulse response expðjkcy2=xcÞ in the y-direction.

3.Identify x as ct=2 or t as 2x=c.

322

COMPUTERIZED IMAGING TECHNIQUES I: SYNTHETIC APERTURE RADAR

Figure 17.6. Block diagram of SAR reconstruction with Fresnel approximation [Courtesy of Soumekh, 1999].

The block diagram for the algorithm is shown in Figure 17.6. When the approximations discussed above are not valid, more sophisticated image reconstruction algorithms are used, as discussed in the next section.

17.10ALGORITHMS FOR DIGITAL IMAGE RECONSTRUCTION

The goal of SAR imaging is to obtain the reflectivity function uðx; yÞ of the target area from the measured signal sðt; yÞ. There have been developed several digital reconstruction algorithms for this purpose, namely, the spatial frequency interpolation algorithm, the range stacking algorithm, the time domain correlation algorithm, and the backprojection algorithm [Soumekh]. The spatial frequency interpolation algorithm is briefly discussed below.

17.10.1Spatial Frequency Interpolation

The 2-D Fourier transform of sðt; yÞ yields Sðf ; fyÞ, which can be written as Sðk; kyÞ where k ¼ 2pf =c and ky ¼ 2pfy.

After matched filtering, it is necessary to map k to kx so that the reconstructed image can be obtained after 2-D inverse Fourier transform. The relationship between

kx, k, and ky is given by
kx2 ¼ 4k2 ky2	ð17:10-1Þ

When k and ky are sampled in equal intervals in order to use the DFT, kx is sampled in nonequal intervals. This is shown in Figure 17.7. Then, it is necessary to use interpolation in order to sample both kx and ky in a rectangle lattice. Interpolation

ALGORITHMS FOR DIGITAL IMAGE RECONSTRUCTION

323

Figure 17.7. SAR spatial frequency sampling for discrete data [Courtesy of Soumekh, 1999].

is not necessary only when kx ’ 2k, as discussed in the case of narrow bandwidth and narrow beamwidth approximation in Section 17.9.

Interpolation should be done with a lowpass signal. Note that the radar range swath is x 2 ½xc x0; xc þ x0&, y 2 ½yc y0; yc þ y0& where ðxc; ycÞ is the center of the target region. Since ðxc; ycÞ is not (0, 0), Sðkx; kyÞ is a bandpass signal with fast variations. In order to convert it to a lowpass signal, the following computation is performed:

S0ðkx; kyÞ ¼ Sðkx; kyÞe jðkxxcþkyycÞ

ð17:10-2Þ

Next S0ðkx; kyÞ is interpolated so that the sampled values of kx lie in regular intervals. The interpolated S0ðkx; kyÞ is inverse Fourier transformed in order to obtain the reconstructed image. The block diagram of the algorithm is shown in Figure 17.8.

324

COMPUTERIZED IMAGING TECHNIQUES I: SYNTHETIC APERTURE RADAR

Figure 17.8. Block diagram of SAR digital reconstruction algorithm with spatial frequency domain interpolation [Courtesy of Soumekh, 1999].

Interpolation is to be done with unevenly spaced data in the kx and k be sampled as

ky ¼ m ky ¼ kym

k ¼ n k ¼ kn

Then, the sampled values of kx are given by

k ¼ ½4n2ð kÞ2 m2ð k Þ2&1

xmn y

2

direction. Let ky

ð17:10-3Þ

ð17:10-4Þ

The function S0ðkx; kyÞ is to be interpolated for regularly sampled values of kx. This is done as follows [Soumekh, 1988]:

X

S0ðkx; kymÞ ¼ JmnS0ðkxmn; kymÞhðkx kxmnÞjkx kxmnj N kx ð17:10-5Þ n

where N is the number of sampled points along both positive and negative kx directions, kx is the desired sampling interval in the kx direction, hð Þ is the interpolating function to be discussed below, and Jm is the Jacobian given by

		d	½4k2 kym2	1		4kn		ð17:10-6Þ
	Jmn ¼			&2	¼
		dw				c½4kn2 kym2 &
	where k ¼ w=c. kx is chosen as
			kx ¼ i kx ¼ kxi					ð17:10-7Þ

ALGORITHMS FOR DIGITAL IMAGE RECONSTRUCTION

325

The choice of the function hð Þ is dictated by the sampling theorem in digital signal processing and communications. If a signal Vðkx; kyÞ is sampled as Vðn kx; m kyÞ it can be interpolated as

Vðkx; kyÞ ¼	X	Vðn kx; n kyÞ sin c			kx	n	ð17:10-8Þ
	n				kx
where
		sin cðvÞ ¼	sin pv				ð17:10-9Þ
			pv

The sin c function has infinitely long tails. In order to avoid this problem in practise, it is replaced by its truncated version by using a window function wðkxÞ in the form

kx	wðkxÞ	ð17:10-10Þ
hðkxÞ ¼ sin c kx

wðkxÞ can be chosen in a number of ways as discussed in Section 14.2. For example, the Hamming window is given by

	8			kx
		0:54 þ 0:46 cos		p		jkxj N kx
w kx				N kx				17:10-11
ð	:	0				otherwise	ð		Þ
	Þ ¼ <

18

Computerized Imaging II: Image

Reconstruction from Projections

18.1INTRODUCTION

In the second part of computerized imaging involving Fourier-related transforms, image reconstruction from projections including tomography is discussed. The fundamental transform for this purpose is the Radon transform. In this chapter, the Radon transform and its inverse are first described in detail, followed by imaging algorithms used in tomography and other related areas.

Computed tomography (CT) is mostly used as a medical imaging technique in which an area of the subject’s body that is not externally visible is investigated. A 3-D image of the object is obtained from a large series of 2-D x-ray measurements. An example of CT image is shown in Figure 18.1. CT is also used in other fields such as nondestructive materials testing.

This chapter consists of eight sections. The Radon transform is introduced in Section 18.2. The projection slice theorem that shows how the 1-D Fourier transforms of projections are equivalent to the slices in the 2-D Fourier transform plane of the image is discussed in Section 18.3. The inverse Radon transform (IRT) is covered in Section 18.4. The properties of the Radon transform are described in Section 18.5.

The remaining sections are on the reconstruction algorithms. Section 18.6 illustrates the sampling issues involved for the reconstruction of a 2-D signal from its projections. Section 18.7 covers the Fourier reconstruction algorithm, and Section 18.8 describes the filtered backprojection algorithm.

18.2THE RADON TRANSFORM

Consider Figure 18.2. The x y axis are rotated by y to give the axis u v. The relationship between ðx; yÞ and ðu; vÞ is given by a plane rotation of

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

Copyright # 2007 John Wiley & Sons, Inc.

326

THE RADON TRANSFORM

327

Figure 18.1. The CT image of the head showing cerebellum, temporal lobe and sinuses [Courtesy of Wikipedia].

y degrees:

y	¼	sin y cos y	v	ð		Þ
x		cos y sin y	u		18:2-1

The Radon transform pðu; yÞ of a signal gðx; yÞ shown in Figure 18.1 is the line integral of gðx; yÞ parallel to the v-axis at the distance u on the u-axis which makes the angle 0 y < p with the x-axis:

1ð

pðu; yÞ ¼

gðx; yÞdv

1

1ð

ð18:2-2Þ

¼gðu cos y v sin y; u sin y þ v cos yÞdv

1

Figure 18.2. Rotation of Cartesian coordinates for line integration along the v-axis.