===== Lab 12: Algebraic reconstruction technique =====

===== Radon transform =====

[[http://en.wikipedia.org/wiki/Radon_transform|Radon transform]] is a transform that is used for example in computed tomography. Revise the equations describing the Radon transform in [[http://cmp.felk.cvut.cz/cmp/courses/ZSL2/cviceni/ct1/rovnice.pdf|Radon transform equations]]. 

In the previous lab we used the Radon transform in a continuous space and we discritized the equation in order to apply it to digital images. Another form of the Radon transform for digital images is in Equation 1. This form is more suitable for the reconstruction of a sinogram by ART. 

$$
J(p,\theta) = R[f(x,y)] = \sum_q f(x',y') \cdot h \approx \sum_q \sum_{x,y} f(x,y) \cdot \zeta (x,y,q,p,\theta) ~~~~~~~~~ (1)
$$

The function $\zeta$ is defined as

$$
\zeta (x,y,q,p,\theta) = \begin{cases} 1 & x = \textrm{round}(x') \\  & y = \textrm{round}(y') \\ 0 & \textrm{otherwise} \end{cases} ~~~~~~~~~ (2)
$$

where $x′,y′$ are the computed coordinates in the original image for an angle of rotation $\theta$ and the coordinates $p, q$ are given by

$$
x'=p\cos\theta-q\sin\theta ~~~~;~~~~ y'=p\sin\theta+q\cos\theta ~~~~~~~~~ (3)
$$

and where $q=ih$ for $i \in \mathbb{Z}$ and $q$ is from interval $q \in [ -\frac{\sqrt{2}N}{2}; \frac{\sqrt{2}N}{2} ]$ (assuming $h = 1$ for simplicity) for an image of size $N \times N$ . Equation 3 describes the relationship between the $x,y$ coordinates in the original image and the rotated coordinates $p,q$.

===== Algebraic reconstruction technique =====

Slides from [[https://cw.fel.cvut.cz/wiki/_media/courses/zsl/ct2as.pdf|lecture]] (slides 18+ deal with ART). A more detailed description is in an article [1] or in Chapter [[https://engineering.purdue.edu/~malcolm/pct/CTI_Ch07.pdf|Chapter 7 of Principles of Computerized Tomographic Imaging]].

Algebraic reconstruction technique works by iteratively applying corrections to the pixels in the reconstructed image so to minimize the difference between the original sinogram and the sinogram of the reconstructed image.

First, we choose an initial estimate of the reconstructed image or just initialize the reconstructed image with zeros $f_0(x, y) = 0$ of size $N \times N$ pixels.
Then in each iteration of the algorithm we update the reconstructed image for each $(p, \theta)$ as follows:

  - Compute the Radon transform for each $\theta$ and translation $p$ (see Equation 1) $$
\hat{J}(p,\theta) = R [ \hat{f}(x,y) ]. ~~~~~~~~~ (4) $$
  - The difference $J(p,\theta)−\hat{J}(p,\theta)$ between the computed sinogram $J^*(p,\theta)$ and the measured sinogram $J(p,\theta)$ divided by weight function $w(p,\theta)$ indicating how many pixels were measured for each translation $p$ of this projection.
Update the reconstructed image by adding the average error for each q in the projection $(p, \theta)$ $$
f_{k+1} (x,y) = f_k (x,y) + \sum_{p,q,\theta} \zeta (x,y,q,p,\theta) \cdot \frac{J(p,\theta)−\hat{J}(p,\theta)}{w(p,\theta)} ~~~~~~~~~ (5)
$$ the weight function $w$ for image of ones is defined as $$
w(p,\theta) = R[1] (p,\theta) = N_q (p,\theta) \cdot h ~~~~~~~~~ (6)
$$ where $N_q$ is the number of valid values for q for given projection $(p, \theta)$.
Efficient implementation can use different formulation $$
\forall q: ~~~~ f(x,y) \leftarrow f(x,y) + \frac{J(p,\theta)−\hat{J}(p,\theta)}{w(p,\theta)} ~~~~~~~~~ (7)
$$

=== Learning rate ===

The learning rate is a tuning parameter that controls the step size in the iterative update of the image reconstruction. The updated formula with thelearning rate $\alpha$ is the following: $$
\forall q: ~~~~ f(x,y) \leftarrow f(x,y) + \alpha \cdot \frac{J(p,\theta)−\hat{J}(p,\theta)}{w(p,\theta)} ~~~~~~~~~ (8)
$$
 increasing learning rate $ \alpha > 1 $ leads to faster converge, but the algorithm may skip the optimal value and oscillate around it or even become unstable. Decreasing learning rate $0 > \alpha < 1$ leads to slower convergence.

=== Gordon ART ===

The update rule in the form presented in Equation (7) is the Simple ART algorithm. A simple modification of the update rule called the Gordon ART (see [1]) can yield an algorithm with a slightly slower rate of convergence, but more stable and more precise. For Gordon ART consider the update rule in the from:$$\forall q: ~~~~ f(x,y) \leftarrow f(x,y) + \frac{J(p,\theta)−\hat{J}(p,\theta)}{w(p,\theta)^2} ~~~~~~~~~ (9)$$
===== Homework =====

For this exercise we will use data from {{ :courses:zsl:artdata.mat | download here}}. The archive contains both images (''imageCT'') and sinograms (''rawCT''); the phantom dataset contains two additional files with added noise. 

**HW 1.1** [0.5 pt] Visualize the weight function for the reconstruction of the image ''phantom''. The weight function is a sinogram $R[f^*(p,q)]$ of an image of the same size as the coordinate system of the reconstructed image $f^*(p, q) = 1; \forall p, q$, where all the pixel values are 1, see Equation (6). Note that the weights are computed in the larger auxiliary coordinate system of $(p, q)$ and not in the coordinate system of the image to be reconstructed $(x,y)$. Visualize the weights and put the image into the report.

**HW 1.2** [1 pt] Implement function ''image=myARTiter(radonim,theta,nIter)'' which will perform ''nIter'' iterations of algebraic reconstruction technique of the measured sinogram ''radonim'' $[m \times N_\theta]$ for projection angles ''theta'' $[1 \times N_\theta]$.
   - Initialize the reconstructed image as $f_0(x, y) = 0$.
   - Choose the first $\theta$.
   - Compute $J_k(p,\theta)$ and compare it to the input data $J(p,\theta)$.
   - Update all pixels in the reconstructed image lying on projection lines for projection $(p, \theta)$ following the Equation 5
   - Repeat steps 2-4 for all projection angles $\theta$.
   - Perform the total number of iterations ''nIter''
   
Submit your implementation of the ''myARTiter'' alongside the report. 

    Hint - when weighing the update (dividing by) zero or near zero values, the results is ''NaN'' value, which will propagate with each iteration


**HW 1.3** [1.5 pts] Reconstruct images from sinograms ''phantom'', ''phantom_noise1'' and ''head2'' using algebraic reconstruction technique (ART). Put the reconstructed images after 1, 10 and 100 iterations of the algorithm for each sinogram into the report.

**HW 1.4** [1 pt] Estimate the reconstruction errors of the ART as the sum of squared differences (SSD) between the original images and the reconstructed images. $$ SSD = \sum_{\forall x} \sum_{\forall y} (\hat{f}(x,y) - f(x,y))^2 ~~~~~~~~~ (9)
$$ where $f(x,y)$ is the original image and $\hat{f}(x,y)$ is the reconstructed image. Visualize the relationship between the reconstruction error and the number of iterations of the ART used for the reconstruction.

**HW 1.5** [1 pt] Reconstruct the sinograms also with the naïve backprojection and filtered backprojection (use the methods from the previous lab). Estimate the reconstruction error between the original image and the reconstructed image using SSD and compare the results to the ART-reconstructed images.

**BONUS 1.6** [1pt] For the ART and the filtered backprojection examine the relationship between the reconstruction error (SSD) and the magnitude of the added Gaussian noise (sinograms ''phantom_noise1.mat'', ''phantom_noise2.mat'' and ''phantom_noise4.mat'') and the number of projection angles (see bonus exercise in the previous lab).

===== References =====

[1] P.M.V. Subbarao, P. Munshi, and K. Muralidhar. Performance of iterative tomographic algorithms applied to non-destructive evaluation with limited data. //NDT and E International//, 30(6):359 – 370, 1997. {{ :courses:zsl:subbaro_et_al_1996.pdf |PDF}}