V návaznosti na předchozí cvičení a téme rekonstrukce MR obrazových dat se v první části domácího cvičení budeme zabývat artefakty v rekonstruovaných MR snímcích, které jsou způsobené různými nepřesnostmi v k-prostoru. V druhé čáasti pak přijde na řadu technika zrychleného snímání pomocí více detekčních cívek – tzv. parallel imaging.

Update 18/05: původní data pro poslední část obsajují chybu, rekonstrukce nevychází úplně čistě. Pokud již máte odevzdáno, nemusíte zprávu předělávat, jinak pracujte s verzí parallel_acq_fixed.mat.

Budeme vycházet z $K$-matice uložené v axbrain.mat. Nejprve si zobrazte jak matici $K$-matici (nejsnáze lze přes $\log(|K|))$ tak i rekonstruovaný obraz (aplikujte inverzní FFT na matici $K$ a vykreslete absolutní hodnoty). Ve všech následujících částech budeme měnit matici $K$, rekonsturovat obrazu ze změnené matice a oboje vykreslovat, není na škodu si všechno exportovat do funkce.

Bodování: Části (1-3) jsou každá za [1 b]. Část (4) je pak ohodnocena max. [2 body].

Pro ilustraci toho, jak si odpovídají oblasti k-prostoru a výsledného obrázku budeme aplikovat low-pass a high-pass filtr. Vytvořte si masku (datový typ logical), která bude obsahovat kruh vycentrovaný kolem středové souřadnice, radius si zvolte v rozpětí [50, 80]px. Aplikujte masku na matici K tak, aby při low-pass zůstaly zachovány nízké frekvence, při high-pass filtru naopak nízké vyfiltrujte. Vstupní matice je již seřazená (aplikován fftshift) a nízké frekvence jsou uprostřed. Vložte výsleky do zprávy a okomentujte, jak se liší od původního.

Simulujte selhání vyčítání signálu v jednom bodě. Nastavte jeden pixel matice K na vysokou hodnotu (např. maximum původní matice K) a proveďte rekonstrukci. Proveďte alespoň pro dvě různé pozice (blízko a daleko od středu).

Průběh gradientů si v ideálním případě vykreslujeme jako schodový 0/1-signál, ve skutečnosti ale každé zapnutí/vypnutí gradientu chvíli trvá. Při snímací sekvenci EPI (echo-planar imaging), která vzorkuje matici K řádek po řádku – liché zleva a sudé zprava, se efekt zapínání všech nutných gradientů může projevit zpožděním echa a tedy předčasným vyčítáním signálu. Tím jak se K matice snímá cik-cak se jednotlivé řádky (a hlavně jejich středy) vůči sobě posunou. Nasimulujeme to tak, že liché řádky původní matice K posuneme o 2px doprava (vynechané pozice nastavíme na hodnotu 0), sudé řádky analogicky doleva.

http://mriquestions.com/nyquist-n2-ghosts.html

V klinické praxi se běžně setkáme se snímáním pomocí většího počtu snímacích cívek rozmístěných blízko snímaného objektu. Každá cívka má nejlepší signál ve svém okolí, skládáním výstupů více cívek získáme výsledný MR snímek. Při snaze vylepšít rozlišení obrázku zaplatíme za každou snímanou řádku matice K, více cívek samo o sobě tedy pomůže zlepšít SNR, ale snímání neurychlí. Snaha o zrychlení snímání při zachování dobrého rozlišení vedla k technikám parelelního snímání. Když totiž věnujeme více energie do rekonstrukce obrazu, tak nám postačí nasnímat pouze polovinu K-prostoru a ušetřit tedy polovinu času.

Budeme uvažovat paralelní cívky a podíváme na princip rekonstrukční metody Sensitivity Encoding for Fast MRI (SENSE). Ke každé z cívek máme mapu citlivostí, tedy matice $S_1, S_2$ a dále pak podél osy $y$ poloviční K-prostor $K_1, K_2$. Při SENSE nejprve provedeme rekonstrukci částečných snímků $I_j = \|FFT^{-1}(K_j)\|, j=1, 2$ a to i za cenu artefaktu překlopení krajních polovin snímku.

Pozn.:Obrázek používá fázové kódování podél osy $x$, v našem případě (fázové kódování je podél $y$) bude situace otočená o 90 stupňů.

Pro rekonstrukci výsledného obrázku použijeme toho, že každý bod v částečných obrázcích obsahuje váženou kombinaci dvou hodnot, tedy platí

$$ I_1 = X(A) * W_1(A) + X(B) * W_1(B) $$

$$ I_2 = X(A) * W_2(A) + X(B) * W_2(B) $$

A to pro dvojice bodů $A + N/2 = B$ (tedy posun o polovinu výsledné velikosti). V průběhu rekonstrukce tedy řešíme soustavu rovnic pro každý bod rekonstruovaného obrázku $X$.

Data Balíček parallel_acq.mat parallel_acq_fixed.mat obsahuje nasnímané K matice pro cívku 1 a 2 (k_coil_1 a k_coil_2 a mapu sensitivit MapW.

Update 18/05: původní data obsajují chybu, rekonstrukce pak nevychází úplně čistě. Pokud již máte odevzdáno, nemusíte zprávu předělávat, jinak pracujte s verzí parallel_acq_fixed.mat.