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

332 COMPUTERIZED IMAGING II: IMAGE RECONSTRUCTION FROM PROJECTIONS

Property 5: Bounded signal

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

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,

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

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

336 COMPUTERIZED IMAGING II: IMAGE RECONSTRUCTION FROM PROJECTIONS

or

XN 1

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

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

The corresponding impulse response is given by