===== Lab 10 : Rekonstrukce obrazu (MRI) =====

V předchozím cvičení jsme se zabývali závislostí signálu jednoduché sekvence s opakovými RF pulsy na parameterech TR/TE. Průběh signálu jsme uvažovali pro jednu látku, bez prostorového umístění. Ostatně při této sekvenci bychom měřili **jeden** signál -- celkové echo celého excitovaného objemu. Abychom získali 3D obraz, potřebujeme měřené echo //prostorově zakódovat//. To děláme pomocí gradientů, kterými se budeme zabývat v tomto cvičení.


**Frekvenční gradient** (frequency-encoding gradient) Pomocí cívek na vnějším plášti MR skeneru lze vytvořit //gradient frekvenčního kódování// $G_f$, což je doplňkové magnetické pole, které se přičítá k hlavnímu (statickému) poli $B_0$ (princip superpozice). Tento gradient začíná (ve smyslu osy $x$) na levé straně a lineárně roste směrem k pravé. Indukce (efektivní) celkového pole $B(x)$ pole tedy roste podél osy $x$ a tím pádem roste i spin a frekvence echa (//Larmorova rovnost// říká: $f = \gamma B(x)$. S timto kódováním jsme se setkali v předchozím cvičení, nastavovali jsme sílu přidaného gradientu tak, abychom daným pulsem excitovali pouze omezenou oblast osy $x$. Frekvenční gradient je aktivní během echa (resp. vyčítání signálu).

**Fázový gradient** (phase-encoding gradient) Opět použijeme lineární gradient $G_\phi$, kterým dočasně změníme rezonanční frekvenci ($\omega_1$ = $\gamma \cdot (B_0 + B_1) $). Po vypnutí gradientního pole se protony vrací zpátky ke své //původní// rezonanční frekvenci $\omega_0$, ale jejich fáze bude vůči výchozímu stavu posunuta. Velikost posunu je závislosti na síle a délce trvání $G_y$. Fázové kódování se zpravidla aplikuje podél osy $y$. Gradient fázového posunu se zapne a vypne před vytvořením echa.

**K-prostor** Je datová matice nad souřadnicemi $(k_x, k_y)$, někdy také //dočasný obrazový prostor//, kterou plníme vzorkováním signálu echa. Použitím frekvenčního a fázového kódování a definicí $k_x = -\gamma G_f t$ a $k_y = -\gamma G_\phi t$ získáváme souřadnici v reálném obrazu $x$ jako konjugovanou proměnnou k $k_x$, a obdobně pro $y$ a $k_y$. Díky tomu pak probíhá rekonstrukce MRI obrazu aplikací inverzní 2D-FFT.


Další informace jsou opět pěkně zpracovány na Q-MRI: [[http://mriquestions.com/how-to-locate-signals.html|prostorové kódování]], přehledová stránka //[[http://mriquestions.com/hellipmaking-an-image.html|making an image]]//, frekvenční a [[http://mriquestions.com/what-is-phase-encoding.html|fázové kódování]], [[http://mriquestions.com/k-space-basic.html|k-prostor]].

==== Domácí cvičení ====

  - **[3 pt]** Naprogramujte funkci //{{courses:zsl:encodingmri.m?linkonly|[k_space,t]=encodingMRI(pho,T_1,T_2,T_r,T_e)}}// která dostane jako vstupní parametrické mapy $\rho$, $T_1$ a $T_2$, a vrací matici //K// (matici k-prostoru). Jedna řádka matice //K// obsahuje signál echo nasnímaný (navzorkovaný) pro jeden fázový gradient (podrobněji v [[#Detaily implementace]]). {{ :courses:zsl:mri_recon_constants.png?600 |}}
  - **[2 pt]** Pomocí skriptu {{ :courses:a6m33zsl:get_nmr_phantom.m | get_nmr_phantom}} si vygenerujte parametrické mapy pro $T_1$, $T_2$ a $\rho$ (PD). Pro kombinace parametrů $TE, TR$ uvedené v tabulce níže:
    - spočítejte matici $K$
    - vykreslete matici $K$ (zobrazte amplitudu a použijte logaritmus pro převedení na stupně šedi)
    - proveďte rekonstrukci obrazu pomocí inverzní FFT(2D) a také vykreslete.  
    - Jaké kombinace parametrů z tabulky vedou k $T_1$-, $T_2$- a $\rho$-váženému obrazu?

/*  - **[1 pt]** [BONUS] Select one of the TE/TR settings, and add Gaussian noise to the frequency sample values prior to computing each row of the $k$-space (more details in the next section) and compare the resulting reconstructed image to the non-noisy reconstruction [1 pt]

*/

^       ^    A     ^    B    ^    C     ^
^ $TR$  |  2500 ms |  2500 ms    |  400 ms  |   
^ $TE$  |  90 ms  |  10 ms  |  90 ms  |
**<fc #ff1100>Parametrické mapy T1, T2 obsahují hodnoty v [ms], s časy TE a TR je rovněž potřeba pracovat v [ms]</fc>** 

  
==== Detaily implementace ====

Na vstupu máme tři matice pro parametry $\rho, T_1, T_2$, které získáme zavoláním funkce <code = matlab>[P_pd, P_t1, P_t2] = get_nmr_phantom(n);
</code>. V dalším popisu budu pro odlišení řádek a sloupců používat $N, M$, ale generátor fantomu má pouze jeden parameter //n//, tedy platí $M = N = n$ a parametrické mapy jsou čtvercové.

Z parametrických map a z parametrů sekvence //TE// a //TR// pak na základě zjednodušené rovnice 

$$U = \rho \cdot (1 - e^{-\frac{TR}{T_1}}) \cdot e^{-\frac{TE}{T_2}}$$

můžeme určit amplitudu měřeného signálu. Zde si musíme uvědomit, že echo $S(\omega, \phi)$, které lze popsat rovnicí

$$ S = U * \cos[\omega \cdot t + \phi] + \sqrt{-1} * sin [\omega \cdot t + \phi] $$

je souhrnný signál za celou excitovanou 2D-vrstvu. Tedy každý prostorový bod $(x, y)$ ovlivňuje hodnoty více (všech) bodů v k-prostoru $(k_x, k_y)$. Pro fixní fázový posun naplníme vzorkováním signálu jednu řádku matice $k$-prostoru, opakováním měření pro různé fázové posuny pak postupně naplníme všechny řádky.

{{:courses:zsl:labs_2020:mri_kspace.jpg?400|Zdroj: http://www.mri-q.com/data-for-k-space.html}}


Prvním krokem bude diskretizace rovnice výše v proměnných $U$, $t$, $\omega$ a $\phi$, samozřejmě s využitím souvislostí mezi fází $\phi_m$ a prostorovu osou $y$ v $U$ a mezi frekvencí $\omega_m$ a osou $x$. Mějme na vstupu parametrické mapy o velikosti $N \times M$.  Pak již můžeme určit:

pro **frekvenční/fázové kódování**
   * vektor diskrétních fází $\phi_n \in \langle 0, 2\pi)$, $n=1,..,N$ pro každý řádek obrázku.
   * vektor frekvencí pro kódování sloupců $\omega_m = 100 \cdot m$ [Hz] pro každý sloupec $m=1,...,M$.

pro **diskretizaci signálu**
   * vzorkovací frekvenci (minimálně dvojnásobek maximální frekvence, abychom se nedostali přes //Nyquistův vzorkovací limit//) $\Delta T = \frac{2\pi}{\max(\omega)}$
   * vektor vzorkovacích časů $t_s = (s-1) * \Delta T, \quad s=1,...,M$

Vzorkování signálu v nastavených samplovacích bodech $t_s$:

$$ K(r, t_s) = \sum_{i, j}^{N, M} U(i, j) \cos[\omega_j t_s + r \phi_i] + \sqrt{-1} \sin[...]$$

opakujeme pro každou řádku matice $K$. Nezapomeňte, že echo je vždy celkové echo z kompletní vrstvy, proto jedna hodnota K(r, t_s) závisí na součtu přes $N$ a $M$.

/*
In the bonus task, distort the values of $\omega_j$ by generating zero-mean Gaussian noise $\varepsilon_j(r)$ for each repetition to get $\omega_j^{(n)} = \omega_j + \varepsilon_j(r)$ and use it instead of $\omega_j$ when filling the $K$-space.
*/