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. 使用するソフトウェア
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"
著者情報： 斎藤 勇哉
順天堂大学医学部 大学院医学研究科 放射線診断学講座所属脳MRI 画像解析が専門であり、テーマは①神経変性疾患の機序解明、②医用人工知能の開発、③多施設データのハーモナイゼーション、④速読が脳に与える影響や学習効果、⑤SNS解析を用いたマーケティング戦略の改善。
医療分野に関わらず、自然言語処理・スクレイピング・データ分析・Web アプリ開発を得意とし、企業や他大学の研究を支援。
主な使用言語は、Python、Shell Script、MATLAB、HTML、CSS
いつもありがとうございます。詳細な解説ありがとうございます。
私はFSLを用いたProbabilistic tractographyの初心者です。
DWI元画像→NifTI変換→BET処理（脳マスク画像の生成）までは進めましたが次のEPI歪み補正についてつまづいております。脳マスクをFSLeyesで見る限りはすでに歪みがないように見えるのですが、その場合でもEPI歪み補正は必要なのでしょうか。
topup補正を実行する前の準備段階として、位相エンコード情報の抽出や読み出し時間の抽出などを行って必要なファイルを用意する必要があると思いますが、このあたりの処理がかなり難しい印象でした。（正直どのファイルからどのように抽出すれば正しいファイルが作成できるかがわかりません。）
もしEPI歪み補正をスキップできるならばしたいですが、肉眼的に歪みがなくてもtopupの処理は必要なのか教えていただきたいです。topup処理をスキップする場合には次にeddy current およびmotion 補正に移行すると思われます。よろしくお願いいたします。
ご質問ありがとうございます。
まず、確認したいのですが、DWIはB0画像をAP方向、PA方向と撮像されていらっしゃいますか？
それをしていなければ、そもそもTOPUPができません。そのままeddy currentに進むことになるかと思います。
もし、2方向撮像されているならば、TOPUPは是非挑戦していただきたいです。自分が思っている以上に脳が歪んでいたということを認識できます。
ただ、FSLだけだとしんどいと思いますので、MRtrix3を使ったほうがスムーズにできるかもしれません。
根本先生
お忙しい中ご返信ありがとうございます。ご指摘ありがとうございます。抽出の方法自体はわかったのですが、PhaseEncodingDirectionが”j-“を示すjsonファイルが見当たりませんでしたので、改めて元画像の撮像条件を確認しようと思います。
私自身もMRtrix3で前処理をした経験はありますが、eddy current補正に関しては、FSLのeddyよりもdwifslpreprocが優れているということでしょうか？補正処理の仕様自体が異なっているのでしょうか。
休日中にも関わらずご教示くださりありがとうございます。
そのjsonファイルは、dcm2niix を用いて DICOM→NIFTIに変換した時に作成されます。
ちなみに j は PA 方向、 j- はAP方向であることが多いと思います。
j 軸は、PA方向で、プラスは Aになるので、PA方向は、プラスにむかってすすむので、j, AP方向は、マイナスの方向に進むので j- となります。
お忙しい中ありがとうございます。cat コマンドで確認してみると”PhaseEncodingPolarityGE”: “Flipped”,と表示されていましたので上下反転しているのかと誤解しておりました。
現状は”j”はありますが、”j-“がないのでAP方向が存在しないということになるかと思われます。撮像条件を確認いたします。
LRおよびRL方向のファイルはありそうですが、それではEPI補正は不可能でしょうか？
度々ありがとうございます。
お忙しい中ありがとうございます。
catコマンドで見てみますと、”PhaseEncodingPolarityGE”: “Flipped”,の表示があり、反転していると誤解しておりました。
そうなると現状はAP方向のファイルが見つかっていないということになるかと思いますが、もう一度撮像条件を確認いたします。
LR及びRL方向のファイルはありますが、これらでEPI歪み補正を行うことは可能なのでしょうか。一般的ではないとは思われますが。
度々ありがとうございます。
LR方向およびRL方向のファイルがあるとすると、その2つを使って、歪み補正をすることになります。その場合はPA方向の画像はおそらくないでしょう。すみません、たいていはPAかAP方向で撮像されている場合が多いのでそうかなと思ったところでした。
なお、LR、RL方向であれば、phase encoding direction は i か i- となるのかなと思いますがご確認ください。
ピングバック： 【FSL】XTRACTを用いたトラクトグラフィー
先日のABiSではお世話になりありがとうございました。Lin4neuro20.04でのDWI前処理について質問です。dwidenoise, mrdegibbsの後のdwifslpreproc補正ですが、手元にある位相エンコードAP１方向のみのdata(dwi_den_unr.mif)で下記コマンド,dwifslpreproc dwi_den_unr.mif dwi_den_unr_preproc.mif -pe_dir AP -rpe_none -eddy_options “–slm=linear” -readout_time 0.047を実行すると、-eddy_optionsのエラーが出ます。dwifslpreproc -help を見ると、（中略） 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.と出ます。
根本先生の「Ubuntu 20.04/18.04環境でeddy_cuda10.2( in FSL 6.0.5.x),PyTorch, Tensorflow2を使えるようにCUDA10.2,11.0,11.5をセットアップする方法」を実行してみましたが、invida-smi確認時に、CUDA driberが無いことを指摘されます。使用PCがMacOSCatalina10.15.7のため、CUDAdriver for Macをインストールしても起動しません。
以上のような環境で、位相エンコードがAP１方向のdataでdwifslpreproc,eddy_optionsを実行する方法をご教示いただけましたら幸いです。宜しくお願いいたします。
eddy_options のエラーではなく、Warningではないでしょうか？
– 残念なことに、現在、VirtualBoxで、GPUは使えません。したがって、eddy_cudaは使えません。
– さらに残念なことに、AppleとNVIDIAは仲違いをしているので、nVIDIAはApple用のGPUドライバを提供していません。
つまり、現行のmacOSでは、eddy_cudaを使用することができません。
なので、eddy_cudaは使えず、helpにあるとおり、eddy_cudaが使えない場合は、eddy_openmpを使用します。
eddy_openmpは環境にもよりますが、1例4−7時間かかります。
上田先生のチュートリアルのコマンドはしっかり動くことを確認しています。
もし、Warningではなく、エラーでプログラムが止まっているようでしたら、そのエラーメッセージをそのまま貼り付けていただけませんか？
よろしくお願いします。