Radon transform is a transform that is used for example in computed tomography. Revise the equations describing the Radon transform in 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$.

Slides from lecture (slides 18+ deal with ART). A more detailed description is in an article [1] or in Chapter 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:

1. Compute the Radon transform for each $\theta$ and translation $p$ (see Equation 1) $$ \hat{J}(p,\theta) = R [ \hat{f}(x,y) ]. ~~~~~~~~~ (4) $$
2. 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)^2} ~~~~~~~~~ (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)^2} ~~~~~~~~~ (7) $$

The learning rate is a tuning parameter that controls the step size in the iterative update of the image reconstruction. The update formula with modified learning 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)^2} ~~~~~~~~~ (8) $$ for learning rate increasing learning rate \alpha leads to faster converges, but the algorithm may skip the optimal value and oscillate around it or become unstable.

For this exercise we will use data from 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 reconstructed image $f^*(p, q) = 1; \forall p, q$, where all the pixel values are 1, see Equation (6). Put an 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]$.

1. Initialize the reconstructed image as $f_0(x, y) = 0$.
2. For each $\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
3. Repeat steps 2-4 for all projection angles $\theta$.
4. Perform the total number of iterations nIter

Hints - 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).

[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.