/*
===== Lab 10 : Diffusion weighted (MRI) =====


Today we look at Diffusion-weighted Imaging - one specific MR modality. With dMRI, we can use other gradients to measure the movement of molecules in the direction given by this gradient. If there is no obstacle at the measuring point, then the signals at different gradients will not differ. This is not the case when we have obstacles in our way, such as axons of nerve connections. Their beams promote diffusion along the beam direction and restrict it in trash across. Thus, a signal from a sufficient number of diffusion gradients can tell us what structure the tissue has at a given point.

One of the basic models is Diffusion Tensor Imaging (DTI), which uses a symmetric tensor (3 × 3 matrix) to determine the main direction of diffusion.
*/



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






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





==== 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. {{ :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>** 

 


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[\omega_j t_s + r \phi_i]$$

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.
*/