Corregistro de DTI usando epi_reg y ANTs
Tomado de aqui
¿Porqué epi_reg?
Cuando se registra el T1 a espacio nativo DTI usando solo ANTS los resultados no son muy buenos. Con epi_reg el registro resulta ser mucho mejor.
Entonces, ¿Porqué ANTS?
La matriz resultado de epi_reg se usa con flirt o con fnirt pero realmente los mejores resultados en loas atlas se alcanzarian utilizando antsApplyTransforms y utilizando la matriz resultado de epi_reg.
tl;dr
#!/bin/sh
Usage() {
echo ""
echo "Usage: dti_proc.sh <study> <file_name> <T1> <name> <output_dir>"
echo ""
echo "You must have FSL installed in order to run this script"
echo ""
exit 1
}
[ "$5" = "" ] && Usage
debug=1
study=$1
shift
in=`${FSLDIR}/bin/remove_ext $1`
shift
t1=`${FSLDIR}/bin/remove_ext $1`
shift
b_in=$1
shift
out=`${FSLDIR}/bin/remove_ext $1`
shift
td=${out}'/.tmp_'${b_in}
p2a=$(dirname ${in}.nii.gz)'/'${b_in}'_p2a_b0'
if [ ! -d "$td" ]; then
mkdir $td
fi
if [ ! -d "$out" ]; then
mkdir $out
fi
echo "DTI preproccessing begins on ${b_in} ..."
echo [`date`]
echo
${FSLDIR}/bin/fslroi $in ${td}/a2p_b0 0 1
if [[ -f ${p2a}.nii.gz ]]; then
echo "Preparing topup"
${FSLDIR}/bin/imcp ${p2a} ${td}/
${FSLDIR}/bin/fslmerge -t ${td}/a2p_p2a_b0 ${td}/a2p_b0 ${td}/${b_in}_p2a_b0
echo "Running topup"
${FSLDIR}/bin/topup --imain=${td}/a2p_p2a_b0 --datain=${out}/../acqparams.txt --out=${td}/topup_results --iout=${td}/hifi_b0
${FSLDIR}/bin/fslmaths ${td}/hifi_b0 -Tmean ${td}/hifi_b0
else
echo "No P>>A so will continue uncorrected"
${FSLDIR}/bin/fslmaths ${td}/a2p_b0 -Tmean ${td}/hifi_b0
echo ${in} >> ${out}/../dti_uncorrected_subjects.txt
fi
${FSLDIR}/bin/bet ${td}/hifi_b0 ${td}/hifi_b0_brain -m
echo "Copying files for ${b_in} to ${td}/"
echo
cp ${in}.bval ${td}/bvals
cp ${in}.bvec ${td}/bvecs
echo "Doing correction on ${td}/${b_in}_hifi"
echo [`date`]
echo
${FSLDIR}/bin/eddy_cuda --imain=${in} --mask=${td}/hifi_b0_brain --acqp=${out}/../acqparams.txt --index=${out}/../dti_index.txt --bvecs=${td}/bvecs --bvals=${td}/bvals --topup=${td}/topup_results --out=${td}/data
echo "Running dtifit on ${td}/data"
echo [`date`]
echo
${FSLDIR}/bin/dtifit --data=${td}/data --out=${td}/dti --mask=${td}/hifi_b0_brain --bvecs=${td}/bvecs --bvals=${td}/bvals
echo "I will copy all output files to ${out}/${b_in}_XXXXX"
echo
for x in ${td}/dti*; do ${FSLDIR}/bin/imcp ${x} ${out}/${b_in}_$(basename $x); done;
${FSLDIR}/bin/imcp ${td}/hifi_b0_brain ${out}/${b_in}_dti_brain_mask
${FSLDIR}/bin/imcp ${td}/data ${out}/${b_in}_dti_data
echo "I need the T1 image"
${FSLDIR}/bin/fslreorient2std ${t1} ${out}/${b_in}_t1_reoriented
echo "I get the extracted brain"
${FREESURFER_HOME}/bin/mri_vol2vol --mov ${SUBJECTS_DIR}/${study}_${b_in}/mri/brain.mgz --targ ${SUBJECTS_DIR}/${study}_${b_in}/mri/rawavg.mgz --regheader --o ${td}/${b_in}_tmp_brain_in_rawavg.mgz
${FREESURFER_HOME}/bin/mri_convert --in_type mgz --out_type nii ${td}/${b_in}_tmp_brain_in_rawavg.mgz ${td}/${b_in}_tmp_brain.nii.gz
${FSLDIR}/bin/fslreorient2std ${td}/${b_in}_tmp_brain ${out}/${b_in}_brain
echo "Calculating transformation from MNI to T1"
echo [`date`]
echo
antsRegistrationSyN.sh -d 3 -f ${out}/${b_in}_t1_reoriented.nii.gz -m ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz -o ${td}/${b_in}_mni_t1_warp
antsApplyTransforms -d 3 -i ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz -r ${out}/${b_in}_t1_reoriented.nii.gz -t ${td}/${b_in}_mni_t1_warp1Warp.nii.gz -t ${td}/${b_in}_mni_t1_warp0GenericAffine.mat -o ${td}/${b_in}_mni_t1_warped.nii.gz
echo "Calculating transformation from T1 to B0"
echo [`date`]
echo
${FSLDIR}/bin/epi_reg --epi=${out}/${b_in}_dti_data --t1=${out}/${b_in}_t1_reoriented --t1brain=${out}/${b_in}_brain --out=${td}/${b_in}_tmp_diff2std
${FSLDIR}/bin/convert_xfm -omat ${td}/${b_in}_tmp_std2diff.mat -inverse ${td}/${b_in}_tmp_diff2std.mat
c3d_affine_tool -ref ${td}/hifi_b0.nii.gz -src ${out}/${b_in}_t1_reoriented.nii.gz ${td}/${b_in}_tmp_std2diff.mat -fsl2ras -oitk ${td}/${b_in}_epi_reg_ANTS.mat
antsApplyTransforms -d 3 -i ${out}/${b_in}_t1_reoriented.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_t1_b0.nii.gz
antsApplyTransforms -d 3 -i ${td}/${b_in}_mni_t1_warped.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_mni_to_b0.nii.gz
echo "Applying transformations to JHU Atlas (ICBM-labels)"
antsApplyTransforms -d 3 -i ${FSLDIR}/data/atlases/JHU/JHU-ICBM-labels-1mm.nii.gz -r ${out}/${b_in}_t1_reoriented.nii.gz -t ${td}/${b_in}_mni_t1_warp1Warp.nii.gz -t ${td}/${b_in}_mni_t1_warp0GenericAffine.mat -o ${td}/${b_in}_JHU_labels_tmp.nii.gz
antsApplyTransforms -d 3 -i ${td}/${b_in}_JHU_labels_tmp.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_JHU_labels.nii.gz
echo "Applying transformations to JHU Atlas (ICBM-tracts)"
antsApplyTransforms -d 3 -i ${FSLDIR}/data/atlases/JHU/JHU-ICBM-tracts-maxprob-thr25-1mm.nii.gz -r ${out}/${b_in}_t1_reoriented.nii.gz -t ${td}/${b_in}_mni_t1_warp1Warp.nii.gz -t ${td}/${b_in}_mni_t1_warp0GenericAffine.mat -o ${td}/${b_in}_JHU_tracts_tmp.nii.gz
antsApplyTransforms -d 3 -i ${td}/${b_in}_JHU_tracts_tmp.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_JHU_tracts.nii.gz
${FSLDIR}/bin/imcp ${td}/hifi_b0 ${out}/${b_in}_hifi_b0
echo "Game over :-P"
echo [`date`]
exit 0
Correcion con topup y preprocesamiento
Esto se hace del modo usual,
${FSLDIR}/bin/fslroi $in ${td}/a2p_b0 0 1
if [[ -f ${p2a}.nii.gz ]]; then
echo "Preparing topup"
${FSLDIR}/bin/imcp ${p2a} ${td}/
${FSLDIR}/bin/fslmerge -t ${td}/a2p_p2a_b0 ${td}/a2p_b0 ${td}/${b_in}_p2a_b0
echo "Running topup"
${FSLDIR}/bin/topup --imain=${td}/a2p_p2a_b0 --datain=${out}/../acqparams.txt --out=${td}/topup_results --iout=${td}/hifi_b0
${FSLDIR}/bin/fslmaths ${td}/hifi_b0 -Tmean ${td}/hifi_b0
else
echo "No P>>A so will continue uncorrected"
${FSLDIR}/bin/fslmaths ${td}/a2p_b0 -Tmean ${td}/hifi_b0
echo ${in} >> ${out}/../dti_uncorrected_subjects.txt
fi
${FSLDIR}/bin/bet ${td}/hifi_b0 ${td}/hifi_b0_brain -m
Solo que para correr topup se necesita un scan en direccion contraria (P»A). Si este no existe no se puede realizar la correcion. Para poder continuar se toma el B0 como si estuviera corregido pero se anota en un archivo por separado los que no estan corregidos. Por aquello de comprobar que ha pasado.
Luego viene la correcion con eddy. Para esto uso eddy_cuda que aunque me limita el numero de maquinas a usar es muchas veces mas rapido.
cp ${in}.bval ${td}/bvals
cp ${in}.bvec ${td}/bvecs
${FSLDIR}/bin/eddy_cuda --imain=${in} --mask=${td}/hifi_b0_brain --acqp=${out}/../acqparams.txt --index=${out}/../dti_index.txt --bvecs=${td}/bvecs --bvals=${td}/bvals --topup=${td}/topup_results --out=${td}/data
Ahora si, preprocesamiento con dtifit,
${FSLDIR}/bin/dtifit --data=${td}/data --out=${td}/dti --mask=${td}/hifi_b0_brain --bvecs=${td}/bvecs --bvals=${td}/bvals
Registro con epi_reg
Primero voy a tomar todos los archivos necesarios. Para empezar, los resultados del preprocesamiento
for x in ${td}/dti*; do ${FSLDIR}/bin/imcp ${x} ${out}/${b_in}_$(basename $x); done;
${FSLDIR}/bin/imcp ${td}/hifi_b0_brain ${out}/${b_in}_dti_brain_mask
${FSLDIR}/bin/imcp ${td}/data ${out}/${b_in}_dti_data
Ahora el T1,
${FSLDIR}/bin/fslreorient2std ${t1} ${out}/${b_in}_t1_reoriented
Tambien voy a necesitar el cerebro extraido pero el lugar de hacer un bet voy a tomar el que sale de Freesurfer
${FREESURFER_HOME}/bin/mri_vol2vol --mov ${SUBJECTS_DIR}/${study}_${b_in}/mri/brain.mgz --targ ${SUBJECTS_DIR}/${study}_${b_in}/mri/rawavg.mgz --regheader --o ${td}/${b_in}_tmp_brain_in_rawavg.mgz
${FREESURFER_HOME}/bin/mri_convert --in_type mgz --out_type nii ${td}/${b_in}_tmp_brain_in_rawavg.mgz ${td}/${b_in}_tmp_brain.nii.gz
${FSLDIR}/bin/fslreorient2std ${td}/${b_in}_tmp_brain ${out}/${b_in}_brain
Ahora, el template MNI se lleva a T1 con ANTs sin problemas,
antsRegistrationSyN.sh -d 3 -f ${out}/${b_in}_t1_reoriented.nii.gz -m ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz -o ${td}/${b_in}_mni_t1_warp
antsApplyTransforms -d 3 -i ${FSLDIR}/data/standard/MNI152_T1_1mm.nii.gz -r ${out}/${b_in}_t1_reoriented.nii.gz -t ${td}/${b_in}_mni_t1_warp1Warp.nii.gz -t ${td}/${b_in}_mni_t1_warp0GenericAffine.mat -o ${td}/${b_in}_mni_t1_warped.nii.gz
Pero para ir al espacio nativo DTI voy a usar epi_reg. Ojo que tendre que hacer el registro de T1 a DTI nativo, invertir l matriz y convertir la matriz inversa. Esta ya la puedo usar con ANTs.
${FSLDIR}/bin/epi_reg --epi=${out}/${b_in}_dti_data --t1=${out}/${b_in}_t1_reoriented --t1brain=${out}/${b_in}_brain --out=${td}/${b_in}_tmp_diff2std
${FSLDIR}/bin/convert_xfm -omat ${td}/${b_in}_tmp_std2diff.mat -inverse ${td}/${b_in}_tmp_diff2std.mat
c3d_affine_tool -ref ${td}/hifi_b0.nii.gz -src ${out}/${b_in}_t1_reoriented.nii.gz ${td}/${b_in}_tmp_std2diff.mat -fsl2ras -oitk ${td}/${b_in}_epi_reg_ANTS.mat
antsApplyTransforms -d 3 -i ${out}/${b_in}_t1_reoriented.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_t1_b0.nii.gz
antsApplyTransforms -d 3 -i ${td}/${b_in}_mni_t1_warped.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_mni_to_b0.nii.gz
Atlas
Ya solo queda aplicar los las transfomaciones calculadas a los atlas de JHU y copiar el B0 al directorio de trabajo para hacer el QC.
antsApplyTransforms -d 3 -i ${FSLDIR}/data/atlases/JHU/JHU-ICBM-labels-1mm.nii.gz -r ${out}/${b_in}_t1_reoriented.nii.gz -t ${td}/${b_in}_mni_t1_warp1Warp.nii.gz -t ${td}/${b_in}_mni_t1_warp0GenericAffine.mat -o ${td}/${b_in}_JHU_labels_tmp.nii.gz
antsApplyTransforms -d 3 -i ${td}/${b_in}_JHU_labels_tmp.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_JHU_labels.nii.gz
antsApplyTransforms -d 3 -i ${FSLDIR}/data/atlases/JHU/JHU-ICBM-tracts-maxprob-thr25-1mm.nii.gz -r ${out}/${b_in}_t1_reoriented.nii.gz -t ${td}/${b_in}_mni_t1_warp1Warp.nii.gz -t ${td}/${b_in}_mni_t1_warp0GenericAffine.mat -o ${td}/${b_in}_JHU_tracts_tmp.nii.gz
antsApplyTransforms -d 3 -i ${td}/${b_in}_JHU_tracts_tmp.nii.gz -r ${td}/hifi_b0.nii.gz -t ${td}/${b_in}_epi_reg_ANTS.mat -o ${out}/${b_in}_JHU_tracts.nii.gz
${FSLDIR}/bin/imcp ${td}/hifi_b0 ${out}/${b_in}_hifi_b0
QC
El Qc se realiza basicamente comprobando la superposicion del B0 con la plantilla MNI trasladada al espacio B0. Esta es, desde lejos, la mejor aproximacion al registro a DTI que he logrado.