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

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

1. Stáhněte si obrázky z CT a jejich Radonovy obrazy (sinogramy). Dále si stáhněte sinogramy fantomu s přidaným šumem.
2. 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.
3. 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.
4. 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.
5. 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ů.

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

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é $\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)$.
3. 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) $$

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 $\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$.
4. 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)$
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.