Topics:

• Math of CT images – the forward projection

The homework has four consecutive parts. When completed, you will have your own implementation of the forward projection ( Radon transform) used in CT imaging. A set of projections with varying direction $\theta$ is acquired. The object Cartesian coordinate system (grid) is denoted $x, y$. We shall also define a rotated coordinate system $p$, $q$, chosen so that projection corresponds to integration along the $p$ axis, creating a 1D function of $q$ for each angle $\theta$. Illustration of both coordinate systems is given below:

(Implementing our own Radon transform will help us to understand the inverse transform, which is the key concept for image reconstruction from measured CT projections. This will be the topic of the next lab.)

Homework tasks

1. [0.5 pts] The relation between the two coordinate systems is given by \begin{align} p &=& x \cdot \cos \theta + y \cdot \sin \theta \\ q &=& -x \cdot \sin \theta + y \cdot \cos \theta \end{align} Find the smallest range $[-p_m, p_m] \times [-q_m, q_m]$ needed so that all points from $(x, y) \in [-x_m, x_m] \times [-y_m, y_m]$ fall into this range. Provide either a graphical or an analytical solution.
2. [0.5 pts] Discretize the Radon function \begin{equation} J(\theta,p) = \int_{-\infty}^\infty f\left ( p\cdot\cos(\theta) - q\cdot\sin(\theta) ~ , ~ p\cdot\sin(\theta) + q\cdot\cos(\theta) \right ) dq \end{equation} i.e. approximate the integral by a sum.
3. [3 pts] Implement the function

   function fwrad = myRadon(img, theta)

   which will return the forward Radon transform fwrad of the input image img of dimension $[m \times n]$, given the array of projections angles theta of length $k$ (angles are provided in degrees). Size of the resulting sinogram fwrad is $[(2p_m + 1) \times k]$. Submit the code into BRUTE.
4. [1 pt] Transform the Shepp-Logan phantom (matlab function phantom) of size 256 px for projection angles $\theta \in \{0^{\circ}, 1^{\circ}, 2^{\circ}, \ldots, 179^{\circ}\}$ with your function. Show the input and the resulting sinogram. Pay attention to labelling the axes correctly.

Hints

• presentation
• set the image coordinate system such that $[x, y] = [0, 0]$ lies in the centre of the image, this leads to the spatial range $[-x_m, x_m] \times [-y_m, y_m]$
• compute the maximal values $p_m$ and $q_m$
• use the discrete form of (3) and compute $J(\theta, p_h)$
  • think of points $p_h, q_h$ as of grid rotated by (a fixed) $\theta$
  • go over all $(p_h, q_h)$ on a grid, transforming each to a corresponding point $(x_h, y_h)$ in the input image
  • as $(x_h, y_h)$ will in general not coincide with pixel centers, use interpolation ( interp2 function) to get $img(x_h, y_h)$
  • compute the numerical value of $J(\theta, p_h)$ and store it into the fwrad sinogram.