Shepp-Loganův fantom |
Raw-data z X-Ray CT |
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. |
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:
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ě:
Ú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