Search
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.
imgCT
rawCT
phantom
phantom_noise1
head2
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 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).
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)$:
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]$.
image=myARTiter(radonim,theta,nIter)
nIter
radonim
theta
[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.