Cvičení 5: Filtrovaná zpětná projekce (CT)

Zpětná Radonova projekce se využívá k rekonstrukci dat naměřených z CT, více si přečtěte na wiki CT a o samotné Radonově transformci.

1. Stáhněte si archív s daty a kódy pro toto cvičení.
   • phantom.mat obsahuje výsledný obrázek fantomu phantom o velikosti <latex>[m\times n]</latex>
   • radonim.mat obsahuje Radonovu projekci fantomu imgCT o velikosti <latex>[2p_{max}\times \theta_{max}]</latex>
   • designFilter.m je funkce pro navržení/vytvoření filtru
2. Naprogramujte funkci image=myIRadon(radonim,theta,filter) která provede zpětnou projekci, nebo filtrovanou zpětnou projekci Radonovy matice radonim <latex>[m\times n]</latex> pro úhly theta <latex>[1\times n]</latex>. Parametr filter je struktura obsahující parametry filtru. V poli
   • .type nechť je název filtru ('ram-lak', 'shepp-logan', 'hamming'),
   • .len je délka filtru,
   • .d kde <latex>d \in\langle 0,1\rangle</latex> je zlomová frekvence filtru.
     [3b.]
3. Radonův obraz byl získán pro úhly <latex>0,1\ldots 179^\circ</latex> a celá p s krokem 1. Pomocí vámi napsané funkce rekonstruujte data Radonova obrazu ze souboru radonim.mat pro všechny tři typy filtrů a hodnoty parametru <latex> d=\{1,0.9,0.7\}</latex>. Dále proveďte i zpětnou projekci bez filtru (celkem dohromady deset rekonstrukcí). Výsledné obrázky vložte do zprávy a popište. [1b.]
4. Spočítejte rekonstrukční chybu pro všechny rekonstrukce jako sumu čtverců rozdílu původního phantom.mat a rekonstruovaného obrazu a vyberte nejlepší rekonstrukci. Do zprávy ke každé vaší rekonstrukci vložte i vizualizaci rekonstrukční chyby. [1b.]
5. Bonus: Rekonstruujte obrázek s méně hustým vzorkováním úhlu. Použijte postupně každý 2, 5, 10, a 20. vzorek (další 4 rekonstrukce) Výsledky spolu s krátkým hodnocením vložte do zprávy. [+1b.]

Poznámka: Pokud nebude poslední parametr zadán (použijte funkci exist) funkce vrátí pouhou zpětnou projekci bez filtrace. Úhly zadávejte ve stupních. Dodržte rozhraní funkce dané šablonou.

V tomto cvičení se seznámíte s rekonstrukcí obrazu v počítačové tomografii, zpětnou projekcí a filtrovanou zpětnou projekcí. Radonův obraz funkce (obrázku) f(x,y) budeme značit <latex>J(\theta,p)</latex>.

Naivním způsobem jak rekonstruovat obrázek <latex> f(x,y)</latex> je zpětná projekce. Přibližná rekonstrukce <latex> \hat{f}(x,y)</latex> se vypočítá jako <latex> \hat{f}(x,y)=\int_0^\pi J(\theta,p) d \theta </latex> vztah mezi souřadnými soustavami je v následujících rovnicích.

<latex> p=x\cos\theta+y\sin\theta ~~~~;~~~~ q=-x\sin\theta+y\cos\theta</latex>

<latex> x=p\cos\theta-q\sin\theta ~~~~~;~~~ y=p\sin\theta+q\cos\theta</latex>

Přesnou rekonstrukci, tedy inverzní Radonovu transformaci lze provézt pomocí filtrované zpětné projekce. Označme <latex>S(\theta,\omega)</latex> 1D Fourierovu transformaci Radonova prostoru <latex> J(\theta,p)</latex> po sloupcích, tedy přes proměnnou p tak že

<latex> S(\theta,\omega)=\int_{-\infty}^{\infty}J(\theta,p)e^{-2j\pi\omega p} dp </latex>

Rekonstrukci f(x,y) vypočítáme jako zpětnou projekci filtrovaného Radonova obrazu <latex>J^*(\theta,p)</latex>

<latex> f(x,y) =\int_0^\pi J^*(\theta,p) d \theta</latex>, kde <latex>J^*(\theta,p)</latex> je <latex>J_\theta(p)</latex> po filtraci takzvaným “ramp” filtrem $H = |\omega|$, $h=\mathcal{F}^{-1}\{|\omega|\}$:

$J^*(\theta,p)=\mathcal{F}^{-1}\{|\omega|\mathcal{F}\{J_\theta(p)\}\} = \mathcal{F}^{-1}\{|\omega|\}*J_\theta(p)$,

kde <latex>\mathcal{F}</latex> je symbol pro Fourierovu transformaci a <latex>\ast</latex> operátor konvoluce, zde podle 'p'.

Ramp filtr sice teoreticky vede na ideální rekonstrukci, ale má problematické praktické vlastnosti. Protože se jeho zesílení zvyšuje s frekvencí, zesiluje šum. Proto se v praxi používá v kombinaci s okny, která omezí jeho frekvenční rozsah, například Hammingovo okno <latex> \omega_h(n) =0.54-0.46\cos{\left(\frac{2\pi n}{N-1}\right)} </latex> nebo Shepp-Loganovo (sinc) okno <latex> \omega_s(n) =\frac{\sin{\left(\frac{n}{N}\right)}}{\frac{n}{N}} </latex> kde N je délka filtru.

1. Nastavme souřadnou soustavu obrázku tak, aby bod [0,0] ležel uprostřed.
2. Odhadněte velikost rekonstruovaného obrázku.
3. Obdobně jako v předcházejícím cvičení inicializujte prázdný obrázek <latex> f(x,y)=0</latex>.
4. Skládejte dohromady jednotlivé projekce p pro zvolené úhly <latex>\theta</latex>
   • Pro zvolené <latex>\theta</latex>, vypočítejte souřadnici p všech bodů v obrázku <latex>p=x\cos{(\theta)}+y\sin{(\theta})</latex>.
   • V příslušném sloupci <latex>J(\theta,p)</latex> nalezněte hodnoty pro všechny body v obrázku. Protože hodnoty p z předchozího kroku nejsou celočíselné je nutné interpolovat, tentokrát jen v jednom rozměru, <latex>\theta</latex> je celé číslo, použijte funkci interp1. Druhou možností (rychlejší, ale méně přesnou) je místo interpolace použití zaokrouhlování round.
   • Nalezené hodnoty přičtěte k <latex> f(x,y)</latex>
   • Opakujte předešlé tři úkoly z tohoto bodu.
5. Přidejte do kódu filtraci.
   • Pomocí funkce H=designFilter(typ,len,d) vytvořte požadovaný filtr. Máte na výběr typ
     • 'ram-lak' - filtr bez okna,
     • 'shepp-logan' - filtr násobený sinc oknem a
     • 'hamming' - filtr násobený Hammingovým oknem.
   • Filtrujte Radonův obraz fantomu ve frekvenční oblasti (po sloupcích). Funkce designFilter vrací filtry sudé délky, proto je třeba před filtrací přidat řádek k datům s lichým počtem řádků.
   • Po zpětné Fourierově transformaci pracujte jen s reálnou částí signálu a nezapomeňte odebrat přidaný řádek.
6. Relativní rekonstrukční chybu počítejte jako <latex> R=\sqrt{ \frac{1}{\|x\| \cdot \|y\|} \cdot \frac{ \sum_{\forall x, \forall y}(\hat{f}(x,y)-f(x,y))^2}{\sum_{\forall x, \forall y} f(x,y)^2} }</latex>