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

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 expression [7]:

$$f^{n+1} = f^n + \alpha^n d^n \tag{25}$$

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) that are recalculated in a conjugated manner for each iteration and this to optimize the speed of convergence [6].

$$\alpha^{\prime \prime} = \frac{\left\| R^{\prime} p - R^{\prime} R f^{\prime \prime} \right\|^2}{(R^{\prime} p - R^{\prime} R f^{\prime \prime})^{\prime} R (R t \text{ p } - R t \text{ R fm})} \tag{26}$$

p1 p2

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

<sup>1</sup> / <sup>4</sup> <sup>3</sup> / <sup>4</sup> (E.1)

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

103

<sup>2</sup> (E.3)

*<sup>n</sup>* <sup>+</sup> (*pk* <sup>−</sup> *pk <sup>n</sup>*) / *Rii* (E.5)

*n updating. The results obtained are shown n the following matrix:*

*n from the equation of P1(Eq.E.3(2)) for each*

(E.2)

(E.4)

f 1 f 2

*For each projected line, the projector matrix is modeled by considering a value proportional the area of the pixel covered by the ray (projection line).For example for p1 the projection ray cover approximately ¾ of the area of the first pixel (f1) and ¼ of the area of the second pixel (p2). Therefore, the projection operator (matrix) can is given by: <sup>R</sup>* <sup>=</sup> <sup>3</sup> / <sup>4</sup> <sup>1</sup> / <sup>4</sup>

*Supposing that the object function in known (f1=2, f2=3). In this case the projection values are equal to:*

<sup>4</sup> *<sup>f</sup>* <sup>1</sup> <sup>+</sup> <sup>1</sup>

<sup>4</sup> *<sup>f</sup>* <sup>1</sup> <sup>+</sup> <sup>3</sup>

<sup>4</sup> *<sup>f</sup>* <sup>2</sup><sup>=</sup> <sup>9</sup> <sup>4</sup> (1)

<sup>4</sup> *<sup>f</sup>* <sup>2</sup><sup>=</sup> <sup>11</sup> <sup>4</sup> (2)

> *Rki* ∑*<sup>j</sup> Rkj*

*During the iteration process just one ray equation (p1 or p2) is used per iteration alternatively. Results of the iteration and*

*<sup>n</sup>* 2.70 3.26 2.25 2.45 2.20 2.25 2.06 2.10 2.02 2.02 2.02 2.02 2.00

*<sup>n</sup>* 0.90 2.58 2.24 2.84 2.75 2.90 2.84 2.96 2.93 2.98 2.97 2.99 2.98

*<sup>n</sup>* / 3.09 / 2.54 / 2.41 / 2.31 / 2.27 / 2.26 /

*<sup>n</sup>* 1.35 / 2.24 / 2.61 / 2.64 / 2.70 / 2.73 / 2.73

*We can easily verify that this algorithm gives a solution which approximate well the real one (2,3) after 12 iterations*

*We try know to compute the object function data (f1, f2) using some practical iterative algorithms and to compare the values obtained to real ones (f1=2, f2=3). For all of the following algorithms, the initial iteration conditions are:*

*<sup>n</sup>* <sup>+</sup> (*pk* <sup>−</sup> *pk <sup>n</sup>*)

Iteration(n) 1 2 3 4 5 6 7 8 9 10 11 12 13

*<sup>p</sup>*<sup>1</sup> <sup>=</sup><sup>3</sup>

*<sup>p</sup>*<sup>2</sup> <sup>=</sup><sup>1</sup>

*f i <sup>n</sup>*+1 <sup>=</sup> *<sup>f</sup> <sup>i</sup>*

> *f i <sup>n</sup>*+1 <sup>=</sup> *<sup>f</sup> <sup>i</sup>*

*n is calculated from the equation of P1(Eq.E.3(1)) and f1*

*n and p2*

*Fig. E.1. Object (2 pixels) and projection model*

*f1 0=f2*

*0=0="/> p1*

*0= p2 0=0.*

f1

f2

p1

p1

*2. SIRT algorithm with Jacobi method*

*iteration by taking into consideration p1*

*In this method, the iteration updating is given by:*

*(2.02, 2.99)*

*In this method f1*

*1. Additive ART algorithm with Kackzmarz method In this method, the iteration updating is given by:*

*p1, p2 updating are shown on the following matrix.*

$$\begin{aligned} \text{Attention.1}: \qquad &d^1 = \mathcal{R}^t p - \mathcal{R}^t \mathcal{R} f^0\\ \text{Attention.2}: \quad &d^n = (\mathcal{R}^t p - \mathcal{R}^t \mathcal{R} f^n) + b^n d^{n-1} \end{aligned} \tag{27}$$

Through this process, the error between the measured projections and those calculated is minimized progressively. This error evaluation formula is given by:

$$e = \left\| p - Rf \right\|^2\tag{28}$$

#### **3. ART algorithm**

The formula of the iterative updating of this algorithm is given by the following expression:

$$f\_i^{(n+1)} = f\_i^{(n)} + \mathcal{R}\_{ji} \frac{p\_j - \mathcal{R}\_j f^{(n)}}{\left\| \mathcal{R}\_j \right\|^2} \tag{29}$$

The correction in this algorithm can be additive (e=pk-pk n) or multiplicative (e= pk/pk n).

#### **4. SIRT algorithm**

The formula of the iterative updating of this algorithm is given by the following expression:

$$f\_i^{(n+1)} = f\_i^{(n)} + \frac{\sum\_j p\_j}{\sum\_j \sum\_i R\_{ji}} - \frac{\sum\_j R\_j f^{(n)}}{\sum\_j \left\| R\_j \right\|^2} \tag{30}$$

#### *Practical example 1*

*To illustrate the difference between the iterative algorithms, let consider the example of figure (E.1). In this example we consider an object (image) composed by 2 pixels which are projected as shown on this figure.*

#### *Fig. E.1. Object (2 pixels) and projection model*

that are recalculated in a conjugated manner for each iteration and this to optimize the speed

( ) (Rt p - Rt R fn)

*t t nt*


.1 :

minimized progressively. This error evaluation formula is given by:

( 1) ( )

The correction in this algorithm can be additive (e=pk-pk

( 1) ( )

*consider an object (image) composed by 2 pixels which are projected as shown on this figure.*

*i i*

*R p R Rf R*

*iteration d R p R Rf iteration d R p R Rf b d* -

*n*

a

*t tn*

*R p R Rf*

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

.2 : ( )

2

1 0

= - =- +

Through this process, the error between the measured projections and those calculated is

The formula of the iterative updating of this algorithm is given by the following expression:

The formula of the iterative updating of this algorithm is given by the following expression:

å å

*Practical example 1 To illustrate the difference between the iterative algorithms, let consider the example of figure (E.1). In this example we*

*ji j i <sup>J</sup> <sup>j</sup>*

*p Rf*

*j j n n j j*

*f f <sup>R</sup> <sup>R</sup>* <sup>+</sup> =+ -

*n n j j i i ji*

*p Rf f fR*

( )

2 *n*

*j*

*R*

*t t n t t n nn*

1

<sup>2</sup> *e p Rf* = - (28)

<sup>+</sup> - = + (29)

( )

2 *n*

å å <sup>å</sup> (30)

n) or multiplicative (e= pk/pk

n).

(26)

(27)

of convergence [6].

**3. ART algorithm**

**4. SIRT algorithm**

*For each projected line, the projector matrix is modeled by considering a value proportional the area of the pixel covered by the ray (projection line).For example for p1 the projection ray cover approximately ¾ of the area of the first pixel (f1) and ¼ of the area of the second pixel (p2). Therefore, the projection operator (matrix) can is given by:*

$$\mathcal{R} = \begin{bmatrix} \mathbb{B} \ \mathbb{A} & \mathbb{1} \ \mathbb{A} & \mathbb{1} \ \mathbb{A} \\ \mathbb{1} \ \mathbb{A} & \mathbb{B} \ \mathbb{A} \end{bmatrix} \tag{E.1}$$

*Supposing that the object function in known (f1=2, f2=3). In this case the projection values are equal to:*

$$\begin{aligned} \rho\_1 &= \frac{3}{4} \int\_4 \, f\_1 + \frac{1}{4} \int\_4 \, f\_2 = \frac{9}{4} \{1\} \\ \rho\_2 &= \frac{1}{4} \int\_4 \, f\_1 + \frac{3}{4} \int\_4 \, f\_2 = \frac{11}{4} \{2\} \end{aligned} \tag{5.2}$$

*We try know to compute the object function data (f1, f2) using some practical iterative algorithms and to compare the values obtained to real ones (f1=2, f2=3). For all of the following algorithms, the initial iteration conditions are: f1 0=f2 0=0="/> p1 0= p2 0=0.*

#### *1. Additive ART algorithm with Kackzmarz method*

*In this method, the iteration updating is given by:*

$$f\_{l}^{n+1} = f\_{l}^{n} + (\mathfrak{p}\_{k} - \mathfrak{p}\_{k}\,^{n})^{\mathcal{R}\_{kl}} \Big| \sum\_{j} \,^{}\_{k\dot{j}} \mathcal{R}\_{kj} \,^{2} \tag{\text{E.3}}$$

*During the iteration process just one ray equation (p1 or p2) is used per iteration alternatively. Results of the iteration and p1, p2 updating are shown on the following matrix.*


*We can easily verify that this algorithm gives a solution which approximate well the real one (2,3) after 12 iterations (2.02, 2.99)*

#### *2. SIRT algorithm with Jacobi method*

*In this method, the iteration updating is given by:*

$$\{f\_i^{n+1} = f\_i^n + (\mathfrak{p}\_k - \mathfrak{p}\_k \iota^n) \mid \mathcal{R}\_k\} \tag{\mathbb{E}.5}$$

*In this method f1 n is calculated from the equation of P1(Eq.E.3(1)) and f1 n from the equation of P1(Eq.E.3(2)) for each iteration by taking into consideration p1 n and p2 n updating. The results obtained are shown n the following matrix:*


introduction to integral computed tomography. We then go over to projection and inversion algorithms. We give a detailed analysis of the filtered back-projection algorithm in the light of the sampling theorem. We also describe the convergence properties of iterative algorithms.

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

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

105

Good reconstructions from the interlaced lattice can also be obtained by using the direct algebraic reconstruction algorithm, or by increasing the amount of data through the interpo‐ lation according to the sampling theorem. As we have seen, the interpolation step can introduce significant errors in certain cases. It has also been shown that the interpolation can be avoided by choosing the points x where the reconstruction is computed on a polar grid rather than on a square Cartesian grid, and interchanging the order of the two summations. This algorithm should work well for the interlaced lattice and is particularly beneficial in case of the fan-beam sampling geometry, since the method also avoids the homogeneous approx‐

The cone and fan beam scanning are the standard scanning modes in present day's clinical practice. The methods described above assumed that the geometry of the acquisition was a parallel geometry, like that of first generation systems. In the case a conical geometry (or fan),

**1.** The first is to ignore the discrepancy. The error induced by this approximation is consid‐ ered negligible if the beam angle is low (typically below 15 degrees). This method is applicable to systems of second generation, but more to the following where the beam

**2.** The second method is to rearrange the data in parallel projections, mainly by interpolation. **3.** The last method is to reformulate the problem completely. It becomes apparent that the projection theorem cannot be generalized, which does not have direct inversion method. The projection theorem can be adapted to different geometries, and then the algorithm is the same as for a parallel geometry, with the same calculations. Finally, the filtered back projection formulas can be corrected and result in a slightly different algorithms. The algebraic methods have some fewer additional problems, but the flexibility in the choice of basic functions and thus the coefficients of the matrix R, allows these methods to be

Finally, I want just to add that the purpose of this chapter is to give an introduction to the studied topic and treat some related aspects in more detail. The reader interested in a broader overview of the field, its relation to various branches of pure and applied mathematics, and

its development over the years may wish to consult the appropriate bibliography.

Department of Physics, Faculty of Science, University of Ferhat Abbas-Sétif, Algeria

must cover the whole section of the patient to remove the translation.

We shortly mention Fourier based algorithms used in CT.

three methods are possible:

adapted to different geometries.

**Author details**

Faycal Kharfi

imation, whose influence on the reconstruction is difficult to estimate.

*With this algorithm and method the best solution is obtained after 9 iterations.*

#### *3. SIRT algorithm with Gauss-Seidel method*

*This method is similar to the Jacobi method and use the same updating function (E.5) but just for the first iteration which allows the determination of f1 1 value. After that f2 1 is estimated from p2 equation considering the obtained value f1 1 which itself will be used to determinate f1 2 from p1 equation and so on (by matching alternatively between p1 and p2) until obtaining the best results. The results obtained are shown on the following matrix:*

$$\begin{array}{ccccccccc}\text{Attention}(n) & 1 & 2 & 3 & 4\\\text{ $f$ }\_{1}n & 3 & 1.18 & 1.90 & 1.98\\\text{ $f$ }\_{2}n & 2.66 & 3.27 & 3.03 & 3.00\\\end{array} \tag{5.7}$$

*We remark that this algorithm and method gives good results just after 4 iterations; so it is the very fast one in term of convergence.*

#### *End of practical example 1*

#### **2.3. Iteration stopping rules**

The iterative process for image reconstruction must be stopped at the end of the conver‐ gence phase before it starts diverging. This can be done using some regularization meth‐ ods. Indeed, well selected regularization parameter can controls the amount of stabilization imposed on the solution. In iterative methods one can use the stopping index as regulariza‐ tion parameter. When an iterative method is employed for image reconstruction, the user can also study on-line adequate visualizations of the iterates as soon as they are computed, and simply halt the iteration when the approximations reach the desired quality. This may actually be the most appropriate stopping rule in many practical applications, but it requires a good intuitive imagination of what to expect. If this is not the case, a computer can give us some aid to determine the optimal approximation. The stopping rules are divided into two categories: rules which are based on knowledge of the norm of the errors, and rules which do notrequire such information. If the error norm is known within reasonable accuracy, the perhaps most well known stopping rule is the discrepancy principle Morozov [8]. Examples of the second category of methods are the L-curve criterion [9], and the general‐ ized cross-validation criterion [10].

### **3. Conclusions**

After a substantial effort, major breakthroughs have been achieved in the last fourteen years in the mathematical modeling of CT. The aim of this chapter is to survey this progress and to describe the relevant models, mathematical problems and reconstruction procedures used in CT. We give a summary on the mathematics of computed tomography. We start with a short introduction to integral computed tomography. We then go over to projection and inversion algorithms. We give a detailed analysis of the filtered back-projection algorithm in the light of the sampling theorem. We also describe the convergence properties of iterative algorithms. We shortly mention Fourier based algorithms used in CT.

Good reconstructions from the interlaced lattice can also be obtained by using the direct algebraic reconstruction algorithm, or by increasing the amount of data through the interpo‐ lation according to the sampling theorem. As we have seen, the interpolation step can introduce significant errors in certain cases. It has also been shown that the interpolation can be avoided by choosing the points x where the reconstruction is computed on a polar grid rather than on a square Cartesian grid, and interchanging the order of the two summations. This algorithm should work well for the interlaced lattice and is particularly beneficial in case of the fan-beam sampling geometry, since the method also avoids the homogeneous approx‐ imation, whose influence on the reconstruction is difficult to estimate.

The cone and fan beam scanning are the standard scanning modes in present day's clinical practice. The methods described above assumed that the geometry of the acquisition was a parallel geometry, like that of first generation systems. In the case a conical geometry (or fan), three methods are possible:


Finally, I want just to add that the purpose of this chapter is to give an introduction to the studied topic and treat some related aspects in more detail. The reader interested in a broader overview of the field, its relation to various branches of pure and applied mathematics, and its development over the years may wish to consult the appropriate bibliography.
