===== Cvičení 4 : Počítačová tomografie (CT), přímá projekce =====

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á [[http://en.wikipedia.org/wiki/Radon_transform|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ě.

==== Zadání ====

  - 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.]
  - Diskretizujte rovnici (1). [0.5b.]
  - Naprogramujte funkci //{{courses:a6m33zsl:myradon.m?linkonly|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.]
  - Pomocí vámi napsané funkce transformujte Shepp-Loganův fantom, který si můžete vytvořit funkcí {{https://www.mathworks.com/help/images/ref/phantom.html|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.]


==== Podrobný popis ====

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

{{ :courses:a6m33zsl:radonproj.png?400 | Radonova transformace pro projekci p, theta}}
 

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>


==== Možný postup ====

  - Je dobré nastavit souřadnice obrázku tak, aby bod ''[0, 0]'' ležel uprostřed obrázku
  - Určete minimální a maximální hodnoty pro <latex>p, q</latex> (1. krok zadání)
  - Zvolte diskretizaci pro <latex>p, q</latex> v rozsahu <latex> [\min p, \max p] </latex>, respektive <latex> [\min q, \max q] </latex>
  - 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>
    - 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>
    - 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 {{http://www.mathworks.com/help/matlab/ref/interp2.html|interp2}}
    - spočtěte numerickou hodnotu <latex>J(\theta,p_h)</latex>