【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.)
【FSL】fslstatsを使ったROI解析
1. 目的
2. fslstats
2.1. binary maskを使う場合
2.2. index maskを使う場合
3. 使用例
3.1. フォルダ構造
3.2. コード
3.3. 結果
【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)にレジストレーション
【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解析では、次のような処理をする。
- データの準備
- 解析パラメータの設定
gbss_1_T1wpreproc:T1WIの前処理gbss_2_DWIwpreproc:拡散定量値の前処理gbss_3_skelpreproc:拡散定量値をスケルトンに投影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. 解析パラメータの設定
解析に用いるパラメータを設定する。
PEDIR、ECHOSPACINGはepi_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の前処理内容は、次の通り。
- 脳頭蓋除去
- 灰白質をセグメント
- 灰白質を標準脳(MNI152)に位置合わせ
- 標準空間上の各被験者の灰白質を平均化
以上の処理内容を実行するために、関数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:拡散定量値の前処理
拡散定量値の前処理の過程は、次の通り。
- b0をT1WIに位置合わせ
- b0 to T1WIとT1WI to 標準脳のWarp Fieldを使って、拡散定量値マップを標準空間に移動
- 標準空間上にあるすべての被験者の拡散定量値マップをひとつの4D volume dataにまとめる
- 平均灰白質マスク
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_
以下は、標準空間上における灰白質領域のFA画像である。

2.5. gbss_3_skelpreproc:拡散定量値をスケルトンに投影
拡散定量値をスケルトンに投影するまでの手順は、次の通り。
- 平均灰白質画像からスケルトンを生成
- 灰白質スケルトンをしきい値処理をしてスケルトンマスクを生成
- スケルトンの距離マップ(distance map)を生成
- 各拡散定量値マップを、スケルトンに投影
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
【QSM】QSM解析 ~MEDIの使い方~
1. 目的
2. 準備
2.1. MEDI toolboxにパスを通す
2.2. QSM解析実行のコード準備
2.3. Write_DICOM.mを修正
2.4. QSM_3Dを集める
3. 実行
3.1. 作業フォルダ(analyzeフォルダ)へ移動
3.2. runme.mを実行
3.3. Outputファイルの確認
4. 複数の被験者データをまとめて処理したい場合
4.1. run_sequential.m
1. 目的
Cornell大学が開発したQSM software (MEDI)を使って解析

(Cornell MRI Research Lab. QSMより、URL: http://pre.weill.cornell.edu/mri/pages/qsm.html)
2. 準備
####MEDI toolboxのダウンロード
こちらからMEDI_toolboxをダウンロード
MEDI_toolboxの中身はこのような感じ。

2.1. MEDI toolboxにパスを通す
MATLABを起動し、パスの設定からサブフォルダも追加をクリックしfunctionsフォルダを選択した後、保存。

2.2. QSM解析実行のコード準備
QSM解析を実行するためのコードは以下。
以下をrunme.m(※任意の名前でOK)で保存し、作業フォルダとなるanalyzeフォルダ(※任意の名前でOK)に保存する。
clear all;clc;close all;
% Set path
% MEDI_set_path
% STEP 1: Import data
[iField,voxel_size,matrix_size,CF,delta_TE,TE,B0_dir,files]=Read_DICOM('AXL_QSM');
%%%% Generate the Magnitude image %%%%
iMag = sqrt(sum(abs(iField).^2,4));
%%%%% provide a Mask here if possible %%%%%%
if (~exist('Mask','var'))
% Mask = genMask(iField, voxel_size);
Mask = BET(iMag,matrix_size,voxel_size);
end
%%%%% provide a noise_level here if possible %%%%%%
if (~exist('noise_level','var'))
noise_level = calfieldnoise(iField, Mask);
end
%%%%% normalize signal intensity by noise to get SNR %%%
iField = iField/noise_level;
% STEP 2a: Field Map Estimation
%%%%%Estimate the frequency offset in each of the voxel using a
%%%%%complex fitting %%%%
[iFreq_raw N_std] = Fit_ppm_complex(iField);
% STEP 2b: Spatial phase unwrapping %%%%
iFreq = unwrapPhase(iMag, iFreq_raw, matrix_size);
% STEP 2c: Background Field Removal
%%%% Background field removal %%%%
[RDF shim] = PDF(iFreq, N_std, Mask,matrix_size,voxel_size, B0_dir);
% CSF Mask for zero referencing
R2s=arlo(TE,abs(iField));
Mask_CSF = extract_CSF(R2s, Mask, voxel_size);
% STEP 3: Dipole Inversion
save RDF.mat RDF iFreq iFreq_raw iMag N_std Mask matrix_size...
voxel_size delta_TE CF B0_dir Mask_CSF;
%%%% run MEDI %%%%%
QSM = MEDI_L1('lambda',1000, 'smv', 5, 'merit');
%%% save to DICOM
% ignore warnings...
Write_DICOM(QSM,files,'QSM')

2.3. Write_DICOM.mを修正
解析後のQSMの保存には、Write_DICOM.mを使用する。(runme.mを実行すれば自動で実行される。)
修正前のファイルでは、dicomwriteでエラーが発生する。
そのため、MEDI_toolbox/functions/Write_DICOM.mを以下のコードにそっくりそのまま書き換える。
やっていることは、dicomwriteのCreateModeをCopyからCreateへチェンジ。
function Write_DICOM(M,files,outdir,opts)
%WRITE_DICOM Summary of this function goes here
% Detailed explanation goes here
defopts.SeriesDescription = 'QSM';
defopts.SeriesInstanceUID = [];
defopts.SeriesNumber = [];
defopts.Convert = @(x) convert(x);
defopts.Window = 0.500;
defopts.Level = 0;
defopts.flag3D = 0;
% defopts.EchoNumber = [];
% defopts.EchoTime = 0.0;
if nargin<4; opts=struct; end
deffields=fieldnames(defopts);
for i=1:length(deffields)
if ~isfield(opts, deffields{i})
opts.(deffields{i})=defopts.(deffields{i});
end
end
if isfield(files,'M')
filenames=files.M;
elseif isfield(files,'R')
filenames=files.R;
elseif isfield(files, 'info3D')
opts.flag3D=1;
else
error('No filenames (M nor R) found.');
end
flag_signed=min(M(:))<0;
if ~opts.flag3D
if size(M,3) ~= size(filenames,1)
error([num2str(size(filenames,1)) ' filenames given for ' num2str(size(M,3)) ' slices.']);
end
end
if isempty(opts.SeriesInstanceUID)
opts.SeriesInstanceUID=dicomuid;
end
progress='';
mkdir(outdir);
warning('off','images:dicom_add_attr:invalidAttribChar');
load_flag=1;
insert_flag=~opts.flag3D;
for slice=1:size(M,3)
if load_flag
if opts.flag3D
filename=files.info3D.filename;
else
filename=filenames{slice,end};
end
info = dicominfo(filename);
info.SeriesDescription = opts.SeriesDescription;
info.SeriesInstanceUID = opts.SeriesInstanceUID;
if isempty(opts.SeriesNumber)
opts.SeriesNumber=info.SeriesNumber*100;
end
info.SeriesNumber = opts.SeriesNumber;
info.SOPInstanceUID = dicomuid;
info.InstanceNumber = slice;
if opts.flag3D
load_flag=0;
end
end
if opts.flag3D
item=files.info3D.slice2index{slice};
% info.PerFrameFunctionalGroupsSequence.(item).PlanePositionSequence.Item_1.ImagePositionPatient;
info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.SOPInstanceUID = dicomuid;
end
im = M(:,:,slice);
if (isfield(info, 'SmallestImagePixelValue'))
info.SmallestImagePixelValue=opts.Convert(min(im(:)));
end
if (isfield(info, 'LargestImagePixelValue'))
info.LargestImagePixelValue=opts.Convert(max(im(:)));
end
if (isfield(info, 'RescaleIntercept'))
info.RescaleIntercept=0;
end
if (isfield(info, 'RescaleSlope'))
info.RescaleSlope=1;
end
info.WindowCenter=opts.Convert(opts.Level);
info.WindowWidth=opts.Convert(opts.Window);
% if opts.flag3D
% info.PerFrameFunctionalGroupsSequence.Item_1.PixelValueTransformationSequence.Item_1.RescaleIntercept=0;
% info.PerFrameFunctionalGroupsSequence.Item_1.Private_2005_140f.Item_1.RescaleIntercept=0;
% info.PerFrameFunctionalGroupsSequence.Item_1.PixelValueTransformationSequence.Item_1.RescaleSlope=1;
% info.PerFrameFunctionalGroupsSequence.Item_1.Private_2005_140f.Item_1.RescaleSlope=1;
% info.PerFrameFunctionalGroupsSequence.Item_1.FrameVOILUTSequence.Item_1.WindowCenter=opts.Convert(opts.Level);
% info.PerFrameFunctionalGroupsSequence.Item_1.FrameVOILUTSequence.Item_1.WindowWidth=opts.Convert(opts.Window);
% end
info.SamplesPerPixel=1;
info.BitsAllocated=16;
info.BitsStored=16;
info.HighBit=15;
info.PixelRepresentation=flag_signed;
if size(M,3)==slice
insert_flag=1;
end
if insert_flag
outfile=fullfile(outdir,sprintf('IM%05d.dcm', slice));
print_progress(outfile);
if opts.flag3D
f=fieldnames(info.PerFrameFunctionalGroupsSequence);
f=setdiff(f,files.info3D.slice2index,'stable');
for i=1:length(f)
info.PerFrameFunctionalGroupsSequence=rmfield(info.PerFrameFunctionalGroupsSequence, f{i});
end
for i=1:length(files.info3D.slice2index)
item=files.info3D.slice2index{i};
info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.InstanceNumber=1;
info.PerFrameFunctionalGroupsSequence.(item).PixelValueTransformationSequence.Item_1.RescaleIntercept=0;
info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.RescaleIntercept=0;
info.PerFrameFunctionalGroupsSequence.(item).PixelValueTransformationSequence.Item_1.RescaleSlope=1;
info.PerFrameFunctionalGroupsSequence.(item).Private_2005_140f.Item_1.RescaleSlope=1;
info.PerFrameFunctionalGroupsSequence.(item).FrameVOILUTSequence.Item_1.WindowCenter=opts.Convert(opts.Level);
info.PerFrameFunctionalGroupsSequence.(item).FrameVOILUTSequence.Item_1.WindowWidth=opts.Convert(opts.Window);
end
info.NumberOfFrames=length(files.info3D.slice2index);
sz=size(M);
M=reshape(opts.Convert(M), sz(1), sz(2), 1, []);
M=permute(M, [2 1 3 4]);
if isfield(files, 'slices_added')
if files.slices_added
warning('Removing empty slice at bottom of volume');
M=M(:,:,1:end-1);
end
end
%dicomwrite(M,outfile, ...
% 'CreateMode', 'copy', 'WritePrivate', true, info);
dicomwrite(M,outfile, ...
'WritePrivate', true, info);
else
if isfield(files, 'slices_added')
if files.slices_added
warning('Removing empty slice at bottom of volume');
M=M(:,:,1:end-1);
end
end
%dicomwrite(opts.Convert(M(:,:,slice))',outfile, ...
% 'CreateMode', 'copy', 'WritePrivate', true, info);
%
dicomwrite(opts.Convert(M(:,:,slice))',outfile, ...
'WritePrivate', true, info);
end
end
end
fprintf('\n');
function print_progress(arg)
num=length(progress);
num=num-numel(regexp(progress, '\\\\'));
for ii=1:num; fprintf('\b'); end
progress=['Writing file ' arg];
progress=regexprep(progress,'\','\\\\');
fprintf(progress);
end
function y=convert(x)
if flag_signed
y=int16(x*1000);
else
y=uint16(x*1000);
end
end
end
2.4. QSM_3Dを集める
各フォルダには、強度画像と位相画像が入っている。
TEを7回変えて撮像している。
国際脳QSMの撮像は、1回撮像に128 slicesなので強度画像と位相画像はそれぞれ、896枚(=128 slices ×7 phase)


(Cornell MRI Research Lab. QSMより、URL: http://pre.weill.cornell.edu/mri/pages/qsm.html)
QSM_3Dにある強度画像と位相画像をanalyzeフォルダのrawdata(※必ずrawdata)にまとめて保存。


3. 実行
3.1. 作業フォルダ(analyzeフォルダ)へ移動
MATLABを起動し、analyzeフォルダへ移動。

3.2. runme.mを実行
runme.mをMATLABへDrag & Dropし、QSM解析を実行。
3.3. Outputファイルの確認
runme.mを実行後、
- QSM
- results
フォルダが生成される。

QSMフォルダにQSM画像がDICOM形式で保存される。


(Cornell MRI Research Lab. QSMより、URL: http://pre.weill.cornell.edu/mri/pages/qsm.html)
4. 複数の被験者データをまとめて処理したい場合
被験者ごとにフォルダを作成し一つのフォルダにまとめます。
さらに、後で紹介するまとめてMEDIを実行するスクリプトも入れておきます。

各被験者フォルダには強度画像と位相画像のDICOMがまとめて入ったrawdataフォルダがあります。

以下のスクリプトを実行すると、すべての被験者のQSM mapが計算できます。
4.1. run_sequential.m
clear all;clc;close all;
% change direct to study folder
selpath = uigetdir('Select the folder including all subject');
cd(selpath);
% list folder
folderlist = dir(selpath);
folderlist = folderlist(~ismember({folderlist.name}, {'.', '..'}));
folderlist = folderlist([folderlist.isdir]);
% run MEDI
for i = 1:length(folderlist)
disp(['processing ' folderlist(i).name ' ...'])
main(folderlist(i).name, selpath);
end
% define function of MEDI processings
function main(folder, basepath)
% initialize
clearvars -except selpath folderlist folder basepath i
% move to subject folder
cd(folder)
% STEP 1: Import data
[iField,voxel_size,matrix_size,CF,delta_TE,TE,B0_dir,files]=Read_DICOM('rawdata');
%%%% Generate the Magnitude image %%%%
iMag = sqrt(sum(abs(iField).^2,4));
%%%%% provide a Mask here if possible %%%%%%
if (~exist('Mask','var'))
% Mask = genMask(iField, voxel_size);
Mask = BET(iMag,matrix_size,voxel_size);
end
%%%%% provide a noise_level here if possible %%%%%%
if (~exist('noise_level','var'))
noise_level = calfieldnoise(iField, Mask);
end
%%%%% normalize signal intensity by noise to get SNR %%%
iField = iField/noise_level;
% STEP 2a: Field Map Estimation
%%%%%Estimate the frequency offset in each of the voxel using a
%%%%%complex fitting %%%%
[iFreq_raw N_std] = Fit_ppm_complex(iField);
% STEP 2b: Spatial phase unwrapping %%%%
iFreq = unwrapPhase(iMag, iFreq_raw, matrix_size);
% STEP 2c: Background Field Removal
%%%% Background field removal %%%%
[RDF shim] = PDF(iFreq, N_std, Mask,matrix_size,voxel_size, B0_dir);
% CSF Mask for zero referencing
R2s=arlo(TE,abs(iField));
Mask_CSF = extract_CSF(R2s, Mask, voxel_size);
% STEP 3: Dipole Inversion
save RDF.mat RDF iFreq iFreq_raw iMag N_std Mask matrix_size...
voxel_size delta_TE CF B0_dir Mask_CSF;
%%%% run MEDI %%%%%
QSM = MEDI_L1('lambda',1000, 'smv', 5, 'merit');
%%% save to DICOM
% ignore warnings...
Write_DICOM(QSM,files,'QSM')
% back to study folder
cd(basepath)
end
【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の解析手順は、次の通り。
- ファイルの準備
- 脳解剖学的標準化
- 平滑化
- 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解析では、次のような処理をする。
- 拡散MRIからFA画像を計算
tbss_1_preproc:TBSSの準備tbss_2_reg:FAを標準空間に非線形位置合わせtbss_3_postreg:平均FA画像を生成し、FAスケルトンを生成tbss_4_prestats:被験者ごとのFA画像を平均スケルトンに投影tbss_non_FA:FA画像以外の定量値をスケルトンに投影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 :出力画像
- -d :デザインマトリックス
- -t :デザインコントラスト
- -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}<em>skeletonised \
-o stats/tbss</em>${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解析では、次のような処理をする。
- 3D-T1WIの準備
fslvbm_1_bet: 3D-T1WIの脳頭蓋除去fslvbm_2_template: 被験者の脳から灰白質テンプレート画像を生成fslvbm_3_proc: すべての被験者の灰白質を灰白質テンプレート画像に位置合わせし、ボクセルをモジュレーション後、平滑化randomise:GLMと並べ替え検定(permutation test)
2.1. 1. 3D-T1WIの準備
各被験者の3D-T1WIを準備する。
ここでは、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)いることを想定する。
. ├── Con0001_T1w.nii.gz ├── Con0002_T1w.nii.gz ├── Con0003_T1w.nii.gz ... └── Pat0010_T1w.nii.gz
2.2. 2. fslvbm_1_bet: 3D-T1WIの脳頭蓋除去
ファイルの準備ができたら、fslvbm_1_betコマンドを実行して脳頭蓋除去をする。3D-T1WIに首を含まない場合、-bオプションを指定し、首を含む場合は-Nオプションを指定する。
ここでは、3D-T1WIに首を含んでいたとして、-Nオプションを指定する。
fslvbm_1_bet -N # 首を含む場合 # fslvbm_1_bet -b # 首を含まない場合
処理が完了すると「strucフォルダ」が生成され、その中に頭蓋除去後の画像が保存される。頭蓋除去後の画像は「*_brain.nii.gz」というファイル名で保存される。頭蓋除去前後の画像一覧をみるには、「struc/slicesdir/index.html」を開く。下地が頭蓋除去前の画像であり、赤い線が頭蓋除去後の画像である。

2.3. 3. fslvbm_2_template: 被験者の脳から灰白質テンプレート画像を生成
次に、脳の解剖学的標準化で用いるターゲット画像(標準脳)を、頭蓋除去済みの被験者脳から生成する。
まず最初に、標準脳を生成するために用いる被験者の脳画像リスト(template_list)を作成する。このとき、バイアスがかからないように健常者と患者の人数が同じになるようにリストを作る。
ここでは、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)の合計20名を「template_list」に記載した。
cat template_list
Con0001_T1w.nii.gz
Con0001_T1w.nii.gz
Con0001_T1w.nii.gz
...
Pat0010_T1w.nii.gz
template_listは、3D-T1WI(*_T1w.nii.gz)と同じ階層に配置する。
. ├── template_list # ここに配置 ├── Con0001_T1w.nii.gz ├── Con0002_T1w.nii.gz ├── Con0003_T1w.nii.gz ... └── Pat0010_T1w.nii.gz
標準脳を生成するために用いる被験者の脳画像リスト(template_list)が準備できたら、頭蓋除去済みの脳画像を、灰白質・白質・脳脊髄液にセグメントする。次に、全ての被験者の灰白質(*_struc_GM.nii.gz)を元々ある灰白質標準脳(ICBM-152)にアフィン変換で位置合わせし、平均化する。このようにして得られた平均灰白質画像を、x軸に対して反転し再度平均化して、左右対称の平均灰白質画像(template_GM_init.nii.gz)を生成する。既に作成した「template_list」に記載されている灰白質画像を、先ほど作成した平均灰白質画像に位置合わせをし、平均化、さらに左右反転後に再度平均化する。最後に、元々ある灰白質標準脳(ICBM-152)に非線形位置合わせをし、2x2x2 mm^3の灰白質テンプレートを標準空間に生成する。
fslvbm_2_templateコマンドにはオプションがあり、各被験者の灰白質画像を元々ある灰白質標準脳(ICBM-152)に位置合わせする際に、アフィン変換を用いる場合は-aオプションを、非線形変換を用いる場合は-nオプションを指定する。
ここでは、fslvbm_2_templateに-aオプションを指定して実行する。
fslvbm_2_template -a # Affine registration # fslvbm_2_template -n # non-linear registration
処理が完了するとstrucフォルダに灰白質テンプレート(template_GM_4D.nii.gz)が生成される。

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

2.5. randomise:GLMと並べ替え検定(permutation test)
ガウシアンフィルタで平滑化した灰白質画像を用いて、健常群と患者群の群間比較をする。
まず、GLMのデザインマトリックス(計画行列)とコントラストを設定する。
今回は、健常者10名(ID: Con0001~Con0010)と患者10名(ID: Pat0001~Pat0010)のデータがある。すべての被験者の3D-T1WIをlsコマンドでみると、先に健常者10名の3D-T1WI、次に患者10名の3D-T1WIが並んでいることが分かる。
ls |grep nii
Con0001_T1w.nii.gz
Con0002_T1w.nii.gz
Con0003_T1w.nii.gz
...
Pat0010_T1w.nii.gz
次に、GLMのデザインマトリックス(計画行列)とコントラストを決める設定ファイルを生成する。design_ttest2 <出力ファイル> <健常者数> <患者数>でコマンドを実行。
design_ttest2 stats/design 10 10
statsフォルダに、デザインマトリックス(design.mat)とコントラスト(design.con)が生成される。
デザインマトリックス(design.mat)の中身を確認。
/Matrixの一列目は健常者データであるかどうか、二列目は患者データであるかを0, 1で表している。行の順番は、被験者ファイルの順番(昇順)に対応する。したがって、これらは対応があるようにしておかなければならない。
cat stats/design.mat
/NumWaves 2
/NumPoints 20
/PPheights 1 1
/Matrix
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
コントラスト(design.con)の中身を確認してみる。
/Matrix一列目は健常者の偏回帰係数、二列目は患者の偏回帰係数に対するもので、行は別々のコントラストである。この場合、一行目は健常者>患者の検定、二行目は健常者<患者の検定に相当する。
cat stats/design.con
/NumWaves 2
/NumContrasts 2
/PPheights 1 1
/Matrix
1 -1
-1 1
デザインマトリックス(計画行列)とコントラストの確認ができたら、randomiseコマンド使ってGLMと並べ替え検定(permutation test)を実行する。
randomiseコマンドの各オプションは、次の通り。
- -i:入力画像
- -m:マスク画像
- -o :出力画像
- -o :デザインマトリックス
- -o :デザインコントラスト
- -n:並べ替え検定の数
- –T2:2D最適化を用いたTFCE
- -x:voxel-wiseのcorrected P値マップ
- –uncorrp:un-corrected P値マップ
- -R:統計値マップ
fsleyes $FSLDIR/data/standard/MNI152_T1_2mm \
stats/fslvbm_tfce_corrp_tstat1 \
-cm Red-Yellow -dr 0.95,1
次のようなファイルが、生成される。

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

3. おまけ
検定数が多い場合、有意差があったかどうかをすべて確認するのは大変である。そこで、一目で有意差があるかどうかを判断できるように、各検定ごとのP値マップの最大値を自動計測し、テキスト(corrected_P_report.txt)としてまとめる。
for PMAP in $(ls stats/ | grep tfce_corrp); do
PMAX=$(fslstats stats/${PMAP} -R | cut -d " " -f2)
echo ${PMAP} >>stats/tmp1.txt
echo ${PMAX} >>stats/tmp2.txt
done
paste stats/tmp* >stats/tmp_corrected_P_report.txt
echo -e "$(cat stats/tmp_corrected_P_report.txt)\n\n\n$(cat stats/tmp_corrected_P_report.txt | sort -r -n -k 2)" \
>stats/corrected_P_report.txt
rm stats/tmp*
上のコマンドを実行すると、statsフォルダに各検定とそのP値マップの最大値が記された「corrected_P_report.txt」が出力される。
検定結果を、ファイル名でソート(上段)したものと、P値でソートしたもの(下段)に分けて保存している。
cat stats/corrected_P_report.txt
fslvbm_tfce_corrp_tstat1.nii.gz 0.982391
fslvbm_tfce_corrp_tstat2.nii.gz 0.993181
fslvbm_tfce_corrp_tstat2.nii.gz 0.993181
fslvbm_tfce_corrp_tstat1.nii.gz 0.982391
【DIPY】DIPYを用いた拡散尖度イメージング: DKI

1. 目的
2. 準備
2.1. DIPYのインストール
2.2. 使用データ
2.3. 前処理
3. 拡散尖度イメージング(DKI)
3.1. 必要なパッケージをインポート
3.2. 画像およびMPG軸情報の読み込み
3.3. マスク画像の生成
3.4. モデルフィッティング
3.5. 拡散定量値の計算
3.6. NIfTI形式で保存
3.7. 結果
4. おまけ
1. 目的
- DIPYを用いた拡散尖度イメージング: DKI
2. 準備
2.1. DIPYのインストール
pip3 install dipy
2.2. 使用データ
データを次のフォルダ構造で用意する。
Study/
└── Subject
├── DWI.nii.gz # 拡散MRI
├── DWI_mask.nii.gz # 拡散MRIマスク画像
├── bvals # b-values
└── bvecs # b-vectors
2.3. 前処理
DKI(Diffusion Kurtosis Imaging)前に、拡散MRIの前処理をする。
- 拡散MRIのノイズ除去(Software: MRtrix, DIPY)
- ギブズのリンギングアーチファクト(Gibbs ringing)の除去(Software: MRtrix, DIPY)
- 拡散MRIのバイアス(信号ムラ)補正(Software: MRtrix)
- 拡散MRIの前処理 ~歪み・頭の動き・渦電流の補正(Software: FSL, MRtrix)
3. 拡散尖度イメージング(DKI)
Pythonで以下のコマンドを実行。
3.1. 必要なパッケージをインポート
import numpy as np import matplotlib.pyplot as plt import dipy.reconst.dki as dki from dipy.core.gradients import gradient_table from dipy.io.gradients import read_bvals_bvecs from dipy.io.image import load_nifti, save_nifti from dipy.segment.mask import median_otsu
3.2. 画像およびMPG軸情報の読み込み
DWI_FILE = 'DWI.nii.gz' BVALS_FILE = 'bvals' BVECS_FILE = 'bvecs' data, affine = load_nifti(DWI_FILE) bvals, bvecs = read_bvals_bvecs(BVALS_FILE, BVECS_FILE) gtab = gradient_table(bvals, bvecs)
3.3. マスク画像の生成
median_otsu関数を用いて、b=0画像からマスク画像を生成する。vol_idxには、b0 volumeのvolume indexを渡す。
maskdata, mask = median_otsu(
data, vol_idx=np.where(bvals == 0)[0]) # , dilate=3
3.4. モデルフィッティング
以下のコマンドで、DKIのモデルフィッティングを実行。
dkimodel = dki.DiffusionKurtosisModel(gtab) dkifit = dkimodel.fit(maskdata)
3.5. 拡散定量値の計算
モデルフィッティングができたら、拡散定量値を算出する。
MK = dkifit.mk(0, 3) AK = dkifit.ak(0, 3) RK = dkifit.rk(0, 3)
脳周囲の背景では、フィッティングミスをしてnanとなる場合があるため、そのようなnanを0に置き換える。
MK[np.isnan(MK)] = 0 AK[np.isnan(AK)] = 0 RK[np.isnan(RK)] = 0
3.6. NIfTI形式で保存
save_nifti関数で、画像をNIfTI形式で保存する。
save_nifti('DWI_masked.nii.gz', maskdata.astype(np.float32), affine)
save_nifti('DWI_mask.nii.gz', mask.astype(np.float32), affine)
save_nifti('MK.nii.gz', MK.astype(np.float32), affine)
save_nifti('AK.nii.gz', AK.astype(np.float32), affine)
save_nifti('RK.nii.gz', RK.astype(np.float32), affine)
3.7. 結果
DKIによって算出された定量値画像は、以下の通り。

4. おまけ
DKIでもFA, MD, AD, RDを算出することができる。
# 拡散定量値を算出
FA = dkifit.fa
MD = dkifit.md
AD = dkifit.ad
RD = dkifit.rd
# nanを0に置換
FA[np.isnan(FA)] = 0
MD[np.isnan(MD)] = 0
AD[np.isnan(AD)] = 0
RD[np.isnan(RD)] = 0
# NIfTI形式で保存
save_nifti('FA.nii.gz', FA.astype(np.float32), affine)
save_nifti('MD.nii.gz', MD.astype(np.float32), affine)
save_nifti('AD.nii.gz', AD.astype(np.float32), affine)
save_nifti('RD.nii.gz', RD.astype(np.float32), affine)
以下は、FA, MD, AD, RDをDTI(下段)とDKI(上段)で比較した図である。

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

1. 目的
2. 準備
2.1. インストール
2.2. 使用データ
2.3. 前処理
3. 神経突起イメージング(NODDI)
3.1. AMICOのセットアップ
3.2. データの読み込み
3.3. 応答関数(response function)の計算
3.4. モデルフィッティング
3.5. 結果
1. 目的
- AMICOを用いた神経突起イメージング: NODDI
2. 準備
2.1. インストール
Pythonを使って、AMICOを用いた神経突起イメージング(NODDI)をするために、以下のPythonパッケージをインストールする。
pip3 install dmri-amico
2.2. 使用データ
データを次のフォルダ構造で用意する。
Study/
└── Subject
├── DWI.nii.gz # 拡散MRI
├── DWI_mask.nii.gz # 拡散MRIマスク画像
├── bvals # b-values
└── bvecs # b-vectors
2.3. 前処理
NODDI前に、拡散MRIの前処理をする。
- 拡散MRIのノイズ除去(Software: MRtrix, DIPY)
- ギブズのリンギングアーチファクト(Gibbs ringing)の除去(Software: MRtrix, DIPY)
- 拡散MRIのバイアス(信号ムラ)補正(Software: MRtrix)
- 拡散MRIの前処理 ~歪み・頭の動き・渦電流の補正(Software: FSL, MRtrix)
3. 神経突起イメージング(NODDI)
Pythonで以下のコマンドを実行。
3.1. AMICOのセットアップ
今回使用するファイル等の変数設定をする。
STUDY_DIR='Study' SUBJECT_DIR='Subject' DWI_FILE = 'DWI.nii.gz' DWIMASK_FILE = 'DWI_mask.nii.gz' BVALS_FILE = 'bvals' BVECS_FILE = 'bvecs'
次に、使用するamicoパッケージのをインポートし、セットアップと初期化をする。
import amico amico.core.setup()
3.2. データの読み込み
AMICOを実行するために、Study/Subjectフォルダ(PATH)を指定する。
ae = amico.Evaluation(STUDY_DIR, SUBJECT_DIR)
MPG軸情報(bvals/bvecs)の情報が入ったschemeファイルを生成する。
amico.util.fsl2scheme("{}/{}/{}".format(STUDY_DIR,SUBJECT_DIR,BVALS_FILE), "{}/{}/{}".format(STUDY_DIR,SUBJECT_DIR,BVECS_FILE),schemeFilename = "{}/{}/NODDI_protocol.scheme".format(STUDY_DIR,SUBJECT_DIR))
-> Writing scheme file to [ Study/Subject/NODDI_protocol.scheme ]
'Study/Subject/NODDI_protocol.scheme'
画像を読み込む。
ae.load_data(dwi_filename = DWI_FILE, scheme_filename = 'NODDI_protocol.scheme', mask_filename = DWIMASK_FILE, b0_thr = 0)
-> Loading data:
* DWI signal
- dim = 130 x 130 x 82 x 129
- pixdim = 1.769 x 1.769 x 1.800
* Acquisition scheme
- 129 samples, 2 shells
- 1 @ b=0 , 64 @ b=1000.0 , 64 @ b=2000.0
* Binary mask
- dim = 130 x 130 x 82
- pixdim = 1.769 x 1.769 x 1.800
- voxels = 282878
[ 4.4 seconds ]
-> Preprocessing:
* Normalizing to b0... [ min=-3.28, mean=0.25, max=22.86 ]
* Keeping all b0 volume(s)
[ 1.1 seconds ]
3.3. 応答関数(response function)の計算
NODDIモデルを設定して、応答関数(response function)を計算する。計算が完了するとkernelファイルが生成される。
ae.set_model('NODDI')
ae.generate_kernels()
-> Creating LUT for "NODDI" model:
[==================================================] 100.0%
[ 373.3 seconds ]
作成したkernelファイルを読み込む。
ae.load_kernels()
-> Resampling LUT for subject "Subject":
[==================================================] 100.0%
[ 112.8 seconds ]
3.4. モデルフィッティング
NODDIのモデルフィッティングを開始する。
ae.fit()
-> Fitting "NODDI" model to 282878 voxels:
[==================================================] 100.0%
[ 02h 52m 07s ]
最後に、結果をNIfTIフォーマットで保存する。
ae.save_results()
-> Saving output to "AMICO/NODDI/*":
- configuration [OK]
- FIT_dir.nii.gz [OK]
- FIT_ICVF.nii.gz [OK]
- FIT_OD.nii.gz [OK]
- FIT_ISOVF.nii.gz [OK]
[ DONE ]
3.5. 結果
結果は、「Study/Subject/AMICO/NODDI/」フォルダにある。
Study/Subject/AMICO/NODDI/ ├── FIT_ICVF.nii.gz ├── FIT_ISOVF.nii.gz ├── FIT_OD.nii.gz ├── FIT_dir.nii.gz └── config.pickle
画像はこちら。

【DIPY】DIPYを用いた拡散テンソルイメージング: DTI
1. 目的
2. 準備
2.1. DIPYのインストール
2.2. 使用データ
2.3. 前処理
3. 拡散テンソルイメージング(DTI)
3.1. 必要なパッケージをインポート
3.2. 画像およびMPG軸情報の読み込み
3.3. マスク画像の生成
3.4. モデルフィッティング
3.5. 拡散定量値の計算
3.6. NIfTI形式で保存
3.7. 結果
1. 目的
- DIPYを用いた拡散テンソルイメージング: DTI
2. 準備
2.1. DIPYのインストール
pip3 install dipy
2.2. 使用データ
データを次のフォルダ構造で用意する。
Study/
└── Subject
├── DWI.nii.gz # 拡散MRI
├── DWI_mask.nii.gz # 拡散MRIマスク画像
├── bvals # b-values
└── bvecs # b-vectors
2.3. 前処理
DTI(Diffusion Tensor Imaging)前に、拡散MRIの前処理をする。
- 拡散MRIのノイズ除去(Software: MRtrix, DIPY)
- ギブズのリンギングアーチファクト(Gibbs ringing)の除去(Software: MRtrix, DIPY)
- 拡散MRIのバイアス(信号ムラ)補正(Software: MRtrix)
- 拡散MRIの前処理 ~歪み・頭の動き・渦電流の補正(Software: FSL, MRtrix)
3. 拡散テンソルイメージング(DTI)
Pythonで以下のコマンドを実行。
3.1. 必要なパッケージをインポート
from dipy.segment.mask import median_otsu import numpy as np from dipy.io.image import load_nifti, save_nifti from dipy.io.gradients import read_bvals_bvecs from dipy.core.gradients import gradient_table import dipy.reconst.dti as dti
3.2. 画像およびMPG軸情報の読み込み
DWI_FILE = 'DWI.nii.gz' BVALS_FILE = 'bvals' BVECS_FILE = 'bvecs' data, affine = load_nifti(DWI_FILE) bvals, bvecs = read_bvals_bvecs(BVALS_FILE, BVECS_FILE) gtab = gradient_table(bvals, bvecs)
3.3. マスク画像の生成
median_otsu関数を用いて、b=0画像からマスク画像を生成する。vol_idxには、b0 volumeのvolume indexを渡す。
maskdata, mask = median_otsu(data, vol_idx=np.where(bvals == 0)[0])
3.4. モデルフィッティング
以下のコマンドで、DTIのモデルフィッティングを実行。
tenmodel = dti.TensorModel(gtab) tenfit = tenmodel.fit(maskdata)
3.5. 拡散定量値の計算
モデルフィッティングができたら、拡散定量値を算出する。
FA = tenfit.fa MD = tenfit.md AD = tenfit.ad RD = tenfit.rd colour_FA = tenfit.color_fa
脳周囲の背景では、フィッティングミスをしてnanとなる場合があるため、そのようなnanを0に置き換える。
FA[np.isnan(FA)] = 0 MD[np.isnan(MD)] = 0 AD[np.isnan(AD)] = 0 RD[np.isnan(RD)] = 0 colour_FA[np.isnan(colour_FA)] = 0
3.6. NIfTI形式で保存
save_nifti関数で、画像をNIfTI形式で保存する。
save_nifti('DWI_masked.nii.gz', maskdata.astype(np.float32), affine)
save_nifti('DWI_mask.nii.gz', mask.astype(np.float32), affine)
save_nifti('FA.nii.gz', FA.astype(np.float32), affine)
save_nifti('MD.nii.gz', MD.astype(np.float32), affine)
save_nifti('AD.nii.gz', AD.astype(np.float32), affine)
save_nifti('RD.nii.gz', RD.astype(np.float32), affine)
save_nifti('colour_FA.nii.gz', colour_FA.astype(np.float32), affine)
3.7. 結果
DTIによって算出された定量値画像は、以下の通り。

【MRtrix】 MRtrixを用いた拡散テンソルイメージング: DTI
1. 目的
2. コマンド
3. 使用例
3.1. 前準備
3.2. テンソルの推定(コマンド:dwi2tensor)
3.3. 拡散定量値の算出(コマンド:tensor2metric)
1. 目的
- MRtrixを用いた拡散テンソルイメージング: DTI
2. コマンド
MRtrixを用いて、拡散テンソルイメージング(DTI)をするには、dwi2tensorとtensor2metricコマンドを用いる。
dwi2tensorは拡散MRI画像からテンソルを推定するコマンドで、tensor2metricは推定したテンソルから拡散定量値を算出するコマンドである。
dwi2tensorのヘルプは、次の通り。
クリックして展開
SYNOPSIS
Diffusion (kurtosis) tensor estimation
USAGE
dwi2tensor [ options ] dwi dt
dwi the input dwi image.
dt the output dt image.
DESCRIPTION
By default, the diffusion tensor (and optionally its kurtosis) is fitted
to the log-signal in two steps: firstly, using weighted least-squares
(WLS) with weights based on the empirical signal intensities; secondly, by
further iterated weighted least-squares (IWLS) with weights determined by
the signal predictions from the previous iteration (by default, 2
iterations will be performed). This behaviour can be altered in two ways:
* The -ols option will cause the first fitting step to be performed using
ordinary least-squares (OLS); that is, all measurements contribute equally
to the fit, instead of the default behaviour of weighting based on the
empirical signal intensities.
* The -iter option controls the number of iterations of the IWLS
prodedure. If this is set to zero, then the output model parameters will
be those resulting from the first fitting step only: either WLS by
default, or OLS if the -ols option is used in conjunction with -iter 0.
The tensor coefficients are stored in the output image as follows:
volumes 0-5: D11, D22, D33, D12, D13, D23
If diffusion kurtosis is estimated using the -dkt option, these are stored
as follows:
volumes 0-2: W1111, W2222, W3333
volumes 3-8: W1112, W1113, W1222, W1333, W2223, W2333
volumes 9-11: W1122, W1133, W2233
volumes 12-14: W1123, W1223, W1233
OPTIONS
-ols
perform initial fit using an ordinary least-squares (OLS) fit (see
Description).
-mask image
only perform computation within the specified binary brain mask image.
-b0 image
the output b0 image.
-dkt image
the output dkt image.
-iter integer
number of iterative reweightings for IWLS algorithm (default: 2) (see
Description).
-predicted_signal image
the predicted dwi image.
DW gradient table import options
-grad file
Provide the diffusion-weighted gradient scheme used in the acquisition in
a text file. This should be supplied as a 4xN text file with each line is
in the format [ X Y Z b ], where [ X Y Z ] describe the direction of the
applied gradient, and b gives the b-value in units of s/mm^2. If a
diffusion gradient scheme is present in the input image header, the data
provided with this option will be instead used.
-fslgrad bvecs bvals
Provide the diffusion-weighted gradient scheme used in the acquisition in
FSL bvecs/bvals format files. If a diffusion gradient scheme is present in
the input image header, the data provided with this option will be instead
used.
Standard options
-info
display information messages.
-quiet
do not display information messages or progress status; alternatively,
this can be achieved by setting the MRTRIX_QUIET environment variable to a
non-empty string.
-debug
display debugging messages.
-force
force overwrite of output files (caution: using the same file as input and
output might cause unexpected behaviour).
-nthreads number
use this number of threads in multi-threaded applications (set to 0 to
disable multi-threading).
-config key value (multiple uses permitted)
temporarily set the value of an MRtrix config file entry.
-help
display this information page and exit.
-version
display version information and exit.
tensor2metricのヘルプは、次の通り。
クリックして展開
SYNOPSIS
Generate maps of tensor-derived parameters
USAGE
tensor2metric [ options ] tensor
tensor the input tensor image.
OPTIONS
-adc image
compute the mean apparent diffusion coefficient (ADC) of the diffusion
tensor. (sometimes also referred to as the mean diffusivity (MD))
-fa image
compute the fractional anisotropy (FA) of the diffusion tensor.
-ad image
compute the axial diffusivity (AD) of the diffusion tensor. (equivalent to
the principal eigenvalue)
-rd image
compute the radial diffusivity (RD) of the diffusion tensor. (equivalent
to the mean of the two non-principal eigenvalues)
-cl image
compute the linearity metric of the diffusion tensor. (one of the three
Westin shape metrics)
-cp image
compute the planarity metric of the diffusion tensor. (one of the three
Westin shape metrics)
-cs image
compute the sphericity metric of the diffusion tensor. (one of the three
Westin shape metrics)
-value image
compute the selected eigenvalue(s) of the diffusion tensor.
-vector image
compute the selected eigenvector(s) of the diffusion tensor.
-num sequence
specify the desired eigenvalue/eigenvector(s). Note that several
eigenvalues can be specified as a number sequence. For example, '1,3'
specifies the principal (1) and minor (3) eigenvalues/eigenvectors
(default = 1).
-modulate choice
specify how to modulate the magnitude of the eigenvectors. Valid choices
are: none, FA, eigval (default = FA).
-mask image
only perform computation within the specified binary brain mask image.
Standard options
-info
display information messages.
-quiet
do not display information messages or progress status; alternatively,
this can be achieved by setting the MRTRIX_QUIET environment variable to a
non-empty string.
-debug
display debugging messages.
-force
force overwrite of output files (caution: using the same file as input and
output might cause unexpected behaviour).
-nthreads number
use this number of threads in multi-threaded applications (set to 0 to
disable multi-threading).
-config key value (multiple uses permitted)
temporarily set the value of an MRtrix config file entry.
-help
display this information page and exit.
-version
display version information and exit.
DTI拡散定量値(FA, MD, AD, RD, カラーFA)を計算するための基本的な使い方は、以下の通り。
dwi2tensor <入力画像> <出力画像> tensor2metric -fa <出力画像> -adc <出力画像> -ad <出力画像> -rd <出力画像> -vec <出力画像> tensor.mif
3. 使用例
3.1. 前準備
まず、次のファイルを用意する。
. ├── DWI.nii.gz # 拡散MRI ├── DWI_mask.nii.gz ├── bvals # b-values ├── bvecs # b-vectors └── headers.json # ヘッダー情報の入ったJSONファイル
こちらの記事を参考に、拡散MRI(DWI.nii.gz)とそのMPG軸情報(bvecs, bvals)とヘッダー情報(headers.json)をまとめて、MIF形式(DWI.mif)に変換する。
mrconvert -fslgrad bvecs bvals -json_import headers.json DWI.nii.gz DWI.mif
3.2. テンソルの推定(コマンド:dwi2tensor)
ファイルの用意ができたら、dwi2tensorを次のように実行する
dwi2tensor DWI.mif tensor.mif
mrinfoを使って「tensor.mif」の情報を確認すると、6 volumesのデータであることが分かる。
mrinfo tensor.mif
************************************************
Image name: "tensor.mif"
************************************************
Dimensions: 130 x 130 x 82 x 6
Voxel size: 1.76923 x 1.76923 x 1.8 x 1
Data strides: [ -1 2 3 4 ]
Format: MRtrix
Data type: 32 bit float (little endian)
Intensity scaling: offset = 0, multiplier = 1
Transform: 1 0 0 -109
-0 1 0 -103.7
-0 0 1 -58.57
それぞれのボリュームは、各方向の拡散係数に相当する。
The tensor coefficients are stored in the output image as follows:
volumes 0-5: D11, D22, D33, D12, D13, D23
3.3. 拡散定量値の算出(コマンド:tensor2metric)
先程推定した、「tensor.mif」を使って拡散定量値を算出する。
tensor2metric -fa FA.mif -adc MD.mif -ad AD.mif -rd RD.mif -vec color_FA.mif tensor.mif
DTIの各拡散定量値画像は、以下。

【FSL】FSLを用いた拡散テンソルイメージング: DTI
1. 目的
- FSLを用いた拡散テンソルイメージング: DTI
2. コマンド
FSLを用いて、拡散テンソルイメージング(DTI)をするには、dtifitコマンドを用いる。
dtifitのヘルプは、次の通り。
クリックして展開
Usage: dtifit -k <filename> dtifit --verbose Compulsory arguments (You MUST set one or more of): -k,--data dti data file -o,--out Output basename -m,--mask Bet binary mask file -r,--bvecs b vectors file -b,--bvals b values file Optional arguments (You may optionally specify one or more of): -V,--verbose switch on diagnostic messages -h,--help display this message --cni Input confound regressors --sse Output sum of squared errors -w,--wls Fit the tensor with weighted least squares --kurt Output mean kurtosis map (for multi-shell data) --kurtdir Output parallel/perpendicular kurtosis maps (for multi-shell data) --littlebit Only process small area of brain --save_tensor Save the elements of the tensor -z,--zmin min z -Z,--zmax max z -y,--ymin min y -Y,--ymax max y -x,--xmin min x -X,--xmax max x --gradnonlin Gradient Nonlinearity Tensor file
基本的な使い方は、以下の通り。
dtifit -k <DWI images> -o <OUTPUT> -m <MASK> -r <b-vectors file> -b <b-values file>
3. 使用例
まず、次のファイルを用意する。
. ├── DWI.nii.gz # 拡散MRI ├── DWI_mask.nii.gz # 拡散MRIのマスク画像 ├── bvals # b-values └── bvecs # b-vectors
ファイルの用意ができたら、dtifitを次のように実行する
dtifit -k DWI.nii.gz -o output -m DWI_mask.nii.gz -r bvecs -b bvals
処理が終わると次のファイルが生成される。
- output_V1 – 1st eigenvector (固有ベクトル)
- output_V2 – 2nd eigenvector
- output_V3 – 3rd eigenvector
- output_L1 – 1st eigenvalue (固有値)
- output_L2 – 2nd eigenvalue
- output_L3 – 3rd eigenvalue
- output_MD – mean diffusivity
- output_FA – fractional anisotropy (isotropic ~ 0; stick-like ~1)
- output_MO – mode of the anisotropy (oblate ~ -1; isotropic ~ 0; prolate ~ 1)
- output_S0 – raw T2 signal with no diffusion weighting
FA, MD, L1の画像は、以下。

【MRtrix】MRtrixを用いた拡散MRIの前処理 ~歪み・頭の動き・渦電流の補正~
1. 目的
2. コマンド
3. 使用例
3.1. 前準備
3.2. 歪み補正と頭の動き補正
1. 目的
- MRtrixを用いた拡散MRIの前処理(歪み・頭の動き・渦電流の補正)
2. コマンド
MRtrixを用いて、拡散MRIの歪み・頭の動き・渦電流を補正するには、dwifslpreprocを用いる(古いMRtrixバージョンではdwipreproc)。
dwipreprocは、FSLのtopupとeddyを用いるので、前もってFSLをインストールしておく必要がある。
dwifslpreprocのヘルプは、次の通り。
クリックして展開
SYNOPSIS
Perform diffusion image pre-processing using FSL's eddy tool; including
inhomogeneity distortion correction using FSL's topup tool if possible
USAGE
dwifslpreproc [ options ] input output
input The input DWI series to be corrected
output The output corrected image series
DESCRIPTION
This script is intended to provide convenience of use of the FSL software
tools topup and eddy for performing DWI pre-processing, by encapsulating
some of the surrounding image data and metadata processing steps. It is
intended to simply these processing steps for most commonly-used DWI
acquisition strategies, whilst also providing support for some more exotic
acquisitions. The "example usage" section demonstrates the ways in which
the script can be used based on the (compulsory) -rpe_* command-line
options.
The "-topup_options" and "-eddy_options" command-line options allow the
user to pass desired command-line options directly to the FSL commands
topup and eddy. The available options for those commands may vary between
versions of FSL; users can interrogate such by querying the help pages of
the installed software, and/or the FSL online documentation: (topup)
https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup/TopupUsersGuide ; (eddy)
https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy/UsersGuide
The script will attempt to run the CUDA version of eddy; if this does not
succeed for any reason, or is not present on the system, the CPU version
will be attempted instead. By default, the CUDA eddy binary found that
indicates compilation against the most recent version of CUDA will be
attempted; this can be over-ridden by providing a soft-link "eddy_cuda"
within your path that links to the binary you wish to be executed.
Note that this script does not perform any explicit registration between
images provided to topup via the -se_epi option, and the DWI volumes
provided to eddy. In some instances (motion between acquisitions) this can
result in erroneous application of the inhomogeneity field during
distortion correction. Use of the -align_seepi option is advocated in this
scenario, which ensures that the first volume in the series provided to
eddy is also the first volume in the series provided to eddy, guaranteeing
alignment. But a prerequisite for this approach is that the image contrast
within the images provided to the -se_epi option must match the b=0 volumes
present within the input DWI series: this means equivalent TE, TR and flip
angle (note that differences in multi-band factors between two acquisitions
may lead to differences in TR).
EXAMPLE USAGES
A basic DWI acquisition, where all image volumes are acquired in a single
protocol with fixed phase encoding:
$ dwifslpreproc DWI_in.mif DWI_out.mif -rpe_none -pe_dir ap -readout_time 0.55
Due to use of a single fixed phase encoding, no EPI distortion correction
can be applied in this case.
DWIs all acquired with a single fixed phase encoding; but additionally a
pair of b=0 images with reversed phase encoding to estimate the
inhomogeneity field:
$ mrcat b0_ap.mif b0_pa.mif b0_pair.mif -axis 3; dwifslpreproc DWI_in.mif DWI_out.mif -rpe_pair -se_epi b0_pair.mif -pe_dir ap -readout_time 0.72 -align_seepi
Here the two individual b=0 volumes are concatenated into a single 4D image
series, and this is provided to the script via the -se_epi option. Note
that with the -rpe_pair option used here, which indicates that the SE-EPI
image series contains one or more pairs of b=0 images with reversed phase
encoding, the FIRST HALF of the volumes in the SE-EPI series must possess
the same phase encoding as the input DWI series, while the second half are
assumed to contain the opposite phase encoding direction but identical
total readout time. Use of the -align_seepi option is advocated as long as
its use is valid (more information in the Description section).
All DWI directions & b-values are acquired twice, with the phase encoding
direction of the second acquisition protocol being reversed with respect to
the first:
$ mrcat DWI_lr.mif DWI_rl.mif DWI_all.mif -axis 3; dwifslpreproc DWI_all.mif DWI_out.mif -rpe_all -pe_dir lr -readout_time 0.66
Here the two acquisition protocols are concatenated into a single DWI
series containing all acquired volumes. The direction indicated via the
-pe_dir option should be the direction of phase encoding used in
acquisition of the FIRST HALF of volumes in the input DWI series; ie. the
first of the two files that was provided to the mrcat command. In this
usage scenario, the output DWI series will contain the same number of image
volumes as ONE of the acquired DWI series (ie. half of the number in the
concatenated series); this is because the script will identify pairs of
volumes that possess the same diffusion sensitisation but reversed phase
encoding, and perform explicit recombination of those volume pairs in such
a way that image contrast in regions of inhomogeneity is determined from
the stretched rather than the compressed image.
Any acquisition scheme that does not fall into one of the example usages
above:
$ mrcat DWI_*.mif DWI_all.mif -axis 3; mrcat b0_*.mif b0_all.mif -axis 3; dwifslpreproc DWI_all.mif DWI_out.mif -rpe_header -se_epi b0_all.mif -align_seepi
With this usage, the relevant phase encoding information is determined
entirely based on the contents of the relevant image headers, and
dwifslpreproc prepares all metadata for the executed FSL commands
accordingly. This can therefore be used if the particular DWI acquisition
strategy used does not correspond to one of the simple examples as
described in the prior examples. This usage is predicated on the headers of
the input files containing appropriately-named key-value fields such that
MRtrix3 tools identify them as such. In some cases, conversion from DICOM
using MRtrix3 commands will automatically extract and embed this
information; however this is not true for all scanner vendors and/or
software versions. In the latter case it may be possible to manually
provide these metadata; either using the -json_import command-line option
of dwifslpreproc, or the -json_import or one of the -import_pe_* command-
line options of MRtrix3's mrconvert command (and saving in .mif format)
prior to running dwifslpreproc.
OPTIONS
-pe_dir PE
Manually specify the phase encoding direction of the input series; can be a
signed axis number (e.g. -0, 1, +2), an axis designator (e.g. RL, PA, IS),
or NIfTI axis codes (e.g. i-, j, k)
-readout_time time
Manually specify the total readout time of the input series (in seconds)
-se_epi image
Provide an additional image series consisting of spin-echo EPI images,
which is to be used exclusively by topup for estimating the inhomogeneity
field (i.e. it will not form part of the output image series)
-align_seepi
Achieve alignment between the SE-EPI images used for inhomogeneity field
estimation, and the DWIs (more information in Description section)
-json_import file
Import image header information from an associated JSON file (may be
necessary to determine phase encoding information)
-topup_options " TopupOptions"
Manually provide additional command-line options to the topup command
(provide a string within quotation marks that contains at least one space,
even if only passing a single command-line option to topup)
-eddy_options " EddyOptions"
Manually provide additional command-line options to the eddy command
(provide a string within quotation marks that contains at least one space,
even if only passing a single command-line option to eddy)
-eddy_mask image
Provide a processing mask to use for eddy, instead of having dwifslpreproc
generate one internally using dwi2mask
-eddy_slspec file
Provide a file containing slice groupings for eddy's slice-to-volume
registration
-eddyqc_text directory
Copy the various text-based statistical outputs generated by eddy, and the
output of eddy_qc (if installed), into an output directory
-eddyqc_all directory
Copy ALL outputs generated by eddy (including images), and the output of
eddy_qc (if installed), into an output directory
Options for specifying the acquisition phase-encoding design; note that one of the -rpe_* options MUST be provided
-rpe_none
Specify that no reversed phase-encoding image data is being provided; eddy
will perform eddy current and motion correction only
-rpe_pair
Specify that a set of images (typically b=0 volumes) will be provided for
use in inhomogeneity field estimation only (using the -se_epi option)
-rpe_all
Specify that ALL DWIs have been acquired with opposing phase-encoding
-rpe_header
Specify that the phase-encoding information can be found in the image
header(s), and that this is the information that the script should use
Options for importing the diffusion gradient table
-grad GRAD
Provide the diffusion gradient table in MRtrix format
-fslgrad bvecs bvals
Provide the diffusion gradient table in FSL bvecs/bvals format
Options for exporting the diffusion gradient table
-export_grad_mrtrix grad
Export the final gradient table in MRtrix format
-export_grad_fsl bvecs bvals
Export the final gradient table in FSL bvecs/bvals format
Additional standard options for Python scripts
-nocleanup
do not delete intermediate files during script execution, and do not delete
scratch directory at script completion.
-scratch /path/to/scratch/
manually specify the path in which to generate the scratch directory.
-continue <ScratchDir> <LastFile>
continue the script from a previous execution; must provide the scratch
directory path, and the name of the last successfully-generated file.
Standard options
-info
display information messages.
-quiet
do not display information messages or progress status. Alternatively, this
can be achieved by setting the MRTRIX_QUIET environment variable to a non-
empty string.
-debug
display debugging messages.
-force
force overwrite of output files.
-nthreads number
use this number of threads in multi-threaded applications (set to 0 to
disable multi-threading).
-config key value (multiple uses permitted)
temporarily set the value of an MRtrix config file entry.
-help
display this information page and exit.
-version
display version information and exit.
基本的な使い方は、以下の通り。
dwifslpreproc <入力画像> <出力画像> [オプション]
3. 使用例
3.1. 前準備
位相エンコード方向を、APとPAそれぞれで撮像したデータがあったとする。DICOM形式からNIfTI形式に変換する方法は、以下の記事を参考にするとよい。
. ├── DWI_AP.nii.gz # DW images (PE: AP) ├── DWI_PA.nii.gz # DW images (PE: PA) ├── bvals_AP # b-values (PE: AP) ├── bvals_PA # b-values (PE: PA) ├── bvecs_AP # b-vectors (PE: AP) ├── bvecs_PA # b-vectors (PE: PA) ├── headers_AP.json # DICOM headers (PE: AP) └── headers_PA.json # DICOM headers (PE: PA)
まず、こちらの記事を参考に、拡散MRI(DWI.nii.gz)とそのMPG軸情報(bvecs, bvals)とヘッダー情報(headers.json)をまとめて、MIF形式(DWI.mif)に変換する。
mrconvert -fslgrad bvecs_AP bvals_AP -json_import headers_AP.json DWI_AP.nii.gz DWI_AP.mif # PE: AP mrconvert -fslgrad bvecs_PA bvals_PA -json_import headers_PA.json DWI_PA.nii.gz DWI_PA.mif # PE: PA
次に、mrcatを使ってDWI_AP.mifとDWI_PA.mifをひとつの画像(DWI_all.mif)にまとめる。
オプションの-axis 3は、4次元目のt軸(Volume)方向にまとめるという意味である(MRtrixではAxisを0から数える [i.e., x: 0, y: 1, z: 2, t: 3])。mrcatの詳細は、こちら。
mrcat DWI_AP.mif DWI_PA.mif DWI_all.mif -axis 3
次に、dwiextractを用いて、b=0のみを抽出する。dwiextractの詳細は、こちら。
dwiextract -bzero DWI_all.mif DWI_b0.mif
3.2. 歪み補正と頭の動き補正
歪み補正と頭の動き補正をするために、次のコマンドを実行する。
ここで使用した、各オプションは以下。
-rpe_header:位相エンコード情報を読み込む-se_epi:b=0(spin-echo EPI images)を指定-align_seepi:磁場の不均一性場の推定で用いられる、SE-EPI画像とDWIの間の位置合わせを実行
dwifslpreproc DWI_all.mif DWI_preproc.mif -rpe_header -se_epi DWI_b0.mif -align_seepi
歪み補正後の画像は、以下。

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


【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. サンプルコード










