Learning Tomography: ART (Part I)

Algebraic Reconstruction Techniques (ART) is another family of algorithms to compute reconstructions. Here is a first approach to understand the principle of these programs.


1- Principle

In a reconstruction, the aim is to try to determine the pixel values of the reconstructed image (or volume) from various projections. For example, in the 2x2 image of Fig.1, we have four unknowns x₀₀, x₁₀, x₀₁, and x₁₁ and three projections along horizontal, 45°, and vertical directions.

Fig. 1: Three projections along horizontal, 45°, and vertical directions, respectively of a 2x2 image.

 From the given three projections, we are able to build 7 equations:

$\{\begin{array}{c}{x}_{00}+x_{10}=7\\ x_{01}+x_{11}=9\\ x_{10}=6\\ x_{00}+x_{11}=3\\ x_{01}=7\\ x_{10}+x_{11}=8\\ x_{00}+x_{01}=8\end{array}$

Note: In this specific case, that's possible to solve these equations by hand from equations #3 and #5 (the corners of the diagonal projection), but this isn't the purpose of this post (actually, this is used by the Mojette transform, the subject of another post?).

These equations can also be written with a matrix as follows ...

$\left(\begin{array}{c}7\\ 9\\ 6\\ 3\\ 7\\ 8\\ 8\end{array}\right)\text{ = }\left(\begin{array}{cccc}1& 1& 0& 0\\ 0& 0& 1& 1\\ 0& 1& 0& 0\\ 1& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 1\\ 1& 0& 1& 0\end{array}\right)\text{ x }\left(\begin{array}{c}{x}_{00}\\ x_{10}\\ x_{01}\\ x_{11}\end{array}\right)$

p    =          A          .     x

where A is a matrix containing the paths of the various rays traversing the image to yield the projection pixels.
Thus, solving this system can be done by...

x = A^-1 . p

Due to the size of this system (in real projects, more than 100 projections of 2048x2048), it is solved using iterative techniques.

The Kaczmarz method

According to Kaczmarz, it is possible to solve this problem with an iterative process consisting of calculating the difference between the projection and the corresponding projection calculated from the reconstructed - temporary - image and re-injecting this value in the image.
Let's look at the first projection. In Fig. 2, the 2D-reconstruction image is initialized with zero. Then we compute the projection from the image along the horizontal direction yielding [0|0]. Then, we compute the difference 7 - 0 = 7 and divide by the number of pixels (here, 2) yielding 7 / 2 = 3.5. Finally, we add this value (3.5) to the previous pixel values of the row (0 + 3.5 = 3.5).

Fig.2: First projection

If we continue with the diagonal projection,the difference is equal to
(6 - 3.5)/ 1 = 2.5. Then, we update the reconstruction image by adding 2.5 to 3.5 (the previous pixel value) yielding a new pixel value of 6, and we continue with the other pixels as shown in Fig. 3...

Fig.3: Second diagonal projection. The calculations are from left to right: (6 - 3.5)/ 1 = 2.5, ( 3 - (3.5+4.5))/2 = -2.5, and (7 - 4.5)/1 = 2.5. These values are added to the previous pixel values.

 Finally, the last projection doesn't change anything, we stop the iteration process .... and have reconstructed our 2x2 image.


Next step, implement this algorithm... [Link]

Other crazybiocomputing posts

Further readings are available in ...

• Learning Tomography Series  [Link]
• Image Processing TOC [Link]