【FSL】Voxel-BasedAnalysis: VBA



1. 目的
2. VBAとは
2.1. ファイルの準備
2.2. 脳解剖学的標準化
2.3. 平滑化
2.4. GLMと並べ替え検定(permutation test)
3. おまけ


1. 目的

  • Voxel-BasedAnalysis: VBA

2. VBAとは

Voxel-Based Analysis(VBA)は、脳構造解析手法の一つであり、現在では脳科学の分野において幅広く用いられている。VBAは、観測者が関心領域を設定して解析する古典的なマニュアル計測とは異なり、自動処理によって全脳をボクセルごとに解析するため、評価の客観性が高い。

VBAの解析手順は、次の通り。

  1. ファイルの準備
  2. 脳解剖学的標準化
  3. 平滑化
  4. GLMと並べ替え検定(permutation test)

2.1. ファイルの準備

VBA解析をしたファイルを準備する。

ここでは、拡散MRIから計算したFA画像を準備する。健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)いることを想定する。

.
├── Con0001_FA.nii.gz
├── Con0002_FA.nii.gz
├── Con0003_FA.nii.gz
...
└── Pat0010_FA.nii.gz

2.2. 脳解剖学的標準化

非線形変換を用いて、全ての被験者のFA画像を標準FA画像(FMRIB58_FA)に位置合わせする。

以下のコマンドを実行する。

mkdir preprocessing
for SUBJ in $(ls | grep _FA.nii.gz | cut -d . -f1); do
    # 位置合わせ
    ## Affine transform
    flirt -in ${SUBJ}.nii.gz -ref ${FSLDIR}/data/standard/FMRIB58_FA_1mm.nii.gz \
        -dof 12 -omat preprocessing/${SUBJ}_indiv2std.mat
    ## non-linear transform
    fnirt --in=${SUBJ}.nii.gz --aff=preprocessing/${SUBJ}_indiv2std.mat --config=FA_2_FMRIB58_1mm.cnf \
        --cout=preprocessing/${SUBJ}_warp_indiv2std.nii.gz --iout=preprocessing/${SUBJ}_instd.nii.gz
done

位置合わせ方法の詳細は、こちらの記事を参考に。

https://www.nemotos.net/?p=4513

2.3. 平滑化

非線形変換を用いた位置合わせで、脳解剖学的標準化をしたとしても、ボクセル単位でみると完全に標準化できているわけではない。そこで、そのようなエラーを抑えるために平滑化処理をする。

まず、標準FA画像に位置合わせをした、被験者ごとのFA画像をまとめて4D-FA画像(all_FA.nii.gz)とし、4D-FA画像の平均画像からマスク画像を生成する。

mkdir stats
fslmerge -t stats/all_FA.nii.gz preprocessing/*_instd.nii.gz
fslmaths stats/all_FA.nii.gz -Tmean -bin mask stats/all_FA_mask.nii.gz

次に、ガウシアンフィルタを用いた平滑化処理をする。

ここでは、よく用いられるσ=3 voxel(FWHM=約7mm)のガウシアンフィルタで平滑化する。

fslmaths stats/all_FA.nii.gz -fmean -s 3 stats/all_FA_smoothed_s3.nii.gz

2.4. GLMと並べ替え検定(permutation test)

ガウシアンフィルタで平滑化した灰白質画像を用いて、健常群と患者群の群間比較をする。

まず、GLMのデザインマトリックス(計画行列)とコントラストを設定する。

今回は、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)のデータがある。すべての被験者のFAをlsコマンドでみると、先に健常者10名のFA、次に患者10名のFAが並んでいることが分かる。

ls |grep nii
Con0001_FA.nii.gz
Con0002_FA.nii.gz
Con0003_FA.nii.gz
...
Pat0010_FA.nii.gz

次に、GLMのデザインマトリックス(計画行列)とコントラストを決める設定ファイルを生成する。design_ttest2 <出力ファイル> <健常者数> <患者数>でコマンドを実行。

design_ttest2 stats/design 10 10

statsフォルダに、デザインマトリックス(design.mat)とコントラスト(design.con)が生成される。

デザインマトリックス(design.mat)の中身を確認。

/Matrixの一列目は健常者データであるかどうか、二列目は患者データであるかを0, 1で表している。行の順番は、被験者フォルダにあるファイルの順番(昇順)に対応する。したがって、これらは対応があるようにしておかなければならない。

cat stats/design.mat
/NumWaves 2
/NumPoints 20
/PPheights 1 1
/Matrix
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1

コントラスト(design.con)の中身を確認してみる。

/Matrix一列目は健常者の偏回帰係数、二列目は患者の偏回帰係数に対するもので、行は別々のコントラストである。この場合、一行目は健常者>患者の検定、二行目は健常者<患者の検定に相当する。

cat stats/design.con
/NumWaves 2
/NumContrasts 2
/PPheights 1 1
/Matrix
1 -1
-1 1

デザインマトリックス(計画行列)とコントラストの確認ができたら、randomiseコマンド使ってGLMと並べ替え検定(permutation test)を実行する。

randomiseコマンドの各オプションは、次の通り。

  • -i:入力画像
  • -m:マスク画像
  • -o :出力画像
  • -o :デザインマトリックス
  • -o :デザインコントラスト
  • -n:並べ替え検定の数
  • –T2:2D最適化を用いたTFCE
  • -x:voxel-wiseのcorrected P値マップ
  • –uncorrp:un-corrected P値マップ
  • -R:統計値マップ
design_ttest2 stats/design 10 10

randomise -i stats/all_FA_smoothed_s3.nii.gz \  
    -m stats/all_FA_mask.nii.gz \
    -o stats/fslvba \
    -d stats/design.mat \
    -t stats/design.con \
    -n 10000 --T2 -V -x --uncorrp -R

次のようなファイルが、生成される。

TFCEを用いた”健常群>患者群”の検定で、FWE補正をされたP値マップ(fslvba_tfce_corrp_tstat1.nii.gz)を確認する。ここで、得られたP値マップは1-P値のマップであることに注意する。つまり、P<.05を有意とするのであれば、P値マップで0.95-1.00の値を見ればよい

fsleyes $FSLDIR/data/standard/MNI152_T1_2mm \
    stats/fslvba_tfce_corrp_tstat1 \
    -cm Red-Yellow -dr 0.95 1

3. おまけ

検定数が多い場合、有意差があったかどうかをすべて確認するのは大変である。そこで、一目で有意差があるかどうかを判断できるように、各検定ごとのP値マップの最大値を自動計測し、テキスト(corrected_P_report.txt)としてまとめる。

for PMAP in $(ls stats/ | grep tfce_corrp); do
    PMAX=$(fslstats stats/${PMAP} -R | cut -d " " -f2)
    echo ${PMAP} >>stats/tmp1.txt
    echo ${PMAX} >>stats/tmp2.txt
done
paste stats/tmp* >stats/tmp_corrected_P_report.txt
echo -e "$(cat stats/tmp_corrected_P_report.txt)\n\n\n$(cat stats/tmp_corrected_P_report.txt | sort -r -n -k 2)" \
    >stats/corrected_P_report.txt
rm stats/tmp*

上のコマンドを実行すると、statsフォルダに各検定とそのP値マップの最大値が記された「corrected_P_report.txt」が出力される。

検定結果を、ファイル名でソート(上段)したものと、P値でソートしたもの(下段)に分けて保存している。

cat stats/corrected_P_report.txt
fslvba_tfce_corrp_tstat1.nii.gz 0.982391
fslvba_tfce_corrp_tstat2.nii.gz 0.993181


fslvba_tfce_corrp_tstat2.nii.gz 0.993181
fslvba_tfce_corrp_tstat1.nii.gz 0.982391

【FSL】Tract-Based Spatial Statistics: TBSS



1. 目的
2. TBSSとは
2.1. 拡散MRIからFA画像を計算
2.2. tbss_1_preproc:TBSSの準備
2.3. tbss_2_reg:FAを標準空間に非線形位置合わせ
2.4. tbss_3_postreg:平均FA画像を生成し、FAスケルトンを生成
2.5. tbss_4_prestats:被験者ごとのFA画像を平均スケルトンに投影
2.6. tbss_non_FA:FA画像以外の定量値をスケルトンに投影
2.7. randomise:スケルトンに投影された定量値画像を入力したGLMと並べ替え検定(permutation test)
3. おまけ


1. 目的

  • Tract-Based Spatial Statistics: TBSS

2. TBSSとは

Tract-Based Spatial Statistics (TBSS) は、白質路の統計解析をするための手法。

神経線維束の中心線(skeleton)に定量値を投影する。通常の脳画像の統計解析では、脳構造の個人差を除外するために空間的「平滑化」を用いる。しかし、平滑化の程度に原則がなく、平滑化をかけては情報があいまいになり、MRIの高空間分解能を生かせないという問題がある。一方、TBSSでは、神経線維束の中心線と思われるところにskeletonを生成し、そこに個人ごとの定量値を投影するという手法をとる。これにより、平滑化せずに群間比較をすることができるため、平滑化による問題を回避できるという利点がある。

TBSS解析では、次のような処理をする。

  1. 拡散MRIからFA画像を計算
  2. tbss_1_preproc:TBSSの準備
  3. tbss_2_reg:FAを標準空間に非線形位置合わせ
  4. tbss_3_postreg:平均FA画像を生成し、FAスケルトンを生成
  5. tbss_4_prestats:被験者ごとのFA画像を平均スケルトンに投影
  6. tbss_non_FA:FA画像以外の定量値をスケルトンに投影
  7. randomise:スケルトンに投影された定量値画像を入力したGLMと並べ替え検定(permutation test)

2.1. 拡散MRIからFA画像を計算

拡散MRIからFA画像を計算し、被験者ごとのFA画像を準備する。

ここでは、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)いることを想定する。

.
├── Con0001_FA.nii.gz
├── Con0002_FA.nii.gz
├── Con0003_FA.nii.gz
...
└── Pat0010_FA.nii.gz

その他の定量値画像も解析に含めたい場合は、FA画像と合わせて次のように用意する。

ここでは、MDとFWという定量値画像を用意した。この時、定量値画像ごとのフォルダを作成し、各被験者の定量値画像を格納するが、ファイル名はFA画像と全く同じにすること(例. FA画像: Con0001_FA.nii.gz, MD画像: MD/Con0001_FA.nii.gz)。

.
├── Con0001_FA.nii.gz
├── Con0002_FA.nii.gz
├── Con0003_FA.nii.gz
...
└── Pat0010_FA.nii.gz
├── FW
│   ├── Con0001_FA.nii.gz
│   ├── Con0002_FA.nii.gz
│   ├── Con0003_FA.nii.gz
│   ...
│   └── Pat0010_FA.nii.gz
└── MD
    ├── Con0001_FA.nii.gz
    ├── Con0002_FA.nii.gz
    ├── Con0003_FA.nii.gz
    ...
    └── Pat0010_FA.nii.gz

2.2. tbss_1_preproc:TBSSの準備

ファイルの準備ができたら、tbss_1_preprocコマンドでTBSSの準備をする。

tbss_1_preproc *_FA.nii.gz

処理が終わると、「FAフォルダ」と「origdataフォルダ」を生成し、それぞれにFA画像が格納される。FAフォルダには、さらにFAのマスク画像が生成される。

.
├── FA
│   ├── Con0001_FA.nii.gz
│   ├── Con0001_FA_FA_mask.nii.gz
│   ├── Con0002_FA.nii.gz
│   ├── Con0002_FA_FA_mask.nii.gz
│   ├── Con0003_FA.nii.gz
│   ├── Con0003_FA_FA_mask.nii.gz
│   ...
│   └── Pat0010_FA.nii.gz
│   ├── Pat0010_FA._FA_mask.nii.gz
├── FW
│   ├── Con0001_FA.nii.gz
│   ├── Con0002_FA.nii.gz
│   ├── Con0003_FA.nii.gz
│   ...
│   └── Pat0010_FA.nii.gz
├── MD
│   ├── Con0001_FA.nii.gz
│   ├── Con0002_FA.nii.gz
│   ├── Con0003_FA.nii.gz
│   ...
│   └── Pat0010_FA.nii.gz
└── origdata 
    ├── Con0001_FA.nii.gz
    ├── Con0002_FA.nii.gz
    ├── Con0003_FA.nii.gz
    ...
    └── Pat0010_FA.nii.gz

「FA/slicesdir/index.html」をブラウザ(e.g., Chrome)で開くことで、各被験者のFA画像一覧をみることができる。

2.3. tbss_2_reg:FAを標準空間に非線形位置合わせ

tbss_1_preprocコマンドで、全ての被験者のFA画像を1x1x1mmの標準空間に非線形的な位置合わせする。通常は、-Tオプションで標準空間にある標準FA画像(FMRIB58_FA)に位置合わせをするが、-tオプションを用いて任意の画像に位置合わせすることもできる(推奨)。また、-nオプションでは、被験者の中で最も位置合わせ先としてふさわしいFA画像を見つけ出し、そのFA画像にすべての被験者のFA画像を位置合わせすることができる。

ここでは、TBSSで推奨されている-Tオプションを指定しすべての被験者FA画像を標準FA画像(FMRIB58_FA)に位置合わせをする。

tbss_2_reg -T

処理が完了すると、FAフォルダに結果が保存される。

「FA/*_to_target_warp.nii.gz」が、標準FA画像に位置合わせするためのwarp fieldである。

2.4. tbss_3_postreg:平均FA画像を生成し、FAスケルトンを生成

先程生成した標準FA画像に位置合わせするためのwarp fieldを用いて、各被験者のFA画像を標準空間(MNI152)に移動させる。その後、平均FA画像を生成し、その平均FA画像からFAスケルトンを生成する。tbss_3_postregでは、これらの処理を-Sオプションで実行することができが、代わりに標準FA画像(FMRIB58_FA mean)とそのスケルトンを用いたい場合は-Tを指定する。

ここでは、TBSSで推奨されている-Sオプションを指定する。

tbss_3_postreg -S

処理後の画像は、「statsフォルダ」に格納される。

stats/
├── all_FA.nii.gz  # 標準空間上における各被験者のFA画像
├── mean_FA.nii.gz  # 平均FA画像
├── mean_FA_mask.nii.gz  # 平均FA画像のマスク
└── mean_FA_skeleton.nii.gz  # 平均FA画像から生成したスケルトン画像

2.5. tbss_4_prestats:被験者ごとのFA画像を平均スケルトンに投影

tbss_4_prestatsコマンドでは、まず平均FA画像から生成したスケルトン画像(mean_FA_skeleton.nii.gz)をしきい値処理(通常 0.2)をし、スケルトンのバイナリーマスク画像を生成する。次に、このスケルトンマスクからの距離マップ(distance map)が計算され、この距離マップを参考に、被験者ごとのFA画像をスケルトン画像に格納(投影)する。

tbss_4_prestats 0.2

statsフォルダに、新たに次のファイルが生成される。

stats/
├── all_FA_skeletonised.nii.gz  # スケルトンに投影されたすべての被験者のFA画像
├── mean_FA_skeleton_mask.nii.gz  # スケルトンマスク
├── mean_FA_skeleton_mask_dst.nii.gz  # スケルトンマスクからの距離マップ(distance map)
└── thresh.txt  # バイナリースケルトンマスク画像を作る際のしきい値

2.6. tbss_non_FA:FA画像以外の定量値をスケルトンに投影

tbss_non_FAコマンドで、FA画像以外の定量値をスケルトンに投影する。このとき、標準空間への移動やスケルトンを生成するためのパラメータは、FA画像で使ったものが適用される。

NONFA_LIST=$(ls -F | grep / | cut -d / -f 1 | grep -v stats| grep -v origdata)

for MAP in ${NONFA_LIST}; do
    tbss_non_FA ${MAP}
done

2.7. randomise:スケルトンに投影された定量値画像を入力したGLMと並べ替え検定(permutation test)

まず、GLMのデザインマトリックス(計画行列)とコントラストを設定する。

今回は、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)のデータがある。「origdataフォルダ」をみると、先に健常者10名のFA画像、次に患者10名のFA画像が並んでいることが分かる。

ls -1 origdata
Con0001_FA.nii.gz
Con0002_FA.nii.gz
Con0003_FA.nii.gz
...
Pat0010_FA.nii.gz

次に、GLMのデザインマトリックス(計画行列)とコントラストを決める設定ファイルを生成する。design_ttest2 <出力ファイル> <健常者数> <患者数>でコマンドを実行。

design_ttest2 stats/design 10 10

statsフォルダに、デザインマトリックス(design.mat)とコントラスト(design.con)が生成される。

デザインマトリックス(design.mat)の中身を確認。

/Matrixの一列目は健常者データであるかどうか、二列目は患者データであるかを0, 1で表している。行の順番は、origdataフォルダにあるファイルの順番(昇順)に対応する。したがって、これらは対応があるようにしておかなければならない。

cat stats/design.mat
/NumWaves 2
/NumPoints 20
/PPheights 1 1
/Matrix
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1

コントラスト(design.con)の中身を確認してみる。

/Matrix一列目は健常者の偏回帰係数、二列目は患者の偏回帰係数に対するもので、行は別々のコントラストである。この場合、一行目は健常者>患者の検定、二行目は健常者<患者の検定に相当する。

cat stats/design.con
/NumWaves 2
/NumContrasts 2
/PPheights 1 1
/Matrix
1 -1
-1 1

デザインマトリックス(計画行列)とコントラストの確認ができたら、randomiseコマンド使ってGLMと並べ替え検定(permutation test)を実行する。

randomiseコマンドの各オプションは、次の通り。

  • -i:入力画像
  • -m:マスク画像
  • -o :出力画像
  • -o :デザインマトリックス
  • -o :デザインコントラスト
  • -n:並べ替え検定の数
  • –T2:2D最適化を用いたTFCE
  • -x:voxel-wiseのcorrected P値マップ
  • –uncorrp:un-corrected P値マップ
  • -R:統計値マップ
for MAP in FA ${NONFA_LIST}; do
    randomise -i stats/all_${MAP}_skeletonised \
        -o stats/tbss_${MAP} \
        -m stats/mean_FA_skeleton_mask \
        -d stats/design.mat \
        -t stats/design.con \
        -n 10000 --T2 -V -x --uncorrp -R
done

次のようなファイルが、生成される。

TFCEを用いた”健常群>患者群”の検定で、FWE補正をされたP値マップ(tbss_FA_tfce_corrp_tstat1.nii.gz)を確認する。ここで、得られたP値マップは1-P値のマップであることに注意する。つまり、P<.05を有意とするのであれば、P値マップで0.95-1.00の値を見ればよい

fsleyes ${FSLDIR}/data/standard/FMRIB58_FA_1mm.nii.gz \
    ${FSLDIR}/data/standard/FMRIB58_FA-skeleton_1mm.nii.gz -cm Green \
    stats/tbss_FA_tfce_p_tstat1.nii.gz -cm Red-Yellow -dr 0.95 1

スケルトンは細いため有意差が見づらい場合がある。そのような時、tbss_fillコマンドが役に立つ。

tbss_fillコマンドの基本的な使い方は、以下の通り。

tbss_fill <P値マップ> <しきい値> <平均FA画像> <出力画像>

TFCEを用いた”健常群>患者群”の検定で、FWE補正をされたP値マップ(tbss_FA_tfce_corrp_tstat1.nii.gz)の有意差があった領域のみを0.95でしきい値処理をして抽出し、その領域を膨張させる。

tbss_fill stats/tbss_FA_tfce_p_tstat1.nii.gz 0.95 stats/mean_FA stats/tbss_FA_tfce_p_tstat1_fill.nii.gz

赤く表示されている領域は、健常群が患者群よりも有意(FWE-corrected)にFA値が大きいことを示している。

3.おまけ

大量に定量値があり、それらをすべて検定する場合、有意差があったかどうかをすべて確認するのは大変である。そこで、各定量値画像のP値マップが0.95以上の値を持つかどうかを判定し、有意差があった場合のみ、tbss_fillコマンドを実行する。

for PMAP in $(ls stats/ | grep tfce_corrp); do
    PMAX=$(fslstats stats/${PMAP} -R | cut -d " " -f2)
    echo ${PMAP} >>stats/tmp1.txt
    echo ${PMAX} >>stats/tmp2.txt
    if [ $(echo "${PMAX} > 0.95" | bc) == 1 ]; then
        tbss_fill stats/${PMAP} 0.95 stats/mean_FA stats/${PMAP}_fill
    fi
done
paste stats/tmp* >stats/tmp_corrected_P_report.txt
echo -e "$(cat stats/tmp_corrected_P_report.txt)\n\n\n$(cat stats/tmp_corrected_P_report.txt | sort -r -n -k 2)" \
    >stats/corrected_P_report.txt
rm stats/tmp*

上のコマンドを実行すると、statsフォルダに各検定とそのP値マップの最大値が記された「corrected_P_report.txt」が出力される。

検定結果を、ファイル名でソート(上段)したものと、P値でソートしたもの(下段)に分けて保存している。

cat stats/corrected_P_report.txt
tbss_FA_tfce_corrp_tstat1.nii.gz    0.992000
tbss_FA_tfce_corrp_tstat2.nii.gz    0.416000
tbss_FW_tfce_corrp_tstat1.nii.gz    0.361839
tbss_FW_tfce_corrp_tstat2.nii.gz    0.997261
tbss_MD_tfce_corrp_tstat1.nii.gz    0.389816
tbss_MD_tfce_corrp_tstat2.nii.gz    0.985748


tbss_FW_tfce_corrp_tstat2.nii.gz    0.997261
tbss_FA_tfce_corrp_tstat1.nii.gz    0.992000
tbss_MD_tfce_corrp_tstat2.nii.gz    0.985748
tbss_FA_tfce_corrp_tstat2.nii.gz    0.416000
tbss_MD_tfce_corrp_tstat1.nii.gz    0.389816
tbss_FW_tfce_corrp_tstat1.nii.gz    0.361839

【FSL】FSLを用いたVoxel-Based Morphometry: VBM



1. 目的
2. VBMとは
2.1. 3D-T1WIの準備
2.2. fslvbm_1_bet: 3D-T1WIの脳頭蓋除去
2.3. fslvbm_2_template: 被験者の脳から灰白質テンプレート画像を生成
2.4. fslvbm_3_proc: すべての被験者の灰白質を灰白質テンプレート画像に位置合わせし、ボクセルをモジュレーション後、平滑化
2.5. randomise:GLMと並べ替え検定(permutation test)
3. おまけ


1. 目的

  • FSLを用いたVoxel-Based Morphometry: VBM

2. VBMとは

Voxel-Based Morphometry(VBM)は、脳構造解析手法の一つであり、特に脳容積を対象に解析する。

VBMは、古典的なマニュアルの脳容積計測とは異なり、自動処理によって全脳を客観的に評価することができ、現在では脳科学の分野において幅広く用いられている。

VBM解析では、次のような処理をする。

  1. 3D-T1WIの準備
  2. fslvbm_1_bet: 3D-T1WIの脳頭蓋除去
  3. fslvbm_2_template: 被験者の脳から灰白質テンプレート画像を生成
  4. fslvbm_3_proc: すべての被験者の灰白質を灰白質テンプレート画像に位置合わせし、ボクセルをモジュレーション後、平滑化
  5. randomise:GLMと並べ替え検定(permutation test)

2.1. 1. 3D-T1WIの準備

各被験者の3D-T1WIを準備する。

ここでは、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)いることを想定する。

.
├── Con0001_T1w.nii.gz
├── Con0002_T1w.nii.gz
├── Con0003_T1w.nii.gz
...
└── Pat0010_T1w.nii.gz

2.2. 2. fslvbm_1_bet: 3D-T1WIの脳頭蓋除去

ファイルの準備ができたら、fslvbm_1_betコマンドを実行して脳頭蓋除去をする。3D-T1WIに首を含まない場合、-bオプションを指定し、首を含む場合は-Nオプションを指定する。

ここでは、3D-T1WIに首を含んでいたとして、-Nオプションを指定する。

fslvbm_1_bet -N  # 首を含む場合
# fslvbm_1_bet -b  # 首を含まない場合

処理が完了すると「strucフォルダ」が生成され、その中に頭蓋除去後の画像が保存される。頭蓋除去後の画像は「*_brain.nii.gz」というファイル名で保存される。頭蓋除去前後の画像一覧をみるには、「struc/slicesdir/index.html」を開く。下地が頭蓋除去前の画像であり、赤い線が頭蓋除去後の画像である。

2.3. 3. fslvbm_2_template: 被験者の脳から灰白質テンプレート画像を生成

次に、脳の解剖学的標準化で用いるターゲット画像(標準脳)を、頭蓋除去済みの被験者脳から生成する。

まず最初に、標準脳を生成するために用いる被験者の脳画像リスト(template_list)を作成する。このとき、バイアスがかからないように健常者と患者の人数が同じになるようにリストを作る。

ここでは、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)の合計20名を「template_list」に記載した。

cat template_list
Con0001_T1w.nii.gz
Con0001_T1w.nii.gz
Con0001_T1w.nii.gz
...
Pat0010_T1w.nii.gz

template_listは、3D-T1WI(*_T1w.nii.gz)と同じ階層に配置する。

.
├── template_list   # ここに配置
├── Con0001_T1w.nii.gz
├── Con0002_T1w.nii.gz
├── Con0003_T1w.nii.gz
...
└── Pat0010_T1w.nii.gz

標準脳を生成するために用いる被験者の脳画像リスト(template_list)が準備できたら、頭蓋除去済みの脳画像を、灰白質・白質・脳脊髄液にセグメントする。次に、全ての被験者の灰白質(*_struc_GM.nii.gz)を元々ある灰白質標準脳(ICBM-152)にアフィン変換で位置合わせし、平均化する。このようにして得られた平均灰白質画像を、x軸に対して反転し再度平均化して、左右対称の平均灰白質画像(template_GM_init.nii.gz)を生成する。既に作成した「template_list」に記載されている灰白質画像を、先ほど作成した平均灰白質画像に位置合わせをし、平均化、さらに左右反転後に再度平均化する。最後に、元々ある灰白質標準脳(ICBM-152)に非線形位置合わせをし、2x2x2 mm^3の灰白質テンプレートを標準空間に生成する。

fslvbm_2_templateコマンドにはオプションがあり、各被験者の灰白質画像を元々ある灰白質標準脳(ICBM-152)に位置合わせする際に、アフィン変換を用いる場合は-aオプションを、非線形変換を用いる場合は-nオプションを指定する。

ここでは、fslvbm_2_template-aオプションを指定して実行する。

fslvbm_2_template -a  # Affine registration
# fslvbm_2_template -n  # non-linear registration

処理が完了するとstrucフォルダに灰白質テンプレート(template_GM_4D.nii.gz)が生成される。

2.4. fslvbm_3_proc: すべての被験者の灰白質を灰白質テンプレート画像に位置合わせし、ボクセルをモジュレーション後、平滑化

非線形変換を用いた位置合わせで、各被験者の灰白質画像を作成した灰白質テンプレート(template_GM_4D.nii.gz)に合わせ、statsフォルダに全被験者の灰白質4D画像を作る(GM_merg.nii.gz)。この時、被験者ごとの灰白質容積の違いを反映できるよう、非線形変換の際に収縮・拡大された度合いに応じて、灰白質画像のボクセル値を補正(モジュレーション)する。実際には、warp fieldのヤコビアンを灰白質画像の各ボクセルにかけ合わせる。その後、モジュレーションされた画像(GM_mod_merg.nii.gz)はσ=2, 3, 4mm (およそFWHM=2×2.3=4.6mmからFWHM=9mm)のガウシアンフィルタで平滑化する(GM_mod_merg_s[2,3,4].nii.gz)。

fslvbm_3_proc

以下は、モジュレーションした灰白質画像(GM_mod_merg.nii.gz)とσ=3mmのガウシアンフィルタで平滑化した灰白質画像(GM_mod_merg_s3.nii.gz)である。

2.5. randomise:GLMと並べ替え検定(permutation test)

ガウシアンフィルタで平滑化した灰白質画像を用いて、健常群と患者群の群間比較をする。

まず、GLMのデザインマトリックス(計画行列)とコントラストを設定する。

今回は、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)のデータがある。すべての被験者の3D-T1WIをlsコマンドでみると、先に健常者10名の3D-T1WI、次に患者10名の3D-T1WIが並んでいることが分かる。

ls |grep nii
Con0001_T1w.nii.gz
Con0002_T1w.nii.gz
Con0003_T1w.nii.gz
...
Pat0010_T1w.nii.gz

次に、GLMのデザインマトリックス(計画行列)とコントラストを決める設定ファイルを生成する。design_ttest2 <出力ファイル> <健常者数> <患者数>でコマンドを実行。

design_ttest2 stats/design 10 10

statsフォルダに、デザインマトリックス(design.mat)とコントラスト(design.con)が生成される。

デザインマトリックス(design.mat)の中身を確認。

/Matrixの一列目は健常者データであるかどうか、二列目は患者データであるかを0, 1で表している。行の順番は、被験者ファイルの順番(昇順)に対応する。したがって、これらは対応があるようにしておかなければならない。

cat stats/design.mat
/NumWaves 2
/NumPoints 20
/PPheights 1 1
/Matrix
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1

コントラスト(design.con)の中身を確認してみる。

/Matrix一列目は健常者の偏回帰係数、二列目は患者の偏回帰係数に対するもので、行は別々のコントラストである。この場合、一行目は健常者>患者の検定、二行目は健常者<患者の検定に相当する。

cat stats/design.con
/NumWaves 2
/NumContrasts 2
/PPheights 1 1
/Matrix
1 -1
-1 1

デザインマトリックス(計画行列)とコントラストの確認ができたら、randomiseコマンド使ってGLMと並べ替え検定(permutation test)を実行する。

randomiseコマンドの各オプションは、次の通り。

  • -i:入力画像
  • -m:マスク画像
  • -o :出力画像
  • -o :デザインマトリックス
  • -o :デザインコントラスト
  • -n:並べ替え検定の数
  • –T2:2D最適化を用いたTFCE
  • -x:voxel-wiseのcorrected P値マップ
  • –uncorrp:un-corrected P値マップ
  • -R:統計値マップ
fsleyes $FSLDIR/data/standard/MNI152_T1_2mm \
    stats/fslvbm_tfce_corrp_tstat1 \
    -cm Red-Yellow -dr 0.95,1

次のようなファイルが、生成される。

TFCEを用いた”健常群>患者群”の検定で、FWE補正をされたP値マップ(fslvbm_tfce_corrp_tstat1.nii.gz)を確認する。ここで、得られたP値マップは1-P値のマップであることに注意する。つまり、P<.05を有意とするのであれば、P値マップで0.95-1.00の値を見ればよい

fsleyes $FSLDIR/data/standard/MNI152_T1_2mm \
                 stats/fslvbm_tfce_corrp_tstat1 -cm Red-Yellow

3. おまけ

検定数が多い場合、有意差があったかどうかをすべて確認するのは大変である。そこで、一目で有意差があるかどうかを判断できるように、各検定ごとのP値マップの最大値を自動計測し、テキスト(corrected_P_report.txt)としてまとめる。

for PMAP in $(ls stats/ | grep tfce_corrp); do
    PMAX=$(fslstats stats/${PMAP} -R | cut -d " " -f2)
    echo ${PMAP} >>stats/tmp1.txt
    echo ${PMAX} >>stats/tmp2.txt
done
paste stats/tmp* >stats/tmp_corrected_P_report.txt
echo -e "$(cat stats/tmp_corrected_P_report.txt)\n\n\n$(cat stats/tmp_corrected_P_report.txt | sort -r -n -k 2)" \
    >stats/corrected_P_report.txt
rm stats/tmp*

上のコマンドを実行すると、statsフォルダに各検定とそのP値マップの最大値が記された「corrected_P_report.txt」が出力される。

検定結果を、ファイル名でソート(上段)したものと、P値でソートしたもの(下段)に分けて保存している。

cat stats/corrected_P_report.txt
fslvbm_tfce_corrp_tstat1.nii.gz 0.982391
fslvbm_tfce_corrp_tstat2.nii.gz 0.993181


fslvbm_tfce_corrp_tstat2.nii.gz 0.993181
fslvbm_tfce_corrp_tstat1.nii.gz 0.982391

【DIPY】DIPYを用いた拡散尖度イメージング: DKI



1. 目的
2. 準備
2.1. DIPYのインストール
2.2. 使用データ
2.3. 前処理
3. 拡散尖度イメージング(DKI)
3.1. 必要なパッケージをインポート
3.2. 画像およびMPG軸情報の読み込み
3.3. マスク画像の生成
3.4. モデルフィッティング
3.5. 拡散定量値の計算
3.6. NIfTI形式で保存
3.7. 結果
4. おまけ


1. 目的

  • DIPYを用いた拡散尖度イメージング: DKI

2. 準備

2.1. DIPYのインストール

pip3 install dipy

2.2. 使用データ

データを次のフォルダ構造で用意する。

Study/
└── Subject
    ├── DWI.nii.gz  # 拡散MRI
    ├── DWI_mask.nii.gz  # 拡散MRIマスク画像
    ├── bvals  # b-values
    └── bvecs  # b-vectors

2.3. 前処理

DKI(Diffusion Kurtosis Imaging)前に、拡散MRIの前処理をする。

  • 拡散MRIのノイズ除去(Software: MRtrix, DIPY)
  • ギブズのリンギングアーチファクト(Gibbs ringing)の除去(Software: MRtrix, DIPY)
  • 拡散MRIのバイアス(信号ムラ)補正(Software: MRtrix)
  • 拡散MRIの前処理 ~歪み・頭の動き・渦電流の補正(Software: FSL, MRtrix)

3. 拡散尖度イメージング(DKI)

Pythonで以下のコマンドを実行。

3.1. 必要なパッケージをインポート

import numpy as np
import matplotlib.pyplot as plt
import dipy.reconst.dki as dki
from dipy.core.gradients import gradient_table
from dipy.io.gradients import read_bvals_bvecs
from dipy.io.image import load_nifti, save_nifti
from dipy.segment.mask import median_otsu

3.2. 画像およびMPG軸情報の読み込み

DWI_FILE = 'DWI.nii.gz'
BVALS_FILE = 'bvals'
BVECS_FILE = 'bvecs'

data, affine = load_nifti(DWI_FILE)
bvals, bvecs = read_bvals_bvecs(BVALS_FILE, BVECS_FILE)
gtab = gradient_table(bvals, bvecs)

3.3. マスク画像の生成

median_otsu関数を用いて、b=0画像からマスク画像を生成する。vol_idxには、b0 volumeのvolume indexを渡す。

maskdata, mask = median_otsu(
    data, vol_idx=np.where(bvals == 0)[0])  # , dilate=3

3.4. モデルフィッティング

以下のコマンドで、DKIのモデルフィッティングを実行。

dkimodel = dki.DiffusionKurtosisModel(gtab)
dkifit = dkimodel.fit(maskdata)

3.5. 拡散定量値の計算

モデルフィッティングができたら、拡散定量値を算出する。

MK = dkifit.mk(0, 3)
AK = dkifit.ak(0, 3)
RK = dkifit.rk(0, 3)

脳周囲の背景では、フィッティングミスをしてnanとなる場合があるため、そのようなnanを0に置き換える。

MK[np.isnan(MK)] = 0
AK[np.isnan(AK)] = 0
RK[np.isnan(RK)] = 0

3.6. NIfTI形式で保存

save_nifti関数で、画像をNIfTI形式で保存する。

save_nifti('DWI_masked.nii.gz', maskdata.astype(np.float32), affine)
save_nifti('DWI_mask.nii.gz', mask.astype(np.float32), affine)
save_nifti('MK.nii.gz', MK.astype(np.float32), affine)
save_nifti('AK.nii.gz', AK.astype(np.float32), affine)
save_nifti('RK.nii.gz', RK.astype(np.float32), affine)

3.7. 結果

DKIによって算出された定量値画像は、以下の通り。

4. おまけ

DKIでもFA, MD, AD, RDを算出することができる。

# 拡散定量値を算出
FA = dkifit.fa
MD = dkifit.md
AD = dkifit.ad
RD = dkifit.rd

# nanを0に置換
FA[np.isnan(FA)] = 0
MD[np.isnan(MD)] = 0
AD[np.isnan(AD)] = 0
RD[np.isnan(RD)] = 0

# NIfTI形式で保存
save_nifti('FA.nii.gz', FA.astype(np.float32), affine)
save_nifti('MD.nii.gz', MD.astype(np.float32), affine)
save_nifti('AD.nii.gz', AD.astype(np.float32), affine)
save_nifti('RD.nii.gz', RD.astype(np.float32), affine)

以下は、FA, MD, AD, RDをDTI(下段)とDKI(上段)で比較した図である。

【AMICO】AMICOを用いた神経突起イメージング: NODDI



1. 目的
2. 準備
2.1. インストール
2.2. 使用データ
2.3. 前処理
3. 神経突起イメージング(NODDI)
3.1. AMICOのセットアップ
3.2. データの読み込み
3.3. 応答関数(response function)の計算
3.4. モデルフィッティング
3.5. 結果


1. 目的

  • AMICOを用いた神経突起イメージング: NODDI

2. 準備

2.1. インストール

Pythonを使って、AMICOを用いた神経突起イメージング(NODDI)をするために、以下のPythonパッケージをインストールする。

pip3 install dmri-amico

2.2. 使用データ

データを次のフォルダ構造で用意する。

Study/
└── Subject
    ├── DWI.nii.gz  # 拡散MRI
    ├── DWI_mask.nii.gz  # 拡散MRIマスク画像
    ├── bvals  # b-values
    └── bvecs  # b-vectors

2.3. 前処理

NODDI前に、拡散MRIの前処理をする。

  • 拡散MRIのノイズ除去(Software: MRtrix, DIPY)
  • ギブズのリンギングアーチファクト(Gibbs ringing)の除去(Software: MRtrix, DIPY)
  • 拡散MRIのバイアス(信号ムラ)補正(Software: MRtrix)
  • 拡散MRIの前処理 ~歪み・頭の動き・渦電流の補正(Software: FSL, MRtrix)

3. 神経突起イメージング(NODDI)

Pythonで以下のコマンドを実行。

3.1. AMICOのセットアップ

今回使用するファイル等の変数設定をする。

STUDY_DIR='Study'
SUBJECT_DIR='Subject'
DWI_FILE = 'DWI.nii.gz'
DWIMASK_FILE = 'DWI_mask.nii.gz'
BVALS_FILE = 'bvals'
BVECS_FILE = 'bvecs'

次に、使用するamicoパッケージのをインポートし、セットアップと初期化をする。

import amico
amico.core.setup()

3.2. データの読み込み

AMICOを実行するために、Study/Subjectフォルダ(PATH)を指定する。

ae = amico.Evaluation(STUDY_DIR, SUBJECT_DIR)

MPG軸情報(bvals/bvecs)の情報が入ったschemeファイルを生成する。

amico.util.fsl2scheme("{}/{}/{}".format(STUDY_DIR,SUBJECT_DIR,BVALS_FILE), "{}/{}/{}".format(STUDY_DIR,SUBJECT_DIR,BVECS_FILE),schemeFilename = "{}/{}/NODDI_protocol.scheme".format(STUDY_DIR,SUBJECT_DIR))
-> Writing scheme file to [ Study/Subject/NODDI_protocol.scheme ]
'Study/Subject/NODDI_protocol.scheme'

画像を読み込む。

ae.load_data(dwi_filename = DWI_FILE, scheme_filename = 'NODDI_protocol.scheme', mask_filename = DWIMASK_FILE, b0_thr = 0)
-> Loading data:
    * DWI signal
        - dim    = 130 x 130 x 82 x 129
        - pixdim = 1.769 x 1.769 x 1.800
    * Acquisition scheme
        - 129 samples, 2 shells
        - 1 @ b=0 , 64 @ b=1000.0 , 64 @ b=2000.0 
    * Binary mask
        - dim    = 130 x 130 x 82
        - pixdim = 1.769 x 1.769 x 1.800
        - voxels = 282878
   [ 4.4 seconds ]

-> Preprocessing:
    * Normalizing to b0... [ min=-3.28,  mean=0.25, max=22.86 ]
    * Keeping all b0 volume(s)
   [ 1.1 seconds ]

3.3. 応答関数(response function)の計算

NODDIモデルを設定して、応答関数(response function)を計算する。計算が完了するとkernelファイルが生成される。

ae.set_model('NODDI')
ae.generate_kernels()
-> Creating LUT for "NODDI" model:
   [==================================================] 100.0% 
   [ 373.3 seconds ]

作成したkernelファイルを読み込む。

ae.load_kernels()
-> Resampling LUT for subject "Subject":
   [==================================================] 100.0% 
   [ 112.8 seconds ]

3.4. モデルフィッティング

NODDIのモデルフィッティングを開始する。

ae.fit()
-> Fitting "NODDI" model to 282878 voxels:
   [==================================================] 100.0% 
   [ 02h 52m 07s ]

最後に、結果をNIfTIフォーマットで保存する。

ae.save_results()
-> Saving output to "AMICO/NODDI/*":
    - configuration  [OK]
    - FIT_dir.nii.gz  [OK]
    - FIT_ICVF.nii.gz  [OK]
    - FIT_OD.nii.gz  [OK]
    - FIT_ISOVF.nii.gz  [OK]
   [ DONE ]

3.5. 結果

結果は、「Study/Subject/AMICO/NODDI/」フォルダにある。

Study/Subject/AMICO/NODDI/
├── FIT_ICVF.nii.gz
├── FIT_ISOVF.nii.gz
├── FIT_OD.nii.gz
├── FIT_dir.nii.gz
└── config.pickle

画像はこちら。

【DIPY】DIPYを用いた拡散テンソルイメージング: DTI


1. 目的
2. 準備
2.1. DIPYのインストール
2.2. 使用データ
2.3. 前処理
3. 拡散テンソルイメージング(DTI)
3.1. 必要なパッケージをインポート
3.2. 画像およびMPG軸情報の読み込み
3.3. マスク画像の生成
3.4. モデルフィッティング
3.5. 拡散定量値の計算
3.6. NIfTI形式で保存
3.7. 結果


1. 目的

  • DIPYを用いた拡散テンソルイメージング: DTI

2. 準備

2.1. DIPYのインストール

pip3 install dipy

2.2. 使用データ

データを次のフォルダ構造で用意する。

Study/
└── Subject
    ├── DWI.nii.gz  # 拡散MRI
    ├── DWI_mask.nii.gz  # 拡散MRIマスク画像
    ├── bvals  # b-values
    └── bvecs  # b-vectors

2.3. 前処理

DTI(Diffusion Tensor Imaging)前に、拡散MRIの前処理をする。

  • 拡散MRIのノイズ除去(Software: MRtrix, DIPY)
  • ギブズのリンギングアーチファクト(Gibbs ringing)の除去(Software: MRtrix, DIPY)
  • 拡散MRIのバイアス(信号ムラ)補正(Software: MRtrix)
  • 拡散MRIの前処理 ~歪み・頭の動き・渦電流の補正(Software: FSL, MRtrix)

3. 拡散テンソルイメージング(DTI)

Pythonで以下のコマンドを実行。

3.1. 必要なパッケージをインポート

from dipy.segment.mask import median_otsu
import numpy as np
from dipy.io.image import load_nifti, save_nifti
from dipy.io.gradients import read_bvals_bvecs
from dipy.core.gradients import gradient_table
import dipy.reconst.dti as dti

3.2. 画像およびMPG軸情報の読み込み

DWI_FILE = 'DWI.nii.gz'
BVALS_FILE = 'bvals'
BVECS_FILE = 'bvecs'

data, affine = load_nifti(DWI_FILE)
bvals, bvecs = read_bvals_bvecs(BVALS_FILE, BVECS_FILE)
gtab = gradient_table(bvals, bvecs)

3.3. マスク画像の生成

median_otsu関数を用いて、b=0画像からマスク画像を生成する。vol_idxには、b0 volumeのvolume indexを渡す。

maskdata, mask = median_otsu(data, vol_idx=np.where(bvals == 0)[0])

3.4. モデルフィッティング

以下のコマンドで、DTIのモデルフィッティングを実行。

tenmodel = dti.TensorModel(gtab)
tenfit = tenmodel.fit(maskdata)

3.5. 拡散定量値の計算

モデルフィッティングができたら、拡散定量値を算出する。

FA = tenfit.fa
MD = tenfit.md
AD = tenfit.ad
RD = tenfit.rd
colour_FA = tenfit.color_fa

脳周囲の背景では、フィッティングミスをしてnanとなる場合があるため、そのようなnanを0に置き換える。

FA[np.isnan(FA)] = 0
MD[np.isnan(MD)] = 0
AD[np.isnan(AD)] = 0
RD[np.isnan(RD)] = 0
colour_FA[np.isnan(colour_FA)] = 0

3.6. NIfTI形式で保存

save_nifti関数で、画像をNIfTI形式で保存する。

save_nifti('DWI_masked.nii.gz', maskdata.astype(np.float32), affine)
save_nifti('DWI_mask.nii.gz', mask.astype(np.float32), affine)
save_nifti('FA.nii.gz', FA.astype(np.float32), affine)
save_nifti('MD.nii.gz', MD.astype(np.float32), affine)
save_nifti('AD.nii.gz', AD.astype(np.float32), affine)
save_nifti('RD.nii.gz', RD.astype(np.float32), affine)
save_nifti('colour_FA.nii.gz', colour_FA.astype(np.float32), affine)

3.7. 結果

DTIによって算出された定量値画像は、以下の通り。

【MRtrix】 MRtrixを用いた拡散テンソルイメージング: DTI


1. 目的
2. コマンド
3. 使用例
3.1. 前準備
3.2. テンソルの推定(コマンド:dwi2tensor
3.3. 拡散定量値の算出(コマンド:tensor2metric


1. 目的

  • MRtrixを用いた拡散テンソルイメージング: DTI

2. コマンド

MRtrixを用いて、拡散テンソルイメージング(DTI)をするには、dwi2tensortensor2metricコマンドを用いる。

dwi2tensor拡散MRI画像からテンソルを推定するコマンドで、tensor2metric推定したテンソルから拡散定量値を算出するコマンドである。

dwi2tensorのヘルプは、次の通り。

クリックして展開
SYNOPSIS

     Diffusion (kurtosis) tensor estimation

USAGE

     dwi2tensor [ options ] dwi dt

        dwi          the input dwi image.

        dt           the output dt image.


DESCRIPTION

     By default, the diffusion tensor (and optionally its kurtosis) is fitted
     to the log-signal in two steps: firstly, using weighted least-squares
     (WLS) with weights based on the empirical signal intensities; secondly, by
     further iterated weighted least-squares (IWLS) with weights determined by
     the signal predictions from the previous iteration (by default, 2
     iterations will be performed). This behaviour can be altered in two ways:

     * The -ols option will cause the first fitting step to be performed using
     ordinary least-squares (OLS); that is, all measurements contribute equally
     to the fit, instead of the default behaviour of weighting based on the
     empirical signal intensities.

     * The -iter option controls the number of iterations of the IWLS
     prodedure. If this is set to zero, then the output model parameters will
     be those resulting from the first fitting step only: either WLS by
     default, or OLS if the -ols option is used in conjunction with -iter 0.

     The tensor coefficients are stored in the output image as follows:
     volumes 0-5: D11, D22, D33, D12, D13, D23

     If diffusion kurtosis is estimated using the -dkt option, these are stored
     as follows:
     volumes 0-2: W1111, W2222, W3333
     volumes 3-8: W1112, W1113, W1222, W1333, W2223, W2333
     volumes 9-11: W1122, W1133, W2233
     volumes 12-14: W1123, W1223, W1233

OPTIONS

  -ols
     perform initial fit using an ordinary least-squares (OLS) fit (see
     Description).

  -mask image
     only perform computation within the specified binary brain mask image.

  -b0 image
     the output b0 image.

  -dkt image
     the output dkt image.

  -iter integer
     number of iterative reweightings for IWLS algorithm (default: 2) (see
     Description).

  -predicted_signal image
     the predicted dwi image.

DW gradient table import options

  -grad file
     Provide the diffusion-weighted gradient scheme used in the acquisition in
     a text file. This should be supplied as a 4xN text file with each line is
     in the format [ X Y Z b ], where [ X Y Z ] describe the direction of the
     applied gradient, and b gives the b-value in units of s/mm^2. If a
     diffusion gradient scheme is present in the input image header, the data
     provided with this option will be instead used.

  -fslgrad bvecs bvals
     Provide the diffusion-weighted gradient scheme used in the acquisition in
     FSL bvecs/bvals format files. If a diffusion gradient scheme is present in
     the input image header, the data provided with this option will be instead
     used.

Standard options

  -info
     display information messages.

  -quiet
     do not display information messages or progress status; alternatively,
     this can be achieved by setting the MRTRIX_QUIET environment variable to a
     non-empty string.

  -debug
     display debugging messages.

  -force
     force overwrite of output files (caution: using the same file as input and
     output might cause unexpected behaviour).

  -nthreads number
     use this number of threads in multi-threaded applications (set to 0 to
     disable multi-threading).

  -config key value  (multiple uses permitted)
     temporarily set the value of an MRtrix config file entry.

  -help
     display this information page and exit.

  -version
     display version information and exit.

tensor2metricのヘルプは、次の通り。

クリックして展開
SYNOPSIS

     Generate maps of tensor-derived parameters

USAGE

     tensor2metric [ options ] tensor

        tensor       the input tensor image.


OPTIONS

  -adc image
     compute the mean apparent diffusion coefficient (ADC) of the diffusion
     tensor. (sometimes also referred to as the mean diffusivity (MD))

  -fa image
     compute the fractional anisotropy (FA) of the diffusion tensor.

  -ad image
     compute the axial diffusivity (AD) of the diffusion tensor. (equivalent to
     the principal eigenvalue)

  -rd image
     compute the radial diffusivity (RD) of the diffusion tensor. (equivalent
     to the mean of the two non-principal eigenvalues)

  -cl image
     compute the linearity metric of the diffusion tensor. (one of the three
     Westin shape metrics)

  -cp image
     compute the planarity metric of the diffusion tensor. (one of the three
     Westin shape metrics)

  -cs image
     compute the sphericity metric of the diffusion tensor. (one of the three
     Westin shape metrics)

  -value image
     compute the selected eigenvalue(s) of the diffusion tensor.

  -vector image
     compute the selected eigenvector(s) of the diffusion tensor.

  -num sequence
     specify the desired eigenvalue/eigenvector(s). Note that several
     eigenvalues can be specified as a number sequence. For example, '1,3'
     specifies the principal (1) and minor (3) eigenvalues/eigenvectors
     (default = 1).

  -modulate choice
     specify how to modulate the magnitude of the eigenvectors. Valid choices
     are: none, FA, eigval (default = FA).

  -mask image
     only perform computation within the specified binary brain mask image.

Standard options

  -info
     display information messages.

  -quiet
     do not display information messages or progress status; alternatively,
     this can be achieved by setting the MRTRIX_QUIET environment variable to a
     non-empty string.

  -debug
     display debugging messages.

  -force
     force overwrite of output files (caution: using the same file as input and
     output might cause unexpected behaviour).

  -nthreads number
     use this number of threads in multi-threaded applications (set to 0 to
     disable multi-threading).

  -config key value  (multiple uses permitted)
     temporarily set the value of an MRtrix config file entry.

  -help
     display this information page and exit.

  -version
     display version information and exit.

DTI拡散定量値(FA, MD, AD, RD, カラーFA)を計算するための基本的な使い方は、以下の通り。

dwi2tensor <入力画像> <出力画像>
tensor2metric -fa <出力画像> -adc <出力画像> -ad <出力画像> -rd <出力画像> -vec <出力画像> tensor.mif

3. 使用例

3.1. 前準備

まず、次のファイルを用意する。

.
├── DWI.nii.gz  # 拡散MRI
├── DWI_mask.nii.gz
├── bvals  # b-values
├── bvecs  # b-vectors
└── headers.json  # ヘッダー情報の入ったJSONファイル

こちらの記事を参考に、拡散MRI(DWI.nii.gz)とそのMPG軸情報(bvecs, bvals)とヘッダー情報(headers.json)をまとめて、MIF形式(DWI.mif)に変換する。

mrconvert -fslgrad bvecs bvals -json_import headers.json DWI.nii.gz DWI.mif

3.2. テンソルの推定(コマンド:dwi2tensor

ファイルの用意ができたら、dwi2tensorを次のように実行する

dwi2tensor DWI.mif tensor.mif

mrinfoを使って「tensor.mif」の情報を確認すると、6 volumesのデータであることが分かる。

mrinfo tensor.mif
************************************************
Image name:          "tensor.mif"
************************************************
  Dimensions:        130 x 130 x 82 x 6
  Voxel size:        1.76923 x 1.76923 x 1.8 x 1
  Data strides:      [ -1 2 3 4 ]
  Format:            MRtrix
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                    1           0           0        -109
                               -0           1           0      -103.7
                               -0           0           1      -58.57

それぞれのボリュームは、各方向の拡散係数に相当する。

The tensor coefficients are stored in the output image as follows:
volumes 0-5: D11, D22, D33, D12, D13, D23

3.3. 拡散定量値の算出(コマンド:tensor2metric

先程推定した、「tensor.mif」を使って拡散定量値を算出する。

tensor2metric -fa FA.mif -adc MD.mif -ad AD.mif -rd RD.mif -vec color_FA.mif tensor.mif

DTIの各拡散定量値画像は、以下。

【FSL】FSLを用いた拡散テンソルイメージング: DTI


1. 目的
2. コマンド
3. 使用例


1. 目的

  • FSLを用いた拡散テンソルイメージング: DTI

2. コマンド

FSLを用いて、拡散テンソルイメージング(DTI)をするには、dtifitコマンドを用いる。

dtifitのヘルプは、次の通り。

クリックして展開
Usage: 
dtifit -k <filename>
 dtifit --verbose


Compulsory arguments (You MUST set one or more of):
	-k,--data	dti data file
	-o,--out	Output basename
	-m,--mask	Bet binary mask file
	-r,--bvecs	b vectors file
	-b,--bvals	b values file

Optional arguments (You may optionally specify one or more of):
	-V,--verbose	switch on diagnostic messages
	-h,--help	display this message
	--cni	Input confound regressors
	--sse	Output sum of squared errors
	-w,--wls	Fit the tensor with weighted least squares
	--kurt	Output mean kurtosis map (for multi-shell data)
	--kurtdir	Output  parallel/perpendicular kurtosis maps (for multi-shell data)
	--littlebit	Only process small area of brain
	--save_tensor	Save the elements of the tensor
	-z,--zmin	min z
	-Z,--zmax	max z
	-y,--ymin	min y
	-Y,--ymax	max y
	-x,--xmin	min x
	-X,--xmax	max x
	--gradnonlin	Gradient Nonlinearity Tensor file

基本的な使い方は、以下の通り。

dtifit -k <DWI images> -o <OUTPUT> -m <MASK> -r <b-vectors file> -b <b-values file>

3. 使用例

まず、次のファイルを用意する。

.
├── DWI.nii.gz  # 拡散MRI
├── DWI_mask.nii.gz  # 拡散MRIのマスク画像
├── bvals  # b-values
└── bvecs  # b-vectors

ファイルの用意ができたら、dtifitを次のように実行する

dtifit -k DWI.nii.gz -o output -m DWI_mask.nii.gz -r bvecs -b bvals

処理が終わると次のファイルが生成される。

  • output_V1 – 1st eigenvector (固有ベクトル)
  • output_V2 – 2nd eigenvector
  • output_V3 – 3rd eigenvector
  • output_L1 – 1st eigenvalue (固有値)
  • output_L2 – 2nd eigenvalue
  • output_L3 – 3rd eigenvalue
  • output_MD – mean diffusivity
  • output_FA – fractional anisotropy (isotropic ~ 0; stick-like ~1)
  • output_MO – mode of the anisotropy (oblate ~ -1; isotropic ~ 0; prolate ~ 1)
  • output_S0 – raw T2 signal with no diffusion weighting

FA, MD, L1の画像は、以下。

【MRtrix】MRtrixを用いた拡散MRIの前処理 ~歪み・頭の動き・渦電流の補正~


1. 目的
2. コマンド
3. 使用例
3.1. 前準備
3.2. 歪み補正と頭の動き補正


1. 目的

  • MRtrixを用いた拡散MRIの前処理(歪み・頭の動き・渦電流の補正)

2. コマンド

MRtrixを用いて、拡散MRIの歪み・頭の動き・渦電流を補正するには、dwifslpreprocを用いる(古いMRtrixバージョンではdwipreproc)。

dwipreprocは、FSLtopupeddyを用いるので、前もってFSLをインストールしておく必要がある。

dwifslpreprocのヘルプは、次の通り。

クリックして展開
SYNOPSIS

     Perform diffusion image pre-processing using FSL's eddy tool; including
     inhomogeneity distortion correction using FSL's topup tool if possible

USAGE

     dwifslpreproc [ options ] input output

        input        The input DWI series to be corrected

        output       The output corrected image series

DESCRIPTION

     This script is intended to provide convenience of use of the FSL software
     tools topup and eddy for performing DWI pre-processing, by encapsulating
     some of the surrounding image data and metadata processing steps. It is
     intended to simply these processing steps for most commonly-used DWI
     acquisition strategies, whilst also providing support for some more exotic
     acquisitions. The "example usage" section demonstrates the ways in which
     the script can be used based on the (compulsory) -rpe_* command-line
     options.

     The "-topup_options" and "-eddy_options" command-line options allow the
     user to pass desired command-line options directly to the FSL commands
     topup and eddy. The available options for those commands may vary between
     versions of FSL; users can interrogate such by querying the help pages of
     the installed software, and/or the FSL online documentation: (topup)
     https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup/TopupUsersGuide ; (eddy)
     https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy/UsersGuide

     The script will attempt to run the CUDA version of eddy; if this does not
     succeed for any reason, or is not present on the system, the CPU version
     will be attempted instead. By default, the CUDA eddy binary found that
     indicates compilation against the most recent version of CUDA will be
     attempted; this can be over-ridden by providing a soft-link "eddy_cuda"
     within your path that links to the binary you wish to be executed.

     Note that this script does not perform any explicit registration between
     images provided to topup via the -se_epi option, and the DWI volumes
     provided to eddy. In some instances (motion between acquisitions) this can
     result in erroneous application of the inhomogeneity field during
     distortion correction. Use of the -align_seepi option is advocated in this
     scenario, which ensures that the first volume in the series provided to
     eddy is also the first volume in the series provided to eddy, guaranteeing
     alignment. But a prerequisite for this approach is that the image contrast
     within the images provided to the -se_epi option must match the b=0 volumes
     present within the input DWI series: this means equivalent TE, TR and flip
     angle (note that differences in multi-band factors between two acquisitions
     may lead to differences in TR).

EXAMPLE USAGES

     A basic DWI acquisition, where all image volumes are acquired in a single
     protocol with fixed phase encoding:
       $ dwifslpreproc DWI_in.mif DWI_out.mif -rpe_none -pe_dir ap -readout_time 0.55
     Due to use of a single fixed phase encoding, no EPI distortion correction
     can be applied in this case.

     DWIs all acquired with a single fixed phase encoding; but additionally a
     pair of b=0 images with reversed phase encoding to estimate the
     inhomogeneity field:
       $ mrcat b0_ap.mif b0_pa.mif b0_pair.mif -axis 3; dwifslpreproc DWI_in.mif DWI_out.mif -rpe_pair -se_epi b0_pair.mif -pe_dir ap -readout_time 0.72 -align_seepi
     Here the two individual b=0 volumes are concatenated into a single 4D image
     series, and this is provided to the script via the -se_epi option. Note
     that with the -rpe_pair option used here, which indicates that the SE-EPI
     image series contains one or more pairs of b=0 images with reversed phase
     encoding, the FIRST HALF of the volumes in the SE-EPI series must possess
     the same phase encoding as the input DWI series, while the second half are
     assumed to contain the opposite phase encoding direction but identical
     total readout time. Use of the -align_seepi option is advocated as long as
     its use is valid (more information in the Description section).

     All DWI directions & b-values are acquired twice, with the phase encoding
     direction of the second acquisition protocol being reversed with respect to
     the first:
       $ mrcat DWI_lr.mif DWI_rl.mif DWI_all.mif -axis 3; dwifslpreproc DWI_all.mif DWI_out.mif -rpe_all -pe_dir lr -readout_time 0.66
     Here the two acquisition protocols are concatenated into a single DWI
     series containing all acquired volumes. The direction indicated via the
     -pe_dir option should be the direction of phase encoding used in
     acquisition of the FIRST HALF of volumes in the input DWI series; ie. the
     first of the two files that was provided to the mrcat command. In this
     usage scenario, the output DWI series will contain the same number of image
     volumes as ONE of the acquired DWI series (ie. half of the number in the
     concatenated series); this is because the script will identify pairs of
     volumes that possess the same diffusion sensitisation but reversed phase
     encoding, and perform explicit recombination of those volume pairs in such
     a way that image contrast in regions of inhomogeneity is determined from
     the stretched rather than the compressed image.

     Any acquisition scheme that does not fall into one of the example usages
     above:
       $ mrcat DWI_*.mif DWI_all.mif -axis 3; mrcat b0_*.mif b0_all.mif -axis 3; dwifslpreproc DWI_all.mif DWI_out.mif -rpe_header -se_epi b0_all.mif -align_seepi
     With this usage, the relevant phase encoding information is determined
     entirely based on the contents of the relevant image headers, and
     dwifslpreproc prepares all metadata for the executed FSL commands
     accordingly. This can therefore be used if the particular DWI acquisition
     strategy used does not correspond to one of the simple examples as
     described in the prior examples. This usage is predicated on the headers of
     the input files containing appropriately-named key-value fields such that
     MRtrix3 tools identify them as such. In some cases, conversion from DICOM
     using MRtrix3 commands will automatically extract and embed this
     information; however this is not true for all scanner vendors and/or
     software versions. In the latter case it may be possible to manually
     provide these metadata; either using the -json_import command-line option
     of dwifslpreproc, or the -json_import or one of the -import_pe_* command-
     line options of MRtrix3's mrconvert command (and saving in .mif format)
     prior to running dwifslpreproc.

OPTIONS

  -pe_dir PE
     Manually specify the phase encoding direction of the input series; can be a
     signed axis number (e.g. -0, 1, +2), an axis designator (e.g. RL, PA, IS),
     or NIfTI axis codes (e.g. i-, j, k)

  -readout_time time
     Manually specify the total readout time of the input series (in seconds)

  -se_epi image
     Provide an additional image series consisting of spin-echo EPI images,
     which is to be used exclusively by topup for estimating the inhomogeneity
     field (i.e. it will not form part of the output image series)

  -align_seepi
     Achieve alignment between the SE-EPI images used for inhomogeneity field
     estimation, and the DWIs (more information in Description section)

  -json_import file
     Import image header information from an associated JSON file (may be
     necessary to determine phase encoding information)

  -topup_options " TopupOptions"
     Manually provide additional command-line options to the topup command
     (provide a string within quotation marks that contains at least one space,
     even if only passing a single command-line option to topup)

  -eddy_options " EddyOptions"
     Manually provide additional command-line options to the eddy command
     (provide a string within quotation marks that contains at least one space,
     even if only passing a single command-line option to eddy)

  -eddy_mask image
     Provide a processing mask to use for eddy, instead of having dwifslpreproc
     generate one internally using dwi2mask

  -eddy_slspec file
     Provide a file containing slice groupings for eddy's slice-to-volume
     registration

  -eddyqc_text directory
     Copy the various text-based statistical outputs generated by eddy, and the
     output of eddy_qc (if installed), into an output directory

  -eddyqc_all directory
     Copy ALL outputs generated by eddy (including images), and the output of
     eddy_qc (if installed), into an output directory

Options for specifying the acquisition phase-encoding design; note that one of the -rpe_* options MUST be provided

  -rpe_none
     Specify that no reversed phase-encoding image data is being provided; eddy
     will perform eddy current and motion correction only

  -rpe_pair
     Specify that a set of images (typically b=0 volumes) will be provided for
     use in inhomogeneity field estimation only (using the -se_epi option)

  -rpe_all
     Specify that ALL DWIs have been acquired with opposing phase-encoding

  -rpe_header
     Specify that the phase-encoding information can be found in the image
     header(s), and that this is the information that the script should use

Options for importing the diffusion gradient table

  -grad GRAD
     Provide the diffusion gradient table in MRtrix format

  -fslgrad bvecs bvals
     Provide the diffusion gradient table in FSL bvecs/bvals format

Options for exporting the diffusion gradient table

  -export_grad_mrtrix grad
     Export the final gradient table in MRtrix format

  -export_grad_fsl bvecs bvals
     Export the final gradient table in FSL bvecs/bvals format

Additional standard options for Python scripts

  -nocleanup
     do not delete intermediate files during script execution, and do not delete
     scratch directory at script completion.

  -scratch /path/to/scratch/
     manually specify the path in which to generate the scratch directory.

  -continue <ScratchDir> <LastFile>
     continue the script from a previous execution; must provide the scratch
     directory path, and the name of the last successfully-generated file.

Standard options

  -info
     display information messages.

  -quiet
     do not display information messages or progress status. Alternatively, this
     can be achieved by setting the MRTRIX_QUIET environment variable to a non-
     empty string.

  -debug
     display debugging messages.

  -force
     force overwrite of output files.

  -nthreads number
     use this number of threads in multi-threaded applications (set to 0 to
     disable multi-threading).

  -config key value  (multiple uses permitted)
     temporarily set the value of an MRtrix config file entry.

  -help
     display this information page and exit.

  -version
     display version information and exit.

基本的な使い方は、以下の通り。

dwifslpreproc <入力画像> <出力画像> [オプション]

3. 使用例

3.1. 前準備

位相エンコード方向を、APとPAそれぞれで撮像したデータがあったとする。DICOM形式からNIfTI形式に変換する方法は、以下の記事を参考にするとよい。

.
├── DWI_AP.nii.gz  # DW images (PE: AP)
├── DWI_PA.nii.gz  # DW images (PE: PA)
├── bvals_AP  # b-values (PE: AP)
├── bvals_PA  # b-values (PE: PA)
├── bvecs_AP  # b-vectors (PE: AP)
├── bvecs_PA  # b-vectors (PE: PA)
├── headers_AP.json  # DICOM headers (PE: AP)
└── headers_PA.json  # DICOM headers (PE: PA)

まず、こちらの記事を参考に、拡散MRI(DWI.nii.gz)とそのMPG軸情報(bvecs, bvals)とヘッダー情報(headers.json)をまとめて、MIF形式(DWI.mif)に変換する。

mrconvert -fslgrad bvecs_AP bvals_AP -json_import headers_AP.json DWI_AP.nii.gz DWI_AP.mif  # PE: AP
mrconvert -fslgrad bvecs_PA bvals_PA -json_import headers_PA.json DWI_PA.nii.gz DWI_PA.mif  # PE: PA

次に、mrcatを使ってDWI_AP.mifとDWI_PA.mifをひとつの画像(DWI_all.mif)にまとめる。

オプションの-axis 3は、4次元目のt軸(Volume)方向にまとめるという意味である(MRtrixではAxisを0から数える [i.e., x: 0, y: 1, z: 2, t: 3])。mrcatの詳細は、こちら。

mrcat DWI_AP.mif DWI_PA.mif DWI_all.mif -axis 3

次に、dwiextractを用いて、b=0のみを抽出する。dwiextractの詳細は、こちら。

dwiextract -bzero DWI_all.mif DWI_b0.mif

3.2. 歪み補正と頭の動き補正

歪み補正と頭の動き補正をするために、次のコマンドを実行する。

ここで使用した、各オプションは以下。

  • -rpe_header:位相エンコード情報を読み込む
  • -se_epi:b=0(spin-echo EPI images)を指定
  • -align_seepi:磁場の不均一性場の推定で用いられる、SE-EPI画像とDWIの間の位置合わせを実行

dwifslpreproc DWI_all.mif DWI_preproc.mif -rpe_header -se_epi DWI_b0.mif -align_seepi

歪み補正後の画像は、以下。

頭の動き補正後の画像は、以下。


【MRtrix】MRtrixを用いた5TT(five-tissue-type)画像の生成


1. 目的
2. コマンド
3. 使用例
3.1. FSLアルゴリズムを用いる場合
3.2. FreeSurferアルゴリズムを用いる場合
4. 結果


1. 目的

  • MRtrixを用いた5TT(five-tissue-type)画像の生成

2. コマンド

MRtrixの5ttgenを用いて、次の5つの組織(five-tissue-type: 5TT)画像を生成する。

  1. Cortical grey matter
  2. Sub-cortical grey matter
  3. White matter
  4. CSF
  5. Pathological tissue

5ttgenのヘルプは、次の通り。

クリックして展開
SYNOPSIS

     Generate a 5TT image suitable for ACT

USAGE

     5ttgen [ options ] algorithm ...

        algorithm    Select the algorithm to be used to complete the script operation;
                     additional details and options become available once an
                     algorithm is nominated. Options are: freesurfer, fsl, gif,
                     hsvs

DESCRIPTION

     5ttgen acts as a 'master' script for generating a five-tissue-type (5TT)
     segmented tissue image suitable for use in Anatomically-Constrained
     Tractography (ACT). A range of different algorithms are available for
     completing this task. When using this script, the name of the algorithm to
     be used must appear as the first argument on the command-line after
     '5ttgen'. The subsequent compulsory arguments and options available depend
     on the particular algorithm being invoked.

     Each algorithm available also has its own help page, including necessary
     references; e.g. to see the help page of the 'fsl' algorithm, type '5ttgen
     fsl'.

Options common to all 5ttgen algorithms

  -nocrop
     Do NOT crop the resulting 5TT image to reduce its size (keep the same
     dimensions as the input image)

  -sgm_amyg_hipp
     Represent the amygdalae and hippocampi as sub-cortical grey matter in the
     5TT image

Additional standard options for Python scripts

  -nocleanup
     do not delete intermediate files during script execution, and do not delete
     scratch directory at script completion.

  -scratch /path/to/scratch/
     manually specify the path in which to generate the scratch directory.

  -continue <ScratchDir> <LastFile>
     continue the script from a previous execution; must provide the scratch
     directory path, and the name of the last successfully-generated file.

Standard options

  -info
     display information messages.

  -quiet
     do not display information messages or progress status. Alternatively, this
     can be achieved by setting the MRTRIX_QUIET environment variable to a non-
     empty string.

  -debug
     display debugging messages.

  -force
     force overwrite of output files.

  -nthreads number
     use this number of threads in multi-threaded applications (set to 0 to
     disable multi-threading).

  -config key value  (multiple uses permitted)
     temporarily set the value of an MRtrix config file entry.

  -help
     display this information page and exit.

  -version
     display version information and exit.

基本的な使い方は、以下の通り。5ttgenのアルゴリズムは、freesurfer, fsl, gif, hsvsがあるが、ここではfreesurferとfslのアルゴリズムについて使い方を解説する。

5ttgen [アルゴリズム] <入力画像> <出力画像>

3. 使用例

3.1. FSLアルゴリズムを用いる場合

FSLアルゴリズムを用いる場合、3D-T1WI(T1w.nii.gz)が必要となる。また、オプションとして3D-T2WIも入力することができる。

5ttgen fsl T1w.nii.gz 5tt.nii.gz

3.2. FreeSurferアルゴリズムを用いる場合

FreeSurferアルゴリズムを用いる場合、Freesurferの生成ファイルであるaparc+aseg.mgz(asegとついたファイル)が必要となる。

FreeSurferの使い方は、こちらの記事を参考にするとよい。

aparc+aseg.mgzが準備できたら、以下のコマンドを実行する。

5ttgen freesurfer aparc+aseg.mgz 5tt.nii.gz

4. 結果

5ttgenで生成された画像は、5ボリュームデータであり、各ボリュームと対応する組織は次の通り。

  1. Cortical grey matter
  2. Sub-cortical grey matter
  3. White matter
  4. CSF
  5. Pathological tissue

以下に、FSLとFreeSurferのアルゴリズムを用いて5ttgenした結果(緑)を示す。

【MRtrix】MRtrixを用いた拡散MRIのバイアス(信号ムラ)補正


1. 目的
2. コマンド
3. 使用例
3.1. 前準備
3.2. 拡散MRIのバイアス(信号ムラ)補正


1. 目的

  • MRtrixを用いた拡散MRIのバイアス(信号ムラ)補正

2. コマンド

MRtrixを用いて拡散MRIのバイアス(信号ムラ)補正をするには、dwibiascorrectを使用する。

ここでは、ANTsのN4アルゴリズムを用いたバイアス補正を紹介する。ANTsアルゴリズムを使用する場合は、ANTsを前もってインストールしておく必要がある。

dwibiascorrectのヘルプは、次の通り。

クリックして展開
SYNOPSIS

     Perform B1 field inhomogeneity correction for a DWI volume series

USAGE

     dwibiascorrect [ options ] algorithm ...

        algorithm    Select the algorithm to be used to complete the script operation;
                     additional details and options become available once an
                     algorithm is nominated. Options are: ants, fsl

Options for importing the diffusion gradient table

  -grad GRAD
     Provide the diffusion gradient table in MRtrix format

  -fslgrad bvecs bvals
     Provide the diffusion gradient table in FSL bvecs/bvals format

Options common to all dwibiascorrect algorithms

  -mask image
     Manually provide a mask image for bias field estimation

  -bias image
     Output the estimated bias field

Additional standard options for Python scripts

  -nocleanup
     do not delete intermediate files during script execution, and do not delete
     scratch directory at script completion.

  -scratch /path/to/scratch/
     manually specify the path in which to generate the scratch directory.

  -continue <ScratchDir> <LastFile>
     continue the script from a previous execution; must provide the scratch
     directory path, and the name of the last successfully-generated file.

Standard options

  -info
     display information messages.

  -quiet
     do not display information messages or progress status. Alternatively, this
     can be achieved by setting the MRTRIX_QUIET environment variable to a non-
     empty string.

  -debug
     display debugging messages.

  -force
     force overwrite of output files.

  -nthreads number
     use this number of threads in multi-threaded applications (set to 0 to
     disable multi-threading).

  -config key value  (multiple uses permitted)
     temporarily set the value of an MRtrix config file entry.

  -help
     display this information page and exit.

  -version
     display version information and exit.

基本的な使い方は、以下の通り。

dwibiascorrect ants <入力画像> <出力画像>

3.使用例

3.1.前準備

まず、こちらの記事を参考に、拡散MRI(DWI.nii.gz)とそのMPG軸情報(bvecs, bvals)とヘッダー情報(headers.json)をまとめて、MIF形式(DWI.mif)に変換する。

mrconvert -fslgrad bvecs bvals -json_import headers.json DWI.nii.gz DWI.mif

3.2.拡散MRIのバイアス(信号ムラ)補正

以下のコマンドを実行する。-biasオプションを指定することで、バイアスフィールドを出力することができる。

dwibiascorrect ants DWI.mif DWI_unbiased.mif -bias bias.mif

補正後の画像は、以下。

【DIPY】DIPYを用いた拡散MRIのノイズ除去 ~Denoise~


1. 目的
2. 準備
2.1. DIPYのインストール
2.2. 使用データ
3. 拡散MRIのノイズ除去
3.1. 必要なパッケージをインポート
3.2. 画像およびMPG軸情報の読み込み
3.3. ノイズ除去(デノイズ)
3.4. NIfTI形式で保存
3.5. 結果
4. おまけ


1. 目的

  • DIPYを用いた拡散MRIのノイズ除去 ~Denoise~

2. 準備

2.1. DIPYのインストール

pip3 install dipy

2.2. 使用データ

データを次のフォルダ構造で用意する。

Study/
└── Subject
    ├── DWI.nii.gz  # 拡散MRI
    ├── DWI_mask.nii.gz  # 拡散MRIマスク画像
    ├── bvals  # b-values
    └── bvecs  # b-vectors

3. 拡散MRIのノイズ除去

Pythonで以下のコマンドを実行。

3.1. 必要なパッケージをインポート

import numpy as np
import matplotlib.pyplot as plt
from time import time
from dipy.denoise.localpca import mppca
from dipy.core.gradients import gradient_table
from dipy.io.image import load_nifti, save_nifti
from dipy.io.gradients import read_bvals_bvecs
from dipy.segment.mask import median_otsu

3.2. 画像およびMPG軸情報の読み込み

DWI_FILE = 'DWI.nii.gz'
BVALS_FILE = 'bvals'
BVECS_FILE = 'bvecs'

data, affine = load_nifti(DWI_FILE)
bvals, bvecs = read_bvals_bvecs(BVALS_FILE, BVECS_FILE)
gtab = gradient_table(bvals, bvecs)

3.3. ノイズ除去(デノイズ)

mppca関数を用いて、Marchenko-Pastur PCAを用いたデノイズをする。

denoised_arr = mppca(data, patch_radius=3)

3.4. NIfTI形式で保存

save_nifti関数で、画像をNIfTI形式で保存する。

save_nifti('DWI_denoised.nii.gz', denoised_arr.astype(np.float32), affine)

3.5. 結果

拡散強調像(b=2000 s/mm^2)のデノイズ前後の比較と差分画像は、以下。

実際に、デノイズ前(上段)とデノイズ後(下段)でDTIおよびDKIを計算し、比較してみる。

4. おまけ

Marchenko-Pastur PCAアルゴリズムは、ノイズの標準偏差も推定することができる。ノイズの標準偏差を算出するためには、return_sigmaのフラグを「True」にする。

denoised_arr, sigma = mppca(data, patch_radius=3, return_sigma=True)
save_nifti('DWI_noise_sigma.nii.gz', sigma .astype(np.float32), affine)

ノイズの標準偏差マップは、以下の通り。

脳領域における平均ノイズ標準偏差は、次のようにして計測できる。

mean_sigma = np.mean(sigma[mask])
print(mean_sigma)

推定した平均ノイズ標準偏差を用いて、b=0 (s/mm^2)画像のSNRを算出する。

b0 = denoised_arr[..., 0]
mean_signal = np.mean(b0[mask])
snr = mean_signal / mean_sigma
print(snr)

【MRtrix】MRtrixを用いた拡散MRIのノイズ除去 ~Denoise~


1. 目的
2. コマンド
2.1. 使用例


1. 目的

  • MRtrixを用いた拡散MRIのノイズ除去(Denoise)

2. コマンド

拡散MRIのノイズ除去には、MRtrixdwidenoiseを用いる。dwidenoiseは、Marchenko-Pastur PCAを用いたデノイズである。

拡散MRIのノイズ除去は、前処理の一番最初に実行する必要がある。

dwidenoiseのヘルプは、次の通り。

クリックして展開
SYNOPSIS

     dMRI noise level estimation and denoising using Marchenko-Pastur PCA

USAGE

     dwidenoise [ options ] dwi out

        dwi          the input diffusion-weighted image.

        out          the output denoised DWI image.


DESCRIPTION

     DWI data denoising and noise map estimation by exploiting data redundancy
     in the PCA domain using the prior knowledge that the eigenspectrum of
     random covariance matrices is described by the universal Marchenko-Pastur
     (MP) distribution. Fitting the MP distribution to the spectrum of
     patch-wise signal matrices hence provides an estimator of the noise level
     'sigma', as was first shown in Veraart et al. (2016) and later improved in
     Cordero-Grande et al. (2019). This noise level estimate then determines
     the optimal cut-off for PCA denoising.

     Important note: image denoising must be performed as the first step of the
     image processing pipeline. The routine will fail if interpolation or
     smoothing has been applied to the data prior to denoising.

     Note that this function does not correct for non-Gaussian noise biases
     present in magnitude-reconstructed MRI images. If available, including the
     MRI phase data can reduce such non-Gaussian biases, and the command now
     supports complex input data.

OPTIONS

  -mask image
     Only process voxels within the specified binary brain mask image.

  -extent window
     Set the patch size of the denoising filter. By default, the command will
     select the smallest isotropic patch size that exceeds the number of DW
     images in the input data, e.g., 5x5x5 for data with <= 125 DWI volumes,
     7x7x7 for data with <= 343 DWI volumes, etc.

  -noise level
     The output noise map, i.e., the estimated noise level 'sigma' in the data.
     Note that on complex input data, this will be the total noise level across
     real and imaginary channels, so a scale factor sqrt(2) applies.

  -datatype float32/float64
     Datatype for the eigenvalue decomposition (single or double precision).
     For complex input data, this will select complex float32 or complex
     float64 datatypes.

  -estimator Exp1/Exp2
     Select the noise level estimator (default = Exp2), either: 
     * Exp1: the original estimator used in Veraart et al. (2016), or 
     * Exp2: the improved estimator introduced in Cordero-Grande et al. (2019).

Standard options

  -info
     display information messages.

  -quiet
     do not display information messages or progress status; alternatively,
     this can be achieved by setting the MRTRIX_QUIET environment variable to a
     non-empty string.

  -debug
     display debugging messages.

  -force
     force overwrite of output files (caution: using the same file as input and
     output might cause unexpected behaviour).

  -nthreads number
     use this number of threads in multi-threaded applications (set to 0 to
     disable multi-threading).

  -config key value  (multiple uses permitted)
     temporarily set the value of an MRtrix config file entry.

  -help
     display this information page and exit.

  -version
     display version information and exit.

基本的な使い方は、次の通り。

dwidenoise <入力画像> <出力画像>

2.1. 使用例

前処理する前の拡散MRI(DWI.nii.gz)に、dwidenoiseを実行する。

dwidenoise DWI.nii.gz DWI_denoised.nii.gz

処理後の画像は、以下。

【FSL】FSLを用いた構造MRIと拡散MRIの位置合わせ ~Boundary-Based Registration: BBR~


1. 目的
2. Boundary-Based Registration (BBR)とは
3. コマンド
4. 目的使用例
4.1. 目的前準備
4.2. 目的Boundary-Based Registration: BBB


1. 目的

2. Boundary-Based Registration (BBR)とは

Doug Greve氏によって開発されたEPI画像(拡散MRIや機能的MRI)用の位置合わせツールで、EPI画像と構造MRI画像(例:T1WI)との位置合わせで、灰白質・白質境界(Boundary)を頼りに位置合わせをする。以下のウェブサイトに、詳細な説明と分かりやすい図がある。

FLIRT_BBR

3. コマンド

Boundary-Based Registration (BBR) に基づいた構造MRIと拡散MRIの位置合わせをFSLで実装するには、epi_regを用いる。

epi_regのヘルプは、以下の通り。

クリックして展開
Usage: epi_reg [options] --epi=<EPI image> --t1=<wholehead T1 image> --t1brain=<brain extracted T1 image> --out=<output name>
 
Optional arguments
  --fmap=<image>         : fieldmap image (in rad/s)
  --fmapmag=<image>      : fieldmap magnitude image - wholehead extracted
  --fmapmagbrain=<image> : fieldmap magnitude image - brain extracted
  --gdc=<image>          : Gradient-distortion corection warpfield
  --wmseg=<image>        : white matter segmentation of T1 image
  --echospacing=<val>    : Effective EPI echo spacing (sometimes called dwell time) - in seconds
  --pedir=<dir>          : phase encoding direction, dir = x/y/z/-x/-y/-z
  --weight=<image>       : weighting image (in T1 space)
  --nofmapreg            : do not perform registration of fmap to T1 (use if fmap already registered) 
  --noclean              : do not clean up intermediate files
  -v                     : verbose output
  -h                     : display this help message
 
e.g.:  epi_reg --epi=example_func --t1=struct --t1brain=struct_brain --out=epi2struct --fmap=fmap_rads --fmapmag=fmap_mag --fmapmagbrain=fmap_mag_brain --echospacing=0.0005 --pedir=-y
 
Note that if parallel acceleration is used in the EPI acquisition then the *effective* echo spacing is the actual echo spacing between acquired lines in k-space divided by the acceleration factor.

基本的な使い方は、以下の通り。

epi_reg --epi=<b=0画像(spin-echo EPI)> --t1=<頭蓋除去前の3D-T1WI> --t1brain=<頭蓋除去後の3D-T1WI> --echospacing=<echo spacing時間(sec)> --pedir=<位相エンコード方向> --out=<出力画像>

4. 目的使用例

4.1 目的前準備

次のファイルを準備する。

.
├── CSF_GM_WM_seg.nii.gz:FASTで作成したCSF, GM, WMのセグメント
├── DWI.nii.gz:拡散MRI画像(b≠0)
├── DWI_b0.nii.gz:拡散MRIのb=0画像(spin-echo EPI)
├── T1_skull_stripped.nii.gz:頭蓋除去後の3D-T1WI
└── T1w.nii.gz:頭蓋除去前の3D-T1WI

頭蓋除去とFASTを用いたCSF, GM, WMのセグメントは、以下の記事を参考にするとよい。

また、拡散MRIからb値ごとに画像を抽出する方法は、以下の記事を参考にするとよい。

4.2. 目的Boundary-Based Registration: BBB

ここで、構造MRI(3D-T1WI)空間にあるCSF, GM, WMのセグメント(CSF_GM_WM_seg.nii.gz)を、拡散MRIの空間に移動させることを考える。

そのために、まずBBBを使って拡散MRIのb=0画像(spin-echo EPI)を、構造MRI(3D-T1WI)に位置合わせする。

この時、echo spacing時間( sec)--echospacingと位相エンコード方向--pedirを前もって調べておく必要がある。これらの情報は、dcm2niixを用いてDICOM形式からNIfTI形式に変換する際に出力されるJSONファイルに記載されている。dcm2niixの使い方は、以下の記事を参考にするとよい。

準備ができたら、以下のようにepi_regコマンドを実行して、拡散MRIから構造MRIへ位置合わせする変換行列を生成する。

epi_reg --epi=DWI_b0.nii.gz --t1=T1w.nii.gz --t1brain=T1_skull_stripped.nii.gz --echospacing=0.0380544 --pedir=-y --out=DWI2T1w

次に、先ほどの変換行列(拡散MRI空間→構造MRI空間)の逆変換行列(構造MRI空間→拡散MRI空間)を生成する。

convert_xfm -omat T1w2DWI.mat -inverse DWI2T1w.mat

逆変換行列(構造MRI空間→拡散MRI空間)を構造MRI空間にいるCSF, GM, WMのセグメント(CSF_GM_WM_seg.nii.gz)に適用することで、セグメントを拡散MRI空間に移動させる。このとき、-applyxfm-interp nearestneighbourとしていることに注意する。FLIRTの使い方の詳細はこちら

flirt -in CSF_GM_WM_seg.nii.gz -ref DWI_b0.nii.gz -out CSF_GM_WM_seg_DWIspace.nii.gz -init T1w2DWI.mat -applyxfm -interp nearestneighbour

処理前後のセグメントの様子をみてみる。

fsleyes DWI.nii.gz \
CSF_GM_WM_seg.nii.gz -cm random \
CSF_GM_WM_seg_DWIspace.nii.gz -cm random

軸位断で後頭葉付近を拡大してみる。

【MRtrix】拡散MRIからb値ごとに画像を抽出


1. 目的
2. コマンド
3.使用例
3.1.前準備
3.2.b=0のみを抽出
3.3.b≠0を抽出
3.4.b値ごとに抽出


1. 目的

  • 拡散MRIからb値ごとに画像を抽出

2. コマンド

拡散MRIからb値ごとに画像を抽出するには、MRtrixdwiextractを用いる。

dwiextractのヘルプは、次の通り。

クリックして展開
SYNOPSIS

     Extract diffusion-weighted volumes, b=0 volumes, or certain shells from a
     DWI dataset

USAGE

     dwiextract [ options ] input output

        input        the input DW image.

        output       the output image (diffusion-weighted volumes by default).


EXAMPLE USAGES

     Calculate the mean b=0 image from a 4D DWI series:
       $ dwiextract dwi.mif - -bzero | mrmath - mean mean_bzero.mif -axis 3
     The dwiextract command extracts all volumes for which the b-value is
     (approximately) zero; the resulting 4D image can then be provided to the
     mrmath command to calculate the mean intensity across volumes for each
     voxel.

OPTIONS

  -bzero
     Output b=0 volumes (instead of the diffusion weighted volumes, if
     -singleshell is not specified).

  -no_bzero
     Output only non b=0 volumes (default, if -singleshell is not specified).

  -singleshell
     Force a single-shell (single non b=0 shell) output. This will include b=0
     volumes, if present. Use with -bzero to enforce presence of b=0 volumes
     (error if not present) or with -no_bzero to exclude them.

DW gradient table import options

  -grad file
     Provide the diffusion-weighted gradient scheme used in the acquisition in
     a text file. This should be supplied as a 4xN text file with each line is
     in the format [ X Y Z b ], where [ X Y Z ] describe the direction of the
     applied gradient, and b gives the b-value in units of s/mm^2. If a
     diffusion gradient scheme is present in the input image header, the data
     provided with this option will be instead used.

  -fslgrad bvecs bvals
     Provide the diffusion-weighted gradient scheme used in the acquisition in
     FSL bvecs/bvals format files. If a diffusion gradient scheme is present in
     the input image header, the data provided with this option will be instead
     used.

DW shell selection options

  -shells b-values
     specify one or more b-values to use during processing, as a
     comma-separated list of the desired approximate b-values (b-values are
     clustered to allow for small deviations). Note that some commands are
     incompatible with multiple b-values, and will report an error if more than
     one b-value is provided. 
     WARNING: note that, even though the b=0 volumes are never referred to as
     shells in the literature, they still have to be explicitly included in the
     list of b-values as provided to the -shell option! Several algorithms
     which include the b=0 volumes in their computations may otherwise return
     an undesired result.

DW gradient table export options

  -export_grad_mrtrix path
     export the diffusion-weighted gradient table to file in MRtrix format

  -export_grad_fsl bvecs_path bvals_path
     export the diffusion-weighted gradient table to files in FSL (bvecs /
     bvals) format

Options for importing phase-encode tables

  -import_pe_table file
     import a phase-encoding table from file

  -import_pe_eddy config indices
     import phase-encoding information from an EDDY-style config / index file
     pair

Options for selecting volumes based on phase-encoding

  -pe desc
     select volumes with a particular phase encoding; this can be three
     comma-separated values (for i,j,k components of vector direction) or four
     (direction & total readout time)

Stride options

  -strides spec
     specify the strides of the output data in memory; either as a
     comma-separated list of (signed) integers, or as a template image from
     which the strides shall be extracted and used. The actual strides produced
     will depend on whether the output image format can support it.

Standard options

  -info
     display information messages.

  -quiet
     do not display information messages or progress status; alternatively,
     this can be achieved by setting the MRTRIX_QUIET environment variable to a
     non-empty string.

  -debug
     display debugging messages.

  -force
     force overwrite of output files (caution: using the same file as input and
     output might cause unexpected behaviour).

  -nthreads number
     use this number of threads in multi-threaded applications (set to 0 to
     disable multi-threading).

  -config key value  (multiple uses permitted)
     temporarily set the value of an MRtrix config file entry.

  -help
     display this information page and exit.

  -version
     display version information and exit.

基本的な使い方は、以下の通り。

dwiextract -bzero <入力画像> <出力画像>  # b=0のみを抽出
dwiextract -no_bzero <入力画像> <出力画像>  # b=0以外の拡散強調像を抽出
dwiextract -singleshell <入力画像> <出力画像>  # b=0以外の拡散強調像を抽出

3. 使用例

3.1. 前準備

まず、こちらの記事を参考に、拡散MRI(DWI.nii.gz)とそのMPG軸情報(bvecs, bvals)とヘッダー情報(headers.json)をまとめて、MIF形式(DWI.mif)に変換する。

mrconvert -fslgrad bvecs bvals -json_import headers.json DWI.nii.gz DWI.mif

ここで使用する拡散MRI(DWI.mif)は、b=0が1枚、b=1000が64枚、b=2000が64枚で構成されている(全部で129 volumes)。

mrinfo DWI.mif  |grep Dimensions
Dimensions:        130 x 130 x 82 x 129

3.2. b=0のみを抽出

オプション-bzeroを指定する。

dwiextract -bzero DWI.mif DWI_b0.mif

b=0の画像のみ抽出される。

mrinfo DWI_b0.mif  |grep Dimensions
Dimensions:        130 x 130 x 82 x 1

3.3. b≠0を抽出

オプション-no_bzeroを指定する。

dwiextract -no_bzero DWI.mif DWI_nonb0.mif

b≠0の画像のみ抽出される。

mrinfo DWI_nonb0.mif  |grep Dimensions
Dimensions:        130 x 130 x 82 x 128

3.4. b値ごとに抽出

オプション-singleshellを指定する。

例えば、b=1000のみを抽出する場合、以下のようになる。

dwiextract -shells 1000 DWI.mif DWI_b1000.mif

b=1000の画像のみ抽出される。

mrinfo DWI_b1000.mif  |grep Dimensions
Dimensions:        130 x 130 x 82 x 64

【MRtrix】MRtrixを用いた拡散MRIのマスク画像の作成


1. 目的
2. コマンド
3. 使用例
3.1. 前準備
3.2. 拡散MRIのマスク画像の作成


1. 目的

  • MRtrixを用いた拡散MRIのマスク画像の作成

2. コマンド

MRtrixを用いて拡散MRIのマスク画像の作成するには、dwi2maskを使用する。

dwi2maskのヘルプは、次の通り。

クリックして展開
SYNOPSIS

     Generates a whole brain mask from a DWI image

USAGE

     dwi2mask [ options ] input output

        input        the input DWI image containing volumes that are both
                     diffusion weighted and b=0

        output       the output whole-brain mask image


DESCRIPTION

     All diffusion weighted and b=0 volumes are used to obtain a mask that
     includes both brain tissue and CSF.

     In a second step peninsula-like extensions, where the peninsula itself is
     wider than the bridge connecting it to the mask, are removed. This may
     help removing artefacts and non-brain parts, e.g. eyes, from the mask.

OPTIONS

  -clean_scale value
     the maximum scale used to cut bridges. A certain maximum scale cuts
     bridges up to a width (in voxels) of 2x the provided scale. Setting this
     to 0 disables the mask cleaning step. (Default: 2)

DW gradient table import options

  -grad file
     Provide the diffusion-weighted gradient scheme used in the acquisition in
     a text file. This should be supplied as a 4xN text file with each line is
     in the format [ X Y Z b ], where [ X Y Z ] describe the direction of the
     applied gradient, and b gives the b-value in units of s/mm^2. If a
     diffusion gradient scheme is present in the input image header, the data
     provided with this option will be instead used.

  -fslgrad bvecs bvals
     Provide the diffusion-weighted gradient scheme used in the acquisition in
     FSL bvecs/bvals format files. If a diffusion gradient scheme is present in
     the input image header, the data provided with this option will be instead
     used.

基本的な使い方は、以下の通り。

dwi2mask <入力画像> <出力画像>

3. 使用例

3.1. 前準備

まず、こちらの記事を参考に、拡散MRI(DWI.nii.gz)とそのMPG軸情報(bvecs, bvals)とヘッダー情報(headers.json)をまとめて、MIF形式(DWI.mif)に変換する。

mrconvert -fslgrad bvecs bvals -json_import headers.json DWI.nii.gz DWI.mif

3.2. 拡散MRIのマスク画像の作成

以下のコマンドを実行する。

dwi2mask DWI.mif DWI_mask.mif

拡散MRIとマスク画像(緑)を重ね合わせてみる。

Ubuntu 20.04 / 18.04 環境で eddy_cuda10.2 (in FSL 6.0.5.x), PyTorch, Tensorflow 2 を使えるようにCUDA 10.2, 11.0, 11.5をセットアップする方法

私のメインマシンは Lin4Neuro 18.04 ですが、そろそろ Lin4Neuro 20.04 への移行を考えています。

今、実験機には NVIDIA GeForce RTX 2070 が備え付けられています。
これを使って、FSL 6.0.5 の eddy をGPUが使えるように設定し、なおかつ、Tensorflow, Pytorch といった Deep Learning のフレームワークも使えるようにしたいと思います。

FSL 6.0.5 にはデフォルトで CUDA 10.2 に対応した eddy_cuda10.2 が配布されています。なので、CUDA 10.2を入れることにします。

なお、これは Ubuntu 18.04 でも全く問題なくできることがわかりましたので、タイトルを変更しました。

続きを読む

【MRtrix】MRtrixを用いた解像度の変更 ~Upsampling~


1. 目的
2. コマンド
3. 使用例
3.1. ボクセルサイズを指定(オプション:-voxel)
3.2. スケールを指定(オプション:-scale))
3.3. ボクセルサイズを指定(オプション:-voxel))
3.4. 目的の解像度を持つ画像を指定(オプション:-template))


1. 目的

  • MRtrixを用いたアップサンプリング(Upsampling)

2. コマンド

MRtrixのmrgridを用いる。

mrgridのヘルプは、次の通り。

クリックして展開
SYNOPSIS

     Modify the grid of an image without interpolation (cropping or padding) or
     by regridding to an image grid with modified orientation, location and or
     resolution. The image content remains in place in real world coordinates.

USAGE

     mrgrid [ options ] input operation output

        input        input image to be regridded.

        operation    the operation to be performed, one of: regrid, crop, pad.

        output       the output image.


DESCRIPTION

     - regrid: This operation performs changes of the voxel grid that require
     interpolation of the image such as changing the resolution or location and
     orientation of the voxel grid. If the image is down-sampled, the
     appropriate smoothing is automatically applied using Gaussian smoothing
     unless nearest neighbour interpolation is selected or oversample is
     changed explicitly. The resolution can only be changed for spatial
     dimensions. 

     - crop: The image extent after cropping, can be specified either manually
     for each axis dimensions, or via a mask or reference image. The image can
     be cropped to the extent of a mask. This is useful for axially-acquired
     brain images, where the image size can be reduced by a factor of 2 by
     removing the empty space on either side of the brain. Note that cropping
     does not extend the image beyond the original FOV unless explicitly
     specified (via -crop_unbound or negative -axis extent).

     - pad: Analogously to cropping, padding increases the FOV of an image
     without image interpolation. Pad and crop can be performed simultaneously
     by specifying signed specifier argument values to the -axis option.

     This command encapsulates and extends the functionality of the superseded
     commands 'mrpad', 'mrcrop' and 'mrresize'. Note the difference in -axis
     convention used for 'mrcrop' and 'mrpad' (see -axis option description).

EXAMPLE USAGES

     Crop and pad the first axis:
       $ mrgrid in.mif crop -axis 0 10,-5 out.mif
     This removes 10 voxels on the lower and pads with 5 on the upper bound,
     which is equivalent to padding with the negated specifier (mrgrid in.mif
     pad -axis 0 -10,5 out.mif).

     Right-pad the image to the number of voxels of a reference image:
       $ mrgrid in.mif pad -as ref.mif -all_axes -axis 3 0,0 out.mif -fill nan
     This pads the image on the upper bound of all axes except for the volume
     dimension. The headers of in.mif and ref.mif are ignored and the output
     image uses NAN values to fill in voxels outside the original range of
     in.mif.

     Regrid and interpolate to match the voxel grid of a reference image:
       $ mrgrid in.mif regrid -template ref.mif -scale 1,1,0.5 out.mif -fill nan
     The -template instructs to regrid in.mif to match the voxel grid of
     ref.mif (voxel size, grid orientation and voxel centres). The -scale
     option overwrites the voxel scaling factor yielding voxel sizes in the
     third dimension that are twice as coarse as those of the template image.

Regridding options (involves image interpolation, applied to spatial axes only)

  -template image
     match the input image grid (voxel spacing, image size, header
     transformation) to that of a reference image. The image resolution
     relative to the template image can be changed with one of -size, -voxel,
     -scale.

  -size dims
     define the size (number of voxels) in each spatial dimension for the
     output image. This should be specified as a comma-separated list.

  -voxel size
     define the new voxel size for the output image. This can be specified
     either as a single value to be used for all spatial dimensions, or as a
     comma-separated list of the size for each voxel dimension.

  -scale factor
     scale the image resolution by the supplied factor. This can be specified
     either as a single value to be used for all dimensions, or as a
     comma-separated list of scale factors for each dimension.

  -interp method
     set the interpolation method to use when reslicing (choices: nearest,
     linear, cubic, sinc. Default: cubic).

  -oversample factor
     set the amount of over-sampling (in the target space) to perform when
     regridding. This is particularly relevant when downsamping a
     high-resolution image to a low-resolution image, to avoid aliasing
     artefacts. This can consist of a single integer, or a comma-separated list
     of 3 integers if different oversampling factors are desired along the
     different axes. Default is determined from ratio of voxel dimensions
     (disabled for nearest-neighbour interpolation).

Pad and crop options (no image interpolation is performed, header transformation is adjusted)

  -as reference image
     pad or crop the input image on the upper bound to match the specified
     reference image grid. This operation ignores differences in image
     transformation between input and reference image.

  -uniform number
     pad or crop the input image by a uniform number of voxels on all sides

  -mask image
     crop the input image according to the spatial extent of a mask image. The
     mask must share a common voxel grid with the input image but differences
     in image transformations are ignored. Note that even though only 3
     dimensions are cropped when using a mask, the bounds are computed by
     checking the extent for all dimensions. Note that by default a gap of 1
     voxel is left at all edges of the image to allow valid trilinear
     interpolation. This gap can be modified with the -uniform option but by
     default it does not extend beyond the FOV unless -crop_unbound is used.

  -crop_unbound
     Allow padding beyond the original FOV when cropping.

  -axis index spec  (multiple uses permitted)
     pad or crop the input image along the provided axis (defined by index).
     The specifier argument defines the number of voxels added or removed on
     the lower or upper end of the axis (-axis index delta_lower,delta_upper)
     or acts as a voxel selection range (-axis index start:stop). In both
     modes, values are relative to the input image (overriding all other
     extent-specifying options). Negative delta specifier values trigger the
     inverse operation (pad instead of crop and vice versa) and negative range
     specifier trigger padding. Note that the deprecated commands 'mrcrop' and
     'mrpad' used range-based and delta-based -axis indices, respectively.

  -all_axes
     Crop or pad all, not just spatial axes.

General options

  -fill number
     Use number as the out of bounds value. nan, inf and -inf are valid
     arguments. (Default: 0.0)

Stride options

  -strides spec
     specify the strides of the output data in memory; either as a
     comma-separated list of (signed) integers, or as a template image from
     which the strides shall be extracted and used. The actual strides produced
     will depend on whether the output image format can support it.

Data type options

  -datatype spec
     specify output image data type. Valid choices are: float32, float32le,
     float32be, float64, float64le, float64be, int64, uint64, int64le,
     uint64le, int64be, uint64be, int32, uint32, int32le, uint32le, int32be,
     uint32be, int16, uint16, int16le, uint16le, int16be, uint16be, cfloat32,
     cfloat32le, cfloat32be, cfloat64, cfloat64le, cfloat64be, int8, uint8,
     bit.

Standard options

  -info
     display information messages.

  -quiet
     do not display information messages or progress status; alternatively,
     this can be achieved by setting the MRTRIX_QUIET environment variable to a
     non-empty string.

  -debug
     display debugging messages.

  -force
     force overwrite of output files (caution: using the same file as input and
     output might cause unexpected behaviour).

  -nthreads number
     use this number of threads in multi-threaded applications (set to 0 to
     disable multi-threading).

  -config key value  (multiple uses permitted)
     temporarily set the value of an MRtrix config file entry.

  -help
     display this information page and exit.

  -version
     display version information and exit.

解像度の変更する場合の基本的な使い方は、以下の通り。

mrgrid <入力画像> regrid -voxel <値> <出力画像>  # ボクセルサイズを指定
mrgrid <入力画像> regrid -scale <値> <出力画像>  # スケールを指定
mrgrid <入力画像> regrid -template <目的の解像度を持つ画像> <出力画像>  # 目的の解像度を持つ画像を指定

3. 使用例

3D-T1WI(T1w.nii.gz)の解像度を変更する。

3D-T1WI(T1w.nii.gz)の解像度を確認してみる。

mrinfo T1w.nii.gz
************************************************
Image name:          "T1w.nii.gz"
************************************************
  Dimensions:        192 x 256 x 256
  Voxel size:        0.9 x 0.9375 x 0.9375
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         signed 16 bit integer (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -82.89
                         -0.01788      0.9946      0.1023      -113.6
                         0.001492     -0.1023      0.9948      -114.6
  comments:          6.0.3:b862cdd5

3.1. ボクセルサイズを指定(オプション:-voxel

-voxelオプションを用いて、以下のコマンドを実行。

ボクセルサイズを1mm isotropicにする。

mrgrid T1w.nii.gz regrid -voxel 1 T1w_1mm_iso.nii.gz

解像度を確認してみる。ボクセルサイズが1 x 1 x 1(1mm iso)になっている

mrinfo T1w_1mm_iso.nii.gz
************************************************
Image name:          "T1w_1mm_iso.nii.gz"
************************************************
  Dimensions:        173 x 240 x 240
  Voxel size:        1 x 1 x 1
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -82.94
                         -0.01788      0.9946      0.1023      -113.6
                         0.001492     -0.1023      0.9948      -114.5
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225

3.2. スケールを指定(オプション:-scale

-scaleオプションを用いて、以下のコマンドを実行。

スケールを2にして、解像度を2倍にする。

mrgrid T1w.nii.gz regrid -scale 2 T1w_scale2.nii.gz

解像度を確認してみる。解像度が173 x 240 x 240からになっている。

mrinfo T1w_scale2.nii.gz
************************************************
Image name:          "T1w_scale2.nii.gz"
************************************************
  Dimensions:        384 x 512 x 512
  Voxel size:        0.45 x 0.46875 x 0.46875
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -83.12
                         -0.01788      0.9946      0.1023      -113.9
                         0.001492     -0.1023      0.9948      -114.8
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225

3.3. ボクセルサイズを指定(オプション:-voxel

-voxelオプションを用いて、以下のコマンドを実行。

ボクセルサイズを1mm isotropicにする。

mrgrid T1w.nii.gz regrid -voxel 1 T1w_1mm_iso.nii.gz

解像度を確認してみる。ボクセルサイズが1 x 1 x 1(1mm iso)になっている。

mrinfo T1w_1mm_iso.nii.gz
************************************************
Image name:          "T1w_1mm_iso.nii.gz"
************************************************
  Dimensions:        173 x 240 x 240
  Voxel size:        1 x 1 x 1
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -82.94
                         -0.01788      0.9946      0.1023      -113.6
                         0.001492     -0.1023      0.9948      -114.5
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225

3.4. 目的の解像度を持つ画像を指定(オプション:-template

標準脳(MNI152)の3D-T1WI(MNI152_T1_2mm.nii.gz)と同じ解像度にする。標準脳(MNI152_T1_2mm.nii.gz)の解像度は以下。

mrinfo MNI152_T1_2mm.nii.gz
************************************************
Image name:          "MNI152_T1_2mm.nii.gz"
************************************************
  Dimensions:        91 x 109 x 91
  Voxel size:        2 x 2 x 2
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         signed 16 bit integer (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                    1           0           0         -90
                               -0           1           0        -126
                               -0           0           1         -72
  comments:          FSL5.0

個人脳(T1w.nii.gz)を標準脳(MNI152_T1_2mm.nii.gz)の解像度に合わせるには、-templateオプションを用いて、以下のコマンドを実行。

mrgrid T1w.nii.gz regrid -template MNI152_T1_2mm.nii.gz T1w_MNIreso.nii.gz

解像度を確認してみる。解像度が標準脳(MNI152_T1_2mm.nii.gz)と同じになっている。

mrinfo T1w_MNIreso.nii.gz
************************************************
Image name:          "T1w_MNIreso.nii.gz"
************************************************
  Dimensions:        91 x 109 x 91
  Voxel size:        2 x 2 x 2
  Data strides:      [ -1 2 3 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:                    1           0           0         -90
                               -0           1           0        -126
                               -0           0           1         -72
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225

【FSL/MRtrix】4D画像から3D画像を抽出


1. 目的
2. FSLを用いる場合
2.1. コマンド
2.2. 使用例
3. MRtrixを用いる場合
3.1. コマンド
3.2. 使用例


1. 目的

  • 4D画像から3D画像を抽出

2. FSLを用いる場合

2.1. コマンド

FSLfslroiコマンドを用いる。

fslroiのヘルプは、次の通り。

クリックして展開
Usage: fslroi <input> <output> <xmin> <xsize> <ymin> <ysize> <zmin> <zsize>
       fslroi <input> <output> <tmin> <tsize>

       fslroi <input> <output> <xmin> <xsize> <ymin> <ysize> <zmin> <zsize> <tmin> <tsize>
Note: indexing (in both time and space) starts with 0 not 1! Inputting -1 for a size will set it to the full image extent for that dimension.

4D画像から3D画像を抽出する際の、基本的な使い方は以下の通り。

fslroi <入力画像> <出力画像> <Volume Index> <Volume Indexから残したいVolume数>

2.2. 使用例

例えば、5ttgen等で作成した以下のような5つの組織画像(5tt.nii.gz)が4D画像となっている場合。

Pathological tissue(Volume 4th)を取り除くには、次のようにコマンドを実行する。FSLではVolumeのIndexを0から数える。つまり、1番目のVolumeのIndexは0となる。以下のコードを翻訳すると、「Volume Index0番から数えて4 Volumesまでを残す」ということになる。

fslroi 5tt.nii.gz 4tt.nii.gz 0 4

fslhdコマンドを用いて、ボリューム数を確認すると、処理前で5 Volumesだったのが処理後に4 Volumesになっていることが分かる。使い方の詳細は、こちらの記事を参考に。

fslhd 5tt.nii.gz |grep ^dim4
fslhd 4tt.nii.gz |grep ^dim4
dim4        5  # 5tt.nii.gz
dim4        4  # 4tt.nii.gz

3. MRtrixを用いる場合

3.1. コマンド

MRtrixmrconvertコマンドを用いる。

mrconvertのヘルプは、次の通り。

クリックして展開
SYNOPSIS

     Perform conversion between different file types and optionally extract a
     subset of the input image

USAGE

     mrconvert [ options ] input output

        input        the input image.

        output       the output image.


DESCRIPTION

     If used correctly, this program can be a very useful workhorse. In
     addition to converting images between different formats, it can be used to
     extract specific studies from a data set, extract a specific region of
     interest, or flip the images. Some of the possible operations are
     described in more detail below.

     Note that for both the -coord and -axes options, indexing starts from 0
     rather than 1. E.g. -coord 3 <#> selects volumes (the fourth dimension)
     from the series; -axes 0,1,2 includes only the three spatial axes in the
     output image.

     Additionally, for the second input to the -coord option and the -axes
     option, you can use any valid number sequence in the selection, as well as
     the 'end' keyword (see the main documentation for details); this can be
     particularly useful to select multiple coordinates.

     The -vox option is used to change the size of the voxels in the output
     image as reported in the image header; note however that this does not
     re-sample the image based on a new voxel size (that is done using the
     mrresize command).

     By default, the intensity scaling parameters in the input image header are
     passed through to the output image header when writing to an integer
     image, and reset to 0,1 (i.e. no scaling) for floating-point and binary
     images. Note that the -scaling option will therefore have no effect for
     floating-point or binary output images.

     The -axes option specifies which axes from the input image will be used to
     form the output image. This allows the permutation, omission, or addition
     of axes into the output image. The axes should be supplied as a
     comma-separated list of axis indices. If an axis from the input image is
     to be omitted from the output image, it must either already have a size of
     1, or a single coordinate along that axis must be selected by the user by
     using the -coord option. Examples are provided further below.

     The -bvalue_scaling option controls an aspect of the import of diffusion
     gradient tables. When the input diffusion-weighting direction vectors have
     norms that differ substantially from unity, the b-values will be scaled by
     the square of their corresponding vector norm (this is how multi-shell
     acquisitions are frequently achieved on scanner platforms). However in
     some rare instances, the b-values may be correct, despite the vectors not
     being of unit norm (or conversely, the b-values may need to be rescaled
     even though the vectors are close to unit norm). This option allows the
     user to control this operation and override MRrtix3's automatic detection.

EXAMPLE USAGES

     Extract the first volume from a 4D image, and make the output a 3D image:
       $ mrconvert in.mif -coord 3 0 -axes 0,1,2 out.mif
     The -coord 3 0 option extracts, from axis number 3 (which is the fourth
     axis since counting begins from 0; this is the axis that steps across
     image volumes), only coordinate number 0 (i.e. the first volume). The
     -axes 0,1,2 ensures that only the first three axes (i.e. the spatial axes)
     are retained; if this option were not used in this example, then image
     out.mif would be a 4D image, but it would only consist of a single volume,
     and mrinfo would report its size along the fourth axis as 1.

     Extract slice number 24 along the AP direction:
       $ mrconvert volume.mif slice.mif -coord 1 24
     MRtrix3 uses a RAS (Right-Anterior-Superior) axis convention, and
     internally reorients images upon loading in order to conform to this as
     far as possible. So for non-exotic data, axis 1 should correspond
     (approximately) to the anterior-posterior direction.

     Extract only every other volume from a 4D image:
       $ mrconvert all.mif every_other.mif -coord 3 1:2:end
     This example demonstrates two features: Use of the colon syntax to
     conveniently specify a number sequence (in the format 'start:step:stop');
     and use of the 'end' keyword to generate this sequence up to the size of
     the input image along that axis (i.e. the number of volumes).

     Alter the image header to report a new isotropic voxel size:
       $ mrconvert in.mif isotropic.mif -vox 1.25
     By providing a single value to the -vox option only, the specified value
     is used to set the voxel size in mm for all three spatial axes in the
     output image.

     Alter the image header to report a new anisotropic voxel size:
       $ mrconvert in.mif anisotropic.mif -vox 1,,3.5
     This example will change the reported voxel size along the first and third
     axes (ideally left-right and inferior-superior) to 1.0mm and 3.5mm
     respectively, and leave the voxel size along the second axis (ideally
     anterior-posterior) unchanged.

     Turn a single-volume 4D image into a 3D image:
       $ mrconvert 4D.mif 3D.mif -axes 0,1,2
     Sometimes in the process of extracting or calculating a single 3D volume
     from a 4D image series, the size of the image reported by mrinfo will be
     "X x Y x Z x 1", indicating that the resulting image is in fact also 4D,
     it just happens to contain only one volume. This example demonstrates how
     to convert this into a genuine 3D image (i.e. mrinfo will report the size
     as "X x Y x Z".

     Insert an axis of size 1 into the image:
       $ mrconvert XYZD.mif XYZ1D.mif -axes 0,1,2,-1,3
     This example uses the value -1 as a flag to indicate to mrconvert where a
     new axis of unity size is to be inserted. In this particular example, the
     input image has four axes: the spatial axes X, Y and Z, and some form of
     data D is stored across the fourth axis (i.e. volumes). Due to insertion
     of a new axis, the output image is 5D: the three spatial axes (XYZ), a
     single volume (the size of the output image along the fourth axis will be
     1), and data D will be stored as volume groups along the fifth axis of the
     image.

     Manually reset the data scaling parameters stored within the image header
     to defaults:
       $ mrconvert with_scaling.mif without_scaling.mif -scaling 0.0,1.0
     This command-line option alters the parameters stored within the image
     header that provide a linear mapping from raw intensity values stored in
     the image data to some other scale. Where the raw data stored in a
     particular voxel is I, the value within that voxel is interpreted as:
     value = offset + (scale x I).  To adjust this scaling, the relevant
     parameters must be provided as a comma-separated 2-vector of
     floating-point values, in the format "offset,scale" (no quotation marks).
     This particular example sets the offset to zero and the scale to one,
     which equates to no rescaling of the raw intensity data.

Options for manipulating fundamental image properties

  -coord axis selection  (multiple uses permitted)
     retain data from the input image only at the coordinates specified in the
     selection along the specified axis. The selection argument expects a
     number sequence, which can also include the 'end' keyword.

  -vox sizes
     change the voxel dimensions reported in the output image header

  -axes axes
     specify the axes from the input image that will be used to form the output
     image

  -scaling values
     specify the data scaling parameters used to rescale the intensity values

Options for handling JSON (JavaScript Object Notation) files

  -json_import file
     import data from a JSON file into header key-value pairs

  -json_export file
     export data from an image header key-value pairs into a JSON file

Options to modify generic header entries

  -clear_property key  (multiple uses permitted)
     remove the specified key from the image header altogether.

  -set_property key value  (multiple uses permitted)
     set the value of the specified key in the image header.

  -append_property key value  (multiple uses permitted)
     append the given value to the specified key in the image header (this adds
     the value specified as a new line in the header value).

  -copy_properties source
     clear all generic properties and replace with the properties from the
     image / file specified.

Stride options

  -strides spec
     specify the strides of the output data in memory; either as a
     comma-separated list of (signed) integers, or as a template image from
     which the strides shall be extracted and used. The actual strides produced
     will depend on whether the output image format can support it.

Data type options

  -datatype spec
     specify output image data type. Valid choices are: float32, float32le,
     float32be, float64, float64le, float64be, int64, uint64, int64le,
     uint64le, int64be, uint64be, int32, uint32, int32le, uint32le, int32be,
     uint32be, int16, uint16, int16le, uint16le, int16be, uint16be, cfloat32,
     cfloat32le, cfloat32be, cfloat64, cfloat64le, cfloat64be, int8, uint8,
     bit.

DW gradient table import options

  -grad file
     Provide the diffusion-weighted gradient scheme used in the acquisition in
     a text file. This should be supplied as a 4xN text file with each line is
     in the format [ X Y Z b ], where [ X Y Z ] describe the direction of the
     applied gradient, and b gives the b-value in units of s/mm^2. If a
     diffusion gradient scheme is present in the input image header, the data
     provided with this option will be instead used.

  -fslgrad bvecs bvals
     Provide the diffusion-weighted gradient scheme used in the acquisition in
     FSL bvecs/bvals format files. If a diffusion gradient scheme is present in
     the input image header, the data provided with this option will be instead
     used.

  -bvalue_scaling mode
     enable or disable scaling of diffusion b-values by the square of the
     corresponding DW gradient norm (see Desciption). Valid choices are yes/no,
     true/false, 0/1 (default: automatic).

DW gradient table export options

  -export_grad_mrtrix path
     export the diffusion-weighted gradient table to file in MRtrix format

  -export_grad_fsl bvecs_path bvals_path
     export the diffusion-weighted gradient table to files in FSL (bvecs /
     bvals) format

Options for importing phase-encode tables

  -import_pe_table file
     import a phase-encoding table from file

  -import_pe_eddy config indices
     import phase-encoding information from an EDDY-style config / index file
     pair

Options for exporting phase-encode tables

  -export_pe_table file
     export phase-encoding table to file

  -export_pe_eddy config indices
     export phase-encoding information to an EDDY-style config / index file
     pair

Standard options

  -info
     display information messages.

  -quiet
     do not display information messages or progress status; alternatively,
     this can be achieved by setting the MRTRIX_QUIET environment variable to a
     non-empty string.

  -debug
     display debugging messages.

  -force
     force overwrite of output files (caution: using the same file as input and
     output might cause unexpected behaviour).

  -nthreads number
     use this number of threads in multi-threaded applications (set to 0 to
     disable multi-threading).

  -config key value  (multiple uses permitted)
     temporarily set the value of an MRtrix config file entry.

  -help
     display this information page and exit.

  -version
     display version information and exit.

4D画像から3D画像を抽出する際の、基本的な使い方は以下の通り。

mrconvert <入力画像> <出力画像> -coord <軸番号> <残したいボリューム数>

3.2. 使用例

例えば、5ttgen等で作成した以下のような5つの組織画像(5tt.nii.gz)が4D画像となっている場合。

Pathological tissue(Volume 4th)を取り除くには、次のようにコマンドを実行する。MRtrixでもFSLと同様に、VolumeのIndexを0から数える。つまり、1番目のVolumeのIndexは0となる。また軸番号は、x, y, z, tの順番に0, 1, 2, 3であり、Volume数を操作するには、t軸(-coord 3)を操作することになる。以下のコードを翻訳すると、「Volume Index0番からVolume Index3番までを残す」ということになる。

mrconvert 5tt.nii.gz 4tt.nii.gz -coord 3 0:3

mrinfoコマンドを用いて、ボリューム数を確認すると、処理前で5 Volumesだったのが処理後に4 Volumesになっていることが分かる。使い方の詳細は、こちらの記事を参考に。

mrinfo 5tt.nii.gz 4tt.nii.gz
************************************************
Image name:          "5tt.nii.gz"
************************************************
  Dimensions:        168 x 185 x 169 x 5
  Voxel size:        0.9 x 0.9375 x 0.9375 x ?
  Data strides:      [ 1 2 3 4 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -70.81
                         -0.01788      0.9946      0.1023       -88.1
                         0.001492     -0.1023      0.9948      -56.89
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225
************************************************
Image name:          "4tt.nii.gz"
************************************************
  Dimensions:        168 x 185 x 169 x 4
  Voxel size:        0.9 x 0.9375 x 0.9375 x ?
  Data strides:      [ 1 2 3 4 ]
  Format:            NIfTI-1.1 (GZip compressed)
  Data type:         32 bit float (little endian)
  Intensity scaling: offset = 0, multiplier = 1
  Transform:               0.9998     0.01794   0.0003439      -70.81
                         -0.01788      0.9946      0.1023       -88.1
                         0.001492     -0.1023      0.9948      -56.89
  comments:          6.0.3:b862cdd5
  mrtrix_version:    3.0.0-40-g3e1ed225