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

328 COMPUTERIZED IMAGING II: IMAGE RECONSTRUCTION FROM PROJECTIONS

pðu; yÞ is also called the ray-sum or the ray integral of gðx; yÞ at the angle y. A set of ray integrals forms a projection. Below y will be implicitly treated by writing pðu; yÞ as pðuÞ.

A projection taken along a set of parallel rays as shown in Figure 18.2 is called a parallel projection. If the set of rays emanates from a point source, the resulting projection is called a fan-beam projection. Only the case of parallel projections is considered in this chapter.

The Radon transform can also be expressed in terms of the Dirac-delta function. u can be written as atr where a ¼ ½cos y sin y&t, r ¼ ½x y&t. Then, Eq. (18.2-1 ) is the

same as

ð

pðu; yÞ ¼ gðrÞdðu atrÞdr

ðð	ð		Þ
1		18:2-3

¼gðx; yÞdðu cos y x sin y yÞdxdy

1

The Radon transform has many applications in areas such as computerized tomography, geophysical and multidimensional signal processing. The most common problem is to reconstruct gðx; yÞ from its projections at a finite number of the values of y. The reconstruction is exact if projections for all y are known. This is due to the theorem given in the next section.

18.3THE PROJECTION SLICE THEOREM

The 1-D FT of pðuÞ denoted by Pðf Þ is equal to the central slice at angle y of the 2-D FT of gðx; yÞ denoted by Gðfx; fyÞ:

Pð f Þ ¼ Gð f cos y; f sin yÞ

ð18:3-1Þ

The slice of frequency components in the 2-D frequency plane is visualized in Figure 18.3.

Figure 18.3. The slice of frequencies used in the projection slice theorem.

THE PROJECTION SLICE THEOREM

329

Proof:

1ð

Pðf Þ ¼ pðuÞe j2pftudu

1

ðð1

¼gðu cos y v sin y; u sin y þ v cos yÞe j2pfudvdu

1

In the unrotated coordinate system, this is the same as

ðð1

Pð f Þ ¼

gðx; yÞe j2pð fx cos yþfy sin yÞdxdy

1

¼ Gðf cos y; f sin yÞ

EXAMPLE 18.1 Find the Radon transform of gðx; yÞ ¼ e x2 y2 :

Solution:

ðð1

pðu; yÞ ¼ e x2 y2 dðu x cos y y sin yÞdxdy

1

Since the transformation ðx; yÞ ! ðu; vÞ as given by Eq. (18.2-1) is orthonormal, x2 þ y2 ¼ u2 þ v2. Hence,

1ð

p	u;	yÞ ¼	e u2 e v2 dv
ð

1

¼ ppe u2

EXAMPLE 18.2 Find the Radon transform of gðx; yÞ if

g

x; y

ð

1

x2

y2Þl 1

x2 þ y2 < 1

ð

Þ ¼

0

otherwise

Solution: gðx; dyÞ is nonzero inside the unit circle. Since x2

þ y2

¼ u2 þ v2, the

p

2

limits of integration with respect to u are

p1 u2. Thus,

p u; y

1

u2

v2 l 1dv

1 u

ð

Þ ¼

ð

½

&

p

1 u2

330 COMPUTERIZED IMAGING II: IMAGE RECONSTRUCTION FROM PROJECTIONS

The following definite integral can be found in [Erdelyi]:

	ð ða2 t2		Þl 1dt ¼						1ð Þ
	a						a2l 1pp l
	a							l þ
									2
Utilizing this result above yields
p u; y	8	ð	1Þ			ð1 u2Þl 1=2			1 < u < 1
	>	pp	l
ð Þ ¼	<
	>	l þ	2
	>
	>
	:

0otherwise

18.4THE INVERSE RADON TRANSFORM

The 2-D FT representation of gðx; yÞ is given by

	gðx; yÞ ¼			ðð Gðfx; fyÞej2pðfxxþfyyÞdfxdfy	ð18:4-1Þ
				1
				1
Converting ð fx; fyÞ to polar coordinates ð f ; yÞ gives
p	1
gðx; yÞ ¼ ð	ð	Gðf cos y; f sin yÞej2pf ðx cos yþy sin yÞjf jdf dy			ð18:4-2Þ
0	1
Since Pðf Þ ¼ Gðf cos y;		f sin yÞ, this is the same as
		p	1	fPðf Þsgnðf Þej2pf ðx cos yþy sin yÞdf dy
gðx; yÞ ¼ ð			ð		ð18:4-3Þ
		0	1

By convolution theorem, with u ¼ x cos y þ y sin y, the first integral is given by

1

e j2pfudf

1

@

1

fP

f

sgn

f

p u

ð

Þ

Þ

¼

2jp @u

Þ

jpu

ð

ð

ð

1

1

ð18:4-4Þ

¼ 2p2

ð

@@ðttÞ u t dt

1

p

1

1

PROPERTIES OF THE RADON TRANSFORM

331

The last result is the Hilbert transform of 21p @@u pðuÞ. Let

		1		@ @ðttÞ u t dt			ð18:4-5Þ
^pðuÞ ¼ 2p2 ð
		1		p	1
		1
^pðuÞ is called the filtered projection. The inverse Radon transform becomes
	p
gðx; yÞ ¼ ð0		^pðuÞdy					ð18:4-6Þ
	p
¼ ð0		^pðx cos y þ y sin yÞdy

This is equivalent to the back-projection of ^pð This means the projections for all y are needed for not band-limited.

uÞ along the angular direction y. perfect reconstruction if gðx; yÞ is

18.5PROPERTIES OF THE RADON TRANSFORM

The Radon transform has a number of properties that are very useful in applications. The most important ones are summarized below. When needed, pðuÞ will be written as pgðuÞ to show that the corresponding signal is gðx; yÞ.

Property 1: Linearity

If hðx; yÞ ¼ g1ðx; yÞ þ g2ðx; yÞ;
phðuÞ ¼ pg1 ðuÞ þ pg2 ðuÞ		ð18:5-1Þ
Property 2: Periodicity
pðu; yÞ ¼ pðu; y þ 2pkÞ; k an integer		ð18:5-2Þ
Property 3: Mass conservation	ð pðuÞdu
ðð gðx; yÞdxdy ¼		ð18:5-3Þ
1	1
1	1
Property 4: Symmetry
pðu; yÞ ¼ pð u; y pÞ		ð18:5-4Þ

332 COMPUTERIZED IMAGING II: IMAGE RECONSTRUCTION FROM PROJECTIONS

Property 5: Bounded signal

If gðx; yÞ ¼ 0 for jxj, jyj> D2 ,
		D	ð18:5-5Þ
pðuÞ ¼ 0 for juj> p2
Property 6: Shift
If hðx; yÞ ¼ xðx x0; y y0Þ,
phðu; yÞ ¼ pgðu u0; yÞ			ð18:5-6Þ
where u0 ¼ x0cosy þ y0siny.
Property 7: Rotation
gðx; yÞ in polar coordinates is gðr; fÞ. If hðr; fÞ ¼ gðr; f þ f0Þ,
phðu; yÞ ¼ pgðu; y þ f0Þ			ð18:5-6Þ
Property 8: Scaling
If hðx; yÞ ¼ gðax; ayÞ, a ¼6 0,
1
phðuÞ ¼		pgðauÞ	ð18:5-7Þ
	jaj

Property 9: Convolution

If hðx; yÞ ¼ g1ðx; yÞ g2ðx; yÞ (where represents 2-D convolution),

phðuÞ ¼ pg1 ðuÞ pg2 ðuÞ

1ð

¼pg1 ðtÞpg2 ðu tÞdt

1 ð18:5-8Þ

1ð

¼pg1 ðu tÞpg2 ðtÞdt

1

This property is useful to implement 2-D linear filters by 1-D linear filters.

18.6RECONSTRUCTION OF A SIGNAL FROM ITS PROJECTIONS

The signal is given by its inverse Radon transform. For computer reconstruction, it is necessary to discretize both y and u. A major consideration is how many samples are needed.

Assume that projections are available at angles y0; y1; y2 . . . yN 1. In order to estimate N, the number of projections, the signal is assumed to be both space-limited

THE FOURIER RECONSTRUCTION METHOD

333

Figure 18.4. Samples in the Fourier domain when all projections are sampled at the same sampling rate.

and band-limited (even though this is not theoretically possible, it is a good approximation in practice). Thus, gðx; yÞ will be assumed to be band-limited such

that Gðfx; fyÞ ¼ 0 for fx2 þ fy2 F02. By the 2-D sampling theorem, gðx; yÞ can be sampled as gð xm; ynÞ without loss of information if x 1=2F0 and

y 1=2F0.

gðx; yÞ will also be assumed to be space-limited such that gðx; yÞ ¼ 0 for x2 þ y2 T02. By the 2-D sampling theorem with the domains reversed, Gðfx; fyÞ can be sampled as Gð fxm; fynÞ without loss of information if fx 1=2T0,fy 1=2T0.

Consider Figure 18.4 in which a polar raster of samples in the Fourier domain is shown [Dudgeon]. On the outermost circle with N projections, the distance between the samples can be estimated as pF0=N. This distance is assumed to satisfy the same criterion as fx and fy. Thus,

pF0		1
							ð18:6-1Þ
	N			2T0
With the sampling interval in		the signal					domain approximated as
¼ x ¼ y 1=2F0, the number of projections N should satisfy
	N		pT0				ð18:6-2Þ

Equation (18.6-2) should be considered as a rule of thumb.

18.7THE FOURIER RECONSTRUCTION METHOD

We will assume that there are N projections, each projection is sampled, and the number of samples per projection is M. A DFT of each sampled projection is

334 COMPUTERIZED IMAGING II: IMAGE RECONSTRUCTION FROM PROJECTIONS

Figure 18.5. The sampling points on a polar lattice to be interpolated.

computed. These DFT values are interpreted as samples of the FT of gðx; yÞ on a regular polar lattice, as shown in Figure 18.5 [Dudgeon]. These samples are interpolated to yield estimates on a rectangular lattice. Finally, the 2-D inverse DFT is computed to obtain gðm; nÞ. Usually, the size of the inverse DFT is chosen as two to three times that of each dimension of gðm; nÞ.

Appropriate windowing as discussed in Section 14.2 before the inverse DFT allows the minimization of the effects of Fourier domain truncation.

The simplest approach to interpolation is zeroth order or linear interpolation. Each desired rectangular sample point ðfx; fyÞ is surrounded by four polar samples. In zeroth order interpolation, the rectangular sample is assigned the value of the nearest polar sample. In linear interpolation, a weighted average of the four nearest polar samples is used. The weighting can be taken inversely proportional to the Euclidian distance between the points.

A simple method is to use bilinear interpolation in terms of the polar coordinates: let ðf ; yÞ be the coordinates of the point at which the interpolated sample is desired; let ðfi; yjÞ, i; j ¼ 1; 2 denote the coordinates of the four nearest polar samples. Then, the bilinear interpolation can be written as

X X

Gðf ; yÞ ¼

2

2

Gðfi; yjÞh1ðf fiÞh2ðy yjÞ

ð18:7-1Þ

i¼1 j¼1

where

jf 0j

h1

f 0

Þ ¼

1

j

f 0

j

f

f

ð

ð18:7-2Þ

jy0j

h

0

Þ ¼

1

0

j

y

y

2

ðy

jy

and f , y are the sampling intervals in polar coordinates.

THE FILTERED-BACKPROJECTION ALGORITHM

335

Since the density of the polar samples is less dense at high frequencies than at low frequencies, the interpolation results are also less accurate at high frequencies. This causes some degradation in the reconstructed image.

18.8THE FILTERED-BACKPROJECTION ALGORITHM

The inverse Radon transform can be written as

p
gðx; yÞ ¼ ð0	^pðx cos y þ y sin yÞdy		ð18:8-1Þ
where
	1
^pðuÞ ¼	ð	jf jPðf Þej2pfudf	ð18:8-2Þ
	1

Pðf Þ is zero outside ½ F0; F0&. Its samples can be approximated at intervals of

f ¼ 2F0=N. pðuÞ is

also

sampled at

intervals of

1=2F0. By

choosing

f ¼ 1=N, Pð fkÞ can be approximated by a DFT of size N as follows:

X

N=2

N

N

P k

P

fk

Þ ¼

p

n

e j2pnk=N ;

2

< k

2

ð

18:8-3

Þ

ð Þ ¼

ð

ð

Þ

n¼N=2þ1

or

X

1

N 1

PðkÞ ¼

pðnÞe j2pnk=N ;

k ¼ 0; 1 . . . ; ðN 1Þ

ð18:8-4Þ

2F0

n¼0

since pð j‘jÞ and Pð jmjÞ can be chosen equal to pðN j‘jÞ and PðN jmjÞ, respectively.

Equation (18.8-2 ) can be written as
^pðuÞ ¼	Fð0	jf jPðf Þej2pfudf	ð18:8-5Þ
	F0

^pðuÞ can be approximated by an inverse DFT by using the same sampling interval as follows:

^pðnÞ ¼ ^pðn Þ ¼ f	N=2	ð18:8-6Þ
	X Pð fkÞjk f jej2pnk=N

k¼N=2þ1

336 COMPUTERIZED IMAGING II: IMAGE RECONSTRUCTION FROM PROJECTIONS

or

XN 1

^pðnÞ ¼ f		P0ðkÞej2pnk=N ;		n ¼ 0; 1; . . . ; ðN 1Þ					ð18:8-7Þ
	k¼0
where
						N
		<
P0	k	8	PðkÞjk f j		0 k 2				ð	18:8-8	Þ
	ð	Þ ¼ > P k N k f			N < k < N
		:		Þ j
		> ð Þjð			2

Better reconstructions are obtained if PðkÞ is multiplied by a windowing function as discussed in Section 14.2. The windowing function deemphasizes the high frequencies at which measurement noise is more serious.

Finally, the reconstructed image is obtained by the numerical integration of Eq. (18.8-1). If there are M projections, gðx; yÞ can be approximated as g0ðx; yÞ given by

	p	X
g0ðx; yÞ ¼	M	^pðx cos yk þ y sin ykÞ	ð18:8-9Þ
		k¼1

u ¼ ðx cos yk þ y sin ykÞ may not correspond to pð nÞ of Eq. (18.8-4). Then, interpolation is necessary again to find pðuÞ given pðn Þ, n ¼ 0; 1; 2 . . . ; ðN 1Þ. Linear interpolation is usually sufficient for this purpose.

The algorithm as described above has certain drawbacks. The use of the DFT with N points means that linear convolution of Eq. (18.4-5) is replaced by a ‘‘partial’’ circular convolution. The word ‘‘partial’’ is used since the discretization of jf j is exact. If the size N DFTs are replaced by size 2N DFTs by zero-padding, the reconstruction improves since circular convolution becomes linear in the desired output terms.

However, best results are achieved if the convolution expressed in the transform domain by Eq. (18.4-5) is first written in the time-domain, and then the whole processing is done by DFTs of zero-padded sequences.

The transfer function in Eq. (18.4-5 ) is given by

H f	f	f F0	ð	18:8-10	Þ
ð	Þ ¼ j0j	otherwisej j

The corresponding impulse response is given by

	hðtÞ ¼	Fð0	Hðf Þej2pftdf				ð18:8-11Þ
		F0
			sin 2ptF0			sin ptF0	2
	¼ 2F02			F02
			2pxF0			pxF0

THE FILTERED-BACKPROJECTION ALGORITHM									337
The discretized impulse response at t ¼ n t is given by
			8	0			n even
			>	F02			n ¼ 0
					4F	0
ð		Þ ¼	<					ð		Þ
h	n t		>						18:8-12
			>
			>
			:
			>		n2p2		n odd
Equation (18.4-5) as a convolution in the time-domain is given by
			1
^pðuÞ ¼			ð	pðtÞhðu tÞdt				ð18:8-13Þ
		1
Its discrete version is given by		X
^pðn Þ ’			1	pð mÞhððn mÞ Þ				ð18:8-14Þ
		m¼ 1
Since pð tmÞ is zero for jmj > N2,
		¼X
		N=2		pð tmÞhððn mÞ tÞ				ð18:8-15Þ
^pðn tÞ ’ t
		m	N=2

This linear convolution can now be computed with DFTs of size 2N after zeropadding both pð Þ and hð Þ. Further improvement is possible by properly windowing the frequency domain results. The whole procedure can be written as follows:

1. Zero-pad pð Þ and hð Þ in the form

8

pðnÞ

0 n

N

2

>

>

>

N

3N

<

ð

Þ ¼

>

ð

Þ

>

2

2

p n

>

0

< n

<

18:8-16

>

n > 3N

>

>

>

2N

> p n

:

> ð

Þ

2

and similarly for hð Þ.

2.Compute the size 2N DFTs of pð Þ and hð Þ.

3.Do the transform domain operations.

4.Window the results of step 3 by a proper window.

5.Compute the size 2N inverse DFT.