Iterativní
algebraická rekonstrukce
Podívejte se znovu na
rovnice
popisující Radonovou transformaci (1) a zpětnou projekci (2).
Projděte si téz slidy z
přednásky a prostudujte si Iterativní rekonstrukci (slide 77+).
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 obsahující obrazovou funkci
J(p,θ) se budeme odkazovat pod názvem imgCT.
Shepp-Loganův
fantom
Raw-data z X-Ray
CT
Iterativní
rekonstrukce
Principem iterativní rekonstrukce je aplikování korekcí na libovolné
počáteční hodnoty pixelů v rekonstruovaném obrázku tak, aby po
provedení Radonové transformace na rekonstruovaný obrázek bylo dosazeno shody
se vstupními daty.
Zvolíme počáteční odhad rekonstruovaného obrázku f
0(x,y) =
0 o velikosti NxN pixelů.
Pro kazdé θ postupně provedeme aktualizaci rekonstruovaných hodnot:
-
Spočítáme Radonovou transformaci pro θ a pro vsechny posuny p - Jk(p,θ)
= R[fk(x,y)] - viz rovnice (1).
-
Porovnáme Jk(p,θ) se skutečnou Radonovou transformací
obrázku zadanou na vstupu J(p,θ).
-
Rozdíl J(p,θ)-Jk(p,θ) vydělíme váhovou funkci
w(p,θ) udávající přes kolik pixelů na paprsku (p,θ) byla
provedena projekce. Takto získáme přibliznou chybu pro kazdý pixel lezící
na paprsku (p,θ).
-
Pro kazdý pixel lezící na paprsku (p,θ) provedeme aktualizaci jeho hodnoty
fk+1(x,y) = fk(x,y) + (J(p,θ)-Jk(p,θ))/w(p,θ)*pi/(2*Nθ).
Krok 1 -
Radonova transformace
Stáhněte si obrázek fantomu.
Proveďte na něm Radonovou transformaci. Pro kontrolu porovnejte
výstup s radon.mat.
Radonova transformace:
-
Zvolte si jedno pevné θ.
-
Pro kazdé p= -Np/2 : Np/2 a q=-Np/2 : Np/2
spočítejte jejich souřadnice v souřadné soustavě snímaného
objektu [x,y]:
-
Xcor = p cos θ - q sin θ
-
Ycor = p sin θ + q cos θ
-
Poznamka: promenne p a q jsou dvou-rozmerne matice souradnic o velikosti NpxNp
(p se meni v jedne dimenzi a q se meni v druhe dimenzi). Predchozi dve
rovnice spocitaji souradnice XY do obrazku (nektere muzou byt mimo rozsah
obrazku). V nasledovnim kroku udelame projekci podle q a tedy vypocteme
Radonovu transformaci J(p,θ) pro jeden uhel θ a vsechny hodnoty p.
-
Pro kazdé p spočtěte sumu podle q: J(p,θ) = Σq
f(p,q,θ) = Σq f(x,y) = Σq f(x',y'), kde x
= x'-Nx/2. Je nutné pouzít interpolaci a při korekci chyb je
nutné pouzít stejnou interpolaci (rozdělit chybu mezi vsechny body f(x,y)
ze kterých se f(x,y) pro necelé x,y dointerpolovala). Pro interpolaci pouzijte
funkci interp2(Xax,Yax,img,Xcor,YCor,method,0), kde Xax je repmat((1:Nx) -
Nx/2,Nx,1), Yax dostanete obdobně, img je vstupní obrázek a XCor, YCor
jsou souřadnice, na kterých chcete znát hodnotu obrazové funkce img (viz
souřadnice spočítané v předchozím kroku).
-
Opakujte postup pro vsechna θ.
-
J(p,θ) je Radonova transformace obrazu f(x',y').
Uploadujte soubor radon.m, do zpravy nevkladejte zdrojovy
text!
Krok 2 - Váhová funkce w
Při zpětném promítnutí chyby
během iterativní rekonstrukce se celková chyba pro jeden paprsek
distribuuje do vsech pixelů lezících na tomto paprsku. Proto je třeba
před odečtením chyby od vsech pixelů chybu vydělit
počtem pixelů, které na paprsku
lezí. Tento
počet jednoduse zjistíme následujícím postupem:
-
Provedeme Radonovou transformaci obrazu f(x',y') obsahující samé 1. w(p,θ)
je pak rovna J(p,θ).
Do zpravy vlozte matici w vykreslenou jako obrazek!
Krok 3 -
Iterativní úprava rekonstruovaných hodnot
Oprava hodnot bude probíhat v několika
iteracích.
-
Inicializujte rekonstruovaný obrázek f0(x,y) = 0.
-
Zvolte jedno pevné θ.
-
Spočítejte Jk(p,θ) a porovnejte se vstupními daty
J(p,θ).
-
Upravte vsechny body v rekonstruovaném obrázku lezící na paprsku (p,θ):
fk+1(x,y) = fk(x,y) + (J(p,θ)-Jk(p,θ))/w(p,θ)*pi/(2*Nθ).
-
Souřadnice bodů [x',y'] máte z předcházejícího kroku. Vzhledem k
tomu, ze souřadnice nemusí být celá čísla, je potřeba pouzít
interpolaci (v kroku 1). V kroku 3 při zápisu je potom potřeba
upravit hodnoty ve vsech pixelech, ze kterých se brala data.
-
Nejprve implementujte program, který bude pouzívat interpolaci pomocí
nejblizsího souseda a v kroku 3 upravte pixel, který vznikne zaokrouhlením
příslusných souřadnic [x',y'].
-
Opakujte pro vsechny θ.
-
Iterujte dokud se bude zmensovat celková chyba.
Můzete pouzít funkci sub2indX(siz,Y,X), která převede vektory
souřadnic Y,X odpovídající y-ové a x-ové souřadnici v matici
velikosti Ny,Nx (siz = [Ny,Nx]) na tzv. lineární index -
pomocí kterého můzete přímo indexovat matici. Funkce sub2indX vyuzívá
fci sub2ind s tím, ze z vektorů souřadnic automaticky odstraní
souřadnice odkazující se na hodnoty mimo rozměr matice.
Uploadujte soubor
iterativniRekonstrukce.m, do zpravy nevkladejte zdrojovy text!
Úkol
Zrekonstrujte obrázek fantomu (pouzijte
bikubickou interpolaci pro Radonovu transformaci). Porovnejte výsledek s
výsledkem z filtrované zpětné projekce (z minulého cvičení) -
spočíteje chybu oproti originálnímu obrázku fantomu.
Proveďte rekonstrukci pro zasuměné vstupní
obrázky a opět porovnejte s výsledky z filtrované zpětné
projekce.
Do zpravy vlozte rekonstruovany obrazek po 1, 10 a 100
iteracich.
Zaznamenajte chybu rekonstrukce ciselne od
poctu iteraci. Pouzijte SSD kriterium (suma kvadratickych odchylek od
originalniho obrazku).
Porovnejte
vysledek zpetne rekonstrukce a algebraicke rekonstrukce (subjektivne, casove
naroky, a pomoci SSD).