====== Cvičení 10 : Diffusion-weighted imaging ======

Dnes se podíváme na //Diffusion-weighted Imaging// - jednu specifickou modalitu MR. Při dMRI můžeme pomocí dalších gradientů měřit pohyb molekul ve směru daném tímto gradientem. Pokud se v místě měření nenachází žádná překážka, pak se signály při odlišných gradientech nebudou lišit. Jinak tomu je, když máme v cestě překážky, například axony nervových propojení. Jejich svazky podporují difůzi podél směru svazku a omezují ji ve smětu napříč. Signál z dostatečného množství difúzních gradientů nám tedy může napovědět, jakou strukturu má tkáň v daném bodě.

Jedním ze základních modelů je Diffusion Tensor Imaging (DTI), které používá symetrický tensor (matice 3x3) k určení hlavního směru difůze.

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

  - Ze souboru {{ :courses:a6m33zsl:diffusion_data.mat |}} načtěte 4D obrázek (dimenze ''X,Y,1,grad'') ''dwimage'' a informace o gradientech ''bvecs'', ''bvals''. Najděte pro ''grad'' > 3 dva obrázky pro různé směry ''grad'', které se od sebe výrazně liší a vykreslete. Proč je signál odlišný a co nám to napovídá o měřené tkáni? [1b]
  - Spočtěte pomocí funkce ''dti_recon_2d'' (dostupné {{ :courses:a6m33zsl:dti_recon_2d.m | zde}}) tensory v každé voxelu. Z tensorů pak pro každý voxel určete hodnoty //fractional anisotropy// (FA), //mean diffusivity//(MD) a //radial diffusivity// (RD) a výsledné FA, MD a RD mapy vykreslete. [1.5b]
  - Spočtěte mapy z předchozího bodu také na {{ :courses:a6m33zsl:diffusion_data_2.mat |}}. Čím se od sebe liší výsledné mapy pro ''diffusion_data'' a ''diffusion_data_2''? Pro vygenerování dat byl pro oba obrázky použit simulátor,který na základě konfigurace //vláken// počítá příslušný (realistický) dMRI signál. Vstupní konfigurace vláken je na následujících obrázcích. Přiřaďte k sobě snímky, konfigurace a tensory (obrázek níže) odůvodněte Vaši volbu.[1.5b]
{{ :courses:a6m33zsl:dti_fiber_combi.png?nolink&600 |}}
  - Vytvořte pro ''diffusion_data'' v každém voxelu vizualizaci hlavního směru difůze (vlastní vektor příslušící k největší vlastní hodnotě) a vykreslete (nejsnáze to jde pomocí funkce [[https://www.mathworks.com/help/matlab/ref/quiver.html|quiver]]). Jsou na snímku místa, kde je informace o hlavním směru spočtená pomocí DTI nedostatečná? [1b]



==== Návod ====

Soubor ''dwimage'' obsahuje pro každý gradient jeden 2D řez velikosti 20x20, gradientní schéma obsahuje 3x nevážený obraz (<latex>b_0</latex>) a 30x gradient o //síle// <latex> b = 1000 </latex>. Celkově má tedy rozměry 20 x 20 x 1 x 33. Odpovídající směry gradientů jsou uložené jako 33x3 matice ''bvecs''. Směry si můžeme vykreslit pomocí
<code matlab>
scatter3(bvecs(1,:), bvecs(2,:), bvecs(3,:), 50, 'filled');
</code>

=== Reconstrukce tensorů ===

//Diffusion tensor// druhého řádu je symetrická matice o velikosti 3x3, která popisuje difůzi do směru <latex> x, y</latex> a <latex>z</latex>. Pokud je difůze //izotropní//, pak je v ideálním případě <latex> 
T = \begin{pmatrix} %
D & 0 & 0 \\
0 & D & 0 \\
0 & 0 & D \\
\end{pmatrix}
</latex>

V //anizotropním// případě (některý ze směrů neumožňuje volnou difůzi) se začínají objevovat členy mimo diagonálu. 

Pro rekonstrukci tensorů musíme nejdříve převést do strukturované formy, abychom mohli následně zavolat funkci ''dti_recon'':

<code matlab>
%% DTI Reconstruction
for i=1:numel(bvals),
    DTIdata(i).VoxelData = single(dwimage(:,:,:,i));
    DTIdata(i).Gradient = bvecs(:,i);
    DTIdata(i).Bvalue = bvals(i);
end

DTens = dti_recon_2d(DTIdata, 223);
</code>

Funkce nám v každém bodě vrátí vektor <latex>M</latex> s koeficienty matice <latex>T</latex>, respektive její vrchní trojúhelníkové části, což je díky symetrii dostačující popis. Tj. dostáváme

<latex>
T = \begin{pmatrix} %
M(1) & M(2) & M(3) \\
. & M(4) & M(5) \\
. & . & M(6) \\
\end{pmatrix}
</latex>

Pomocí vlastních vektorů <latex>\varepsilon_i</latex> a vlastních čísel <latex>\lambda_i</latex> tensoru <latex>T</latex> jsme schopni získat minimální reprezentaci pro elipsoid popisující preferované směry difúze. Pro (ideální) izotropní případ dostáváme sféru, při zesilující anisotropii v jednom směru se formuje //doutník//((Obrázek převzat ze stránky {{http://mriquestions.com/diffusion-tensor.html}} )):
{{ :courses:a6m33zsl:dti_tensor_ellipsoid.png?nolink&800 | Diffusion Tensor ellipsoid }}

Na základě získaných vlastních čísel a vektorů se pak informace o elipsoidu redukuje na skalární hodnoty. Jednou z nejvýznamnějších je //Fractional Anisotropy// (FA), která popisuje, jak silně převládá jeden směr nad ostatními. Rozezí hodnot je 0 - isotropní až 1 - pouze jeden směr. Vzorec pro výpočet
{{ :courses:a6m33zsl:faformula.png?nolink&400 | }}

//Mean Diffusivity //(MD) je střední hodnota vlastních čísel <latex>(\lambda_1 + \lambda_2 + \lambda_3) \cdot 1/3</latex>. Tato hodnota se také označuje jako Apparent Diffusion Coefficient (ADC). //Radial diffusivity // (RD) je průměr vlastních čísel dvou menších vlastních čísel <latex>(\lambda_2 + \lambda_3)/2</latex>. 

Výsledné konfigurace tensorů vypadají následovně. Přiřaďte ke konfiguracím vláken a spočteným FA, RD a MD mapám.
{{:courses:a6m33zsl:dti_tensor_combi.png?nolink&600 |}}

=== Hlavní směr ===
Z reprezentace elipsoidu je zřejmé, že  největší vlastní číslo a související vektor nám ukazují hlavní směr difúze. Pro vizualizaci tedy můžeme v každém bodě spočítat vektor <latex>(u, v)</latex> jako x-ové a y-ové souřadnice vektoru <latex>\lambda_1 \cdot \varepsilon_1</latex> (<latex>z</latex> komponenty můžeme v tomto případě ignorovat) a pomocí funkce ''quiver(x, y, u, v)'' vykreslit na odpovídající souřadnice pixelů <latex>x, y</latex>.