Warning
This page is located in archive.

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.

Stručné zadání

  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 [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
  2. Naprogramujte funkci 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.]
  3. 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.]
  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.

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

  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  f(x,y)=0.
  4. 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 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  f(x,y)
    • 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  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} }
courses/a6m33zsl/lab05_ct_bwd.txt · Last modified: 2018/03/17 10:55 by herinjan