===== Lab 04 : CT Reconstruction, forward projection =====

<note important>**Update 11.03.2020: **
Protože neběží výuka, přidávám sem {{ :courses:zsl:lab04_hints.pdf | PDF Dokument}} s pár komentáři k dnešnímu cvičení. Neváhejte se na mne obrátit mailem, pokud budete mít nějaké nejasnosti či dotazy. Odevzdávací lhůta bude tentokráte o 14 dnů posunuta.</note>

Topics:
  * Math of CT images -- the forward projection

/*
  * [[courses:zsl:labs2020_07_ultrasound | Ultrasound assignment (Lab 7)]]-  acquire data for US assignment if not done already
*/ 
==== HW04 Assignment ====

The homework has four consecutive parts. When completed, you will have your own implementation of the forward projection ([[http://en.wikipedia.org/wiki/Radon_transform| 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:

{{ :courses:a6m33zsl:radonproj.png?400 | Radon transform of f for projection p }}

(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**
  - **[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.
  - **[0.5 pts]** Diskretize 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 pts]** Implement the function<code matlab>
function fwrad = myRadon(img, theta)
</code> 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.
  - **[1 pt]** Transform the Shepp-Logan phantom (matlab function {{https://www.mathworks.com/help/images/ref/phantom.html| 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 axis correctly.

** Hints **
  * 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$ /* and then the discretized points $p_h, q_h$ */
  * 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 ({{https://www.mathworks.com/help/matlab/ref/interp2.html | 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.