====== 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 [[http://en.wikipedia.org/wiki/X-ray_computed_tomography|CT]] a o samotné [[http://en.wikipedia.org/wiki/Radon_transform|Radonově transformci]].
===== Stručné zadání =====
- Stáhněte si {{courses:a6m33zsl:ct_bpdata.zip|archív}} s daty a kódy pro toto cvičení.
* //phantom.mat// obsahuje výsledný obrázek fantomu //phantom// o velikosti [m\times n]
* //radonim.mat// obsahuje Radonovu projekci fantomu //imgCT// o velikosti [2p_{max}\times \theta_{max}]
* //designFilter.m// je funkce pro navržení/vytvoření filtru
- Naprogramujte funkci //{{courses:a6m33zsl:myiradon.m?linkonly|image=myIRadon(radonim,theta,filter)}}// která provede zpětnou projekci, nebo filtrovanou zpětnou projekci Radonovy matice //radonim// [m\times n] pro úhly //theta// [1\times n]. 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 d \in\langle 0,1\rangle je zlomová frekvence filtru. \\ [3b.]
- Radonův obraz byl získán pro úhly 0,1\ldots 179^\circ 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 d=\{1,0.9,0.7\}. 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.]
- 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.]
- **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 [[http://www.mathworks.com/help/distcomp/exist.html|exist]]) funkce vrátí pouhou zpětnou projekci bez filtrace. Úhly zadávejte ve stupních. Dodržte rozhraní funkce dané šablonou. //
===== Podrobný popis =====
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 J(\theta,p).
Naivním způsobem jak rekonstruovat obrázek f(x,y) je zpětná projekce. Přibližná rekonstrukce \hat{f}(x,y) se vypočítá jako
\hat{f}(x,y)=\int_0^\pi J(\theta,p) d \theta
vztah mezi souřadnými soustavami je v následujících rovnicích.
p=x\cos\theta+y\sin\theta ~~~~;~~~~ q=-x\sin\theta+y\cos\theta
x=p\cos\theta-q\sin\theta ~~~~~;~~~ y=p\sin\theta+q\cos\theta
Přesnou rekonstrukci, tedy inverzní Radonovu transformaci lze provézt pomocí filtrované zpětné projekce. Označme S(\theta,\omega) 1D Fourierovu transformaci Radonova prostoru J(\theta,p) po sloupcích, tedy přes proměnnou //p// tak že
S(\theta,\omega)=\int_{-\infty}^{\infty}J(\theta,p)e^{-2j\pi\omega p} dp
Rekonstrukci //f(x,y)// vypočítáme jako zpětnou projekci filtrovaného Radonova obrazu J^*(\theta,p)
f(x,y) =\int_0^\pi J^*(\theta,p) d \theta,
kde J^*(\theta,p) je J_\theta(p) 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 \mathcal{F} je symbol pro Fourierovu transformaci a \ast 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
\omega_h(n) =0.54-0.46\cos{\left(\frac{2\pi n}{N-1}\right)}
nebo Shepp-Loganovo (sinc) okno
\omega_s(n) =\frac{\sin{\left(\frac{n}{N}\right)}}{\frac{n}{N}}
kde //N// je délka filtru.
===== Možný postup =====
- Nastavme souřadnou soustavu obrázku tak, aby bod //[0,0]// ležel uprostřed.
- Odhadněte velikost rekonstruovaného obrázku.
- Obdobně jako v předcházejícím cvičení inicializujte prázdný obrázek f(x,y)=0.
- Skládejte dohromady jednotlivé projekce //p// pro zvolené úhly \theta
* Pro zvolené \theta, vypočítejte souřadnici //p// všech bodů v obrázku p=x\cos{(\theta)}+y\sin{(\theta}).
* V příslušném sloupci J(\theta,p) 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, \theta je celé číslo, použijte funkci //[[http://www.mathworks.com/help/matlab/ref/interp1.html|interp1]]//. Druhou možností (rychlejší, ale méně přesnou) je místo interpolace použití zaokrouhlování //[[http://www.mathworks.com/help/matlab/ref/round.html|round]]//.
* Nalezené hodnoty přičtěte k f(x,y)
* Opakujte předešlé tři úkoly z tohoto bodu.
- 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.
- Relativní rekonstrukční chybu počítejte jako 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} }
{{tag>projection ct medical filter radon}}