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

*1.1.3. Discretization of the analytical methods*

reconstruct a continuous function f in R2

structed on a discrete grid for a finite number of points (f(xi

*k*

q

algebraic method will be presented in the next section.

*i j*

 q = D D = q p

*x ix y jy*

*k P*

approached as follows: for given a set of projection measurements:

continuous domain (R2

with:

categories:

presented.

y in the reconstruction plan.

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

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

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

by the used numerical reconstruction algorithms. The reconstruction problem then will be

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

, , , , ,

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

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

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

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

*l*

*u ld*

), and tried to reconstruct a continuous function f from continuous

,yj

{P (u ),0 l S,0 k P , θk l ££ £ £ } (19)

{f(x ,y ),0 i S,0 j S i j ££ ££ } (20)

<sup>=</sup> = D = D (21)

. Therefore, the object function f(x,y) will be recon‐

)). This is also the limit imposed

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 "pl " 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 each projection line, by the following equation [6]:

$$\mathbf{p} \text{=} \text{R.f} \tag{22}$$

$$
\begin{bmatrix} p\_1 \\ p\_2 \\ \vdots \\ \vdots \\ p\_n \end{bmatrix} = \begin{bmatrix} r\_{11} & \dots & \dots & \dots & r\_{14} \\ \vdots & \ddots & & & \\ \vdots & & \ddots & & \\ \vdots & & & \ddots & \\ r\_{41} & \dots & \dots & \dots & r\_{44} \end{bmatrix} \begin{bmatrix} f\_1 \\ f\_2 \\ \vdots \\ \vdots \\ f\_n \end{bmatrix} \tag{23}
$$

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 iterative process, the back-projection is modelled by a back-projection operator Rt who is none other than the transposed matrix of R. Thus, the reconstruction problem is limited in solving the inverse problem f = Rt .p.

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 frequencies of the image solution. The results of the first iterations are smooth because of the predominance of low frequencies (internal structure) of the object. Subsequently, more iteration are applied more high frequencies (overall shape and background noise) are repre‐ sented. The image produced by iterations approximate gradually the image solution (the algorithm converges) [6]. However, we show that when using such iterative methods and after a number of iterations, the process begins to diverge (under the influence of noise) and the image moves away from the true solution. To overcome this inconvenient, we impose a constraint to the reconstruction process to interrupt the iteration after a certain number of iterations. This is equivalent to use a low pass filter as in the case of filtered back-projection method.

**2. Geometric modeling of the projection operator**: for the determination of the coefficients of the projection matrix R, we must consider the number of projections and their angular distributions and the projection beam geometry which can be parallel or fan. If, for example, the intensity model distribution is that of Dirac, a given pixel with an index k crossed by a ray of an index pl generates a projection coefficient rlk equal to 1; if this is not

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

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

101

**3. Physical modeling of the projection operator**: this model is based on the distance between the position of the object pixel to the detector. A pixel located away from the detector will see its contribution to the projection ray reduced compared to a nearest one. With this modeling aspect the beam attenuation will be included in the resulting reconstructed

There are many algorithms that have been developed that used the iterative method for tomographic reconstruction (see practical example 1). Among these algorithms, the main ones are: the Algebraic Reconstruction Technique (ART), the Simultaneous Iterative Reconstruction Technique (SIRT), and the Iterative Least Squared Technique (ILST). Currently, the most commonly used algorithms are: the Expectation Maximization algorithm (EM) and the

This algorithm was developed and proposed by Lange and Carson. The formula of this

.

Where n is the number of the actual iteration. This algorithm is characterized by the fact that it keeps the number of iterations for each projection. Moreover its multiplicative form gives a

It is an algorithm which is used because it converges very quickly. It is based on a classical descent method. Its formula of iterative updating can be, roughly, given by the following

We can verify that the correction is not multiplicative as for EM algorithm, but it is additive. This formula is characterized by a descent direction d (Eq. 26) and a descent speed α (Eq.27)

*n n nn* <sup>1</sup> *ffd* a

*n*

<sup>+</sup> <sup>=</sup> (24)

<sup>+</sup> = + (25)

1

positivity constraint although it implies a slow convergence.

**2. Conjugate gradient algorithm (CG)**

*n nt*

*<sup>p</sup> f fR R f*

the case, zero value is assigned.

image.

**1. EM algorithm**

expression [7]:

**2.2. Main iterative algorithms**

Conjugate Gradient algorithm (CG).

algorithm is the following [6]:

**Figure 7.** Illustration of projection and back-projection processes.

#### **2.1. Projection modeling**

The projection process is modeled by considering the coefficients of the projection matrix R that generate the acquisition data. Some special geometric and physical considerations are necessary for this modeling. In the iterative reconstruction method the modeling concerns the following points:

**1. Modeling of the pixel (detector) intensity distribution**: it is necessary to specify the real conditions of image pixels projection [6]. It is based on the evaluation of the contribution of each pixel in the corresponding projection line (ray). Knowing that the most accurate model (perfect) and the more complicated to be applied consists in considering square pixel (uniform); simpler models are possible. Among these simple models, there is the model called Dirac model in which all the pixel intensity is concentrated in the centre of the pixel. Thus, the whole intensity of the pixel contrib‐ utes to the projection line (ray) if and only if it passes through the Dexel. There is also another model called the model of concave disc which is considered as a compromise between the two previous models. In this model, the intensity of the pixel is geometri‐ cally limited to a disk included in the pixel and distributed so that its projection is rectangular regardless of the direction of projection [6].

