Radon Transform and Fourier Slice Theorem

With respect to reconstruction, it is important to note that the projec- tions under all the angles ϑ can be regarded as another 2-D represen- tation of the image. One coordinate is the position in the projection profi le, r, the other the angle ϑ (Fig. 8.15). Consequently, we can re- gard the parallel projection as a transformation of the image into an- other 2-D representation. Reconstruction then just means applying the inverse transformation. The critical issue, therefore, is to describe the

226                                                                                                         8 3-D Imaging

 

Figure 8.15: Geometry of a projection beam.

 

tomographic transform mathematically and to investigate whether the inverse transform exists.

A projection beam is characterized by the angle ϑ and the off set r (Fig. 8.15). The angle ϑ is the angle between the projection plane and the x axis. Furthermore, we assume that we slice the 3-Dobject parallel to the xy plane. Then, the scalar product between a vector x on the projection beam and a unit vector

n ¯ = [cos ϑ, sin ϑ ]^T                                                (8.23)

normal to the projection beam is constant and equal to the off set r of the beam

x n ¯ − r = x cos ϑ + y sin ϑ − r = 0.                                      (8.24)

The projected intensity P (r, ϑ ) is given by integration along the projec- tion beam:

 

P (r, ϑ ) = ∫

g( x )d s = ∫ ∫ g( x )δ (x₁ cos ϑ + x₂ sin ϑ − r )d²x. (8.25)

∞

∞

path

− ∞ − ∞

The projective transformation of a 2-D function g( x ) onto P (r, ϑ ) is named after the mathematician Radon as the Radon transform.

To better understand the properties of the Radon transform, we an- alyze it in the Fourier space. The Radon transform can be understood as a special case of a linear shift-invariant fi lter operation, the projec- tion operator. All gray values along the projection beam are added up.

8.6 Depth from Multiple Projections: Tomography                          227

 

±

Therefore the point spread function of the projection operator is a δ line in the direction of the projection beam. In the Fourier domain this convolution operation corresponds to a multiplication with the transfer function, which is a δ line (2-D) or δ plane (3-D) normal to the δ line in the spatial domain (see R5). In this way, the projection operator slices a line or plane out of the spectrum that is perpendicular to the projection beam.

This elementary relation can be computed most easily, without loss of generality, in a rotated coordinate system in which the projection di- rection coincides with the y' axis. Then the r coordinate in P (r, ϑ ) co- incides with the x' coordinate and ϑ becomes zero. In this special case, the Radon transform reduces to an integration along the y' direction:

∞

P (x', 0) =

∫ g(x', y')dy'.                         (8.26)

− ∞

The Fourier transform of the projection function can be written as

∞

Pˆ (k_x', 0) =

∫ P (x', 0) exp(− 2π ik_x' x')dx'.                             (8.27)

− ∞

Replacing P (x', 0) by the defi nition of the Radon transform, Eq. (8.26) yields

 

Pˆ (k

∞   ∞

, 0)   ∫  ∫ g(x', y')dy'  exp(

 

2π ik

 

x')dx'.            (8.28)

x' =

− ∞

 − ∞

   −    _x'

−           =

If we insert the factor exp( 2π i0y') 1 in this double integral, we recognize that the integral is a 2-D Fourier transform of g(x', y') for k_y' = 0:

 

∞

∞

Pˆ (k_x', 0)  =    ∫  ∫  g(x', y') exp(− 2π ik_x' x') exp(− 2π i0y')dx'dy'

− ∞ − ∞

= g(k_x', 0).

 

 

(8.29)

Back transformation into the original coordinate system fi nally yields

Pˆ (q, ϑ ) = gˆ ( k )δ ( k − ( k n ¯ ) n ¯ ),                                        (8.30)

where q is the coordinate in the k space in the direction of ϑ and n ¯ the normal vector introduced in Eq. (8.23). The spectrum of the projection is identical to the spectrum of the original object on a beam normal to the direction of the projection beam. This important result is called the Fourier slice theorem or projection theorem.

228                                                                                                         8 3-D Imaging

 

Filtered Back-Projection

| |

| |

If the projections from all directions are available, the slices of the spec- trum obtained cover the complete spectrum of the object. Inverse Fourier transform then yields the original object. Filtered back-projection uses this approach with a slight modifi cation. If we just added the spectra of the individual projection beams to obtain the complete spectrum of the object, the spectral density for small wave numbers would be too high as the beams are closer to each other for small radii. Thus, we must correct the spectrum with a suitable weighting factor. In the continuous case, the geometry is very easy. The density of the projection beams goes with k − 1. Consequently, the spectra of the projection beams must be multiplied by k . Thus, fi ltered back-projection is a two-step process. First, the individual projections must be fi lled before the reconstruction can be performed by summing up the back-projections.

| |

| |

In the fi rst step, we thus multiply the spectrum of each projection direction by a suitable weighting function wˆ ( k ). Of course, this oper- ation can also be performed as a convolution with the inverse Fourier transform of wˆ ( k ), w(r). Because of this step, the procedure is called the fi ltered back-projection.

In the second step, the back-projection is performed and each pro- jection gives a slice of the spectrum. Adding up all the fi ltered spectra yields the complete spectrum. As the Fourier transform is a linear op- eration, we can add up the fi ltered projections in the space domain. In the space domain, each fi ltered projection contains the part of the ob- ject that is constant in the direction of the projection beam. Thus, we can back-project the corresponding gray value of the fi ltered projection along the direction of the projection beam and add it up to the contri- butions from the other projection beams.

After this illustrative description of the principle of the fi ltered back- projection algorithm we derive the method for the continuous case. We start with the Fourier transform of the object and write the inverse Fourier transformation in polar coordinates (q, ϑ ) in order to make use of the Fourier slice theorem

∫ ∫

2π ∞

g( x ) =     qgˆ (q, ϑ ) exp[iq(x₁ cos ϑ + x₂ sin ϑ )]dqdθ.                               (8.31)

0 0

 

In this formula, the spectrum is already multiplied by the wave number,

− ∞   ∞

q. The integration boundaries, however, are not yet correct to be applied to the Fourier slice theorem (Eq. (8.30)). The coordinate, q, should run from to and ϑ only from 0 to π. In Eq. (8.31), we integrate only over half a beam from the origin to infi nity. We can compose a full beam from two half beams at the angles ϑ and ϑ +π. Thus, we split the integral

8.6 Depth from Multiple Projections: Tomography                          229

in Eq. (8.31) into two over the angle ranges [0, π [ and [π, 2π [ and obtain

π ∞

g( x )  =  ∫   ∫ qgˆ (q, ϑ ) exp[iq(x₁ cos ϑ + x₂ sin ϑ )]dqdϑ

0 0

∞

π

+  ∫   ∫ qgˆ (− q, ϑ ') exp[− iq(x₁ cos ϑ ' + x₂ sin ϑ ')]dqdϑ '

0 0

using the following identities:

ϑ ' = ϑ +π, gˆ (− r, ϑ ) = gˆ (r, ϑ '), cos(ϑ ') = − cos(ϑ ), sin(ϑ ') = − sin(ϑ ).

Now we can recompose the two integrals again, if we substitute q by

−

q in the second integral and replace gˆ (q, ϑ ) by Pˆ (q, ϑ ) because of the Fourier slice theorem Eq. (8.30):

∞

π

g( x ) = ∫ ∫ |q|Pˆ (q, ϑ ) exp[iq(x₁ cos ϑ + x₂ sin ϑ )]dqdϑ.                                    (8.32)

0− ∞

Equation (8.32) gives the inverse Radon transform and is the basis for the fi ltered back-projection algorithm. The inner integral performs the back-projection of a single projection:

P' = F− 1(|q|FP ).                                               (8.33)

| |

F

denotes the 1-DFourier transform operator. P' is the projection func- tion P multiplied in the Fourier space by q . If we perform this operation as a convolution in the space domain, we can formally write

P' = [F− 1(|q|)] ∗ P.                                (8.34)

The outer integral in Eq. (8.32) over the angle ϑ,

∫

π

g( x ) = P'(r, ϑ ) dϑ,                                               (8.35)

0

 

sums up the back-projected and fi ltered projections over all directions and thus forms the reconstructed image. Note that the fi ltered projec- tion profi le P'(r, ϑ ) in Eq. (8.35) must be regarded as a 2-D function to built up a 2-D object g( x ). This means that the projection profi le is projected back into the projection direction.

230                                                                                                         8 3-D Imaging

 

a                                                                         b

c

Figure 8.16: Illustration of the fi ltered back-projection algorithm with a point ob- ject: a projections from diff erent directions; b fi ltering of the projection functions; c back-projection: adding up the fi ltered projections.

 

⇐ Предыдущая37383940414243444546Следующая ⇒

Поделиться: