====== Cvičení 9 : Rekonstrukce MRI obrazu ======

Úkolem cvičení je vytvořit obraz v k-prostoru a poté rekonstruovat <latex>\rho</latex>-, <latex>T_1</latex>- a <latex>T_2</latex>-vážený obraz z NMR (MRI) pomocí Fourierovské rekonstrukce.

/*
Stáhněte si {{courses:a6m33zsl:mri_fantom.mat?linkonly|mat-file}} s 3 maticemi (//pd, t1, t2// odpovídající <latex>\rho</latex>, <latex>T_1</latex>, <latex>T_2</latex>), které budou tvořit náš virtuální fantom. Relaxační doby <latex>T_1</latex> a <latex>T_2</latex> jsou v //ms// (proto je třeba zadávat <latex>T_r</latex> a <latex>T_e</latex> v ms).


Pro hlubší porozumění tématu doporučujeme se seznámit s interaktivními animacemi na následující strance [[http://www.imaios.com/en/e-Courses/e-MRI/MRI-Sequences/Spin-echo|Spin echo]] (login: ZSL; heslo: zsl<lonsky rok>)
*/

Pro rekonstrukci budeme používat virtuální fantom -- tři matice (''pd, t1, t2'' opovídající <latex>\rho</latex>, <latex>T_1</latex>, <latex>T_2</latex>). Fantom si vygenerujte pomocí skriptu {{ :courses:a6m33zsl:get_nmr_phantom.m | get_nmr_phantom}}. 

{{ :courses:a6m33zsl:phantom_properties.png?800 |}}

Šablona pro cvičení, včetně vygenerování fantomu je k dispozici ve skriptu {{ :courses:a6m33zsl:mri_reconstruction.m | mri_reconstruction}}. 

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

  - Implementujte funkci //{{courses:a6m33zsl:encodingmri.m?linkonly|[k_space,t]=encodingMRI(pho,T_1,T_2,T_r,T_e)}}// kodující obrázek NMR, jejímž vstupem bude virtuální fantom (matice <latex>\rho</latex>, <latex>T_1</latex> a <latex>T_2</latex>) a výstupem matice //K// (k-prostor), kde řadky budou odpovídat jednotlivým excitacím (viz. [[#Implementace]]). [2b]
  - Aplikujte vytvořenou funkci na virtuální fantom, nastavte <latex>T_r=5</latex>s, <latex>T_e=1</latex>ms. Zobrazte absolutní hodnotu obrazu v k-prostoru, hodnoty jednotlivých prvků matice <latex>K</latex> mapujte na jasy logaritmicky. [1b]
  - Proveďte výpočet obrazu v k-prostoru a zpětnou rekonstrukci obrazu pomocí inverzní FFT ([[http://www.mathworks.com/help/matlab/ref/ifft2.html|ifft2]]) pro 3 následující kombinace parametrů v tabulce. [2b]
    * Která z kombinací odpovídá <latex>T_1</latex>, <latex>T_2</latex> a která <latex>\rho</latex> váženému obrazu a proč?
    * Uveďte rovnici a vysvětlete jednotlivé případy.
    * Zobrazte všechny tři rekonstruované obrazy, vyřízněte jen tu část, která odpovídá vstupu. Porovnejte s maticemi parametrů <latex>T_1</latex>, <latex>T_2</latex>, <latex>\rho</latex>.

^   ^   1  ^   2  ^   3  ^
^ <latex>T_r</latex>  |  5s  |  5ms  |  5s  |  
^ <latex>T_e</latex>  |  1ms  |  1ms  |  1s  |  

===== Rekonstrukce obrazu =====

MRI nám dává možnost pomocí nastavení parametrů <latex>T_r</latex> a <latex>T_e</latex> vytvářet obrazy vážené podle <latex>\rho</latex>  (hustoty protonů, proton density), nebo podle relaxačních časů <latex>T_1</latex> a <latex>T_2</latex>. Jak se dnes přesvědčíte, nepodaří se nám nikdy získat obraz reprezentující čistě <latex>\rho</latex>, <latex>T_1</latex>, nebo <latex>T_2</latex>. V obraze bude vždy přitomná i jistá složka ostatních parametrů. Vhodnou volbou parametrů <latex>T_e</latex> a <latex>T_r</latex> však můžeme vliv těchto nežádoucích složek potlačit.

Nastudujte si {{:courses:a6m33zsl:mri3.pdf|přednášky}} o MRI.

===== Implementace =====

Implementujte funkci rekonstruující obrázek NMR, jejímž vstupem bude virtuální fantom (matice <latex>\rho</latex>, <latex>T_1</latex> a <latex>T_2</latex>) a výstupem bude matice //K// (k-prostor), kde řádky budou odpovídat jednotlivým excitacím.
  - Pro každý pixel obrázku vypočtěte amplitudu výstupního signálu podle zjednodušeného vzorce: <latex>U=\rho(1-e^{-\frac{T_r}{T_1}})e^{-\frac{T_e}{T_2}} </latex>
  - Zvolte počet excitací (opakování) pro fázové kódování. Počet opakování //M// je roven počtu řádků vstupního obrazu.
  - Vytvořte vektor fází <latex>\varphi\in\langle0,2\pi)</latex>, <latex>\varphi_n=(n-1)\frac{2\pi}{M},~n=1\ldots M</latex>. Každá řádka obrazu je kódována jednou fází <latex>\varphi_n</latex>, každý element vektoru <latex>\boldsymbol{\varphi}</latex> bude odpovídat jedné excitaci.
  - Vytvořte vektor frekvencí <latex>\boldsymbol{\omega}</latex> pro frekvenční kódování sloupců <latex>\omega_m=100m\textrm{ [Hz]},~m=1\ldots S</latex>, kde //S// odpovídá počtu sloupců vstupního obrazu.
  - Zvolte počet časových vzorků na řádku, <latex>N=aS,\textrm{ kde }a\ge2,~a\in \mathbb{N}</latex>.
  - Vypočítejte vzorkovací frekvenci <latex>\Omega=a\max{(\boldsymbol{\omega)}}</latex>.
  - Vypočítejte interval vzorkování <latex>T=\frac{2\pi}{\Omega}</latex>.
  - Vygenerujte vektor vzorkovacích časů <latex>t_s=(s-1)T,~s=1\ldots N</latex>.
  - Spočítejte reálnou část obrazu v k-prostoru. Obraz v k-prostoru je matice o velikosti <latex>M\times N</latex> s elementy: <latex> \mathbf{K}(r,s)=\sum_{i,j=1}^{i=M,j=S}U_{i,j}cos(\omega_j\mathbf{t}(s)+(r-1)\varphi_i)\label{suma}</latex>, kde //i// je index fázového kódování, //j// index frekvenčního kódování, //r// pořadové číslo opakování a <latex>\mathbf{t}</latex> vektor vzorkovacích časů.
  - Pomocí Hilbertovy transformace dopočítejte imaginární složku obrazu v k-prostoru. Víme že výsledkem má být reálný obraz, hledáme tedy imaginární část k-prostoru tak, aby zpětná fourierova transformace celého k-prostoru byla reálná. K tomu poslouží Hilbertova transformace, v Matlabu //hilbert(x)//. Pozor, Matlabská funkce funguje po sloupcích, my potřebujeme Hilbertovu transformaci po řádcích ($\textrm{k_space} = hilbert(\mathbf{K}^\intercal)^\intercal$).

{{ :courses:a6m33zsl:mri_kspace.png?400 |Absolutní hodnota k-prostoru pro T_r = 5s, T_e = 1ms. Jasy mapovány logaritmicky, a = 4}}

Absolutní hodnota k-prostoru pro <latex>T_r = 5</latex>s, <latex>T_e = 1</latex>ms. Jasy mapovány **logaritmicky**, a = 4.


{{tag>medical MRI reconstruction matlab}}