===== Lab 10: Algebraic reconstruction technique =====




/*

Vstupní data pro úlohu jsou v archivu {{ :courses:zsl:artdata.mat | ke stažení zde}}. V archivu jsou jak obrázky (''imgCT'') tak i jejich sinogramy (''rawCT''), pro dataset phantom jsou k dispozici navíc dva soubory s přidaným rušením.

 - **[2 pts]** Rekonstruujte data pro sinogramy ''phantom'', ''phantom_noise1'' a ''head2'' pomocí iterativní algebraické rekontrukční metody. Do zprávy vložte rekonstruované obrázky po 1, 10 a 100 iteracích.
  - **[2 pt]** Spočtěte rekonstrukční chybu iterativní metody jako součet druhých mocnin rozdílu originálního obrázku a vaší rekonstrukce.$$
SSD = \sum_{\forall x} \sum_{\forall y} (\hat{f}(x,y) - f(x,y))^2 ~~~~~~~~~ (1)
$$ kde $f(x,y)$ je originální a $\hat{f}(x,y)$ je rekonstruovaný obraz. Vykreslete jako graf v závislosti na počtu iterací.
  - **[1 pt]** Rekonstruujte všechny obrazy také pomocí přímé zpětné a filtrované zpětné projekce (použijte funkci z předcházejícího cvičení). Opět spočítejte $SSD$ rekonstrukcí oproti originálním obrázkům.
  - BONUS: **[1 pt]** Porovnejte závislost $SSD$ chyby obou metod (filtrované zpětné projekce a algebraické rekonstrukce) na amplitudě přidaného aditivního Gaussovského šumu (data phanom_noise1, 2, 4) a na počtu projekčních úhlů (viz bonus předchozího cvičení).


===== Radonova transformace =====

[[http://en.wikipedia.org/wiki/Radon_transform|Radon transform]] is a transform that is used for example in computed tomography. Podívejte se znovu na [[http://cmp.felk.cvut.cz/cmp/courses/ZSL2/cviceni/ct1/rovnice.pdf|rovnice]] popisující Radonovou transformaci (2). Projděte si též slidy z [[http://cmp.felk.cvut.cz/cmp/courses/33zsl2leto2006/slidy/ct-hozman-jk.pdf|přednášky]] a prostudujte si Iterativní rekonstrukci (slide 43+). Detailnější popis námi implementované ART metody můžete najít v článku [1].

V minulém cvičení na dopřednou Radonovu projekci je uvedena rovnice pro spojitou transformaci. Následující rovnice popisuje dopřednou diskretizovanou Radonovu transformaci

$$
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) ~~~~~~~~~ (2)
$$

kde funkce $\zeta$ je definována

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

kde vypočtené souřadnice $x′,y′$ v původním pro dané otočení $\theta$ a souřadnice $p, q$ jsou dány předpisem

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

a kde $q=ih$ pro $i \in \mathbb{Z}$ a $q$ je z intervalu $q \in [ -\frac{\sqrt{2}N}{2}; \frac{\sqrt{2}N}{2} ]$ (pro jednoduchost v našem případě budeme uvažovat $h = 1$) pro obrázek o velikosti $N \times N$ . Vztah mezi souřadnicemi $x,y$ v původním obrázku a otočenými souřadnicemi $p,q$ je dán předpisem (3).

===== Iterativní rekonstrukce - postup =====

Principem iterativní rekonstrukce je aplikování korekcí na libovolné počáteční hodnoty pixelů v rekonstruovaném obrázku tak, aby rozdíl Radonovy transformace rekonstruovaného obrazu a vstupního sinogramu byl minimální.

Zvolíme počáteční odhad rekonstruovaného obrázku $f_0(x, y) = 0$ o velikosti $N \times N$ pixelů.
V každém kroku iterativního algoritmu provedeme aktualizaci rekonstruovaných hodnot pro všechny $(p, \theta)$:

  - Spočítáme Radonovou transformaci pro všechny $\theta$ a posuny $p$ (viz rovnice (2)) $$
\hat{J}(p,\theta) = R [ \hat{f}(x,y) ]. ~~~~~~~~~ (5) $$
  - Rozdíl vypočtené $J^*(p,\theta)$ a změřené $J(p,\theta)$ Radonovy transformace obrázku $J(p,\theta)−\hat{J}(p,\theta)$ vydělíme váhovou funkci $w(p,\theta)$ udávající přes kolik pixelů na paprsku $(p,\theta)$ byla provedena projekce. Takto získáme průměrnou chybu pro každý pixel ležící na paprsku $(p, \theta)$.
  - Pro každé q provedeme aktualizaci hodnoty pixelů ležících na paprsku $(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)} ~~~~~~~~~ (6)
$$ kde váhová funkce $w$ pro jednotkový obrázek je definována $$
w(p,\theta) = R[1] (p,\theta) = N_q (p,\theta) \cdot h ~~~~~~~~~ (7)
$$ kde $N_q$ je počet platných q pro daný úhel a posunutí $(p, \theta)$. Pro potřeby efektivní implementace můžeme úlohu formulovat i následovně $$
\forall q: ~~~~ f(x,y) \leftarrow f(x,y) + \frac{J(p,\theta)−\hat{J}(p,\theta)}{w(p,\theta)} ~~~~~~~~~ (8)
$$


==== Krok 1 - Radonova transformace ====

Proveďte Radonovu transformaci rekonstrukce obrazu z minulé iterace (5).

==== Krok 2 - Váhová funkce $w$ ====

Při zpětném promítnutí chyby během iterativní rekonstrukce se celková chyba pro jeden paprsek distribuuje do všech pixelů ležících na tomto paprsku. Proto je třeba před odečtením chyby od všech pixelů chybu vydělit počtem pixelů, které na paprsku leží. Tento počet jednoduše zjistíme následujícím postupem: Provedeme Radonovou transformaci obrazu $f^*(x, y) = 1; \forall x, y$ obsahující samé 1. Váhová funkce $w(p,\theta)$ je pak rovna $R[f^*(x,y)]$.

Do zprávy vložte matici $w$ vykreslenou jako obrázek.

==== Krok 3 - Iterativní úprava rekonstruovaných hodnot ====

Naprogramujte funkci
''image=myARTiter(radonim,theta,nIter)''
která provede ''nIter'' kroků iterativní rekonstrukce měřených dat ''radonim'' $[m \times N_\theta]$ pro úhly ''theta'' $[1 \times N_\theta]$.
   - Inicializujte rekonstruovaný obrázek $f_0(x, y) = 0$.
   - Zvolte první $\theta$.
   - Spočítejte $J_k(p,\theta)$ a porovnejte se vstupními daty $J(p,\theta)$.
   - Upravte všechny body v rekonstruovaném obrázku ležící na paprsku $(p, \theta)$ podle rovnice (6)
   - Opakujte body 2-4 pro všechny $\theta$.
   - Iterujte stanovený počet iterací (viz zadání úkolu).

===== Reference =====

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

*/

===== 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/b212/_media/courses/zsl/ct2.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)^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)
$$

=== Learning rate ===

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.
===== 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 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]$.
   - Initialize the reconstructed image as $f_0(x, y) = 0$.
   - 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
   - Repeat steps 2-4 for all projection angles $\theta$.
   - 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).

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