====== 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}}