*1.1.2. 2D and 3D image reconstruction*

*Fig. D.4. Density of the measured values in the frequency domain.*

domain (Fig. D.4) is given by:

*parallel.*

described in the next section.

*For the S values of each projection in spatial domain correspond S measured values in the frequency space (measured* for each line). Therefore, the distance ε between two consecutive values measured on a radial line in the frequency

94 Imaging and Radioanalytical Techniques in Interdisciplinary Research - Fundamentals and Cutting Edge Applications

A sufficient condition to obtain a good reconstruction is to ensure that the worst azimuth resolution (s) in the frequency

DS <sup>⇒</sup>P≈<sup>S</sup> <sup>π</sup>

determined by their dimensions or pixel size. The exploration beam is not perfectly parallel, too; this is another source of errors on the measured data obtained especially when using a reconstruction method which assumes that the beam is

In practice, insufficient number of projection generates artefacts such as those shown in figure (D.5) [3]. The projections of figure (D.5) have a dimension of 64x64 pixels. So the number of rows S is 64. As it has already been demonstrated, a number P of projections around a value of 100 (P ≈ Sπ / 2) is sufficient to produce an acceptable reconstructed image. For the case of this example, it is clear that the reconstruction image obtained from 64 projections is the closest one to the original image. The 2D reconstruction is performed using the Filtered Back-Projection method (FBP) that will be

Thus, the ratio between the number of projections (P) and the number of rows (S) must be in the order of π / 2. An insufficient number of projections may produce a 2D or 3D reconstructed image which present many undesirable artefacts. In practice, most of the tomography detectors cannot measure below the nominal Nyquist resolution

DS (D.21)

<sup>2</sup> (D.22)

ε= 2ωmax <sup>S</sup> <sup>=</sup> <sup>1</sup>

domain (Fig. D.4) is in the same order of radial resolution (ε). This condition can be expressed as follows:

1 2D π <sup>P</sup> <sup>≈</sup> <sup>1</sup> As already mentioned, the inverse transformation (simple inversion) of the Fourier transform Pθ(ω) can be directly used to produce the 2D reconstructed layer of μ (x, y). However, a more efficient way and more elegant has been developed called " Filtered Back-Projection (FBP)" making the reconstruction process less expensive and less complicated when compared to the direct calculation of the inverse Fourier transform [1,2]. The basic idea of this method is derived from the tomography principle and the scanning mode: the object (set of layers) is scanned projection by projection from 0 ° to 180 ° which means that Pθ(t) = Pθ+180(-t) and suggesting the use of a polar coordinates rather than Cartesian square ones. As it was demonstrated the digitization and sampling details by the right selection of P and S are very important in Computed tomography. If the projections are recorded in polar coordinates, the projections Fourier transforms will be discrete values of a polar function. Thus, it seems appropriate to write the 2D Fourier transform of the object function μ (x, y) in polar coordinates. The inverse Fourier transform μ (x, y) of F(u, v) written in Cartesian coordinates is given by:

$$\text{Im}\,\text{(x,y)} = \int\_{-\infty}^{+\infty} \int\_{-\infty}^{+\infty} \text{F}(\mathbf{u}, \mathbf{v}) . \text{e}^{\#\mathbf{2}(\mathbf{u}\mathbf{x} + \mathbf{v}\mathbf{y})} \quad \text{dudv} \tag{13}$$

The same inverse transform in polar coordinates is given by:

$$\mathfrak{u}(\mathbf{x}, \mathbf{y}) = \int\_0^{2\pi} \int\_{-\infty}^{+\infty} \mathbf{F}(\omega, \theta) . \mathbf{e}^{\mathfrak{A}\omega(\mathbf{x}\cos\theta + \mathbf{y}\sin\theta)} \quad \text{od}\,\mathfrak{o}\mathbf{d}\,\Theta \tag{14}$$

The substitution of x cos (θ) + y sin (θ) by t and the application of Fourier Central Slice theorem, allow us to get the following expression:

$$\begin{aligned} \mathsf{h}(\mathsf{x}, \mathsf{y}) &= \Big[\_{0}^{\mathsf{T}} \Big[ \mathsf{F}(\mathsf{a}, \mathsf{θ}) . \mathsf{e}^{\star 2\pi \mathsf{i} \omega \mathsf{t}} \mid \omega \mid \mathsf{d}\omega \Big] \mathsf{d}\theta \\\ &= \Big[\_{0}^{\mathsf{T}} \Big[ \int\_{-\infty}^{+\infty} \mathsf{F}(\mathsf{a}, \mathsf{θ}) . \mathsf{e}^{\star 2\pi \mathsf{i} \omega \mathsf{t}} \mid \omega \mid \mathsf{d}\omega \Big] \mathsf{d}\theta \end{aligned} \tag{15}$$

The integral in brackets can be regarded as the inverse Fourier transform Pθ(t) of Pθ(ω). However, we note that it is multiplied by the function |ω| which plays the role of a special ramp filter function in the frequency domain. So, we define the filtered projection by:

$$\mathbf{Q}\_{\Theta}(\mathbf{t}) = \left| \int\_{-\infty}^{+\infty} \mathbf{P}\_{\Theta}(\omega) \mid \omega \mid \mathbf{e}^{\ast 2\pi i \omega t} \mathbf{d}\omega \right. \tag{16}$$

**Figure 4.** Scanning a 2D object (details: circle and square) and the corresponding projections [7].

layer of the object.

artefacts when the number of projection is not sufficient.

**Figure 5.** Example of projection Sinogram of the object of Fig. 4: the y-axis shows the projections (gray level variation as a function of the pixel's or detector's position) and on the x-axis, these one pixel width projections are arranged from the first projection (θ=0) to the last one (θ=180°). Each sinogram will be used for the reconstruction of one pixel

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples

http://dx.doi.org/10.5772/52351

97

**Figure 6.** Reconstructions obtained for 4 different numbers of projections to illustrate the appearance of star-shaped

Thus, the object function will be simply given by:

$$\mu(\mathbf{x}, \mathbf{y}) = \int\_0^\pi \left[ \mathbb{Q}\_\theta(t) \right] d\theta = \int\_0^\pi \left[ \mathbb{Q}\_\theta(\mathbf{x}\cos\theta + \mathbf{y}\sin\theta) \right] d\theta \tag{17}$$

Knowing previously that a product (multiplication) in Fourier space (frequency domain) corresponds to a convolution of inverse Fourier transforms in real space (spatial domain), the following relation can be written:

$$\text{FT}^{-1}[\text{P}\_{\theta}(\omega) \times \mid \omega \mid \mathbf{I} \mathbf{=P}\_{\theta}(\mathbf{t}) \otimes \text{FT}^{-1}\mathbf{I} \mid \omega \mid \mathbf{I} \tag{18}$$

The convolution operator is denoted by the symbol ⊗. The function |ω| is not a squareintegrable function, thus it has no inverse Fourier transform. However, the inverse transform FT-1[|ω|] can be approximated by different filter response functions (convolution with gains, filter functions). Considering these last remarks and regarding the result of Eq.17, The reconstruction process of μ (x, y) can be performed from the projections Pθ(t) obtained as follows:

<sup>&</sup>quot;*In FBP reconstruction process, each projection Pθ(t) is first convoluted with a specific and suitable filtering function (eg Shepp-Logan, ω/2π sinc(ω)) and the measured values are recorded in (x, y) plans as illustrated in Fig. 4. To control the measurement process of tomography, the projections are arranged in a sinograms*<sup>2</sup> *recorded as shown in Fig. 5. On the sinogram, each vertical line is a projection at an angle θ and represents the variation in the gray level as a function of the pixel's or detector's position. On Fig. 6 are shown 4 different reconstructions. The projections Pθ(t) are converted to grayscale, convoluted with a filter with a specific gain and back-projected onto the entire plan (x, y). The addition of all projections result in the reconstruction of the 2D layer desired. More the number of projections is high more the reconstruction plans is well covered in terms of data and less the star-shaped artefacts are present on the reconstructed image. With such reconstructions, 3D images can be obtained by stacking all the 2D layers and a reconstructed 3D volume data (details of the object) can be extracted from the stack obtained*".

<sup>2</sup> the number of sinograms is equal to the sampling number of S

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples http://dx.doi.org/10.5772/52351 97

**Figure 4.** Scanning a 2D object (details: circle and square) and the corresponding projections [7].

μ(x,y)=*∫* 0 *π ∫* −*∞* +*∞*

> Qθ(t)=*∫* -*∞* +*∞*

p

q

*reconstructed 3D volume data (details of the object) can be extracted from the stack obtained*".

2 the number of sinograms is equal to the sampling number of S

=*∫* 0 *π ∫* −*∞* +*∞*

Thus, the object function will be simply given by:

following relation can be written:

follows:

F(ω,θ).e+2πiωt|*ω* |d*ω* dθ

Pθ(ω)|ω|e+2πiωtdω (16)

éù é <sup>ù</sup> ò ò ëû ë <sup>û</sup> (17)

FT-1 Pθ(ω)<sup>×</sup> <sup>|</sup>ω<sup>|</sup> =Pθ(t) <sup>⊗</sup> FT-1 <sup>|</sup>ω<sup>|</sup> (18)

(15)

F(ω,θ).e+2πiωt|*ω* |d*ω* dθ

96 Imaging and Radioanalytical Techniques in Interdisciplinary Research - Fundamentals and Cutting Edge Applications

ramp filter function in the frequency domain. So, we define the filtered projection by:

<sup>θ</sup> 0 0 μ(x,y)= ( ) dθ= Q (xcosθ+ysinθ) dθ *Q t*

Knowing previously that a product (multiplication) in Fourier space (frequency domain) corresponds to a convolution of inverse Fourier transforms in real space (spatial domain), the

The convolution operator is denoted by the symbol ⊗. The function |ω| is not a squareintegrable function, thus it has no inverse Fourier transform. However, the inverse transform FT-1[|ω|] can be approximated by different filter response functions (convolution with gains, filter functions). Considering these last remarks and regarding the result of Eq.17, The reconstruction process of μ (x, y) can be performed from the projections Pθ(t) obtained as

"*In FBP reconstruction process, each projection Pθ(t) is first convoluted with a specific and suitable filtering function (eg Shepp-Logan, ω/2π sinc(ω)) and the measured values are recorded in (x, y) plans as illustrated in Fig. 4. To control the measurement process of tomography, the projections are arranged in a sinograms*<sup>2</sup> *recorded as shown in Fig. 5. On the sinogram, each vertical line is a projection at an angle θ and represents the variation in the gray level as a function of the pixel's or detector's position. On Fig. 6 are shown 4 different reconstructions. The projections Pθ(t) are converted to grayscale, convoluted with a filter with a specific gain and back-projected onto the entire plan (x, y). The addition of all projections result in the reconstruction of the 2D layer desired. More the number of projections is high more the reconstruction plans is well covered in terms of data and less the star-shaped artefacts are present on the reconstructed image. With such reconstructions, 3D images can be obtained by stacking all the 2D layers and a*

 p

The integral in brackets can be regarded as the inverse Fourier transform Pθ(t) of Pθ(ω). However, we note that it is multiplied by the function |ω| which plays the role of a special

**Figure 5.** Example of projection Sinogram of the object of Fig. 4: the y-axis shows the projections (gray level variation as a function of the pixel's or detector's position) and on the x-axis, these one pixel width projections are arranged from the first projection (θ=0) to the last one (θ=180°). Each sinogram will be used for the reconstruction of one pixel layer of the object.

**Figure 6.** Reconstructions obtained for 4 different numbers of projections to illustrate the appearance of star-shaped artefacts when the number of projection is not sufficient.

### *1.1.3. Discretization of the analytical methods*

In the last sections, many times we have considered the ideal continuous case to model the projection and the reconstruction processes of tomography and to explain the principle of the analytical reconstruction methods such as FBP. Indeed, we have worked in an infinite continuous domain (R2 ), and tried to reconstruct a continuous function f from continuous projections P<sup>θ</sup> defined for all angle θ in the interval [0, π [. These conditions are obviously not possible in practice. Acquisition systems allow just obtaining a number of projections Pθ for a finite number of angles, denoted θk. The limited number of detectors makes these projections sampled and known only at discrete points uk. It would be unrealistic with such data to try to reconstruct a continuous function f in R2 . Therefore, the object function f(x,y) will be recon‐ structed on a discrete grid for a finite number of points (f(xi ,yj )). This is also the limit imposed by the used numerical reconstruction algorithms. The reconstruction problem then will be approached as follows: for given a set of projection measurements:

$$\left\{\mathbf{P}\_{\theta\mathbf{k}}(\mathbf{u}\_{1}), \mathbf{0} \le \mathbf{l} \le \mathbf{S}, \mathbf{0} \le \mathbf{k} \le \mathbf{P}\right\},\tag{19}$$

**2. Algebraic methods of image reconstruction in tomography**

each projection line, by the following equation [6]:

the inverse problem f = Rt

.p.

"pl

The algebraic methods are especially iterative methods. These methods are less used when compared to the easy and well known filtered back-projection method. In this method the reconstruction problem is see differently and no longer refers to the Radon transform. The image of the object consists of a number k of pixels whose values fk are unknown. Similarly, the projections are discrete and formed by a number of l dexels (depth pixels) whose values

Mathematics and Physics of Computed Tomography (CT): Demonstrations and Practical Examples

http://dx.doi.org/10.5772/52351

99

" are known since they correspond to measurements in each line of projection. Reconstruc‐ tion of the object image by the iterative method is based on the following hypothesis: each detected values in a dexel is a linear combination of pixel's values to be reconstructed [6]. The reconstruction problem is formulated by a discrete expression of matrix (p = R.f) describing the projection process. The set of values of the projection lines (dexels) is arranged in projection vector p. All pixels of the image to be reconstructed are also grouped in a image vector f. The coefficients which characterize the contribution of each pixel in each line of projection is determined and stored in a matrix R. The projection of an object is modeled in the case of an original image of n x n pixels, a number n of projection directions and a number n of dexels in

> 1 1 11 14 2 2

KKK

é ù é ù é ù ê ú ê ú ê ú ê ú ê ú ê ú <sup>=</sup> ê ú ê ú ê ú ê ú ê ú ê ú ê ú ë û ë û ë û

M M M O M M M O KKK

*p f r r p f*

M O

41 44

This last equation expresses the fact that what is detected (p) is the result of values (f) of the image to be reconstructed, subject to a projection operation represented by the projection operator R [6]). Through this modeling of the projection process, one looks in practice to find f according to p by solving the inverse problem f = R-1.p. Because of the size of this equations system, the resolution cannot be performed except by successive iterations [7]. Fig. 7 shows an example of projection and back-projection. In the algebraic reconstruction method based on

other than the transposed matrix of R. Thus, the reconstruction problem is limited in solving

The resolution of the inverse problem by iterative methods consists in finding a solution f minimizing the distance d between p and R.f where p and R are known. Here, it is question to start from an arbitrary estimate of the image solution and to proceed schematically to the correction the first estimate basing on a principle of trial and error. Each next estimate is projected and the result obtained is compared to the measured projection. The returned error is used to improve the next estimate. This method leads to build gradually low-to-high

*n n*

*p f r r*

iterative process, the back-projection is modelled by a back-projection operator Rt

.

p=R.f (22)

(23)

who is none

the reconstruction problem is reduced to finding f at any point of a finite discrete grid:

$$\left| \{ \mathbf{f}(\mathbf{x}\_{i}, \mathbf{y}\_{j}) , 0 \le \mathbf{i} \le \mathbf{S} , 0 \le \mathbf{j} \le \mathbf{S} \} \right| \tag{20}$$

with:

$$\begin{aligned} \boldsymbol{\theta}\_{k} &= \mathbf{k} \Delta \boldsymbol{\theta}\_{\prime} & \Delta \boldsymbol{\theta} &= \boldsymbol{\pi} / \mathbf{P}\_{\prime} \\ \mathbf{x}\_{l} &= \mathbf{i} \Delta \mathbf{x}\_{\prime} & \mathbf{y}\_{j} &= \mathbf{j} \Delta \mathbf{y}\_{\prime} \end{aligned} \quad \begin{aligned} \boldsymbol{\mu} &= \boldsymbol{\pi} / \mathbf{P}\_{\prime} \\ \boldsymbol{u}\_{l} &= \mathbf{i} \boldsymbol{\Delta} \boldsymbol{y}\_{\prime} \end{aligned} \tag{21}$$

where Δθ is the sampling step of the rotation angles, P is the number of angles (projections), d is the sampling step on each projection line (ray), Δx and Δy are the sampling step on x and y in the reconstruction plan.

The discrete reconstruction methods according to this definition can be divided into two categories:

1. The first class method consists in the definition of discrete operators and functions equivalent to those defined in the continuous case (Radon transform, back-projection, Fourier transform, etc.) and all the inversion formulas and theorems defined for the analytical methods already presented.

2. The second category is based on a completely different approach: in these methods the projection equation (p=Rf) is directly discredited and a linear equations system is built. The resolution of this system is possible only by iterative methods called algebraic methods. The algebraic method will be presented in the next section.
