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⁰(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 - J^k(p,θ) = R[f^k(x,y)] - viz rovnice (1).
• Porovnáme J^k(p,θ) se skutečnou Radonovou transformací obrázku zadanou na vstupu J(p,θ).
• Rozdíl J(p,θ)-J^k(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 f^k+1(x,y) = f^k(x,y) + (J(p,θ)-J^k(p,θ))/w(p,θ)*pi/(2*N_θ).

• Zvolte si jedno pevné θ.
• Pro kazdé p= -N_p/2 : N_p/2 a q=-N_p/2 : N_p/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 N_pxN_p (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'-N_x/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!

• 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 f⁰(x,y) = 0.
• Zvolte jedno pevné θ.
• Spočítejte J^k(p,θ) a porovnejte se vstupními daty J(p,θ).
• Upravte vsechny body v rekonstruovaném obrázku lezící na paprsku (p,θ):   f^k+1(x,y) = f^k(x,y) + (J(p,θ)-J^k(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 N_y,N_x (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).