Radonova transformace je matematická transformace, který se používá například při počítačové tomografii.

Vstupní data pro úlohu jsou v archivu 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.

1. [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. [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í.
3. [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.
4. 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í).

Podívejte se znovu na rovnice popisující Radonovou transformaci (2). Projděte si též slidy z přednášky a prostudujte si Iterativní rekonstrukci (slide 77+). 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).

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)$:

1. 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) $$
2. 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)$.
3. 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) $$

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

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.

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]$.

1. Inicializujte rekonstruovaný obrázek $f_0(x, y) = 0$.
2. Zvolte první $\theta$.
3. Spočítejte $J_k(p,\theta)$ a porovnejte se vstupními daty $J(p,\theta)$.
4. Upravte všechny body v rekonstruovaném obrázku ležící na paprsku $(p, \theta)$ podle rovnice (6)
5. Opakujte body 2-4 pro všechny $\theta$.
6. Iterujte stanovený počet iterací (viz zadání úkolu).

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