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.

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

Šablona pro cvičení, včetně vygenerování fantomu je k dispozici ve skriptu mri_reconstruction.

1. Implementujte funkci [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]
2. 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]
3. Proveďte výpočet obrazu v k-prostoru a zpětnou rekonstrukci obrazu pomocí inverzní FFT (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>.

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 přednášky o MRI.

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.

1. 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>
2. 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.
3. 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.
4. 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.
5. Zvolte počet časových vzorků na řádku, <latex>N=aS,\textrm{ kde }a\ge2,~a\in \mathbb{N}</latex>.
6. Vypočítejte vzorkovací frekvenci <latex>\Omega=a\max{(\boldsymbol{\omega)}}</latex>.
7. Vypočítejte interval vzorkování <latex>T=\frac{2\pi}{\Omega}</latex>.
8. Vygenerujte vektor vzorkovacích časů <latex>t_s=(s-1)T,~s=1\ldots N</latex>.
9. 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ů.
10. 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$).

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