【初心者向け】基礎&実践プログラミング

初心者がつまづきやすいところ、最短で実力が身につく方法をお伝えします。

【脳画像前処理】拡散MRIの前処理(歪み補正&頭の動き補正)

f:id:AIProgrammer:20201209202006p:plain

動かしながら学ぶ PyTorchプログラミング入門

動かしながら学ぶ PyTorchプログラミング入門

  • 作者:斎藤勇哉
  • 発売日: 2020/11/30
  • メディア: Kindle版

目次

目的

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

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

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

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

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

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

f:id:AIProgrammer:20201209164110p:plain

撮像中の頭の動き

f:id:AIProgrammer:20201209171525g:plain

使用するソフトウェア

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

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の出力ファイルの接頭辞
--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は、次の式で算出できる。


Total readout time = Echo spacing \times ( EPI factor ― 1)

例えば、--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

--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 
画像の歪み補正の効果

f:id:AIProgrammer:20201209165126g:plain

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の出力ファイルの接頭辞

頭の動き補正の効果

f:id:AIProgrammer:20201209171525g:plain f:id:AIProgrammer:20201209204316g:plain

サンプルコード

以下のコードは、前処理前の拡散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"


頑張れ!喝!!の代わりにB!ブックマークを押していただけるとただただうれしいです(^^)! ↓

動かしながら学ぶ PyTorchプログラミング入門

動かしながら学ぶ PyTorchプログラミング入門

  • 作者:斎藤勇哉
  • 発売日: 2020/11/30
  • メディア: Kindle版