Na dnešním cvičení se budeme zabývat principem měření v počítačové tomografii (CT). CT pracuje při měření objektu (v kartézských souřadnicích <latex>x, y</latex>) pomocí sady projekcí pod různými úhly <latex>\theta</latex>. To odpovídá Radonově transformaci vstupního obrázku a dnešním úkolem je naprogramovat tuto transformaci pomocí vlastní funkce myRadon. Při rekonstrukci dat CT z měřeného signálu se pak používá inverzní Radonova transformace, kterou se budeme zabývat příště.

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, 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.]

Stručná ilustrace souřadných systémů pro danou projekci pod úhlem <latex>\theta</latex> je na následujícím obrázku.

Souřadná soustava <latex>(p, q)</latex> vznikne otočením soustavy <latex>(x, y)</latex> o úhel <latex>\theta</latex>. Pro transformaci tedy platí

<latex> \begin{equation*} p = x \cos(\theta) + y \sin(\theta) \qquad \qquad q = -x \sin(\theta) + y \cos(\theta) \end{equation*} </latex>

a pro inverzní transformaci pak platí

<latex> \begin{equation*} x = p \cos(\theta) - q \sin (\theta) \qquad \qquad y = p \sin(\theta) + q \cos(\theta) \end{equation*} </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. Je dobré nastavit souřadnice obrázku tak, aby bod [0, 0] ležel uprostřed obrázku
2. Určete minimální a maximální hodnoty pro <latex>p, q</latex> (1. krok zadání)
3. Zvolte diskretizaci pro <latex>p, q</latex> v rozsahu <latex> [\min p, \max p] </latex>, respektive <latex> [\min q, \max q] </latex>
4. Pro všechny úhly <latex>\theta</latex> a všechny diskretizační body <latex>p_h, q_h</latex> spočtěte <latex>J(\theta,p_h)</latex>
   1. pro uvažované body <latex>p_h, q_h</latex> spočtěte souřadnice <latex>x, y</latex> při otočení o <latex>\theta</latex>
   2. intenzitu obrázku v daném bodě, který nebude ve většině případů ležet přímo ve středu pixelu, dostaneme interpolací (funkce interp2
   3. spočtěte numerickou hodnotu <latex>J(\theta,p_h)</latex>