EPI distortion correction形变矫正, eddy, fieldmap等五种不同方法

在快速采集功能核磁共振（fMRI，functional MRI）和弥散核磁共振（dMRI，diffusion MRI）数据时，EPI成像方式使得数据中存在很大的畸变。这种畸变使得额叶、颞叶位置和真实大脑差异巨大，因而无法进行后续的跨模态配准。例如，将dMRI数据中的b0配准到T1项。
因此，在传统的数据预处理流程中，无论是基于PANDA的临床低精度数据预处理流程，还是基于HCP的HARDI高精度数据预处理流程，形变矫正都是一个关键的步骤，决定着后续的数据分析任务的可靠度。

Fig.1 Distortion correction (eddy) & Coregistration steps

在本文档，将主要介绍EPI数据，也就是fMRI和dMRI的数据的形变矫正（Distortion correction）五种不同的方法。采用这五种方法，基本上可以处理当前在临床和科研场景中绝大多数的数据。本文档将分散于网络各处的文字汇总，并抽取关键的步骤，使得新手可以快速参考上手。

主要使用的软件有ANTs，FSL，MRtrix3

1 topup + eddy

topup+eddy是HCP处理流程中的一个标准步骤。在采集dMRI数据时，需要同时采集两个phase-encoding方向的数据。如Fig.1中的标为AP和PA都图片所示，可以明显发现，AP方向的数据额叶位置被拉长，而PA方向的数据额叶位置被压进去，表现为一个凹陷。

echo "Topuping  $1 data ..."
printf "0 -1 0 0.07176\n0 1 0 0.07176" > $dwi_dir/acqparams.txt
topup --imain=$dwi_dir/ap_pa_b0_mean.nii.gz --datain=$dwi_dir/acqparams.txt --config=b02b0.cnf --out=$dwi_dir/topup_results --iout=$dwi_dir/b0

利用AP和PA数据，topup可以估计出dMRI数据的distortion形变场，为–out的输出，以topup_results为前缀的fieldmap。在上述代码中，0 0.07176是total readout time，它是通过Echo_spacing和EPI-factor计算得到的，具体计算公式为：
然后eddy利用topup生成的fieldmap，作为参数–topup，就可以将被拉长或者凹进去的数据恢复成正常的大脑形状，变得更加圆润和饱满。在下面的代码中，使用的是cuda版本的eddy，能显著加速计算过程。结果如Fig.1所示。

echo "Eddy  $1 data ..."

b_num=`cat bval | wc -w`        
indx=""
for ((iii=1; iii<=$b_num; iii+=1)); 
do 
   indx="$indx 1"; 
done
echo $indx > $dwi_dir/index.txt

eddy_cuda --imain=$dwi_dir/ap.nii.gz --mask=$dwi_dir/mask.nii.gz --acqp=$dwi_dir/acqparams.txt --index=$dwi_dir/index.txt --bvecs=bvec --bvals=bval --topup=$dwi_dir/topup_results --out=$dwi_dir/dwi --repol --data_is_shelled #  --flm=quadratic --mporder=6 --slspec=slspec.txt --s2v_niter=5 --s2v_lambda=1 --s2v_interp=trilinear

在上面生成acqparams.txt的代码中，还需要注意的是0 -1 0中的-1，通常它代表的是在AP方向采集到的数据，而1代表的是PA方向的数据。如果采集到的是其它方向的数据建议阅读FSL的官方文档。
扩展阅读
FSL： topup
FSL：topup-userguide
FSL：eddy
FSL：eddy-userguide

2 fieldmap + eddy

除了可以使用第一种topup方法，eddy还可以使用采集到的fieldmap，做distortion correction。在使用fieldmap前，需要对fieldmap做预处理。在这里，除了dMRI数据以外，还需要同时采集到magnitude数据（mag.nii.gz）和phase数据（phase.nii.gz）。

2.1 对mag做去脑壳

mrconvert mag.nii.gz -coord 3 0 -axes 0,1,2 mag_1volume.nii.gz
bet mag_1volume.nii.gz mag_1volume_bet.nii.gz -R -f 0.2 -m

注：

1. 在此，mrconvert是MRtrix3库函数，本行代码含义为提取第一个volume。在此步骤，使用fslmaths也可以实现相同功能
2. 在去脑壳前，mag.nii.gz有两个volume，如果有较大的运动，则需要对mag.nii.gz的两个volumes做一次刚性配准，对齐大脑。

2.2 基于去过脑壳的mag_1volume_bet.nii.gz数据，对fieldmap进行预处理

fsl_prepare_fieldmap SIEMENS phase.nii.gz mag_1volume_bet.nii.gz fieldmap.nii.gz 7.38 
fslmaths fieldmap.nii.gz -div 6.28318 fieldmap_hz.nii.gz -odt float

注：
1.fsl_prepare_fieldmap仅支持SIEMENS采集的数据，单位为rad/s。其他机器型号的处理方式详见官方文档（https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FUGUE/Guide）
2. 7.38是回拨时间差（echo time difference），从操作员可以获得此项数据
3. 在eddy中，供使用的fieldmap单位为Hz，需要根据情况进行修正：fsl_prepare_fieldmap生成的fieldmap单位为rad/s，需转换成Hz。公式如下：1 rad/s= 1 /2π Hz

3 fieldmap + fugue

3.1 fieldmap进行预处理
基本上和2.2类似，但由于在fugue中fieldmap单位为rad/s，不需要将其单位从rad/s转为Hz。
3.2 利用fugue完成校正

fugue -i hardi_PA.nii.gz --dwell=0.00069 --loadfmap=fieldmap -u fugue_test.nii.gz

注：
dwell是回拨间隔（echo spacing），且以秒为单位，一般数据记录以ms为单位。如果写错，则报如下错信息：
#ERROR:: dwell time should be in seconds but the value of 0.69is unusually large and maybe incorrectly specified in units of milliseconds.

4 fieldmap + epi_reg

4.1 fieldmap进行预处理
epi_reg 也需要rad/s为单位的，所以喝3.1，2.2类似
4.2 利用epi_reg完成校正

epi_reg --echospacing=0.00095 --wmseg=t1_wmseg.nii.gz --fmap=field_maps_rad.nii.gz --fmapmag=mag_single.nii.gz --fmapmagbrain=mag_brain_single.nii.gz --pedir=y  --epi=b0.nii.gz --t1=t1.nii.gz --t1brain=t1_brain.nii.gz --out=dwi_epi_reg
5. 非线性配准

5 非线性配准

To be continued.