Zpetna projekce
(cviceni pro ZSL2)
Projděte si rovnice popisující Radonovou transformaci, zpětnou projekci a filtrovanou zpětnou projekci - jedná se o vybrané rovnice prezentované na přednáskách.

Stáhněte si matici radon.mat, obsahující Radonovu transformaci Shepp-Loganova fantomu (viz. obrázek napravo). Takto vypadají raw-data z X-Ray CT a Vasím úkolem bude zrekonstruovat z těchto dat původní obrázek. Na obrázek z radon.mat se budeme odkazovat pod názvem imgCT.

Shepp-Loganův fantom

Raw-data z X-Ray CT
 
Úkol 1 - Zpetna projekce (Backprojection)

Zpětná projekce je popsána rovnicí (Eq. 2). Napistě program, který numericky spočítá přibliznou rekonstrukci obrazu fantomu - tedy integrál z Jθ (p) podle θ. Jako vstupní data pouzijte obraz imgCT - projekce obrazu fantomu (Jθ(p)) pro θ je 0:Nθ stupňů a p= -NP:NP s krokem odpovídajícím jednomu pixelu v originálním obrázku (NP = 183 a Nθ = 179). Pripominame ze p zavisi na x,y,θ (Eq. 1).

Podukol: Nejdříve spočítejte přibliznou velikost rekonstruovaného obrázku NxN (N spočítáte z NP).

Centrum souřadného systému [p,θ] i systému [x,y] je ve středu rekonstruovaného obrázku. Pro θ=0 je osa p shodná s osou x. Souřadná soutava [x',y'] (odpovídající souřadné soustavě obrázku v Matlabu) je tudíz posunutá oproti [x,y] o 'půl obrázku'.

Těd musime pro pevně zvolenéθ spočítat souřadnici p pro kazdý bod [x',y'] v rekonstruovaném obrázku (viz. Eq.1 - převod souřadnic). P = XB cos θ + YB sin θ. Muzeme to udelat i bez pouziti for-cyklu (navod je v dodatku na konci teto stranky).

Krok 1

Nasím cílem je provést zpětnou projekci z Jθ(p), kde jsou data ulozena v polárních souřadnicích. Proto je třeba před samotnou zpětnou transformací provést převod zpět do kartézských souřadnic.

Krok 1.2:Sestrojte vektory x = (1:N) a y = (N:-1:1)' a upravte je (přičtením vhodných hodnot) tak, aby pro oba vektory x a y index vektoru odpovídal souřadnici v [x',y'] a hodnota na příslusném místě vektoru odpovídala souřadnici v [x,y].

Krok  1.3: Vytvořte matice, kde indexy v matici budou odpovídat souřadnicím v [x',y'] a hodnoty v matici x a y budou odpovídat souřadnicím v [x,y].
  xb = repmat(x,N,1);
  yb = repmat(y,1,N);
Krok 1.4:Zvolte jedno konkrétní θ (napr. θ = π/3) a pro vsechny souřadnice [x',y'] v rekonstruovaném obrázku spočtěte souřadnici p - vyuzijte výse popsaný postup(P = XB cos θ + YB sin θ). Do reportu vlozte matici indexu P pro zvolene θ vykreslenou pomoci 'imagesc' (jako zmenseny obrazek!) spolu se souradnicami.

Souřadný systém v orig.obrázku:


Krok 2
 Při zpětné projekci vezmeme kazdý z paprsků [p,θ] a promítneme ho zpět do rekonstruovaného obrázku. Vsechny tyto projekce sčítáme. Upozornění: Tímto postupem dostaneme pouze přiblizný obraz původního objektu, v zádném případě se nejedná o inverzní Radonovou transformaci!!

Krok 2.1: Ve hlavnim for-cyklu budeme prochazet vsetky uhly θ=0:179. Pro konkrétní θ vybereme jeden stloupec z matice imgCT(:,θ), a spočteme souřadnici P (podle kroku 1 bez for-cyklu) p ro kazdý bod [x',y'] . Dostaneme matici P = XB cos θ + YB sin θ.

Krok 2.2: Hodnoty v matici P jsou souřadnice do vstupní matice imgCT(:,θ ). Nejsou to vsak celá čísla, nelze pomocí nich přímo indexovat v matici reprezentující radonovu transformaci. Proto je nutné pouzít interpolaci! (Matice imgCT je obrázek po radonové transformaci, který jste dostali na začátku.) Pro interpolaci volime jednu z metod 'nearest', 'linear', 'cubic'.

  osaP = (1:size(imgCT,1)) - ceil(size(imgCT,1)/2);
  projContrib_linarr = interp1(osaP,imgCT(:,θ),P(:),'metoda',0);
  projContrib = reshape(projContrib,N,N);
Krok 2.3: Takto dostaneme obrázek, který je zpětnou projekci pro jedno θ. Celý postup opakujeme pro vsechny θ od 0 do 179 a výstupy sčítáme. Výsledný obraz je zpětnou projekcí Radonové transformace obrazu.

Pozor! Můze se stát, ze pro některé hodnoty [x',y'] dostaneme hodnoty p, které jsou mimo rozsah matice imgCT. Pokud zadáme 0 jako poslední parametr funkce interp1, vsechny hodnoty P mimo interval osaP se naplní 0.

Krok 2.4: Po rekonstrukci je třeba normalizovat rekonstruovaný obrázek - img = img*pi/(2*Nθ).

Ukol: Do reportu vlozte vysledek zpetne rekonstrukce obrazku!

Rekonstrukce nebude přílis kvalitni, ale bude mozné v ni původní obraz rozeznat. Pro kvalitnějsí rekonstrukci pouzijeme flitrovanu zpetnou projekci, viz. nasledujici ukol.


Úkol 2 - Filterovana zpetna projekce (Filtered backprojection)

S vyuzítím kódu pro zpětnou projekci implementujte filtrovanou zpětnou projekci (Eq. 5).

Krok 1: Nejdrive J*θ vytvoříte z Jθ pomoci filtrace ve frekvencnim spektru: F-1{|ω| F{Jθ(p)}} - viz. posledni rovnice (Eq. 6). Budeme delat 1D filtraci pro kazdy stloupec matice Jθ (p) (= imgCT).  Pro získání funkce |ω| (filtru ve frekvenční oblasti) pouzijte funkci designFilter(filter,len,d) , která vytvoří filtr pomocí kterého získáte J*θ z Jθ.

  H = designFilter(filter, len, d);

Parametry funkce designFilter:

  • filter - typ filtru - 'ram-lak', 'shepp-logan', 'hamming'
  • len - délka (mělo by to být sudé číslo)
  • d - číslo z intervalu [0,1] - specifikuje, jaké procento nejnizsích frekvencí bude filtrem pouzito při filtraci. 0 - nepouzije se nic. 1 - pouzijí se vsechny frekvence az po Nyquistovu frekvenci.

Upozornění - funkce designFilter vytvoří filtr sudé délky. Pokud má Jθ (=imgCT) lichou délku, je potřeba před Fourierovou transformací Jθ(p) přidat řádek nul a posléze ho odebrat .

Krok 2: Filtr pouzijte ve frekvenční oblasti, tak ze vynásobite Fourierovou transformací Jθ s filtrem H zvlásť pro kazdé p. Pouzijte funkci fft.

  fft(imgCT(:,p)) .* H;

Krok 3: Po zpětné Fourierově transformaci vezměte pouze reálnou slozku. Pouzijte funkci iffta potom funkci real.

Nezapomente odebrat posledni radek z vysledku (viz upozorneni).

Krok 4: Provedte zpetnou rekonstrukci na filtrovane matici J*θ (postup z ukolu 1).

Ukol: Do reportu vlozte vysledek filtrovany zpetne rekonstrukce obrazku!


Úkol 3 - Posuďte kvalitu rekonstrukce

Řesení tohoto úkolu musí být ve zprávě:

  • Porovnejte kvalitu obrázků zrekonstruovaných pomocí zpětné projekce a filtrované zpětné projekce.
  • Spočítejte rekonstrukční chybu jako sumu čtverců rozdílu mezi zrekonstruovaným I'(x,y) a původním obrázkem I(x,y).
    (Error(I(x,y), I'(x,y)) = suma( (I(m,n) - I'(m,n)).^2 ), pre [m,n]=vsetky pixle obrazku)
  • Porovnejte (kvalitativně i kvantitativně) kvalitu rekonstrukce při pouzití různých filtru ('ram-lak', 'shep-logan', 'hamming') a ruzné hodnoty parametru d (1, 0.9 a 0.7). Zvolte jednu nejlepsi hodnotu pro kazdy parametr.

Úkol 4 - Bonusový úkol + 0.5 bodu

Rekonstrujte obrázek s méně hustým vzorkováním úhlu otočení. Radonová transformace byla provedena pro úhly 0 az 179 s krokem 1 stupeň. Pouzijte pro rekonstrukci pouze kazdý 2, 5, 10 a 20 vzorek. Zhodnoťte kvalitu takto zrekonstruovaných obrázků.

Dodatek - vypocet souradnic bez for-cyklu


Následující matice XA a YA na pozici [x,y] obsahují x-ové (v matici XA) a y-ové (YA) souřadnice příslusného bodu v matici:
  XA = repmat(1:5,N,1);
  YA = repmat((1:5)',1,N);
Pokud má souřadná soustava střed ve středu obrázku jsou x-ové a y-ové souřadnice:
  XB = repmat((1:5)-3,N,1);
  YB = repmat((1:5)'-3,1,N);

Pomocí funkce arctg můzeme pro kazdý bod lehce spočítat úhel, který svírá vektor  ([x,y], [0,0]) s osou x. Toto lze vyuzít při transformaci do polárních souřadnic. TOTO JE JEN PŘIKLAD POUZITÍ - NEMUSÍTE POČÍTAT ZÁDNÝ UHEL!!! V nasem případě budeme počítat souřadnici p - viz. níze.
coords = atan(YB./XB) nebo coords = atan2(YB,XB)
 		coords = 
[135
 117
 90
 63
 45;
153
 135
 90
 45
 27;
180
 180
 0
 0
 0;
-153
 -135
 -90
 -45
 -27
-135
 -117
 -90
 63
 -45]

Posledni uprava: 30.3.2009