【FSL】ROI(VOI)を用いた画像計測



1. 目的
2. 計測コマンド
3. 複数人被験者がいる場合
3.1. フォルダ構造
3.2. ソースコード


1. 目的

  • 関心領域ROI(VOI)内の画像の定量値を計測

2. 計測コマンド

FSLのfslstatsを使う。
各オプションは、以下の通り。

  • -K: 関心領域ROI
  • -M: ROI内の平均値の算出
fslstats -K <ROI> <input> -M

3. 複数人被験者がいる場合

複数人の被験者がいる場合、以下のようなフォルダ構造にならって各被験者ごとにファイルを用意し、後述するソースコードを実行することで、すべての被験者の計測結果をresult_mapstats.csvにまとめてくれる。

3.1. フォルダ構造

data
├── sub001
│   ├── basalganglia_roi_in_qsm.nii.gz
│   └── qsm_map.nii.gz
├── sub002
│   ├── basalganglia_roi_in_qsm.nii.gz
│   └── qsm_map.nii.gz
└── sub003
    ├── basalganglia_roi_in_qsm.nii.gz
    └── qsm_map.nii.gz

3.2. ソースコード

すべての被験者のデータが入っている「data」フォルダと同じ階層で以下のコマンドを実行。

# define name
foldername="data"
map="std_qsm_map"
roi="basalganglia_roi_in_qsm"

# make ROI name file
for roiname in \
"Region" \
"Left-Thalamus-Proper" \
"Left-Caudate" \
"Left-Putamen" \
"Left-Pallidum" \
"Brain-Stem/4th/Ventricle" \
"Left-Hippocampus" \
"Left-Amygdala" \
"Left-Accumbens-area" \
"Right-Thalamus-Proper" \
"Right-Caudate" \
"Right-Putamen" \
"Right-Pallidum" \
"Right-Hippocampus" \
"Right-Amygdala" \
"Right-Accumbens-area"
do
echo $roiname >> tmp1
done

# segment basal ganglia and calc volume
cd $foldername
for k in *;do
echo "processing $k"
echo $k > ../tmp2  # make ID name file
# measure mean value from QSM map using basal ganglia ROI
fslstats -K $k/$roi $k/$map -M \
|grep -v missing >> ../tmp2
## save ID and volume ,and accumulate them in result_vol.csv 
cat ../tmp2 > ../tmp3
paste ../tmp1 ../tmp3 > ../result_mapstats.csv
cat ../result_mapstats.csv > ../tmp1
# cd ..
done
cd ..

# remove temporaly files
rm tmp*

【FSL】TBSSで出た有意領域をマスクとして画像定量値を計測



1. 目的
2. ソースコード


1. 目的

2. ソースコード

TBSSの結果が格納されているフォルダ(「stats」フォルダ)で以下のコードを実行。
この例では、健常群と患者群の2群比較をした結果を用いている。

c1: 患者群 > 健常群
c2: 患者群 < 健常群
の片側検定の結果をまとめている。

手順は次の通りです。

  1. corrected p mapを0.95でしきい処理することで有意差のあった領域のマスクを作成
  2. 1.のマスクを用いて各被験者の画像定量値を計測
  3. skeletonisedされた画像(all_[map]skeletonised.nii.gz)とそうでない画像(all[map].nii.gz)のそれぞれを別々に計測
  4. 計測した結果をまとめる(c?_allmap.txt)
.nii.gz)
    # C1 (ASD>HC)
    fslstats -t all_${map}.nii.gz -k tbss${map}_tfce_corrp_tstat1_sigmask.nii.gz -M >> stats_sigmap/original_map/c1_${i}_${map}.txt
    # C2 (ASD<HC)
    fslstats -t all_${map}.nii.gz -k tbss${map}_tfce_corrp_tstat2_sigmask.nii.gz -M >> stats_sigmap/original_map/c2_${i}_${map}.txt

    # calc skeletonised map (all_[map]_skeletonised.nii.gz)
    # C1 (ASD>HC)
    fslstats -t all_${map}_skeletonised.nii.gz -k tbss${map}_tfce_corrp_tstat1_sigmask.nii.gz -M >> stats_sigmap/skeletonised_map/c1_${i}_${map}.txt
    # C2 (ASD<HC)
    fslstats -t all_${map}_skeletonised.nii.gz -k tbss${map}_tfce_corrp_tstat2_sigmask.nii.gz -M >> stats_sigmap/skeletonised_map/c2_${i}_${map}.txt
done

# summary result
paste stats_sigmap/original_map/c1_* > stats_sigmap/original_map/c1_allmap.txt
paste stats_sigmap/original_map/c2_* > stats_sigmap/original_map/c2_allmap.txt
paste stats_sigmap/skeletonised_map/c1_* > stats_sigmap/skeletonised_map/c1_allmap.txt
paste stats_sigmap/skeletonised_map/c2_* > stats_sigmap/skeletonised_map/c2_allmap.txt

【FSL】XTRACTを用いたトラクトグラフィー



1. 目的
2. XTRACT
3. 前処理
4. 必要なファイル
5. 実行
6. 結果
7. XTRACTのトラクトを用いて、定量値のサンプル


1. 目的

  • XTRACTを用いたトラクトグラフィー(Tractography)
  • 主要白質路の抽出

2. XTRACT

XTRACTの論文はこちら

https://www.sciencedirect.com/science/article/pii/S1053811920304092

XTRACTのドキュメントはこちら

https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/XTRACT:embed:cite

XTRACTによって以下の白質路がトラクトグラフィーによって抽出される。

Tract Abbreviation Left/Right? Seeding Strategy
Association Fibres Arcuate Fasciculus AF Yes Reverse-seeding
Frontal Aslant Tract FA Yes Single-ROI
Inferior Longitudinal Fasciculus ILF Yes Reverse-seeding
Inferior Fronto-Occipital fasciculus IFO Yes Reverse-seeding
Middle Longitudinal Fasciculus MdLF Yes Reverse-seeding
Superior Longitudinal Fasciculus I SLF1 Yes Single-ROI
Superior Longitudinal Fasciculus II SLF2 Yes Single-ROI
Superior Longitudinal Fasciculus III SLF3 Yes Single-ROI
Uncinate Fasciculus UF Yes Single-ROI
Vertical Occipital Fasciculus VOF Yes Reverse-seeding
Commissural Fibres Anterior Commissure AC No Reverse-seeding
Forceps Major FMA No Reverse-seeding
Forceps Minor FMI No Reverse-seeding
Middle Cerebellar Peduncle MCP No Reverse-seeding
Limbic Fibres Cingulum subsection: Dorsal CBD Yes Single-ROI
Cingulum subsection: Peri-genual CBP Yes Single-ROI
Cingulum subsection: Temporal CBT Yes Single-ROI
Fornix FX Yes Single-ROI
Projection Fibres Acoustic Radiation AR Yes Reverse-seeding
Anterior Thalamic Radiation ATR Yes Single-ROI
Corticospinal Tract CST Yes Single-ROI
Optic Radiation OR Yes Reverse-seeding
  Superior Thalamic Radiation STR Yes Single-ROI

(Shaun W, et al. NeuroImage. 2020より)

XTRACTコマンドのヘルプは次の通り。

 __  _______ ____      _    ____ _____ 
 \ \/ /_   _|  _ \    / \  / ___|_   _|
  \  /  | | | |_) |  / _ \| |     | |  
  /  \  | | |  _ <  / ___ \ |___  | |  
 /_/\_\ |_| |_| \_\/_/   \_\____| |_|  
 

Usage: 
    xtract -bpx <bedpostX_dir> -out <outputDir> -str <structuresFile> -p <protocolsFolder> [options]
    xtract -bpx <bedpostX_dir> -out <outputDir> -species HUMAN [options]
    xtract -bpx <bedpostX_dir> -out <outputDir> -species MACAQUE [options]

    Compulsory arguments:

       -bpx <folder>                     Path to bedpostx folder
       -out <folder>                     Path to output folder
       
       And EITHER:
       -str <file>                       Structures file (format: <tractName> [samples=1], 1 means 1000, '#' to skip lines)
       -p   <folder>                     Protocols folder (all masks in same standard space)

       Or:
       -species <SPECIES>                One of HUMAN or MACAQUE

    Optional arguments:

       -stdwarp <std2diff> <diff2std>    Standard2diff and Diff2standard transforms (Default=bedpostx_dir/xfms/{standard2diff,diff2standard}) 
       -gpu                              Use GPU version 
       -native                           Run tractography in native (diffusion) space
       -res <mm>                         Output resolution (Default=same as in protocol folders unless '-native' used)
       -ptx_options <options.txt>	 Pass extra probtrackx2 options as a text file to override defaults, e.g. --steplength=0.2 --distthresh=10)

3. 前処理

XTRACTに必要な前処理は次の通り。FDT processing pipelineで処理することが可能。

4. 必要なファイル

必要なファイルは次の通り。

.
├── DTI.bedpostX  # BEDPOSTXの出力フォルダ
└── T1_brain.nii.gz  # BET後のBrain 3D-T1WI

「DTI.bedpostX/xfms/」にDiffusion空間からMNI空間への変換行列があることを確認する。

DTI.bedpostX/xfms/
├── diff2standard.mat
├── diff2standard_warp.nii.gz
├── diff2str.mat
├── eye.mat
├── standard2diff.mat
├── standard2diff_warp.nii.gz
├── standard2str.mat
├── standard2str_warp.nii.gz
├── str2diff.mat
├── str2standard.mat
└── str2standard_warp.nii.gz

5. 実行

XTRACTを実行するには、次のコマンドを実行。

  • GPUを使用する場合:-gpuオプション追加
  • 個人のDiffusion空間で出力する場合:-nativeオプション追加
BPX_DIR='DTI.bedpostX'

xtract -bpx $BPX_DIR -out XTRACT_output -str T1_brain.nii.gz -species HUMAN -stdwarp $BPX_DIR/xfms/standard2diff_warp $BPX_DIR/xfms/diff2standard_warp

6. 結果

結果の可視化には、xtract_viewerを用いるとよい。xtract_viewerのヘルプは次の通り。

 __  _______ ____      _    ____ _____         _                        
 \ \/ /_   _|  _ \    / \  / ___|_   _| __   _(_) _____      _____ _ __ 
  \  /  | | | |_) |  / _ \| |     | |   \ \ / / |/ _ \ \ /\ / / _ \ '__|
  /  \  | | |  _ <  / ___ \ |___  | |    \ V /| |  __/\ V  V /  __/ |   
 /_/\_\ |_| |_| \_\/_/   \_\____| |_|     \_/ |_|\___| \_/\_/ \___|_|                                                                           
                                                             

Usage:
    xtract_viewer -dir <xtractDir> [options]

    Compulsory arguments:

       -dir FOLDER                       Path to XTRACT output folder

    Optional arguments:

       -str STRUCTURE,STRUCTURE,...      Structures (comma separated (default = display all that is found in input folder)

       -thr NUMBER NUMBER                The lower and upper thresholds applied to the tracts for viewing
                                         Default = 0.001 0.1

       -brain                            The brain image to use for the background overlay - must be in the same space as tracts.
                                         Default is the FSL_HCP065_FA map

今回の例では、次のようなコマンドを実行する。引数として、XTRACTの出力フォルダを指定する。個人空間でXTRACTをした場合:-brainで個人脳を指定。

xtract_viewer -dir XTRACT_output

7. XTRACTのトラクトを用いて、定量値のサンプル

XTRACTのトラクトグラフィーで抽出した、白質ROIを用いて定量値を算出することも可能。コマンドは、xtract_statsを用いる。

xtract_statsのヘルプは次の通り。

__  _______ ____      _    ____ _____    _        _
\ \/ /_   _|  _ \    / \  / ___|_   _|__| |_ __ _| |_ ___
 \  /  | | | |_) |  / _ \| |     | |/ __| __/ _  | __/ __|
 /  \  | | |  _ <  / ___ \ |___  | |\__ \ || (_| | |_\__ \
/_/\_\ |_| |_| \_\/_/   \_\____| |_||___/\__\__ _|\__|___/


Usage:
    xtract_stats -d <dir_basename> -xtract <XTRACT_dir> -w <xtract2diff> [options]

    Compulsory arguments:

       -d <folder_basename>                   Path to microstructure folder and basename of data (e.g. /home/DTI/dti_)
       -xtract <folder>                       Path to XTRACT output folder
       -w <xtract2diff>                       EITHER XTRACT results to diffusion space transform OR 'native' if tracts are already in diffusion space

    Optional arguments:
       -r <reference>                         If not 'native', provide reference image in diffusion space (e.g. /home/DTI/dti_FA)
       -out <path>                            Output filepath (Default <XTRACT_dir>/stats.csv)
       -str <file>                            Structures file (as in XTRACT) (Default is all tracts under <XTRACT_dir>)
       -thr <float>                           Threshold applied to tract probability map (default = 0.001 = 0.1%)

       -meas <list>                           Comma separated list of features to extract (Default = vol,prob,length,FA,MD - assumes DTI folder has been provided)
                                              vol = tract volume, prob = tract probability, length = tract length
                                              Additional metrics must follow file naming conventions. e.g. for dti_L1 use 'L1'

       -keepfiles                             Keep temporary files

計測するための、Diffusion定量値を計算する。DTIの定量値を計算するには、dtifitを用いる。

DTIFITの詳細はこちらをご覧ください。

dtifitのために次のファイルを用意する。

DTI_metrics/
├── bvals  # DWIのGradient Table
├── bvecs  # DWIのGradient Table
├── data.nii.gz  # DWI
└── nodif_brain_mask.nii.gz  # b=0のマスク

DTIの定量値を計算するために、dtifitを用いて、次のようにコマンドを実行する。

dtifit -k DTI_metrics/data -o DTI_metrics/dti -m DTI_metrics/nodif_brain_mask -r DTI_metrics/bvecs -b DTI_metrics/bvals

これで、DTIの定量値であるFAおよびMDが計算されました。

次に、定量値の計測に移ります。今回の例では、標準空間にトラクトがあるため、次のようなコマンドを実行する。デフォルトでは、トラクトの容積・確率・長さ、FAとMDの平均値のみしか計算しないが、-measオプションで計測したい定量値をリストで指定することで計測できる。

BPX_DIR='DTI.bedpostX'
xtract_stats -d DTI_metrics/dti_ -xtract XTRACT_output -w $BPX_DIR/xfms/standard2diff_warp  -r DTI_metrics/dti_FA
tract vol (mm3) median_prob mean_prob stddev_prob median_length (mm) mean_length (mm) stddev_length (mm) median_FA mean_FA stddev_FA median_MD (mm2.s-1) mean_MD (mm2.s-1) stddev_MD (mm2.s-1)
ac 3414.398 0.002678 0.017128 0.045917 28.7 29.74347 17.81838 0.521566 0.523296 0.236441 0.000546 0.000549 0.000086
af_l 20666.68 0.002328 0.003994 0.004256 42.12476 44.20898 23.76673 0.613976 0.577222 0.173791 0.000502 0.00051 0.000074
af_r 24948.77 0.002255 0.003645 0.003219 47.85464 48.74961 25.61298 0.615803 0.592681 0.1526 0.000504 0.00051 0.000188
ar_l 8164.129 0.003856 0.012093 0.022575 28.42037 36.09538 29.55519 0.55656 0.520793 0.225014 0.000536 0.00054 0.000158
ar_r 8293.718 0.002954 0.012143 0.022982 22.85014 32.43816 28.56032 0.533049 0.499066 0.221096 0.000544 0.000546 0.00018
atr_l 18289 0.002699 0.004647 0.004877 57.80597 49.69841 29.79699 0.575311 0.575069 0.165324 0.000503 0.000506 0.000063
atr_r 17821.35 0.002619 0.00464 0.005357 54.94699 49.30058 29.63312 0.557883 0.545817 0.170427 0.000514 0.000517 0.000065
cbd_l 11657.41 0.002613 0.010754 0.028017 98.31573 81.48625 42.73049 0.607537 0.596432 0.165071 0.000514 0.000522 0.000058
cbd_r 8603.606 0.002572 0.012132 0.027985 89.72475 75.92992 42.34245 0.592488 0.580259 0.179246 0.000525 0.000534 0.000068
cbp_l 2197.385 0.003905 0.022907 0.053306 22.45625 25.27225 19.48524 0.545097 0.533311 0.195338 0.000531 0.000542 0.000076
cbp_r 1295.893 0.004301 0.030165 0.056452 22.73876 25.37224 16.77534 0.584728 0.562756 0.176806 0.000554 0.000556 0.000074
cbt_l 3211.562 0.005425 0.01773 0.030143 26.74173 26.77412 18.85376 0.430775 0.435227 0.207792 0.000581 0.000573 0.000085
cbt_r 3256.637 0.005607 0.016774 0.028038 16.79167 22.45104 18.22606 0.382306 0.4093 0.212651 0.000577 0.000572 0.000237
cst_l 17832.62 0.003415 0.007274 0.012186 90.84464 85.81865 32.18107 0.640534 0.621043 0.179487 0.000484 0.000488 0.000102
cst_r 18311.54 0.003422 0.007005 0.010465 87.59891 82.39769 33.99891 0.660138 0.637923 0.176238 0.000489 0.00049 0.000082
fa_l 8981.105 0.00313 0.008087 0.014326 30.88156 35.4233 21.24802 0.551998 0.524526 0.186808 0.000505 0.000513 0.000085
fa_r 11020.73 0.00295 0.007052 0.011075 32.56344 35.19621 21.1158 0.556506 0.530259 0.18864 0.000504 0.000513 0.000133
fma 20244.11 0.002896 0.007776 0.011268 83.8533 83.33609 57.70774 0.681275 0.659728 0.211593 0.000525 0.000559 0.001264
fmi 26627.79 0.002735 0.004226 0.003915 50.42218 53.67333 34.15696 0.616961 0.618517 0.170364 0.000526 0.000536 0.000067
fx_l 7392.227 0.002952 0.015563 0.036943 61.5648 63.38259 41.55957 0.584486 0.591175 0.244622 0.000576 0.002147 0.053695
fx_r 6062.528 0.003308 0.018033 0.046488 79.0278 72.24227 42.73615 0.505915 0.527057 0.251368 0.000612 0.000642 0.001626
ifo_l 24374.07 0.00259 0.006683 0.011182 55.74806 73.31606 51.36399 0.644492 0.632296 0.150666 0.000522 0.000528 0.000061
ifo_r 24531.83 0.002684 0.006509 0.011672 63.05942 78.94564 55.58284 0.629515 0.614532 0.158353 0.000527 0.000536 0.000268
ilf_l 16356.43 0.00263 0.005763 0.008443 51.05518 48.10694 28.54763 0.610046 0.568141 0.206338 0.000554 0.000558 0.000068
ilf_r 16750.83 0.002442 0.004905 0.006768 43.45034 43.53699 24.9466 0.566343 0.532399 0.205195 0.000564 0.000571 0.000078
mcp 20756.83 0.002571 0.005435 0.006748 53.16116 59.9034 37.39916 0.638924 0.606495 0.217408 0.00044 0.000441 0.000085
mdlf_l 22621.79 0.002328 0.004845 0.006854 62.31984 57.58495 34.51005 0.619331 0.583481 0.190562 0.000524 0.001737 0.076146
mdlf_r 26210.85 0.002244 0.004559 0.006158 64.63949 60.43966 34.84705 0.618442 0.584208 0.190205 0.000529 0.000542 0.000303
or_l 15900.05 0.002753 0.005406 0.006641 58.00233 52.80807 25.41879 0.687723 0.664749 0.148922 0.000517 0.000522 0.000067
or_r 19567.99 0.002326 0.004878 0.006523 66.28878 62.15824 29.92579 0.670205 0.647372 0.15949 0.000518 0.000525 0.000071
slf1_l 11668.68 0.002888 0.009406 0.016493 44.12645 45.98651 26.8693 0.634192 0.59082 0.19284 0.000514 0.000522 0.000098
slf1_r 10479.83 0.003252 0.009862 0.019153 41.26475 42.87262 22.49848 0.600393 0.553235 0.204767 0.00053 0.000537 0.000091
slf2_l 11629.24 0.003415 0.012128 0.023267 33.56578 36.14035 24.11797 0.557181 0.528692 0.198405 0.000509 0.000514 0.000085
slf2_r 14666.13 0.00291 0.010663 0.019574 52.40123 47.55139 21.13686 0.591315 0.548816 0.193379 0.000506 0.000512 0.000094
slf3_l 12801.17 0.00283 0.008231 0.016514 42.1875 41.84864 21.03708 0.584062 0.531114 0.189168 0.00052 0.000534 0.000109
slf3_r 21331.53 0.002617 0.00664 0.011388 50.27447 52.50993 26.76177 0.591089 0.549534 0.177987 0.000517 0.000533 0.000283
str_l 17781.91 0.002877 0.004736 0.004635 45.51739 43.66699 24.95468 0.572124 0.557619 0.189184 0.000481 0.000492 0.000123
str_r 17117.06 0.002982 0.004948 0.004979 44.51994 43.55247 24.74726 0.584902 0.560258 0.183142 0.000482 0.000487 0.000082
uf_l 11189.76 0.002795 0.008007 0.01537 51.86394 47.84278 28.20727 0.532622 0.518818 0.178007 0.000562 0.000845 0.010806
uf_r 10981.29 0.002616 0.007184 0.011671 49.07641 46.89259 27.77912 0.560397 0.53906 0.170793 0.000566 0.001967 0.058334
vof_l 7223.198 0.002525 0.004833 0.005789 19.95182 19.58273 9.168275 0.578046 0.540677 0.183096 0.000504 0.00051 0.000065
vof_r 9183.941 0.002322 0.004658 0.006384 27.81706 27.70939 11.97529 0.553992 0.522878 0.201648 0.000505 0.000519 0.00008

【FSL】FLICAを用いたマルチモーダル独立成分分析


1. 目的
2. リポジトリ
3. 必要なデータ
4. 注意点
5. 実行
5.1. run_FLICA.m
6. 結果
6.1. index.html
6.2. コンポーネントごとの結果(All_in_one_page)
6.3. モダリティごとの結果(FA, MD, AD, RD, etc.)
7. GBSSへの適応
7.1. index.html
7.2. コンポーネントごとの結果(All_in_one_page)
7.3. モダリティごとの結果(FA, MD, AD, RD, etc.)


1. 目的

  • FLICAを用いた画像のマルチモーダル独立成分分析(ICA)

2. リポジトリ

FLICAに必要なスクリプトを次のコマンドで取得。今のところプライベートリポジトリにしています。ご興味のある方はコメントしてください。

git clone https://github.com/Hexans/FLICA_corrected.git

3. 必要なデータ

必要なデータは、モダリティごとに全ての被験者の画像がまとめられた4D画像である。例えば、次のようなファイルを用いる。

  • VBM: GM_mod_merg_s?.nii.gz
  • TBSS: all_??_skeletonised.nii.gz
  • FreeSurfer: ?h.thick.fsaverage.mgh and ?h.pial.area.fsaverage.mgh
  • etc.

さらに、共変量(年齢・疾患グループ・性別)との相関もみることができる。共変量ごとに1列に数値をまとめテキストファイル(.txt)として保存する。

  • 年齢: age.txt
  • 疾患グループ: group.txt
  • 性別: sex.txt

4. 注意点

NANボクセルがある場合、FLICA処理でエラーが発生する。

例えば、次のようなファイルがあったとする。

.
├── all_median_FA_skeletonised.nii.gz
├── all_median_ICVF_skeletonised.nii.gz
├── all_median_ISO_skeletonised.nii.gz
├── all_median_MD_skeletonised.nii.gz
├── all_median_MVF_skeletonised.nii.gz
└── all_median_OD_skeletonised.nii.gz

NANボクセルがあるかは、FSLのfslstatsで検出できる。-Mオプションは、non zeroボクセルにおける平均値を算出するオプションである。

Method:

fslstats <INPUT> -M

次のコマンドを実行することで、すべての画像に対して、NANボクセルがありなしをチェックする。

In:

for image in *nii.gz; do
    echo "${image} : $(fslstats $image -M)"
done

以上のコマンドを実行した例が、以下である。この場合、「all_median_ICVF_skeletonised.nii.gz」と「all_median_MVF_skeletonised.nii.gz」にNANボクセルが含まれていることが分かる。

Out:

all_median_FA_skeletonised.nii.gz : 0.157853 
all_median_ICVF_skeletonised.nii.gz : -nan 
all_median_ISO_skeletonised.nii.gz : 0.246872 
all_median_MD_skeletonised.nii.gz : 1.141417 
all_median_MVF_skeletonised.nii.gz : -nan 
all_median_OD_skeletonised.nii.gz : 0.516476

NANボクセルをゼロに置換するには、FSLのfslmathsコマンドを用いる。

Method:

fslmaths <INPUT> -nan <OUTPUT>

すべての画像に対して、NANボクセルがある場合、0に置換する処理をするには、以下のコマンドを実行する。以下の処理が完了すると、「corrected_nan」フォルダに、NANボクセルが含まれない画像が保存される。

mkdir corrected_nan
for image in *nii.gz; do
    fslmaths $image -nan corrected_nan/${image}
done

NANボクセル補正後のデータに、NANが含まれていないかチェックするには、以下のコマンドを実行する。

In:

for image in $(ls corrected_nan); do
    echo "${image} : $(fslstats corrected_nan/$image -M)"
done

「all_median_ICVF_skeletonised.nii.gz」と「all_median_MVF_skeletonised.nii.gz」にあったNANボクセルが、解消されていることが分かる。

Out:

all_median_FA_skeletonised.nii.gz : 0.157853 
all_median_ICVF_skeletonised.nii.gz : 0.422611 
all_median_ISO_skeletonised.nii.gz : 0.246872 
all_median_MD_skeletonised.nii.gz : 1.141417 
all_median_MVF_skeletonised.nii.gz : 0.086699 
all_median_OD_skeletonised.nii.gz : 0.516476

5. 実行

FLICAを実行するには、MATLABから次のコマンドを実行。

ICAのComponent数を被験者数の4分の1以下(Components < subjects/4)に設定し、Overfittingを避ける。また、被験者ごとのノイズ推定を省く(opts.lambda_dims = ”;)ことでもOverfittingを避けることができる。

5.1. run_FLICA.m

%% Set up:
addpath([getenv('FSLDIR') '/etc/matlab/'])
addpath('~/Documents/software/flica')

%% Load data
Yfiles = {'all_FA_skeletonised.nii.gz',
    'all_MD_skeletonised.nii.gz',
    'all_AD_skeletonised.nii.gz',
    'all_RD_skeletonised.nii.gz'};
% NOTE that these should be downsampled to around 20k voxels in the mask,
% per modality... hoping to increase this memory/cpu-related limitation.
[Y, fileinfo] = flica_load(Yfiles);
fileinfo.shortNames = {'FA', 'MD', 'AD', 'RD'};

%% Run FLICA
clear opts
opts.num_components = 5;
opts.maxits = 100;
% opts.lambda_dims = '';  # Reduce subject-noise estimation to avoid overfitting

Morig = flica(Y, opts);
[M, weights] = flica_reorder(Morig);

%% Save results
outdir = [pwd '/result/'];
mkdir(outdir)
flica_save_everything(outdir, M, fileinfo);

%% Produce correlation plots
clear des
des.Group = load('group.txt');
des.Age = load('age.txt');
des.Gender = load('gender.txt');
flica_posthoc_correlations(outdir, des);

実行が成功すると次のように画像を読み込み、計算が始まる。

Loading "all_FA_skeletonised.nii.gz"... 
Converting back to single... Generating mask... Flattening data matrix... Freeing memory... 
Loading "all_MD_skeletonised.nii.gz"... 
Converting back to single... Generating mask... Flattening data matrix... Freeing memory... 
Loading "all_AD_skeletonised.nii.gz"... 
Converting back to single... Generating mask... Flattening data matrix... Freeing memory... 
Loading "all_RD_skeletonised.nii.gz"... 
Converting back to single... Generating mask... Flattening data matrix... Freeing memory... 
Setup... doing PCA... more setup... done (0.72 sec)

Iteration 1... eta H W X.... lambda.... done updates (0.2 sec) F (0.061 sec)No previous lastF
Iteration 2... eta H W X.... lambda.... done updates (0.17 sec) F (0.029 sec), dF = 6.38536e+07
Iteration 3... eta H W X.... lambda.... done updates (0.16 sec) F (0.027 sec), dF = 21559.9
Iteration 4... eta H W X.... lambda.... done updates (0.13 sec) F (0.019 sec), dF = 253.743
...
Iteration 100... eta H W X.... lambda.... done updates (0.13 sec) F (0.019 sec), dF = 253.743

計算が終わったら、bashでƒƒflica_report.shを実行することで、HTML形式で結果をまとめることができる。まとめられた結果(HTML)を見るには「index.html」を開く。

flica_report.sh result  # 引数としてresultフォルダを指定

6. 結果

6.1. index.html

index.htmlの開くと次のようにFLICAの、まとめの結果を見ることができる。

6.2. コンポーネントごとの結果(All_in_one_page)

ICAの各Componentごとの結果は、index.htmlの「All_in_one_page」から確認することができる。

6.3. モダリティごとの結果(FA, MD, AD, RD, etc.)

モダリティごとの結果は、index.htmlの「モダリティ名」をクリックすることで確認できる。以下の図は、FAの結果である。

7. GBSSへの適応

GBSSの結果に対してもFLICAが使えるか、不明であったためできるか試してみた。

特に問題なく、使えそうな感じ。

7.1. index.html

7.2. コンポーネントごとの結果(All_in_one_page)

7.3. モダリティごとの結果(FA, MD, AD, RD, etc.)

【FSL】fslmathsの使ったROI解析


1. 目的
2. fslmaths
2.1. binary maskを使う場合
2.2. index maskを使う場合
3. 使用例
3.1. フォルダ構造
3.2. コード
3.3. 結果


1. 目的

  • MRIデータをROIを用いて計測

2. fslmaths

fslmathsを使った(non zeroボクセルの)平均値および標準偏差の計算方法は次のとおり。

2.1. binary maskを使う場合

[preparation]として-k <binary mask>を指定

# 平均値
fslmaths -k <index mask> <input image> -M

# 標準偏差
fslmaths -k <index mask> <input image> -S

2.2. index maskを使う場合

maskの中には、1以上の離散値を持つ場合もある(例:前頭葉: 1, 側頭葉: 2, 頭頂葉: 3, 後頭葉: 4)。このようなmaskをFSLでは「index mask」と呼ぶ。
[preparation]として-K <index mask>を指定(ラージKであることに注意)。

# 平均値
fslmaths -K <index mask> <input image> -M

# 標準偏差
fslmaths -K <index mask> <input image> -S

3. 使用例

firstというディレクトリにfactor 1.0 ~ 2.0の範囲で0.2刻みでデータ収集し、1つのfactorには5種類の画像がある。

さらに、3種類のindex maskがある(20190616_${map}.nii.gz)。

新たにデータセットを用意する場合、

  • index maskは、「_${map}.nii.gz(map =T1, T2, PD)」の部分だけ揃える
  • 画像データは、「_${i}.nii(i =1, 2, 3, 4, 5))」の部分は揃える

3.1. フォルダ構造

first/
├── 1.0
│   ├── 20190616_135659WIPCS103DSyMRIs901a1009_1.nii
│   ├── 20190616_135659WIPCS103DSyMRIs901a1009_2.nii
│   ├── 20190616_135659WIPCS103DSyMRIs901a1009_3.nii
│   ├── 20190616_135659WIPCS103DSyMRIs901a1009_4.nii
│   └── 20190616_135659WIPCS103DSyMRIs901a1009_5.nii
├── 1.2
├── 1.4
...
└──2.0
20190616_T1.nii.gz
20190616_T2.nii.gz
20190616_PD.nii.gz

3.2. コード

以上の準備ができたら、以下のコマンドをfirstフォルダのあるディレクトリと同じ場所で実行する。

計測結果は、「temp」フォルダにまとめられ、最終的なまとめ結果は「temp/summary」にまとめられる。

mkdir -p temp temp/summary
for img in $(seq 1 5);do
    echo "Calc image$img ..."
    for i in $(ls first);do
        for map in T1 T2 PD;do
            echo ${map}_img${img}_${i} > temp/${map}_img${img}_${i}.txt
            fslstats -K *_${map}.nii.gz first/$i/*_${img}.nii -M >> temp/${map}_img${img}_${i}.txt
        done
    done
    # summary
    paste temp/T1_img${img}_*.txt > temp/summary/T1_img${img}.txt
    paste temp/T2_img${img}_*.txt > temp/summary/T2_img${img}.txt
    paste temp/PD_img${img}_*.txt > temp/summary/PD_img${img}.txt
done

3.3. 結果

「temp/summary」フォルダの中身を確認すると次のようなファイルが生成されている。

これらの結果は、各index maskおよび各factorにおける5種類の画像ごとに計測結果をまとめたものである。

$ ls temp/summary/
PD_img1.txt  PD_img4.txt  T1_img2.txt  T1_img5.txt  T2_img3.txt
PD_img2.txt  PD_img5.txt  T1_img3.txt  T2_img1.txt  T2_img4.txt
PD_img3.txt  T1_img1.txt  T1_img4.txt  T2_img2.txt  T2_img5.txt

「PD_img1.txt」の中身は次の通り。

PD_1_1.0	PD_1_1.2	PD_1_1.4	PD_1_1.6	PD_1_1.8	PD_1_10	PD_1_2.0	PD_1_2.2	PD_1_2.4	PD_1_2.6	PD_1_2.8	PD_1_3.0	PD_1_3.5	PD_1_4.0	PD_1_4.5	PD_1_5.0	PD_1_5.5	PD_1_6.0	PD_1_6.5	PD_1_7.0	PD_1_7.5	PD_1_8.0	PD_1_8.5	PD_1_9.0	PD_1_9.5
56.071172 	53.571421 	57.104340 	65.576766 	56.634142 	164.695774 	59.164714 	66.530548 	64.210779 	72.075169 	88.882530 	68.767341 	78.006755 	82.362270 	102.866073 	95.148962 	157.649307 	135.169850 		136.416867 	209.898826 	143.974304 	206.080797 	144.660493 	152.182925 
59.219002 	59.634794 	67.469257 	57.990326 	68.822425 	135.478640 	61.384532 	70.464787 	71.207407 	73.836757 	85.293909 	69.715787 	85.544729 	98.812038 	120.032300 	113.378021 	137.899509 	162.424688 		186.341014 	182.449928 	168.822476 	226.915020 	168.971066 	156.942100 
65.593938 	73.217330 	78.069650 	72.319718 	72.329499 	229.141061 	72.960663 	85.006750 	71.576984 	86.298830 	101.787668 	68.924592 	111.201372 	124.327963 	136.776286 	160.195031 	187.189908 	202.490824 		222.796888 	278.856697 	198.264230 	250.119406 	227.780093 	328.350670 
60.599368 	73.082315 	75.607594 	78.199978 	77.303263 	264.188672 	82.484480 	87.064013 	90.335936 	86.806134 	99.869636 	91.716010 	100.842814 	124.390488 	152.835757 	158.187163 	190.057344 	231.426025 		255.339445 	233.320259 	210.198541 	268.658566 	277.634037 	304.397152 
62.596507 	73.224809 	73.435417 	74.447392 	82.949692 	200.789283 	95.077401 	76.603784 	91.650407 	86.794202 	97.143067 	96.004177 	91.693638 	114.771571 	147.396567 	128.985783 	156.368281 	180.671239 		219.410221 	212.203824 	191.402947 	209.449536 	208.555533 	208.911982 
70.904297 	74.822798 	82.031894 	81.959065 	81.159454 	204.884733 	84.416149 	84.009042 	105.231751 	102.618848 	109.251246 	87.436376 	87.707241 	103.725383 	124.369558 	148.823854 	178.081095 	170.385173 		216.685645 	324.138544 	221.210790 	270.552265 	221.176516 	209.839571 
73.786061 	76.172245 	78.538106 	88.747996 	75.945391 	210.691364 	96.988982 	93.988295 	92.378013 	94.479654 	102.801537 	105.931549 	100.643631 	119.402497 	133.575660 	158.669452 	168.026760 	208.149685 		216.036306 	200.196084 	221.846031 	254.239667 	222.523944 	208.206735 
87.495523 	88.813652 	87.880272 	85.964579 	98.281155 	349.093039 	98.692713 	94.235417 	101.215252 	112.179749 	109.359223 	92.050666 	103.862347 	126.171527 	137.595066 	170.217829 	185.742386 	242.152668 		271.305374 	222.041778 	273.485643 	306.670516 	315.758490 	264.002075 
98.739737 	100.720432 	97.328787 	99.053358 	97.533872 	316.762011 	103.369259 	113.277178 	99.181052 	109.115135 	118.181639 	109.618875 	118.613545 	125.297607 	148.418087 	162.936641 	185.542536 	241.243367 		219.280341 	239.735891 	205.617819 	270.702777 	229.903270 	235.547812 
103.277802 	106.846910 	104.091932 	106.542414 	102.750478 	207.280336 	102.150977 	113.678573 	109.015869 	102.060726 	112.479080 	106.674224 	115.947559 	116.607501 	144.210974 	143.748899 	149.864118 	209.414898 		207.632018 	212.082947 	221.947801 	243.473957 	262.798043 	187.095868 
126.149411 	128.146098 	123.403289 	129.219018 	120.576404 	307.939330 	127.342593 	130.534624 	136.748147 	135.015205 	130.025685 	127.302044 	129.562001 	181.892504 	175.129838 	214.058127 	226.497845 	261.011736 		293.147569 	369.290946 	318.991315 	334.730611 	391.036790 	375.232291 
142.132261 	136.881755 	140.397574 	139.790297 	137.662219 	372.306719 	144.546539 	141.836003 	157.273254 	146.306116 	172.242885 	156.343952 	157.901801 	186.644997 	198.501073 	221.079046 	239.737134 	273.049238 		337.570792 	359.821785 	327.865160 	385.884041 	403.470958 	392.444883 
153.414811 	149.560298 	153.284503 	153.066674 	161.247593 	389.665428 	167.865698 	162.732419 	163.929323 	170.265862 	175.904991 	159.839045 	172.798252 	191.501041 	201.477961 	243.534578 	247.816996 	296.709427 		355.380479 	385.017319 	365.828550 	409.998897 	380.333908 	410.494708 
123.627007 	125.429929 	125.827541 	128.563233 	124.580074 	305.157273 	133.259714 	136.405932 	140.542948 	138.670648 	136.025625 	134.278958 	149.233630 	156.115582 	180.495762 	207.980796 	201.333914 	247.976671 		304.193507 	306.642640 	293.506436 	327.611026 	335.249888 	354.644789 

これをエクセルに貼り付けると、次のようになる。行はindex maskのindexに対応し、列は各factorに対応する。factorのネーミングの仕方が悪く「PD_1_10」が6列目に出てきてしまっていることに注意。

PD_1_1.0 PD_1_1.2 PD_1_1.4 PD_1_1.6 PD_1_1.8 PD_1_10 PD_1_2.0 PD_1_2.2 PD_1_2.4 PD_1_2.6 PD_1_2.8 PD_1_3.0 PD_1_3.5 PD_1_4.0 PD_1_4.5 PD_1_5.0 PD_1_5.5 PD_1_6.0 PD_1_6.5 PD_1_7.0 PD_1_7.5 PD_1_8.0 PD_1_8.5 PD_1_9.0 PD_1_9.5
56.07117 53.57142 57.10434 65.57677 56.63414 164.6958 59.16471 66.53055 64.21078 72.07517 88.88253 68.76734 78.00676 82.36227 102.8661 95.14896 157.6493 135.1699 102.8661 136.4169 209.8988 143.9743 206.0808 144.6605 152.1829
59.219 59.63479 67.46926 57.99033 68.82243 135.4786 61.38453 70.46479 71.20741 73.83676 85.29391 69.71579 85.54473 98.81204 120.0323 113.378 137.8995 162.4247 120.0323 186.341 182.4499 168.8225 226.915 168.9711 156.9421
65.59394 73.21733 78.06965 72.31972 72.3295 229.1411 72.96066 85.00675 71.57698 86.29883 101.7877 68.92459 111.2014 124.328 136.7763 160.195 187.1899 202.4908 136.7763 222.7969 278.8567 198.2642 250.1194 227.7801 328.3507
60.59937 73.08232 75.60759 78.19998 77.30326 264.1887 82.48448 87.06401 90.33594 86.80613 99.86964 91.71601 100.8428 124.3905 152.8358 158.1872 190.0573 231.426 152.8358 255.3394 233.3203 210.1985 268.6586 277.634 304.3972
62.59651 73.22481 73.43542 74.44739 82.94969 200.7893 95.0774 76.60378 91.65041 86.7942 97.14307 96.00418 91.69364 114.7716 147.3966 128.9858 156.3683 180.6712 147.3966 219.4102 212.2038 191.4029 209.4495 208.5555 208.912
70.9043 74.8228 82.03189 81.95907 81.15945 204.8847 84.41615 84.00904 105.2318 102.6188 109.2512 87.43638 87.70724 103.7254 124.3696 148.8239 178.0811 170.3852 124.3696 216.6856 324.1385 221.2108 270.5523 221.1765 209.8396
73.78606 76.17225 78.53811 88.748 75.94539 210.6914 96.98898 93.9883 92.37801 94.47965 102.8015 105.9315 100.6436 119.4025 133.5757 158.6695 168.0268 208.1497 133.5757 216.0363 200.1961 221.846 254.2397 222.5239 208.2067
87.49552 88.81365 87.88027 85.96458 98.28116 349.093 98.69271 94.23542 101.2153 112.1797 109.3592 92.05067 103.8623 126.1715 137.5951 170.2178 185.7424 242.1527 137.5951 271.3054 222.0418 273.4856 306.6705 315.7585 264.0021
98.73974 100.7204 97.32879 99.05336 97.53387 316.762 103.3693 113.2772 99.18105 109.1151 118.1816 109.6189 118.6135 125.2976 148.4181 162.9366 185.5425 241.2434 148.4181 219.2803 239.7359 205.6178 270.7028 229.9033 235.5478
103.2778 106.8469 104.0919 106.5424 102.7505 207.2803 102.151 113.6786 109.0159 102.0607 112.4791 106.6742 115.9476 116.6075 144.211 143.7489 149.8641 209.4149 144.211 207.632 212.0829 221.9478 243.474 262.798 187.0959
126.1494 128.1461 123.4033 129.219 120.5764 307.9393 127.3426 130.5346 136.7481 135.0152 130.0257 127.302 129.562 181.8925 175.1298 214.0581 226.4978 261.0117 175.1298 293.1476 369.2909 318.9913 334.7306 391.0368 375.2323
142.1323 136.8818 140.3976 139.7903 137.6622 372.3067 144.5465 141.836 157.2733 146.3061 172.2429 156.344 157.9018 186.645 198.5011 221.079 239.7371 273.0492 198.5011 337.5708 359.8218 327.8652 385.884 403.471 392.4449
153.4148 149.5603 153.2845 153.0667 161.2476 389.6654 167.8657 162.7324 163.9293 170.2659 175.905 159.839 172.7983 191.501 201.478 243.5346 247.817 296.7094 201.478 355.3805 385.0173 365.8286 409.9989 380.3339 410.4947
123.627 125.4299 125.8275 128.5632 124.5801 305.1573 133.2597 136.4059 140.5429 138.6706 136.0256 134.279 149.2336 156.1156 180.4958 207.9808 201.3339 247.9767 180.4958 304.1935 306.6426 293.5064 327.611 335.2499 354.6448

【FSL】MRI画像のレジストレーション(Registration)と関心領域(ROI)解析



1. 目的
2. JHU White-matter labels & tractography atlasでROI解析する方法
2.1. ディレクトリ構造
2.2. 実行方法
2.3. ソースコード
3. Desikan Killiany AtlasでROI解析する方法
3.1. 前提条件
3.2. ディレクトリ構造
3.3. 実行方法
3.4. ソースコード
4. FreeSurferのwm.seg.nii.gzをDWI空間に位置合わせする方法
4.1. 必要なファイル
4.2. ソースコード
5. 標準空間にあるROIを各被験者脳に位置合わせ
6. bertのwmparcをMNI空間に移動後、個人脳(Perfusion image)にレジストレーション


1. 目的

  • JHU White-matter labels & tractography atlasでROI解析する方法
  • Desikan Killiany AtlasでROI解析
  • FreeSurferのwm.seg.nii.gzをDWI空間に位置合わせする方法

2. JHU White-matter labels & tractography atlasでROI解析する方法

標準空間(template (standard) space)にあるJHU-ICBM-labels-1mm.nii.gz、JHU-ICBM-tracts-maxprob-thr25-1mm.nii.gz をDWI空間に位置合わせしてROI解析するスクリプト。

2.1. ディレクトリ構造

最低限 FA.nii.gz は用意。

SUBJ1/
 + map/
    + FA.nii.gz
    + MD.nii.gz
    + AD.nii.gz
       :
SUBJ2/
 + map/
    + FA.nii.gz
    + MD.nii.gz
    + AD.nii.gz
       :

2.2. 実行方法

cd SUBJ1
bash ../reg_FMRIB_to_FA.sh

2.3. ソースコード

function reg_FMRIB-Atlas_to_FA () {
    flirt -in map/FA -ref ${FSLDIR}/data/standard/FMRIB58_FA_1mm -omat FA_to_FMRIB.mat
    fnirt --in=map/FA --aff=FA_to_FMRIB.mat --config=FA_2_FMRIB58_1mm --cout=warp_FA_to_FMRIB
    invwarp -w warp_FA_to_FMRIB -o warp_FMRIB_to_FA -r map/FA

    applywarp --in=${FSLDIR}/data/atlases/JHU/JHU-ICBM-labels-1mm --ref=map/FA --warp=warp_FMRIB_to_FA --interp=nn --out=JHU-ICBM-labels_to_FA
    applywarp --in=${FSLDIR}/data/atlases/JHU/JHU-ICBM-tracts-maxprob-thr25-1mm --ref=map/FA --warp=warp_FMRIB_to_FA --interp=nn --out=JHU-ICBM-tracts-maxprob-thr25_to_FA
}

function split_JHU-ICBM-labels () {
    OUT_DIR=mask/JHU-ICBM-labels
    mkdir -p ${OUT_DIR}

    fslmaths JHU-ICBM-labels_to_FA -thr  1 -uthr  1 ${OUT_DIR}/01_Middle-cerebellar-peduncle
    fslmaths JHU-ICBM-labels_to_FA -thr  2 -uthr  2 ${OUT_DIR}/02_Pontine-crossing-tract
    fslmaths JHU-ICBM-labels_to_FA -thr  3 -uthr  3 ${OUT_DIR}/03_Genu-of-corpus-callosum
    fslmaths JHU-ICBM-labels_to_FA -thr  4 -uthr  4 ${OUT_DIR}/04_Body-of-corpus-callosum
    fslmaths JHU-ICBM-labels_to_FA -thr  5 -uthr  5 ${OUT_DIR}/05_Splenium-of-corpus-callosum
    fslmaths JHU-ICBM-labels_to_FA -thr  6 -uthr  6 ${OUT_DIR}/06_Fornix
    fslmaths JHU-ICBM-labels_to_FA -thr  7 -uthr  7 ${OUT_DIR}/07_Corticospinal-tract_R
    fslmaths JHU-ICBM-labels_to_FA -thr  8 -uthr  8 ${OUT_DIR}/08_Corticospinal-tract_L
    fslmaths JHU-ICBM-labels_to_FA -thr  9 -uthr  9 ${OUT_DIR}/09_Medial-lemniscus_R
    fslmaths JHU-ICBM-labels_to_FA -thr 10 -uthr 10 ${OUT_DIR}/10_Medial-lemniscus_L
    fslmaths JHU-ICBM-labels_to_FA -thr 11 -uthr 11 ${OUT_DIR}/11_Inferior-cerebellar-peduncle_R
    fslmaths JHU-ICBM-labels_to_FA -thr 12 -uthr 12 ${OUT_DIR}/12_Inferior-cerebellar-peduncle_L
    fslmaths JHU-ICBM-labels_to_FA -thr 13 -uthr 13 ${OUT_DIR}/13_Superior-cerebellar-peduncle_R
    fslmaths JHU-ICBM-labels_to_FA -thr 14 -uthr 14 ${OUT_DIR}/14_Superior-cerebellar-peduncle_L
    fslmaths JHU-ICBM-labels_to_FA -thr 15 -uthr 15 ${OUT_DIR}/15_Cerebral-peduncle_R
    fslmaths JHU-ICBM-labels_to_FA -thr 16 -uthr 16 ${OUT_DIR}/16_Cerebral-peduncle_L
    fslmaths JHU-ICBM-labels_to_FA -thr 17 -uthr 17 ${OUT_DIR}/17_Anterior-limb-of-internal-capsule_R
    fslmaths JHU-ICBM-labels_to_FA -thr 18 -uthr 18 ${OUT_DIR}/18_Anterior-limb-of-internal-capsule_L
    fslmaths JHU-ICBM-labels_to_FA -thr 19 -uthr 19 ${OUT_DIR}/19_Posterior-limb-of-internal-capsule_R
    fslmaths JHU-ICBM-labels_to_FA -thr 20 -uthr 20 ${OUT_DIR}/20_Posterior-limb-of-internal-capsule_L
    fslmaths JHU-ICBM-labels_to_FA -thr 21 -uthr 21 ${OUT_DIR}/21_Retrolenticular-part-of-internal-capsule_R
    fslmaths JHU-ICBM-labels_to_FA -thr 22 -uthr 22 ${OUT_DIR}/22_Retrolenticular-part-of-internal-capsule_L
    fslmaths JHU-ICBM-labels_to_FA -thr 23 -uthr 23 ${OUT_DIR}/23_Anterior-corona-radiata_R
    fslmaths JHU-ICBM-labels_to_FA -thr 24 -uthr 24 ${OUT_DIR}/24_Anterior-corona-radiata_L
    fslmaths JHU-ICBM-labels_to_FA -thr 25 -uthr 25 ${OUT_DIR}/25_Superior-corona-radiata_R
    fslmaths JHU-ICBM-labels_to_FA -thr 26 -uthr 26 ${OUT_DIR}/26_Superior-corona-radiata_L
    fslmaths JHU-ICBM-labels_to_FA -thr 27 -uthr 27 ${OUT_DIR}/27_Posterior-corona-radiata_R
    fslmaths JHU-ICBM-labels_to_FA -thr 28 -uthr 28 ${OUT_DIR}/28_Posterior-corona-radiata_L
    fslmaths JHU-ICBM-labels_to_FA -thr 29 -uthr 29 ${OUT_DIR}/29_Posterior-thalamic-radiation_R
    fslmaths JHU-ICBM-labels_to_FA -thr 30 -uthr 30 ${OUT_DIR}/30_Posterior-thalamic-radiation_L
    fslmaths JHU-ICBM-labels_to_FA -thr 31 -uthr 31 ${OUT_DIR}/31_Sagittal-stratum_R
    fslmaths JHU-ICBM-labels_to_FA -thr 32 -uthr 32 ${OUT_DIR}/32_Sagittal-stratum_L
    fslmaths JHU-ICBM-labels_to_FA -thr 33 -uthr 33 ${OUT_DIR}/33_External-capsule_R
    fslmaths JHU-ICBM-labels_to_FA -thr 34 -uthr 34 ${OUT_DIR}/34_External-capsule_L
    fslmaths JHU-ICBM-labels_to_FA -thr 35 -uthr 35 ${OUT_DIR}/35_Cingulum-cingulate-gyrus_R
    fslmaths JHU-ICBM-labels_to_FA -thr 36 -uthr 36 ${OUT_DIR}/36_Cingulum-cingulate-gyrus_L
    fslmaths JHU-ICBM-labels_to_FA -thr 37 -uthr 37 ${OUT_DIR}/37_Cingulum-hippocampus_R
    fslmaths JHU-ICBM-labels_to_FA -thr 38 -uthr 38 ${OUT_DIR}/38_Cingulum-hippocampus_L
    fslmaths JHU-ICBM-labels_to_FA -thr 39 -uthr 39 ${OUT_DIR}/39_Fornix-Stria-terminalis_R
    fslmaths JHU-ICBM-labels_to_FA -thr 40 -uthr 40 ${OUT_DIR}/40_Fornix-Stria-terminalis_L
    fslmaths JHU-ICBM-labels_to_FA -thr 41 -uthr 41 ${OUT_DIR}/41_Superior-longitudinal-fasciculus_R
    fslmaths JHU-ICBM-labels_to_FA -thr 42 -uthr 42 ${OUT_DIR}/42_Superior-longitudinal-fasciculus_L
    fslmaths JHU-ICBM-labels_to_FA -thr 43 -uthr 43 ${OUT_DIR}/43_Superior-fronto-occipital-fasciculus_R
    fslmaths JHU-ICBM-labels_to_FA -thr 44 -uthr 44 ${OUT_DIR}/44_Superior-fronto-occipital-fasciculus_L
    fslmaths JHU-ICBM-labels_to_FA -thr 45 -uthr 45 ${OUT_DIR}/45_Uncinate-fasciculus_R
    fslmaths JHU-ICBM-labels_to_FA -thr 46 -uthr 46 ${OUT_DIR}/46_Uncinate-fasciculus_L
    fslmaths JHU-ICBM-labels_to_FA -thr 47 -uthr 47 ${OUT_DIR}/47_Tapetum_R
    fslmaths JHU-ICBM-labels_to_FA -thr 48 -uthr 48 ${OUT_DIR}/48_Tapetum_L

    fsladd ${OUT_DIR}/51_Corticospinal-tract ${OUT_DIR}/07_Corticospinal-tract_R ${OUT_DIR}/08_Corticospinal-tract_L > /dev/null
    fsladd ${OUT_DIR}/52_Medial-lemniscus ${OUT_DIR}/09_Medial-lemniscus_R ${OUT_DIR}/10_Medial-lemniscus_L > /dev/null
    fsladd ${OUT_DIR}/53_Inferior-cerebellar-peduncle ${OUT_DIR}/11_Inferior-cerebellar-peduncle_R ${OUT_DIR}/12_Inferior-cerebellar-peduncle_L > /dev/null
    fsladd ${OUT_DIR}/54_Superior-cerebellar-peduncle ${OUT_DIR}/13_Superior-cerebellar-peduncle_R ${OUT_DIR}/14_Superior-cerebellar-peduncle_L > /dev/null
    fsladd ${OUT_DIR}/55_Cerebral-peduncle ${OUT_DIR}/15_Cerebral-peduncle_R ${OUT_DIR}/16_Cerebral-peduncle_L > /dev/null
    fsladd ${OUT_DIR}/56_Anterior-limb-of-internal-capsule ${OUT_DIR}/17_Anterior-limb-of-internal-capsule_R ${OUT_DIR}/18_Anterior-limb-of-internal-capsule_L > /dev/null
    fsladd ${OUT_DIR}/57_Posterior-limb-of-internal-capsule ${OUT_DIR}/19_Posterior-limb-of-internal-capsule_R ${OUT_DIR}/20_Posterior-limb-of-internal-capsule_L > /dev/null
    fsladd ${OUT_DIR}/58_Retrolenticular-part-of-internal-capsule ${OUT_DIR}/21_Retrolenticular-part-of-internal-capsule_R ${OUT_DIR}/22_Retrolenticular-part-of-internal-capsule_L > /dev/null
    fsladd ${OUT_DIR}/59_Anterior-corona-radiata ${OUT_DIR}/23_Anterior-corona-radiata_R ${OUT_DIR}/24_Anterior-corona-radiata_L > /dev/null
    fsladd ${OUT_DIR}/60_Superior-corona-radiata ${OUT_DIR}/25_Superior-corona-radiata_R ${OUT_DIR}/26_Superior-corona-radiata_L > /dev/null
    fsladd ${OUT_DIR}/61_Posterior-corona-radiata ${OUT_DIR}/27_Posterior-corona-radiata_R ${OUT_DIR}/28_Posterior-corona-radiata_L > /dev/null
    fsladd ${OUT_DIR}/62_Posterior-thalamic-radiation ${OUT_DIR}/29_Posterior-thalamic-radiation_R ${OUT_DIR}/30_Posterior-thalamic-radiation_L > /dev/null
    fsladd ${OUT_DIR}/63_Sagittal-stratum ${OUT_DIR}/31_Sagittal-stratum_R ${OUT_DIR}/32_Sagittal-stratum_L > /dev/null
    fsladd ${OUT_DIR}/64_External-capsule ${OUT_DIR}/33_External-capsule_R ${OUT_DIR}/34_External-capsule_L > /dev/null
    fsladd ${OUT_DIR}/65_Cingulum-cingulate-gyrus ${OUT_DIR}/35_Cingulum-cingulate-gyrus_R ${OUT_DIR}/36_Cingulum-cingulate-gyrus_L > /dev/null
    fsladd ${OUT_DIR}/66_Cingulum-hippocampus ${OUT_DIR}/37_Cingulum-hippocampus_R ${OUT_DIR}/38_Cingulum-hippocampus_L > /dev/null
    fsladd ${OUT_DIR}/67_Fornix-Stria-terminalis ${OUT_DIR}/39_Fornix-Stria-terminalis_R ${OUT_DIR}/40_Fornix-Stria-terminalis_L > /dev/null
    fsladd ${OUT_DIR}/68_Superior-longitudinal-fasciculus ${OUT_DIR}/41_Superior-longitudinal-fasciculus_R ${OUT_DIR}/42_Superior-longitudinal-fasciculus_L > /dev/null
    fsladd ${OUT_DIR}/69_Superior-fronto-occipital-fasciculus ${OUT_DIR}/43_Superior-fronto-occipital-fasciculus_R ${OUT_DIR}/44_Superior-fronto-occipital-fasciculus_L > /dev/null
    fsladd ${OUT_DIR}/70_Uncinate-fasciculus ${OUT_DIR}/45_Uncinate-fasciculus_R ${OUT_DIR}/46_Uncinate-fasciculus_L > /dev/null
    fsladd ${OUT_DIR}/71_Tapetum ${OUT_DIR}/47_Tapetum_R ${OUT_DIR}/48_Tapetum_L > /dev/null

    fsladd ${OUT_DIR}/81_Corpus-callosum ${OUT_DIR}/03_Genu-of-corpus-callosum ${OUT_DIR}/04_Body-of-corpus-callosum ${OUT_DIR}/05_Splenium-of-corpus-callosum > /dev/null
    fsladd ${OUT_DIR}/82_Internal-capsule ${OUT_DIR}/56_Anterior-limb-of-internal-capsule ${OUT_DIR}/57_Posterior-limb-of-internal-capsule ${OUT_DIR}/58_Retrolenticular-part-of-internal-capsule > /dev/null
    fsladd ${OUT_DIR}/83_Corona-radiata ${OUT_DIR}/59_Anterior-corona-radiata ${OUT_DIR}/60_Superior-corona-radiata ${OUT_DIR}/61_Posterior-corona-radiata > /dev/null
}

function split_JHU-ICBM-tracts-maxprob-thr25 () {
    OUT_DIR=mask/JHU-ICBM-tracts-maxprob-thr25
    mkdir -p ${OUT_DIR}

    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr  1 -uthr  1 ${OUT_DIR}/01_Anterior-thalamic_Radiation_L
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr  2 -uthr  2 ${OUT_DIR}/02_Anterior-thalamic_Radiation_R
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr  3 -uthr  3 ${OUT_DIR}/03_Corticospinal-tract_L
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr  4 -uthr  4 ${OUT_DIR}/04_Corticospinal-tract_R
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr  5 -uthr  5 ${OUT_DIR}/05_Cingulum-cingulate-gyrus_L
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr  6 -uthr  6 ${OUT_DIR}/06_Cingulum-cingulate-gyrus_R
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr  7 -uthr  7 ${OUT_DIR}/07_Cingulum-hippocampus_L
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr  8 -uthr  8 ${OUT_DIR}/08_Cingulum-hippocampus_R
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr  9 -uthr  9 ${OUT_DIR}/09_Forceps-major
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr 10 -uthr 10 ${OUT_DIR}/10_Forceps-minor
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr 11 -uthr 11 ${OUT_DIR}/11_Inferior-fronto-occipital-fasciculus_L
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr 12 -uthr 12 ${OUT_DIR}/12_Inferior-fronto-occipital-fasciculus_R
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr 13 -uthr 13 ${OUT_DIR}/13_Inferior_Longitudinal-fasciculus_L
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr 14 -uthr 14 ${OUT_DIR}/14_Inferior_Longitudinal-fasciculus_R
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr 15 -uthr 15 ${OUT_DIR}/15_Superior_Longitudinal-fasciculus_L
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr 16 -uthr 16 ${OUT_DIR}/16_Superior_Longitudinal-fasciculus_R
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr 17 -uthr 17 ${OUT_DIR}/17_Uncinate-fasciculus_L
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr 18 -uthr 18 ${OUT_DIR}/18_Uncinate-fasciculus_R
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr 19 -uthr 19 ${OUT_DIR}/19_Superior_Longitudinal-fasciculus-temporal-part_L
    fslmaths JHU-ICBM-tracts-maxprob-thr25_to_FA -thr 20 -uthr 20 ${OUT_DIR}/20_Superior_Longitudinal-fasciculus-temporal-part_R

    fsladd ${OUT_DIR}/21_Anterior-thalamic_Radiation ${OUT_DIR}/01_Anterior-thalamic_Radiation_L ${OUT_DIR}/02_Anterior-thalamic_Radiation_R > /dev/null
    fsladd ${OUT_DIR}/22_Corticospinal-tract ${OUT_DIR}/03_Corticospinal-tract_L ${OUT_DIR}/04_Corticospinal-tract_R > /dev/null
    fsladd ${OUT_DIR}/23_Cingulum-cingulate-gyrus ${OUT_DIR}/05_Cingulum-cingulate-gyrus_L ${OUT_DIR}/06_Cingulum-cingulate-gyrus_R > /dev/null
    fsladd ${OUT_DIR}/24_Cingulum-hippocampus ${OUT_DIR}/07_Cingulum-hippocampus_L ${OUT_DIR}/08_Cingulum-hippocampus_R > /dev/null
    fsladd ${OUT_DIR}/25_Forceps ${OUT_DIR}/09_Forceps-major ${OUT_DIR}/10_Forceps-minor > /dev/null
    fsladd ${OUT_DIR}/26_Inferior-fronto-occipital-fasciculus ${OUT_DIR}/11_Inferior-fronto-occipital-fasciculus_L ${OUT_DIR}/12_Inferior-fronto-occipital-fasciculus_R > /dev/null
    fsladd ${OUT_DIR}/27_Inferior_Longitudinal-fasciculus ${OUT_DIR}/13_Inferior_Longitudinal-fasciculus_L ${OUT_DIR}/14_Inferior_Longitudinal-fasciculus_R > /dev/null
    fsladd ${OUT_DIR}/28_Superior_Longitudinal-fasciculus ${OUT_DIR}/15_Superior_Longitudinal-fasciculus_L ${OUT_DIR}/16_Superior_Longitudinal-fasciculus_R > /dev/null
    fsladd ${OUT_DIR}/29_Uncinate-fasciculus ${OUT_DIR}/17_Uncinate-fasciculus_L ${OUT_DIR}/18_Uncinate-fasciculus_R > /dev/null
    fsladd ${OUT_DIR}/30_Superior_Longitudinal-fasciculus-temporal-part ${OUT_DIR}/19_Superior_Longitudinal-fasciculus-temporal-part_L ${OUT_DIR}/20_Superior_Longitudinal-fasciculus-temporal-part_R > /dev/null
}

function output_results () {
    [[ -d results ]] && \rm -rf results
    mkdir -p results

    for MAP in $(\ls -1 map/*.nii.gz | xargs -i basename -s ".nii.gz" {}); do
        for MASK in $(\ls -1 mask/JHU-ICBM-labels/*.nii.gz); do
            fslstats map/${MAP} -k ${MASK} -M >> results/JHU-ICBM-labels_${MAP}_mean.txt &
            fslstats map/${MAP} -k ${MASK} -S >> results/JHU-ICBM-labels_${MAP}_sd.txt &
            wait
        done
    done

    for MAP in $(\ls -1 map/*.nii.gz | xargs -i basename -s ".nii.gz" {}); do
        for MASK in $(\ls -1 mask/JHU-ICBM-tracts-maxprob-thr25/*.nii.gz); do
            fslstats map/${MAP} -k ${MASK} -M >> results/JHU-ICBM-tracts-maxprob-thr25_${MAP}_mean.txt &
            fslstats map/${MAP} -k ${MASK} -S >> results/JHU-ICBM-tracts-maxprob-thr25_${MAP}_sd.txt &
            wait
        done
    done
}

reg_FMRIB-Atlas_to_FA
split_JHU-ICBM-labels
split_JHU-ICBM-tracts-maxprob-thr25
output_results

結果のまとめ方

ls |grep -v temp_ID|tr "\n" "\t" | sed -e 's/^/\t/' | sed -e 's/$/\n/'> temp_ID
ls $(ls |head -n 1)/mask/JHU-ICBM-labels > temp_labels_name
ls $(ls |head -n 1)/mask/JHU-ICBM-tracts-maxprob-thr25 > temp_tracts_name
paste temp_labels_name */results/JHU-ICBM-labels_FA_mean.txt > temp_labels_value
paste temp_tracts_name */results/JHU-ICBM-tracts-maxprob-thr25_FA_mean.txt > temp_tracts_value

cat temp_ID temp_labels_value > labels_result.csv
cat temp_ID temp_tracts_value > tracts_result.csv

rm temp* 

根本先生の解説

斎藤先生の紹介してくれたところを見てくださったらと思うのですが、ポイントは、FSLでは、

invwarp

というコマンドがあるんですね。
この inv は、inverse で、普段の warp は 患者空間 → MNI空間であるものを、invwarp は、MNI空間 → 患者空間に戻します。

じゃ、最初の 患者空間 → MNI空間のパラメータをどうとるのかというところになるわけですが、斎藤先生のスクリプトの2行が肝になります。

flirt -in map/FA -ref ${FSLDIR}/data/standard/FMRIB58_FA_1mm -omat FA_to_FMRIB.mat

これは、患者のFA画像 (map/FA) を入力として、MNI空間に線形変換をし、パラメータを FA_to_FMRIB.mat で出力します。

fnirt --in=map/FA --aff=FA_to_FMRIB.mat --config=FA_2_FMRIB58_1mm --cout=warp_FA_to_FMRIB 

こちらは、入力は同じですが、先程得られた線形変換のパラメータも利用し、
FSLでもともと準備されているFA画像をMNI空間に非線形変換するのに適した設定 (FA_2_FMRIB58_1mm) を利用して
非線形変換を行い、その変換に必要なパラメータを warp_FA_to_FMRIB として出力します。
なので、出力されるものは画像ではなく、変形パラメータになります。

そして、

invwarp -w warp_FA_to_FMRIB -o warp_FMRIB_to_FA -r map/FA

で、非線形変換で得られた warp_FA_to_FMRIB の逆変換を計算して、warp_FMRIB_to_FA として保存します。

そして、最後に、

applywarp --in=${FSLDIR}/data/atlases/JHU/JHU-ICBM-labels-1mm --ref=map/FA --warp=warp_FMRIB_to_FA --interp=nn --out=JHU-ICBM-labels_to_FA
applywarp --in=${FSLDIR}/data/atlases/JHU/JHU-ICBM-tracts-maxprob-thr25-1mm --ref=map/FA --warp=warp_FMRIB_to_FA --interp=nn --out=JHU-ICBM-tracts-maxprob-thr25_to_FA

で、2つのアトラス (JHU-ICBM-labels-1mm と JHU-ICBM-tracts-maxprob-thr25-1mm) に対して、その逆変換したパラメータをあてて、患者空間の(患者用の)アトラスを作成しています。
その後、各アトラスの個々の領域の値を求めています。

3. Desikan Killiany AtlasでROI解析する方法

T1WI空間にあるDesikan Killiany AtlasをDWI空間に位置合わせをしてROI解析するスクリプト。

3.1. 前提条件

「/opt/mrtrix3/share/mrtrix3/labelconvert/fs_default.txt」がある前提で走ります。

3.2. ディレクトリ構造

ディレクトリ構造は次の通り。

SUBJ_N
├── 3DT1.nii.gz
├── aparc+aseg.mgz ... FreeSurfer output
├── b0.nii.gz      ... b0 (ex: `fslroi eddy_result b0 0 1`)
├── wm.seg.mgz     ... FreeSurfer output
└── map
   ├── FA.nii.gz
   ├── MD.nii.gz
   └── uFA.nii.gz

$ ls -1
this_script.sh
SUBJ_1
SUBJ_2
SUBJ_3

3.3. 実行方法

$ cd SUBJ_1
$ bash ../this_script.sh

3.4. ソースコード

#!/bin/bash

# SUBJ_N
# ├── 3DT1.nii.gz
# ├── aparc+aseg.mgz ... FreeSurfer output
# ├── b0.nii.gz      ... b0 (ex: `fslroi eddy_result b0 0 1`)
# ├── wm.seg.mgz     ... FreeSurfer output
# └── map
#     ├── FA.nii.gz
#     ├── MD.nii.gz
#     └── uFA.nii.gz
#
# $ ls -1
# this_script.sh
# SUBJ_1
# SUBJ_2
# SUBJ_3
#  :
# $ cd SUBJ_1
# $ bash ../this_script.sh

function create_nodes () {
    mkdir -p mask/DK

    fslreorient2std 3DT1 T1
    bet T1 T1_brain -B -f 0.1
    mri_label2vol --temp T1.nii.gz --seg aparc+aseg.mgz --regheader aparc+aseg.mgz --o aparc+aseg.nii.gz
    labelconvert aparc+aseg.nii.gz ${FREESURFER_HOME}/FreeSurferColorLUT.txt /opt/mrtrix3/share/mrtrix3/labelconvert/fs_default.txt nodes.nii.gz
    labelsgmfix nodes.nii.gz T1_brain.nii.gz /opt/mrtrix3/share/mrtrix3/labelconvert/fs_default.txt mask/nodes_sgmfixed.nii.gz -premasked
}

function split_nodes () {
    fslmaths mask/nodes_sgmfixed -thr  1 -uthr  1 mask/DK/01-lh-bankssts
    fslmaths mask/nodes_sgmfixed -thr  2 -uthr  2 mask/DK/02-lh-caudalanteriorcingulate
    fslmaths mask/nodes_sgmfixed -thr  3 -uthr  3 mask/DK/03-lh-caudalmiddlefrontal
    fslmaths mask/nodes_sgmfixed -thr  4 -uthr  4 mask/DK/04-lh-cuneus
    fslmaths mask/nodes_sgmfixed -thr  5 -uthr  5 mask/DK/05-lh-entorhinal
    fslmaths mask/nodes_sgmfixed -thr  6 -uthr  6 mask/DK/06-lh-fusiform
    fslmaths mask/nodes_sgmfixed -thr  7 -uthr  7 mask/DK/07-lh-inferiorparietal
    fslmaths mask/nodes_sgmfixed -thr  8 -uthr  8 mask/DK/08-lh-inferiortemporal
    fslmaths mask/nodes_sgmfixed -thr  9 -uthr  9 mask/DK/09-lh-isthmuscingulate
    fslmaths mask/nodes_sgmfixed -thr 10 -uthr 10 mask/DK/10-lh-lateraloccipital
    fslmaths mask/nodes_sgmfixed -thr 11 -uthr 11 mask/DK/11-lh-lateralorbitofrontal
    fslmaths mask/nodes_sgmfixed -thr 12 -uthr 12 mask/DK/12-lh-lingual
    fslmaths mask/nodes_sgmfixed -thr 13 -uthr 13 mask/DK/13-lh-medialorbitofrontal
    fslmaths mask/nodes_sgmfixed -thr 14 -uthr 14 mask/DK/14-lh-middletemporal
    fslmaths mask/nodes_sgmfixed -thr 15 -uthr 15 mask/DK/15-lh-parahippocampal
    fslmaths mask/nodes_sgmfixed -thr 16 -uthr 16 mask/DK/16-lh-paracentral
    fslmaths mask/nodes_sgmfixed -thr 17 -uthr 17 mask/DK/17-lh-parsopercularis
    fslmaths mask/nodes_sgmfixed -thr 18 -uthr 18 mask/DK/18-lh-parsorbitalis
    fslmaths mask/nodes_sgmfixed -thr 19 -uthr 19 mask/DK/19-lh-parstriangularis
    fslmaths mask/nodes_sgmfixed -thr 20 -uthr 20 mask/DK/20-lh-pericalcarine
    fslmaths mask/nodes_sgmfixed -thr 21 -uthr 21 mask/DK/21-lh-postcentral
    fslmaths mask/nodes_sgmfixed -thr 22 -uthr 22 mask/DK/22-lh-posteriorcingulate
    fslmaths mask/nodes_sgmfixed -thr 23 -uthr 23 mask/DK/23-lh-precentral
    fslmaths mask/nodes_sgmfixed -thr 24 -uthr 24 mask/DK/24-lh-precuneus
    fslmaths mask/nodes_sgmfixed -thr 25 -uthr 25 mask/DK/25-lh-rostralanteriorcingulate
    fslmaths mask/nodes_sgmfixed -thr 26 -uthr 26 mask/DK/26-lh-rostralmiddlefrontal
    fslmaths mask/nodes_sgmfixed -thr 27 -uthr 27 mask/DK/27-lh-superiorfrontal
    fslmaths mask/nodes_sgmfixed -thr 28 -uthr 28 mask/DK/28-lh-superiorparietal
    fslmaths mask/nodes_sgmfixed -thr 29 -uthr 29 mask/DK/29-lh-superiortemporal
    fslmaths mask/nodes_sgmfixed -thr 30 -uthr 30 mask/DK/30-lh-supramarginal
    fslmaths mask/nodes_sgmfixed -thr 31 -uthr 31 mask/DK/31-lh-frontalpole
    fslmaths mask/nodes_sgmfixed -thr 32 -uthr 32 mask/DK/32-lh-temporalpole
    fslmaths mask/nodes_sgmfixed -thr 33 -uthr 33 mask/DK/33-lh-transversetemporal
    fslmaths mask/nodes_sgmfixed -thr 34 -uthr 34 mask/DK/34-lh-insula
    fslmaths mask/nodes_sgmfixed -thr 35 -uthr 35 mask/DK/35-Left-Cerebellum-Cortex
    fslmaths mask/nodes_sgmfixed -thr 36 -uthr 36 mask/DK/36-Left-Thalamus-Proper
    fslmaths mask/nodes_sgmfixed -thr 37 -uthr 37 mask/DK/37-Left-Caudate
    fslmaths mask/nodes_sgmfixed -thr 38 -uthr 38 mask/DK/38-Left-Putamen
    fslmaths mask/nodes_sgmfixed -thr 39 -uthr 39 mask/DK/39-Left-Pallidum
    fslmaths mask/nodes_sgmfixed -thr 40 -uthr 40 mask/DK/40-Left-Hippocampus
    fslmaths mask/nodes_sgmfixed -thr 41 -uthr 41 mask/DK/41-Left-Amygdala
    fslmaths mask/nodes_sgmfixed -thr 42 -uthr 42 mask/DK/42-Left-Accumbens-area
    fslmaths mask/nodes_sgmfixed -thr 43 -uthr 43 mask/DK/43-Right-Thalamus-Proper
    fslmaths mask/nodes_sgmfixed -thr 44 -uthr 44 mask/DK/44-Right-Caudate
    fslmaths mask/nodes_sgmfixed -thr 45 -uthr 45 mask/DK/45-Right-Putamen
    fslmaths mask/nodes_sgmfixed -thr 46 -uthr 46 mask/DK/46-Right-Pallidum
    fslmaths mask/nodes_sgmfixed -thr 47 -uthr 47 mask/DK/47-Right-Hippocampus
    fslmaths mask/nodes_sgmfixed -thr 48 -uthr 48 mask/DK/48-Right-Amygdala
    fslmaths mask/nodes_sgmfixed -thr 49 -uthr 49 mask/DK/49-Right-Accumbens-area
    fslmaths mask/nodes_sgmfixed -thr 50 -uthr 50 mask/DK/50-rh-bankssts
    fslmaths mask/nodes_sgmfixed -thr 51 -uthr 51 mask/DK/51-rh-caudalanteriorcingulate
    fslmaths mask/nodes_sgmfixed -thr 52 -uthr 52 mask/DK/52-rh-caudalmiddlefrontal
    fslmaths mask/nodes_sgmfixed -thr 53 -uthr 53 mask/DK/53-rh-cuneus
    fslmaths mask/nodes_sgmfixed -thr 54 -uthr 54 mask/DK/54-rh-entorhinal
    fslmaths mask/nodes_sgmfixed -thr 55 -uthr 55 mask/DK/55-rh-fusiform
    fslmaths mask/nodes_sgmfixed -thr 56 -uthr 56 mask/DK/56-rh-inferiorparietal
    fslmaths mask/nodes_sgmfixed -thr 57 -uthr 57 mask/DK/57-rh-inferiortemporal
    fslmaths mask/nodes_sgmfixed -thr 58 -uthr 58 mask/DK/58-rh-isthmuscingulate
    fslmaths mask/nodes_sgmfixed -thr 59 -uthr 59 mask/DK/59-rh-lateraloccipital
    fslmaths mask/nodes_sgmfixed -thr 60 -uthr 60 mask/DK/60-rh-lateralorbitofrontal
    fslmaths mask/nodes_sgmfixed -thr 61 -uthr 61 mask/DK/61-rh-lingual
    fslmaths mask/nodes_sgmfixed -thr 62 -uthr 62 mask/DK/62-rh-medialorbitofrontal
    fslmaths mask/nodes_sgmfixed -thr 63 -uthr 63 mask/DK/63-rh-middletemporal
    fslmaths mask/nodes_sgmfixed -thr 64 -uthr 64 mask/DK/64-rh-parahippocampal
    fslmaths mask/nodes_sgmfixed -thr 65 -uthr 65 mask/DK/65-rh-paracentral
    fslmaths mask/nodes_sgmfixed -thr 66 -uthr 66 mask/DK/66-rh-parsopercularis
    fslmaths mask/nodes_sgmfixed -thr 67 -uthr 67 mask/DK/67-rh-parsorbitalis
    fslmaths mask/nodes_sgmfixed -thr 68 -uthr 68 mask/DK/68-rh-parstriangularis
    fslmaths mask/nodes_sgmfixed -thr 69 -uthr 69 mask/DK/69-rh-pericalcarine
    fslmaths mask/nodes_sgmfixed -thr 70 -uthr 70 mask/DK/70-rh-postcentral
    fslmaths mask/nodes_sgmfixed -thr 71 -uthr 71 mask/DK/71-rh-posteriorcingulate
    fslmaths mask/nodes_sgmfixed -thr 72 -uthr 72 mask/DK/72-rh-precentral
    fslmaths mask/nodes_sgmfixed -thr 73 -uthr 73 mask/DK/73-rh-precuneus
    fslmaths mask/nodes_sgmfixed -thr 74 -uthr 74 mask/DK/74-rh-rostralanteriorcingulate
    fslmaths mask/nodes_sgmfixed -thr 75 -uthr 75 mask/DK/75-rh-rostralmiddlefrontal
    fslmaths mask/nodes_sgmfixed -thr 76 -uthr 76 mask/DK/76-rh-superiorfrontal
    fslmaths mask/nodes_sgmfixed -thr 77 -uthr 77 mask/DK/77-rh-superiorparietal
    fslmaths mask/nodes_sgmfixed -thr 78 -uthr 78 mask/DK/78-rh-superiortemporal
    fslmaths mask/nodes_sgmfixed -thr 79 -uthr 79 mask/DK/79-rh-supramarginal
    fslmaths mask/nodes_sgmfixed -thr 80 -uthr 80 mask/DK/80-rh-frontalpole
    fslmaths mask/nodes_sgmfixed -thr 81 -uthr 81 mask/DK/81-rh-temporalpole
    fslmaths mask/nodes_sgmfixed -thr 82 -uthr 82 mask/DK/82-rh-transversetemporal
    fslmaths mask/nodes_sgmfixed -thr 83 -uthr 83 mask/DK/83-rh-insula
    fslmaths mask/nodes_sgmfixed -thr 84 -uthr 84 mask/DK/84-Right-Cerebellum-Cortex
}

function conbine_region () {

    #
    # Whole brain
    #

    fsladd mask/wholebrain $(imglob mask/DK/*)

    #
    # Cortical
    #

    # Frontal
    fsladd mask/cortical-frontal \
        mask/DK/27-lh-superiorfrontal \
        mask/DK/76-rh-superiorfrontal \
        mask/DK/26-lh-rostralmiddlefrontal \
        mask/DK/75-rh-rostralmiddlefrontal \
        mask/DK/03-lh-caudalmiddlefrontal \
        mask/DK/52-rh-caudalmiddlefrontal \
        mask/DK/17-lh-parsopercularis \
        mask/DK/66-rh-parsopercularis \
        mask/DK/18-lh-parsorbitalis \
        mask/DK/67-rh-parsorbitalis \
        mask/DK/19-lh-parstriangularis \
        mask/DK/68-rh-parstriangularis \
        mask/DK/11-lh-lateralorbitofrontal \
        mask/DK/13-lh-medialorbitofrontal \
        mask/DK/60-rh-lateralorbitofrontal \
        mask/DK/62-rh-medialorbitofrontal \
        mask/DK/23-lh-precentral \
        mask/DK/72-rh-precentral \
        mask/DK/16-lh-paracentral \
        mask/DK/65-rh-paracentral \
        mask/DK/31-lh-frontalpole \
        mask/DK/80-rh-frontalpole

    # Parietal
    fsladd mask/cortical-parietal \
        mask/DK/28-lh-superiorparietal \
        mask/DK/77-rh-superiorparietal \
        mask/DK/07-lh-inferiorparietal \
        mask/DK/56-rh-inferiorparietal \
        mask/DK/30-lh-supramarginal \
        mask/DK/79-rh-supramarginal \
        mask/DK/21-lh-postcentral \
        mask/DK/70-rh-postcentral \
        mask/DK/24-lh-precuneus \
        mask/DK/73-rh-precuneus

    # Temporal
    fsladd mask/cortical-temporal \
        mask/DK/29-lh-superiortemporal \
        mask/DK/78-rh-superiortemporal \
        mask/DK/14-lh-middletemporal \
        mask/DK/63-rh-middletemporal \
        mask/DK/08-lh-inferiortemporal \
        mask/DK/57-rh-inferiortemporal \
        mask/DK/01-lh-bankssts \
        mask/DK/50-rh-bankssts \
        mask/DK/06-lh-fusiform \
        mask/DK/55-rh-fusiform \
        mask/DK/33-lh-transversetemporal \
        mask/DK/82-rh-transversetemporal \
        mask/DK/05-lh-entorhinal \
        mask/DK/54-rh-entorhinal \
        mask/DK/32-lh-temporalpole \
        mask/DK/81-rh-temporalpole \
        mask/DK/15-lh-parahippocampal \
        mask/DK/64-rh-parahippocampal

    # Occipital
    fsladd mask/cortical-occipital \
        mask/DK/10-lh-lateraloccipital \
        mask/DK/59-rh-lateraloccipital \
        mask/DK/12-lh-lingual \
        mask/DK/61-rh-lingual \
        mask/DK/04-lh-cuneus \
        mask/DK/53-rh-cuneus \
        mask/DK/20-lh-pericalcarine \
        mask/DK/69-rh-pericalcarine

    # Cingulate
    fsladd mask/cortical-cingulate \
        mask/DK/25-lh-rostralanteriorcingulate \
        mask/DK/74-rh-rostralanteriorcingulate \
        mask/DK/02-lh-caudalanteriorcingulate \
        mask/DK/51-rh-caudalanteriorcingulate \
        mask/DK/22-lh-posteriorcingulate \
        mask/DK/71-rh-posteriorcingulate \
        mask/DK/09-lh-isthmuscingulate \
        mask/DK/58-rh-isthmuscingulate

    # Cortical
    fsladd mask/cortical \
        mask/cortical-frontal \
        mask/cortical-parietal \
        mask/cortical-temporal \
        mask/cortical-occipital \
        mask/cortical-cingulate

    #
    # Subcortical
    #

    fsladd mask/subcortical \
        mask/DK/36-Left-Thalamus-Proper \
        mask/DK/37-Left-Caudate \
        mask/DK/38-Left-Putamen \
        mask/DK/39-Left-Pallidum \
        mask/DK/42-Left-Accumbens-area \
        mask/DK/43-Right-Thalamus-Proper \
        mask/DK/44-Right-Caudate \
        mask/DK/45-Right-Putamen \
        mask/DK/46-Right-Pallidum \
        mask/DK/49-Right-Accumbens-area

    #
    # Gray matter
    #

    fsladd mask/graymatter mask/cortical mask/subcortical

    #
    # White matter
    #

    mri_label2vol --temp T1.nii.gz --seg wm.seg.mgz --regheader wm.seg.mgz --o wm.seg.nii.gz
    fslmaths wm.seg -bin mask/whitematter
}

function reg_maps_to_T1 () {
    bet b0 b0_brain -f 0.1 -m
    epi_reg --epi=b0_brain --t1=T1 --t1brain=T1_brain --out=BBR_DWI-to-T1
    for map in $(imglob map/* | sed "s|map/||g" | tr " " "\n"); do
        flirt -in map/${map} -ref T1 -out ${map}_T1 -init BBR_DWI-to-T1.mat -applyxfm
    done
}

function out_results () {
    [[ -d results ]] && rm -fr results
    mkdir results

    for map in $(imglob map/* | sed "s|map/||g" | tr " " "\n"); do
        for mask in $(imglob mask/DK/*); do
            fslstats ${map}_T1 -k ${mask} -M >> results/DK_${map}_mean.txt &
            fslstats ${map}_T1 -k ${mask} -S >> results/DK_${map}_sd.txt &
            wait
        done

        for mask in wholebrain graymatter whitematter cortical subcortical \
                    cortical-frontal cortical-parietal cortical-temporal cortical-occipital cortical-cingulate; do
            fslstats ${map}_T1 -k mask/${mask} -M >> results/DK_${map}_mean.txt &
            fslstats ${map}_T1 -k mask/${mask} -S >> results/DK_${map}_sd.txt &
            wait
        done
    done

}

create_nodes
split_nodes
conbine_region
reg_maps_to_T1
out_results

4. FreeSurferのwm.seg.nii.gzをDWI空間に位置合わせする方法

T1WI空間にあるFreeSurferのwm.seg.nii.gzををDWI空間に位置合わせするスクリプト。

4.1. 必要なファイル

  • dwi.nii.gz (DWI 画像)
  • orig.mgz (T1 画像; T1空間)
  • wm.seg.mgz (White matter ラベル; T1空間)

4.2. ソースコード

# mgz を nii.gz に変換
mrconvert orig.mgz orig.nii.gz
mrconvert wm.seg.mgz wm.seg.nii.gz
fslmaths wm.seg.nii.gz -bin wm.seg.nii.gz
# reorient
fslreorient2std dwi.nii.gz dwi.nii.gz
fslreorient2std orig.nii.gz orig.nii.gz
fslreorient2std wm.seg.nii.gz wm.seg.nii.gz
# epi_reg で DWI→T1 の変換行列を計算し、逆行列 (T1→DWI) も計算
bet orig.nii.gz orig_brain.nii.gz -f 0.1 -B -m
fslroi dwi.nii.gz b0.nii.gz 0 1
epi_reg --epi=b0.nii.gz --t1=orig.nii.gz --t1brain=orig_brain.nii.gz --out=DWI_to_T1
convert_xfm -omat T1_to_DWI.mat -inverse DWI_to_T1.mat
# White matter ラベルを DWI 空間に移動
flirt -in wm.seg.nii.gz -ref b0.nii.gz -out wm.nii.gz -init T1_to_DWI.mat -applyxfm -interp nearestneighbour
# 確認
freeview dwi.nii.gz wm.nii.gz:colormap=jet

5. 標準空間にあるROIを各被験者脳に位置合わせ

# required files
## T1_brain.nii.gz                  : BET T1WI
## T1.nii.gz                        : T1WI
## nodif_brain.nii.gz               : b0 image
## MNI152_T1_1mm_brain.nii.gz       : template image
## SN_labels_LR.nii.gz              : subsegmented mask

# Register SN-labels to Subject space
## MNI to T1WI
flirt -in T1_brain -ref ../MNI152_T1_1mm_brain -omat SubT1_to_MNI.mat
# fnirt --ref=../MNI152_T1_1mm_brain --in=T1_brain --aff=SubT1_to_MNI.mat --config=T1_2_MNI152_2mm --cout=warp_SubT1_to_MNI
fnirt --ref=../MNI152_T1_1mm_brain --in=T1_brain --aff=SubT1_to_MNI.mat --cout=warp_SubT1_to_MNI
invwarp -w warp_SubT1_to_MNI -o warp_MNI_to_SubT1 -r T1_brain

## T1WI to DWI
epi_reg --epi=nodif_brain --t1=T1 --t1brain=T1_brain --out=DWI_to_T1
convert_xfm -omat T1_to_DWI.mat -inverse DWI_to_T1.mat

# warp MNI to Subject T1WI space
# applywarp --in=../SN_labels_L --ref=T1_brain --warp=warp_MNI_to_SubT1 --interp=nn --out=SN_labels_L_to_SubT1
# applywarp --in=../SN_labels_R --ref=T1_brain --warp=warp_MNI_to_SubT1 --interp=nn --out=SN_labels_R_to_SubT1
applywarp --in=../SN_labels_LR --ref=T1_brain --warp=warp_MNI_to_SubT1 --interp=nn --out=SN_labels_LR_to_SubT1
# warp Subject T1WI space to Subject DWI space
# flirt -in SN_labels_L_to_SubT1 -ref nodif_brain -out SN_labels_L_to_SubDWI -init T1_to_DWI.mat -applyxfm -interp nearestneighbour
# flirt -in SN_labels_R_to_SubT1 -ref nodif_brain -out SN_labels_R_to_SubDWI -init T1_to_DWI.mat -applyxfm -interp nearestneighbour
flirt -in SN_labels_LR_to_SubT1 -ref nodif_brain -out SN_labels_LR_to_SubDWI -init T1_to_DWI.mat -applyxfm -interp nearestneighbour

# calc map values
[[ -d results ]] && \rm -rf results
mkdir -p results
MASK='SN_labels_LR_to_SubDWI'
for MAP in $(\ls -1 map/*.nii.gz | xargs -i basename -s ".nii.gz" {}); do
    fslstats -K ${MASK} map/${MAP}  -M >> results/SN-labels_${MAP}_mean.txt &
    # fslstats -K ${MASK} map/${MAP}  -S >> results/SN-labels_${MAP}_sd.txt &
    wait
done

6. bertのwmparcをMNI空間に移動後、個人脳(Perfusion image)にレジストレーション

# Required files
# - brain.mgz: bert's T1WI
# - wmparc.mgz: bert's wmparc
# - <perfusion>.nii.gz: a subjects's perfusion images

# Definition
PERFUSION='E_01_TR05_sms_perfusion_Junten_176_20190212142306_121'
LABEL='std_subj_wmparc'

# 1. Preparation
echo '1. Preparation'
## compress data
find . -name "*nii" |xargs -i gzip {}
## copy files
mkdir preprocessing
cp brain.mgz wmparc.mgz ${PERFUSION}.nii.gz preprocessing
# cp brain.mgz wmparc.mgz ${PERFUSION}.nii.gz std_bert_wmparc_MNI std_bert_brain_MNI preprocessing
cd preprocessing
## convert bert's mgz to nii.gz
mrconvert brain.mgz bert_brain.nii.gz
mrconvert wmparc.mgz bert_wmparc.nii.gz
## abstract first volume from perfusion images
fslroi $PERFUSION perfusion 0 1
## reorient to FMRIB'S standard
fslreorient2std bert_brain std_bert_brain
fslreorient2std bert_wmparc std_bert_wmparc
fslreorient2std perfusion std_perfusion
## BET perfusion image
bet std_perfusion std_perfusion_brain

# 2. bert2MNI
echo '2. bert2MNI'
## flirt
flirt -in std_bert_brain \
      -ref ${FSLDIR}/data/standard/MNI152_T1_1mm_brain \
      -omat bert_to_MNI.mat
## fnirt
fnirt --ref=${FSLDIR}/data/standard/MNI152_T1_1mm_brain \
      --in=std_bert_brain --aff=bert_to_MNI.mat \
      --cout=warp_bert_to_MNI --iout=std_bert_brain_MNI
applywarp --in=std_bert_wmparc \
          --ref=${FSLDIR}/data/standard/MNI152_T1_1mm_brain \
          --warp=warp_bert_to_MNI --interp=nn \
          --out=std_bert_wmparc_MNI


## 3. regster ROIs to subject space (bert2subj)
echo '3. regster ROIs to subject space (bert2subj)'
flirt -in std_perfusion_brain \
      -ref std_bert_brain_MNI \
      -omat subj_to_bert.mat \
      -out std_perfusion_brain_flirt
convert_xfm -omat bert_to_subj.mat  -inverse subj_to_bert.mat
flirt -in std_bert_wmparc_MNI \
      -ref std_perfusion_brain \
      -out std_subj_wmparc \
      -init bert_to_subj.mat \
      -applyxfm \
      -interp nearestneighbour


## 4. abstract Tissues
echo '4. abstract Tissues'
mkdir ROI
## GM
fslmaths $LABEL -thr 800 -uthr 1200  -bin ROI/l_gm
fslmaths $LABEL -thr 1800 -uthr 2200 -bin ROI/r_gm
## WM
fslmaths $LABEL -thr 2800 -uthr 3200 -bin ROI/l_wm
# fslmaths $LABEL -thr 5001 -uthr 5001 -bin l_wm_tmp2
# fsladd l_wm l_wm_tmp*
fslmaths $LABEL -thr 3800 -uthr 4200 -bin ROI/r_wm
# fslmaths $LABEL -thr 5002 -uthr 5002 -bin r_wm_tmp2
# fsladd r_wm r_wm_tmp*
# rm *_wm_tmp*
## Cerebellum
fslmaths $LABEL -thr 7 -uthr 8 -bin ROI/l_cblm
fslmaths $LABEL -thr 46 -uthr 47 -bin ROI/r_cblm
mv ROI ../
cd ..


## 5. check results
fsleyes $PERFUSION \
        ROI/l_gm ROI/r_gm \
        ROI/l_wm ROI/r_wm \
        ROI/l_cblm ROI/r_cblm &

bert脳とwmparcをMNIに移動した画像が用意できていれば、次のコードを実行するとよい。

# Required files
# - brain.mgz: bert's T1WI
# - wmparc.mgz: bert's wmparc
# - <perfusion>.nii.gz: a subjects's perfusion images
# - std_bert_brain_MNI.nii.gz: bert's brain T1WI in MNI
# - std_bert_wmparc_MNI.nii.gz: bert's wmparc in MNI

# Definition
PERFUSION='org_perfusion'
LABEL='std_subj_wmparc'
TEMPLATE='std_bert_brain_MNI'
TEMPLAET_LABEL='std_bert_wmparc_MNI'

# 1. Preparation
echo '1. Preparation'
## compress data
find . -name "*nii" |xargs -i gzip {}
## copy files
mkdir preprocessing
cp ${PERFUSION}.nii.gz ${TEMPLAET_LABEL}.nii.gz ${TEMPLATE}.nii.gz preprocessing
cd preprocessing
## abstract first volume from perfusion images
mcflirt -in $PERFUSION -out ${PERFUSION}_realign
fslmaths ${PERFUSION}_realign -Tmean perfusion
## reorient to FMRIB'S standard
fslreorient2std perfusion std_perfusion
## BET perfusion image
bet std_perfusion std_perfusion_brain

## 2. regster ROIs to subject space (bert2subj)
echo '2. regster ROIs to subject space (bert2subj)'
## subj2bert
# flirt -in std_perfusion_brain \
#       -ref std_bert_brain_MNI \
#       -omat subj_to_bert.mat \
#       -out std_perfusion_brain_flirt
# convert_xfm -omat bert_to_subj.mat  -inverse subj_to_bert.mat

## check std_perfusion_brain_flirt.nii.gz
## if fail subj2bert, use the follwing code:
## bert2subj
flirt -in  std_bert_brain_MNI \
      -ref std_perfusion_brain \
      -omat bert_to_subj.mat \
      -out std_bert_brain_subj

## apply warp to register bert's to subj's wmparc
flirt -in std_bert_wmparc_MNI \
      -ref std_perfusion_brain \
      -out std_subj_wmparc \
      -init bert_to_subj.mat \
      -applyxfm \
      -interp nearestneighbour
# fsleyes std_perfusion_brain.nii.gz std_bert_brain_subj.nii.gz

## 3. abstract Tissues
echo '3. abstract Tissues'
mkdir ROI
## GM
fslmaths $LABEL -thr 800 -uthr 1200 -bin -mul 1 ROI/l_gm
fslmaths $LABEL -thr 1800 -uthr 2200 -bin -mul 2 ROI/r_gm
## WM
fslmaths $LABEL -thr 2800 -uthr 3200 -bin -mul 3 ROI/l_wm
# fslmaths $LABEL -thr 5001 -uthr 5001 -bin l_wm_tmp2
# fsladd l_wm l_wm_tmp*
fslmaths $LABEL -thr 3800 -uthr 4200 -bin -mul 4 ROI/r_wm
# fslmaths $LABEL -thr 5002 -uthr 5002 -bin r_wm_tmp2
# fsladd r_wm r_wm_tmp*
# rm *_wm_tmp*
## Cerebellum
fslmaths $LABEL -thr 7 -uthr 8 -bin -mul 5 ROI/l_cblm
fslmaths $LABEL -thr 46 -uthr 47 -bin -mul 6 ROI/r_cblm
fsladd ROI/subj_WM_GM_CBLM_labels ROI/l_gm ROI/r_gm ROI/l_wm ROI/r_wm ROI/l_cblm ROI/r_cblm
mv ROI ../
cd ..


## 5. check results
# fsleyes $PERFUSION \
#         ROI/l_gm ROI/r_gm \
#         ROI/l_wm ROI/r_wm \
#         ROI/l_cblm ROI/r_cblm &

# for k in *;do
# fsleyes $k/preprocessing/std_perfusion_brain \
#         $k/ROI/l_gm $k/ROI/r_gm \
#         $k/ROI/l_wm $k/ROI/r_wm \
#         $k/ROI/l_cblm $k/ROI/r_cblm 
# done

【FSL】Gray matter-Based Spatial Statistics: GBSS



1. 目的
2. GBSSとは
2.1. データの準備
2.2. 解析パラメータの設定
2.3. gbss_1_T1wpreproc:T1WIの前処理
2.4. gbss_2_DWIwpreproc:拡散定量値の前処理
2.5. gbss_3_skelpreproc:拡散定量値をスケルトンに投影
2.6. randomise:スケルトンに投影された定量値画像を入力したGLMと並べ替え検定(permutation test)
3. おまけ


1. 目的

  • Gray matter-Based Spatial Statistics: GBSS

2. GBSSとは

Gray matter-Based Spatial Statistics(GBSS)は、灰白質の統計解析をするための手法。TBSSの灰白質版。

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

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

  1. データの準備
  2. 解析パラメータの設定
  3. gbss_1_T1wpreproc:T1WIの前処理
  4. gbss_2_DWIwpreproc:拡散定量値の前処理
  5. gbss_3_skelpreproc:拡散定量値をスケルトンに投影
  6. randomise:スケルトンに投影された定量値画像を入力したGLMと並べ替え検定(permutation test)

2.1. データの準備

まず、次のデータを準備する。

  • b0フォルダ:b0 (DWI, b=0 s/mm^2)
  • T1wフォルダ:3D-T1WI
  • 拡散定量値:FA, MD, FW, etc.

フォルダ構造は次のようにする。

ここでは、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)いることを想定する。各フォルダに入れるファイル名は同じにする必要がある(例:FA/Subj001.nii.gz, T1w/Subj001.nii.gzは同じファイル名)。

.
├── FA
│   ├── Con0001.nii.gz
│   ├── Con0002.nii.gz
│   ├── ...
│   └── Pat0010.nii.gz
├── FW
│   ├── Con0001.nii.gz
│   ├── Con0002.nii.gz
│   ├── ...
│   └── Pat0010.nii.gz
├── MD
│   ├── Con0001.nii.gz
│   ├── Con0002.nii.gz
│   ├── ...
│   └── Pat0010.nii.gz
├── T1w
│   ├── Con0001.nii.gz
│   ├── Con0002.nii.gz
│   ├── ...
│   └── Pat0010.nii.gz
└── b0
    ├── Con0001.nii.gz
    ├── Con0002.nii.gz
    ├── ...
    └── Pat0010.nii.gz

2.2. 解析パラメータの設定

解析に用いるパラメータを設定する。

PEDIRECHOSPACINGepi_regコマンドのパラメータであり、それぞれ位相エンコード方向とecho spacing timeを指定する。SKELETON_THRはスケルトン化する際のしきい値である。

MAP_LIST=$(ls | grep -v b0 | grep -v T1w | grep -v stats)
SUBJ_LIST=$(ls T1w | cut -d . -f1)
PEDIR='-y'                        # Phase encoding direction: "epi_reg" config
ECHOSPACING='0.0380544'           # Echo spacing time (s): "epi_reg" config
SKELETON_THR=0.2                  # Skeleton threshold

2.3. gbss_1_T1wpreproc:T1WIの前処理

T1WIの前処理内容は、次の通り。

  1. 脳頭蓋除去
  2. 灰白質をセグメント
  3. 灰白質を標準脳(MNI152)に位置合わせ
  4. 標準空間上の各被験者の灰白質を平均化

以上の処理内容を実行するために、関数gbss_1_T1wpreprocを定義。

function gbss_1_T1wpreproc() {
    echo "T1w preprocessing..."
    # T1w preprocessing
    for SUBJ in ${SUBJ_LIST}; do
        ## Segment GM
        ### Skull-stripping
        bet T1w/${SUBJ} T1w/${SUBJ}_skull_stripped -f 0.3 -R -S -B
        ### Segment T1w into CSF/GM/WM
        fast -t 1 -g -B -b -p -o T1w/${SUBJ}_fast T1w/${SUBJ}_skull_stripped

        ## Register GM to MNI
        ### Register T1w to MNI
        flirt -ref ${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
            -in T1w/${SUBJ}_fast_restore \
            -omat T1w/${SUBJ}_indiv2MNI.mat
        fnirt --in=T1w/${SUBJ}_fast_restore \
            --aff=T1w/${SUBJ}_indiv2MNI.mat \
            --cout=T1w/${SUBJ}_indiv2MNI_warp \
            --config=T1_2_MNI152_2mm
        ### applywarp to GM and move it into MNI
        applywarp --ref=${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
            --in=T1w/${SUBJ}_fast_pve_1 \
            --warp=T1w/${SUBJ}_indiv2MNI_warp \
            --out=T1w/${SUBJ}_fast_pve_1_inMNI &
    done
    ## Create 4D all_GM image
    ### GM_pve1 in MNI are merged into one volume
    mkdir stats
    fslmerge -t stats/all_GM $(imglob T1w/*_fast_pve_1_inMNI.nii.gz)
    ### Create merged GM_pve1 (all_GM) mask and apply it to all_GM
    fslmaths stats/all_GM -max 0 -Tmin -bin stats/mean_GM_mask -odt char
    fslmaths stats/all_GM -Tmean stats/mean_GM
}

関数が定義できたら、実行する。

gbss_1_T1wpreproc

処理が完了すると、statsフォルダにすべての被験者の灰白質画像がまとめられた「all_GM.nii.gz」とそれを平均化した「mean_GM.nii.gz」が生成される。

2.4. gbss_2_DWIwpreproc:拡散定量値の前処理

拡散定量値の前処理の過程は、次の通り。

  1. b0をT1WIに位置合わせ
  2. b0 to T1WIとT1WI to 標準脳のWarp Fieldを使って、拡散定量値マップを標準空間に移動
  3. 標準空間上にあるすべての被験者の拡散定量値マップをひとつの4D volume dataにまとめる
  4. 平均灰白質マスクmean_GM_maskで、3の拡散定量値マップをマスク処理
function gbss_2_DWIwpreproc() {
    echo "Diffusion MAP preprocessing..."
    # Diffusion MAP preprocessing
    for SUBJ in ${SUBJ_LIST}; do
        ## Register b0 to MNI
        ### Register b0 to T1WI
        epi_reg --epi=b0/${SUBJ} \
            --t1=T1w/${SUBJ} --t1brain=T1w/${SUBJ}_skull_stripped \
            --echospacing=${ECHOSPACING} --pedir=${PEDIR} \
            --wmseg=T1w//${SUBJ}_fast_seg_2.nii.gz \
            --out=b0/${SUBJ}_BBR
        for MAP in ${MAP_LIST}; do
            ### applywarp to diffusion map and move it to T1w based on epi_reg
            flirt -in ${MAP}/${SUBJ} \
                -out ${MAP}/${SUBJ}_inT1 \
                -ref T1w/${SUBJ}_skull_stripped \
                -init b0/${SUBJ}_BBR.mat -applyxfm
            ### applywarp to dMAP in T1w and move it to MNI
            applywarp --ref=${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
                --in=${MAP}/${SUBJ}_inT1 \
                --warp=T1w/${SUBJ}_indiv2MNI_warp \
                --out=${MAP}/${SUBJ}_inMNI
        done
    done
    ## 4D all_${dMAP}
    for MAP in ${MAP_LIST}; do
        ### Merge dMAP in MNI into one volume, all_${MAP}
        fslmerge -t stats/all_${MAP} $(imglob ${MAP}/*_inMNI.nii.gz)
        ### Masking dMAP in MNI using mean GM mask
        fslmaths stats/all_${MAP} -mas stats/mean_GM_mask stats/all_${MAP}
    done
}

関数が定義できたら、実行する。

gbss_2_DWIwpreproc

処理が完了すると、拡散定量値の種類に応じてstatsフォルダに「all_.nii.gz」が生成される。

以下は、標準空間上における灰白質領域のFA画像である。

2.5. gbss_3_skelpreproc:拡散定量値をスケルトンに投影

拡散定量値をスケルトンに投影するまでの手順は、次の通り。

  1. 平均灰白質画像からスケルトンを生成
  2. 灰白質スケルトンをしきい値処理をしてスケルトンマスクを生成
  3. スケルトンの距離マップ(distance map)を生成
  4. 各拡散定量値マップを、スケルトンに投影
function gbss_3_skelpreproc() {
    echo "Projecting MAP into skeleton..."
    # Skeletonize dMAP
    ## Create GM skeleton
    tbss_skeleton -i stats/mean_GM -o stats/mean_GM_skeleton
    ## threshold GM skeleton
    fslmaths stats/mean_GM_skeleton -thr ${SKELETON_THR} -bin stats/mean_GM_skeleton_mask
    ## Create skeleton distancemap
    fslmaths stats/mean_GM_mask -mul -1 -add 1 -add stats/mean_GM_skeleton_mask stats/mean_GM_skeleton_mask_dst
    distancemap -i stats/mean_GM_skeleton_mask_dst -o stats/mean_GM_skeleton_mask_dst
    for MAP in ${MAP_LIST}; do
        ## Project masked dMAP in MNI into skeleton
        tbss_skeleton -i stats/mean_GM \
            -p ${SKELETON_THR} \
            stats/mean_GM_skeleton_mask_dst \
            ${FSLDIR}/data/standard/LowerCingulum_1mm \
            stats/all_${MAP} \
            stats/all_${MAP}_skeletonised
    done
}

処理が完了すると、statsフォルダに灰白質のスケルトンmean_GM_skeleton.nii.gz、灰白質のスケルトンの距離マップmean_GM_skeleton_mask_dst.nii.gz、拡散定量値が投影されたスケルトンall_<MAP>_skeletonised.nii.gzが保存される。

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

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

今回は、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)のデータがある。これらのファイルを昇順に並び替えると、先に健常者10名の画像、次に患者10名の画像の順に並ぶ。

Con0001.nii.gz
Con0002.nii.gz
Con0003.nii.gz
...
Pat0010.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:統計値マップ
for MAP in ${MAP_LIST}; do
    randomise -i stats/all_${MAP}_skeletonised \
        -o stats/gbss_${MAP} \
        -m stats/mean_GM_mask \
        -d stats/design.mat \
        -t stats/design.con \
        -n 10000 --T2 -V -x --uncorrp -R
done

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

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

fsleyes ${FSLDIR}/data/standard/MNI152_T1_2mm_brain \
    stats/mean_GM_skeleton -cm Green \
    stats/gbss_FA_tfce_p_tstat1 -cm Red-Yellow -dr 0.95 1

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

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

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

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

tbss_fill stats/gbss_FA_tfce_p_tstat1 0.95 stats/mean_GM stats/gbss_FA_tfce_p_tstat1_fill

赤く表示されている領域は、健常群が患者群よりも有意(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_GM 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
gbss_FA_tfce_corrp_tstat1.nii.gz        0.992000
gbss_FA_tfce_corrp_tstat2.nii.gz        0.416000
gbss_FW_tfce_corrp_tstat1.nii.gz        0.361839
gbss_FW_tfce_corrp_tstat2.nii.gz        0.997261
gbss_MD_tfce_corrp_tstat1.nii.gz    0.389816
gbss_MD_tfce_corrp_tstat2.nii.gz    0.985748


gbss_FW_tfce_corrp_tstat2.nii.gz        0.997261
gbss_FA_tfce_corrp_tstat1.nii.gz        0.992000
gbss_MD_tfce_corrp_tstat2.nii.gz    0.985748
gbss_FA_tfce_corrp_tstat2.nii.gz        0.416000
gbss_MD_tfce_corrp_tstat1.nii.gz    0.389816
tbss_FW_tfce_corrp_tstat1.nii.gz        0.361839

【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.nii.gz)が準備できたら、頭蓋除去済みの脳画像を、灰白質・白質・脳脊髄液にセグメントする。次に、全ての被験者の灰白質(*_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

【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の画像は、以下。

【FSL】FSLを用いた拡散MRIの前処理(歪み補正&頭の動き補正)


1. 目的
2. 脳画像解析における拡散MRIの問題点
2.1. 磁化率効果に起因する画像の幾何学的歪み(位相エンコードをAP方向として撮像した場合)
2.2. 撮像中の頭の動き
3. 使用するソフトウェア
3.1. TOPUP
3.1.1. TOPUP
3.2. –indexについて
3.2.1. TOPUP
3.3. EDDY
3.4. 頭の動き補正の効果
4. サンプルコード


1. 目的

  • 拡散MRIの前処理(歪み補正&頭の動き補正)

2. 脳画像解析における拡散MRIの問題点

拡散MRI (Diffusion MRI) はEcho Planer Imaging (EPI)というシーケンスで撮像されているが、EPIは磁化率効果による影響に敏感であり、画像が幾何学的に歪むという問題がある。

拡散MRIでは、水分子拡散状態を可視化する技術である。水分子の動きをとらえるためには、Motion Probing Gradient (MPG)という水分子の動きを検出するための傾斜磁場を何種類(何軸)も用いる。その結果、脳の拡散MRIでは脳画像というVolumeデータが、MPGの数に応じて複数Volume撮像されることになるが、撮像中に頭が動くことによってVolume間の位置が合わないという問題がある。特に脳画像解析では、Voxel単位という小さなマスの中で解析をするので、撮像中の頭の動きで生じる、Volume間の位置ずれは解析結果に影響する。

そこで、拡散MRIで得られた脳画像に対して、「画像の歪み補正」と「頭の動き補正」をする。

2.1. 磁化率効果に起因する画像の幾何学的歪み(位相エンコードをAP方向として撮像した場合)

2.2. 撮像中の頭の動き

3. 使用するソフトウェア

拡散MRIの前処理にはFSLの、TOPUPおよびEDDY

3.1. TOPUP

TOPUPは、静磁場の不均一(non-zero off-resonance fields)による画像の幾何学的歪みを補正する。

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

topup --imain=<Input image> --datain=<Acquisition_parameters file> --config=b02b0.cnf --out=<Output prefix>

各引数で必要なファイルは、次の通り。

  • –imain: 位相エンコード (PE)が双方向で対になったペア画像 (PEがAP方向の画像 (ex. b=0) とPA方向の画像 (ex. b=0) )。通常b値=0の画像を用いる。
  • –datain: 拡散MRI撮像パラメータ (PE方向およびTotal readout time)が記載されたテキストファイル
  • –config: TOPUPの設定ファイル
  • –out: 出力ファイルの接頭辞

TOPUPによる補正を拡散MRI画像に含まれるすべてのVolumeデータに適応したい場合、applytopupを用いる。

applytopup --imain=<Input image> --datain=<Acquisition_parameters file> --inindex=<Index file> --topup=<TOPUP output prefix> --method=<Resampling method> --out=<Output prefix>
  • –imain: 位相エンコード (PE)が双方向で対になったペア画像 (PEがAP方向の画像 (ex. b=0) とPA方向の画像 (ex. b=0) )。通常b値=0の画像を用いる。
  • –datain: 拡散MRI撮像パラメータ (PE方向およびTotal readout time)が記載されたテキストファイル
  • –inindex: –imainで入力した拡散MRI画像の撮像パラメータが–acqpで指定されたAcquisition_parametersの内どちらのものか、というインデックス
  • –method: 画像のリサンプリン方法
  • –topup: TOPUPの出力ファイルの接頭辞
  • –out: applytopupの出力ファイルの接頭辞
3.1.1. –datain(Acquisition_parameters)について

Acquisition_parametersは、拡散MRI撮像で用いるEPIシーケンスのパラメータである。

必要な情報は次の2つである。

  • 位相エンコード方向 (x, y, z)
  • Total readout time (s)

位相エンコード方向は、x, y, zで表現する。APで撮像した場合は「0 -1 0」、PAで撮像した場合は「0 1 0」となる。

Total readout timeは、次の式で算出できる。

例えば、–imainの入力画像がPEがAPの画像, PAの画像の順番で1 volumeずつマージされている4D画像であり、Echo spacingが0.040であるとすると、Acquisition_parametersは次のように表現する。これをテキストファイル (.txt) として保存して–datainに入力する。

0 -1 0 0.046
0  1 0 0.046

3.2. –indexについて

EDDYでは、拡散MRI画像 (4D画像)のすべてのVolumeに対して補正を実行する。そのため、各Volumeがどのような撮像パラメータ(PE方向およびTotal readout time)で撮像されたものかを入力する必要がある。そこで、–acqpで指定されたAcquisition_parametersをインデックスで用いて参照する。これがインデックスをファイルを入力する理由である。

例えば、64 Volumesの拡散MRI画像があり、1~32 VolumesはPEがAP方向で撮像した画像データ、33~64 VolumesはPEがPA方向で撮像した画像データである場合を考える。

Acquisition_parametersは次のようであるとする。1行目はAP方向の情報、2行目はPA方向の情報である。

0 -1 0 0.046
0  1 0 0.046

この場合のindexファイルは次のようになる。これをテキストファイル (.txt) として保存して–indexに入力する。

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 
3.2.1. 画像の歪み補正の効果

3.3. EDDY

EDDYは、渦電流 (Eddy current) に起因する動的な磁場の不均一や頭の動きを補正する。

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

eddy --imain=<Input image> --mask=<Input brain mask> --acqp=<Acquisition_parameters file> --index=<Index file> --bvecs=<b-vector file> --bvals=<b-values file> --topup=<TOPUP output prefix> --out=<Output prefix>
  • –imain: 拡散MRI画像(すべてのVolume)
  • –mask: 脳実質のマスク画像(Binaryファイル)
  • –acqp: 拡散MRI撮像パラメータ (PE方向およびTotal readout time)が記載されたテキストファイル
  • –index: –imainで入力した拡散MRI画像の撮像パラメータが–acqpで指定されたAcquisition_parametersの内どちらのものか、というインデックス
  • –bvecs: b-vectorが記載されたテキストファイル
  • –bvalues: b-valuesが記載されたテキストファイル
  • –topup: TOPUPの出力ファイルの接頭辞
  • –out: EDDYの出力ファイルの接頭辞

3.4. 頭の動き補正の効果


4. サンプルコード

以下のコードは、前処理前の拡散MRI画像データからTOPUP/EDDYによる補正を実行するものである。

#!/bin/sh

Usage() {
	cat <<EOF

Usage: scr [work patient directory]

EOF
	exit 1
}
# ディレクトリ構造
# Subj001
# ├── DWI_AP.nii.gz: Diffusion MRI (Phase encoding direction: AP)
# ├── b0_PA.nii.gz: b=0 (Phase encoding direction:PA)
# ├── bvals.bval: b-values
# └── bvecs.bvec: b-vectors

# Make pair b0 image using two types of phase encoding directoins (AP/PA)
## b0_APの作成
fslroi DWI_AP b0_AP 0 1
## b0_PAの作成
fslroi b0_PA b0_PA 0 1
fslmerge -t both_b0 b0_AP b0_PA

# TOPUP
## Acquisition_parameters
## EPI factor = 130
## Echo Spacing =0.59
printf "0 -1 0 0.07378\n0 1 0 0.07378" >acquisition_parameters.txt

echo "TOPUP processing..."
topup --imain=both_b0 --datain=acquisition_parameters.txt --config=b02b0.cnf --out=topup_results --fout=field --iout=unwarped_images
echo "TOPUP done"

# Apply TOPUP
echo "Applying TOPUP..."
applytopup --imain=DWI_AP --datain=acquisition_parameters.txt --inindex=1 --topup=topup_results --method=jac --out=applytopup_results
echo "Applying TOPUP done"

# EDDY
fslmaths unwarped_images -Tmean unwarped_images
bet unwarped_images unwarped_images_brain -f 0.2 -m

vol_n=$(fslhd DWI_AP |grep ^dim4|tr -s "\t" " "|cut -d " " -f2)
array=()
for i in $(seq 0 $(( vol_n-1 ))); do
	array+=(1)
done
echo "${array[@]}" >index.txt

echo "EDDY processing..."
eddy --imain=DWI_AP --mask=unwarped_images_brain_mask --acqp=acquisition_parameters.txt --index=index.txt --bvecs=bvecs.bvec --bvals=bvals.bval --topup=topup_results --out=eddy_result_data
echo "EDDY done"

【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した結果(緑)を示す。

【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

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

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 でも全く問題なくできることがわかりましたので、タイトルを変更しました。

続きを読む

【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

【FSL/MRtrix】画像の切り取り・マスキング ~Masking~

目的

  • 画像の切り取り・マスキング ~Masking~

FSLを用いる場合

コマンド

FSLで画像の切り取り・マスキングをするには、fslmaths-masオプションを使用する。

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

クリックして展開
Usage: fslmaths [-dt <datatype>] <first_input> [operations and inputs] <output> [-odt <datatype>]

Datatype information:
 -dt sets the datatype used internally for calculations (default float for all except double images)
 -odt sets the output datatype ( default is float )
 Possible datatypes are: char short int float double input
 "input" will set the datatype to that of the original image

Binary operations:
  (some inputs can be either an image or a number)
 -add   : add following input to current image
 -sub   : subtract following input from current image
 -mul   : multiply current image by following input
 -div   : divide current image by following input
 -rem   : modulus remainder - divide current image by following input and take remainder
 -mas   : use (following image>0) to mask current image
 -thr   : use following number to threshold current image (zero anything below the number)
 -thrp  : use following percentage (0-100) of ROBUST RANGE to threshold current image (zero anything below the number)
 -thrP  : use following percentage (0-100) of ROBUST RANGE of non-zero voxels and threshold below
 -uthr  : use following number to upper-threshold current image (zero anything above the number)
 -uthrp : use following percentage (0-100) of ROBUST RANGE to upper-threshold current image (zero anything above the number)
 -uthrP : use following percentage (0-100) of ROBUST RANGE of non-zero voxels and threshold above
 -max   : take maximum of following input and current image
 -min   : take minimum of following input and current image
 -seed  : seed random number generator with following number
 -restart : replace the current image with input for future processing operations
 -save : save the current working image to the input filename

Basic unary operations:
 -exp   : exponential
 -log   : natural logarithm
 -sin   : sine function
 -cos   : cosine function
 -tan   : tangent function
 -asin  : arc sine function
 -acos  : arc cosine function
 -atan  : arc tangent function
 -sqr   : square
 -sqrt  : square root
 -recip : reciprocal (1/current image)
 -abs   : absolute value
 -bin   : use (current image>0) to binarise
 -binv  : binarise and invert (binarisation and logical inversion)
 -fillh : fill holes in a binary mask (holes are internal - i.e. do not touch the edge of the FOV)
 -fillh26 : fill holes using 26 connectivity
 -index : replace each nonzero voxel with a unique (subject to wrapping) index number
 -grid <value> <spacing> : add a 3D grid of intensity <value> with grid spacing <spacing>
 -edge  : edge strength
 -tfce <H> <E> <connectivity>: enhance with TFCE, e.g. -tfce 2 0.5 6 (maybe change 6 to 26 for skeletons)
 -tfceS <H> <E> <connectivity> <X> <Y> <Z> <tfce_thresh>: show support area for voxel (X,Y,Z)
 -nan   : replace NaNs (improper numbers) with 0
 -nanm  : make NaN (improper number) mask with 1 for NaN voxels, 0 otherwise
 -rand  : add uniform noise (range 0:1)
 -randn : add Gaussian noise (mean=0 sigma=1)
 -inm <mean> :  (-i i ip.c) intensity normalisation (per 3D volume mean)
 -ing <mean> :  (-I i ip.c) intensity normalisation, global 4D mean)
 -range : set the output calmin/max to full data range

Matrix operations:
 -tensor_decomp : convert a 4D (6-timepoint )tensor image into L1,2,3,FA,MD,MO,V1,2,3 (remaining image in pipeline is FA)

Kernel operations (set BEFORE filtering operation if desired):
 -kernel 3D : 3x3x3 box centered on target voxel (set as default kernel)
 -kernel 2D : 3x3x1 box centered on target voxel
 -kernel box    <size>     : all voxels in a cube of width <size> mm centered on target voxel
 -kernel boxv   <size>     : all voxels in a cube of width <size> voxels centered on target voxel, CAUTION: size should be an odd number
 -kernel boxv3  <X> <Y> <Z>: all voxels in a cuboid of dimensions X x Y x Z centered on target voxel, CAUTION: size should be an odd number
 -kernel gauss  <sigma>    : gaussian kernel (sigma in mm, not voxels)
 -kernel sphere <size>     : all voxels in a sphere of radius <size> mm centered on target voxel
 -kernel file   <filename> : use external file as kernel

Spatial Filtering operations: N.B. all options apart from -s use the default kernel or that _previously_ specified by -kernel
 -dilM    : Mean Dilation of non-zero voxels
 -dilD    : Modal Dilation of non-zero voxels
 -dilF    : Maximum filtering of all voxels
 -dilall  : Apply -dilM repeatedly until the entire FOV is covered
 -ero     : Erode by zeroing non-zero voxels when zero voxels found in kernel
 -eroF    : Minimum filtering of all voxels
 -fmedian : Median Filtering 
 -fmean   : Mean filtering, kernel weighted (conventionally used with gauss kernel)
 -fmeanu  : Mean filtering, kernel weighted, un-normalised (gives edge effects)
 -s <sigma> : create a gauss kernel of sigma mm and perform mean filtering
 -subsamp2  : downsamples image by a factor of 2 (keeping new voxels centred on old)
 -subsamp2offc  : downsamples image by a factor of 2 (non-centred)

Dimensionality reduction operations:
  (the "T" can be replaced by X, Y or Z to collapse across a different dimension)
 -Tmean   : mean across time
 -Tstd    : standard deviation across time
 -Tmax    : max across time
 -Tmaxn   : time index of max across time
 -Tmin    : min across time
 -Tmedian : median across time
 -Tperc <percentage> : nth percentile (0-100) of FULL RANGE across time
 -Tar1    : temporal AR(1) coefficient (use -odt float and probably demean first)

Basic statistical operations:
 -pval    : Nonparametric uncorrected P-value, assuming timepoints are the permutations; first timepoint is actual (unpermuted) stats image
 -pval0   : Same as -pval, but treat zeros as missing data
 -cpval   : Same as -pval, but gives FWE corrected P-values
 -ztop    : Convert Z-stat to (uncorrected) P
 -ptoz    : Convert (uncorrected) P to Z
 -rank    : Convert data to ranks (over T dim)
 -ranknorm: Transform to Normal dist via ranks

Multi-argument operations:
 -roi <xmin> <xsize> <ymin> <ysize> <zmin> <zsize> <tmin> <tsize> : zero outside roi (using voxel coordinates). Inputting -1 for a size will set it to the full image extent for that dimension.
 -bptf  <hp_sigma> <lp_sigma> : (-t in ip.c) Bandpass temporal filtering; nonlinear highpass and Gaussian linear lowpass (with sigmas in volumes, not seconds); set either sigma<0 to skip that filter
 -roc <AROC-thresh> <outfile> [4Dnoiseonly] <truth> : take (normally binary) truth and test current image in ROC analysis against truth. <AROC-thresh> is usually 0.05 and is limit of Area-under-ROC measure FP axis. <outfile> is a text file of the ROC curve (triplets of values: FP TP threshold). If the truth image contains negative voxels these get excluded from all calculations. If <AROC-thresh> is positive then the [4Dnoiseonly] option needs to be set, and the FP rate is determined from this noise-only data, and is set to be the fraction of timepoints where any FP (anywhere) is seen, as found in the noise-only 4d-dataset. This is then controlling the FWE rate. If <AROC-thresh> is negative the FP rate is calculated from the zero-value parts of the <truth> image, this time averaging voxelwise FP rate over all timepoints. In both cases the TP rate is the average fraction of truth=positive voxels correctly found.

Combining 4D and 3D images:
 If you apply a Binary operation (one that takes the current image and a new image together), when one is 3D and the other is 4D,
 the 3D image is cloned temporally to match the temporal dimensions of the 4D image.

e.g. fslmaths inputVolume -add inputVolume2 output_volume
     fslmaths inputVolume -add 2.5 output_volume
     fslmaths inputVolume -add 2.5 -mul inputVolume2 output_volume

     fslmaths 4D_inputVolume -Tmean -mul -1 -add 4D_inputVolume demeaned_4D_inputVolume

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

fslmaths <入力画像> -mas <マスク画像> <出力画像>

使用例

頭蓋除去されていないFA画像とマスク画像(緑)を、重ね合わせて表示した画像を以下に示す。

頭蓋除去されていないFA画像(FA.nii.gz)をマスク画像(mask.nii.gz)でマスキングするには、以下のコマンドを実行する。

fslmaths FA.nii.gz -mas mask.nii.gz FA_masked.nii.gz

マスキングして、頭蓋除去したFA画像は以下。

MRtrixを用いる場合

コマンド

MRtrixで画像の切り取り・マスキングをするには、mrcalc-multオプションを使用する。

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

クリックして展開
SYNOPSIS

     Apply generic voxel-wise mathematical operations to images

USAGE

     mrcalc [ options ] operand [ operand ... ]

        operand      an input image, intensity value, or the special keywords
                     'rand' (random number between 0 and 1) or 'randn' (random
                     number from unit std.dev. normal distribution) or the
                     mathematical constants 'e' and 'pi'.


DESCRIPTION

     This command will only compute per-voxel operations. Use 'mrmath' to
     compute summary statistics across images or along image axes.

     This command uses a stack-based syntax, with operators (specified using
     options) operating on the top-most entries (i.e. images or values) in the
     stack. Operands (values or images) are pushed onto the stack in the order
     they appear (as arguments) on the command-line, and operators (specified
     as options) operate on and consume the top-most entries in the stack, and
     push their output as a new entry on the stack.

     As an additional feature, this command will allow images with different
     dimensions to be processed, provided they satisfy the following
     conditions: for each axis, the dimensions match if they are the same size,
     or one of them has size one. In the latter case, the entire image will be
     replicated along that axis. This allows for example a 4D image of size [ X
     Y Z N ] to be added to a 3D image of size [ X Y Z ], as if it consisted of
     N copies of the 3D image along the 4th axis (the missing dimension is
     assumed to have size 1). Another example would a single-voxel 4D image of
     size [ 1 1 1 N ], multiplied by a 3D image of size [ X Y Z ], which would
     allow the creation of a 4D image where each volume consists of the 3D
     image scaled by the corresponding value for that volume in the
     single-voxel image.

EXAMPLE USAGES

     Double the value stored in every voxel:
       $ mrcalc a.mif 2 -mult r.mif
     This performs the operation: r = 2*a  for every voxel a,r in images a.mif
     and r.mif respectively.

     A more complex example:
       $ mrcalc a.mif -neg b.mif -div -exp 9.3 -mult r.mif
     This performs the operation: r = 9.3*exp(-a/b)

     Another complex example:
       $ mrcalc a.mif b.mif -add c.mif d.mif -mult 4.2 -add -div r.mif
     This performs: r = (a+b)/(c*d+4.2).

     Rescale the densities in a SH l=0 image:
       $ mrcalc ODF_CSF.mif 4 pi -mult -sqrt -div ODF_CSF_scaled.mif
     This applies the spherical harmonic basis scaling factor: 1.0/sqrt(4*pi),
     such that a single-tissue voxel containing the same intensities as the
     response function of that tissue should contain the value 1.0.

basic operations

  -abs  (multiple uses permitted)
     |%1| : return absolute value (magnitude) of real or complex number

  -neg  (multiple uses permitted)
     -%1 : negative value

  -add  (multiple uses permitted)
     (%1 + %2) : add values

  -subtract  (multiple uses permitted)
     (%1 - %2) : subtract nth operand from (n-1)th

  -multiply  (multiple uses permitted)
     (%1 * %2) : multiply values

  -divide  (multiple uses permitted)
     (%1 / %2) : divide (n-1)th operand by nth

  -min  (multiple uses permitted)
     min (%1, %2) : smallest of last two operands

  -max  (multiple uses permitted)
     max (%1, %2) : greatest of last two operands

comparison operators

  -lt  (multiple uses permitted)
     (%1 < %2) : less-than operator (true=1, false=0)

  -gt  (multiple uses permitted)
     (%1 > %2) : greater-than operator (true=1, false=0)

  -le  (multiple uses permitted)
     (%1 <= %2) : less-than-or-equal-to operator (true=1, false=0)

  -ge  (multiple uses permitted)
     (%1 >= %2) : greater-than-or-equal-to operator (true=1, false=0)

  -eq  (multiple uses permitted)
     (%1 == %2) : equal-to operator (true=1, false=0)

  -neq  (multiple uses permitted)
     (%1 != %2) : not-equal-to operator (true=1, false=0)

conditional operators

  -if  (multiple uses permitted)
     (%1 ? %2 : %3) : if first operand is true (non-zero), return second
     operand, otherwise return third operand

  -replace  (multiple uses permitted)
     (%1, %2 -> %3) : Wherever first operand is equal to the second operand,
     replace with third operand

power functions

  -sqrt  (multiple uses permitted)
     sqrt (%1) : square root

  -pow  (multiple uses permitted)
     %1^%2 : raise (n-1)th operand to nth power

nearest integer operations

  -round  (multiple uses permitted)
     round (%1) : round to nearest integer

  -ceil  (multiple uses permitted)
     ceil (%1) : round up to nearest integer

  -floor  (multiple uses permitted)
     floor (%1) : round down to nearest integer

logical operators

  -not  (multiple uses permitted)
     !%1 : NOT operator: true (1) if operand is false (i.e. zero)

  -and  (multiple uses permitted)
     (%1 && %2) : AND operator: true (1) if both operands are true (i.e.
     non-zero)

  -or  (multiple uses permitted)
     (%1 || %2) : OR operator: true (1) if either operand is true (i.e.
     non-zero)

  -xor  (multiple uses permitted)
     (%1 ^^ %2) : XOR operator: true (1) if only one of the operands is true
     (i.e. non-zero)

classification functions

  -isnan  (multiple uses permitted)
     isnan (%1) : true (1) if operand is not-a-number (NaN)

  -isinf  (multiple uses permitted)
     isinf (%1) : true (1) if operand is infinite (Inf)

  -finite  (multiple uses permitted)
     finite (%1) : true (1) if operand is finite (i.e. not NaN or Inf)

complex numbers

  -complex  (multiple uses permitted)
     (%1 + %2 i) : create complex number using the last two operands as
     real,imaginary components

  -polar  (multiple uses permitted)
     (%1 /_ %2) : create complex number using the last two operands as
     magnitude,phase components (phase in radians)

  -real  (multiple uses permitted)
     real (%1) : real part of complex number

  -imag  (multiple uses permitted)
     imag (%1) : imaginary part of complex number

  -phase  (multiple uses permitted)
     phase (%1) : phase of complex number (use -abs for magnitude)

  -conj  (multiple uses permitted)
     conj (%1) : complex conjugate

  -proj  (multiple uses permitted)
     proj (%1) : projection onto the Riemann sphere

exponential functions

  -exp  (multiple uses permitted)
     exp (%1) : exponential function

  -log  (multiple uses permitted)
     log (%1) : natural logarithm

  -log10  (multiple uses permitted)
     log10 (%1) : common logarithm

trigonometric functions

  -cos  (multiple uses permitted)
     cos (%1) : cosine

  -sin  (multiple uses permitted)
     sin (%1) : sine

  -tan  (multiple uses permitted)
     tan (%1) : tangent

  -acos  (multiple uses permitted)
     acos (%1) : inverse cosine

  -asin  (multiple uses permitted)
     asin (%1) : inverse sine

  -atan  (multiple uses permitted)
     atan (%1) : inverse tangent

hyperbolic functions

  -cosh  (multiple uses permitted)
     cosh (%1) : hyperbolic cosine

  -sinh  (multiple uses permitted)
     sinh (%1) : hyperbolic sine

  -tanh  (multiple uses permitted)
     tanh (%1) : hyperbolic tangent

  -acosh  (multiple uses permitted)
     acosh (%1) : inverse hyperbolic cosine

  -asinh  (multiple uses permitted)
     asinh (%1) : inverse hyperbolic sine

  -atanh  (multiple uses permitted)
     atanh (%1) : inverse hyperbolic tangent

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.

基本的な使い方は、以下の通り。入力画像とバイナリーマスク画像(二値画像)を掛け算することで、マスキングをする。

mrcalc <入力画像> <バイナリーマスク画像> -mult <出力画像>

使用例

頭蓋除去されていないFA画像とマスク画像(緑)を、重ね合わせて表示した画像を以下に示す。

頭蓋除去されていないFA画像(FA.nii.gz)をマスク画像(mask.nii.gz)でマスキングするには、以下のコマンドを実行する。

mrcalc FA.nii.gz  mask.nii.gz -mult FA_masked.nii.gz

マスキングして、頭蓋除去したFA画像は以下。