Počítačová tomografie (CT) - přímá projekce

Radonova transformace je matematická transformace, který se používá například při počítačové tomografii.

1. Najděte analyticky či graficky maximální absolutní hodnotu <latex>p_m = \max |p|</latex> a <latex>q_m = \max |q|</latex>, kde <latex>p=x\cos\theta+y\sin\theta</latex>, <latex>q=-x\sin\theta+y\cos\theta</latex> a <latex>(x,y)\in [-x_m,x_m] \times [ -y_m, y_m]</latex>. [0.5b.]
2. Diskretizujte rovnici (1). [0.5b.]
3. Naprogramujte funkci imgRadon=myRadon(img,theta) která vrátí přímou Radonovu transformaci imgRadon <latex>[(2p_m+1) \times k]</latex> obrázku image <latex>[m\times n]</latex> pro úhly theta <latex>[1\times k]</latex>. Úhly zadávejte ve stupních. Dodržte šablonu funkce. [3b.]
4. Pomocí vámi napsané funkce transformujte Shepp-Loganův fantom, který si můžete vytvořit funkcí phantom(vel), parametrem funkce je velikost fantomu. Do reportu vložte výsledek spočítaný na fantomu s velikostí 256px pro úhly <latex>0,1\ldots 179^\circ</latex>. Nezapomeňte na správná měřítka a popisy os. [1b.]

Souřadná soustava (p,q) vznikla pootočením souřadné soustavu (x,y) o zvolený úhel <latex>\theta</latex>

<latex> p=x\cos\theta+y\sin\theta ~~~~;~~~~ q=-x\sin\theta+y\cos\theta </latex> (2)

a zpětná transformace je tedy následující

<latex> x=p\cos\theta-q\sin\theta ~~~~;~~~~ y=p\sin\theta+q\cos\theta </latex> (3)

Úhel <latex>\theta</latex> určuje otočení projekce proti ose x a parametr p je vzdálenost projekční přímky od počátku. Vztah soustav x,y a p,q při úhlu <latex>\theta</latex> ilustruje obrázek. Radonova transformace převádí obraz v kartézských souřadnicích x a y do Radonova prostoru o souřadnicích p a <latex>\theta</latex>. Výsledná Radonova projekce ze souřadného systému obrázku (x,y) do Radonova prostoru <latex>(\theta,p)</latex> je dána integrálem

<latex> J(\theta,p) = \int_{-\infty}^\infty f\left ( p\cdot\cos(\theta) - q\cdot\sin(\theta) ~ , ~ p\cdot\sin(\theta) + q\cdot\cos(\theta) \right ) dq </latex> (1)

1. Nastavte si souřadnou soustavu obrázku tak, aby bod [0,0] ležel uprostřed.
2. Ze vztahů odvozených v 1. bodě zadání vypočítejte maximální a minimální hodnoty p a q z transformací (x,y) na (p,q) podle rovnice (2).
3. Zvolte celočíselnou diskretizaci p v rozsahu 'max p', 'min p', a q v rozsahu 'max q', 'min q'..
4. Pro všechny p, q a <latex>\theta</latex> vypočtěte <latex>J(\theta,p)</latex> podle rovnice (1).
   • Pro body vyjádřené v souřadnicích p, q s otočením <latex>\theta</latex> vypočítejte souřadnice x, y.
   • Pomocí interpolace (interp2) najděte intenzity v těchto bodech.
   • Numericky vypočítejte <latex>J(\theta,p)</latex> (integrál diskretizujme).

Poznámka: čím méně for cyklů použijete, tím rychleji vám výpočet poběží.

Doporučujeme odladit kód na malém fantomu, např. <latex>64\times 64</latex>px a teprve pak jej pustit na data <latex>256\times 256</latex>px. Výsledek svého snažení porovnejte s funkcí radon, která je standardně implementovaná v Matlabu a do zprávy vložte rozdílový obrázek.