====== Počítačová tomografie (CT) - Iterativní algebraická rekonstrukce ======

[[http://en.wikipedia.org/wiki/Radon_transform|Radonova transformace]] je matematická transformace, který se používá například při počítačové tomografii.

===== Úkol =====

  - [[http://cmp.felk.cvut.cz/cmp/courses/ZSL2/cviceni/ct2/CTs.zip|Stáhněte si obrázky]] z CT a [[http://cmp.felk.cvut.cz/cmp/courses/ZSL2/cviceni/ct2/RadonCTs.zip|jejich Radonovy obrazy]] (sinogramy). Dále si [[ http://cmp.felk.cvut.cz/cmp/courses/ZSL2/cviceni/ct2/radonNoise.mat|stáhněte sinogramy]] fantomu s přidaným šumem.
  - Rekonstruujte obrazy ze všech dodaných (nezašuměných i zašuměných) sinogramů pomocí iterativní algebraické rekontrukční metody. Do zprávy vložte rekonstruované obrázky po 1, 10 a 100 iteracích.
  - Spočítejte 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.
  - Rekonstruujte všechny obrazy také pomocí 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: Porovnejte závislost $SSD$ chyby obou metod (filtrované zpětné projekce a algebraické rekonstrukce) na amplitudě přidaného aditivního Gaussovského šumu a na počtu projekčních úhlů.


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

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 'Algebraická rekonstrukce' a dále). 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 vztah vypočtených souřadnic $x′,y′$ v původním obrázku 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$ .


===== 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é $\hat{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é $x,y$ provedeme aktualizaci hodnot pixelů ležících na paprsku $(p, \theta)$, kde dva parametry jsou dány a jeden je volitelný $$
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) ~~~~~~~~~ (7)
$$ Pro potřeby efektivní implementace můžeme úlohu formulovat i následovně $$
\forall p,q,\theta: ~~~~ f(x,y) \leftarrow f(x,y) + \zeta (x,y,q,p,\theta) \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. */

==== 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 $\hat{J}(p,\theta)$ a porovnejte se vstupními daty $J(p,\theta)$.
      * Use the matlab function ''interp2(''$x$'', ''$y$'', ''$f(x,y)$'', ''$x'$'', ''$y'$'', 'nearest', 0)'' to calculate $J(p,\theta)$.\\ $x$ is an $N \times N$ matrix where each row is $\{-\frac{N}{2}+1, \cdots, \frac{N}{2} \}$.\\ $y$ is an $N \times N$ matrix where each column is $\{-\frac{N}{2}+1, \cdots, \frac{N}{2} \}$. \\ $x' = p\cos\theta - q\sin\theta$ is an $m \times m$ matrix, $p$ is an $m \times m$ where each row is $\{ $''ceil''$ (-\frac{m}{2}), \cdots, $''ceil''$ (\frac{m}{2}) \}$ and $q$ is an $m \times m$ where each column is $\{ $''ceil''$ (-\frac{m}{2}), \cdots, $''ceil''$ (\frac{m}{2}) \}$. \\ $y' = p\sin\theta + q\cos\theta$ is an $m \times m$ matrix. \\ $f(x,y)$ is the image being reconstructed.
      * ''J = sum( interp2(''$x$'', ''$y$'', ''$f(x,y)$'', ''$x'$'', ''$y'$'', 'nearest', 0), 1)'' will return a vector for the angle $\theta$ where each element corresponds to a different $p$.
   - Upravte všechny body v rekonstruovaném obrázku ležící na paprsku $(p, \theta)$ podle rovnice (8)
      * Obtain a vector $F$ by using the function ''interp1''$(\{$''ceil''$ (-\frac{m}{2}), \cdots, $''ceil''$ (\frac{m}{2})\},\frac{J - \hat{J}}{w(p,\theta)}$,''P(:)'','nearest'$,0)$. \\ where ''P''$~ = x\cos\theta + y\sin\theta$.
      * Obtain a matrix by using the function ''reshape''$(F,N,N)$.
      * Update the reconstructed image $f(x,y) = f(x,y) + \frac{\pi}{2 N_\theta} $''reshape''$(F,N,N)$
   - 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.



{{tag>medical CT projection Radon filter}}