MEDIC susceptibility distortion correction in MicaPipe

[PaulBautin]

What is the current behaviour?

• fMRI susceptibility distortion correction (SDC) is currently run with FSL-TOPUP which estimates a single static B0 non-uniformity map by using the symmetry of the distortion when the PE polarity is reversed.

What is the motivation / use case for changing the behaviour?

• B0 non-uniformity map is not necesseraly static and can be dynamic -- for example with head motion. MEDIC (Multi-Echo DIstortion Correction) can correct B0 non-uniformity distortions for each fMRI frame using the multi-echo phase information.

References:

• Original Paper: https://doi.org/10.1101%2F2023.11.28.568744
• GitHub Repository: https://github.com/vanandrew/warpkit

Before running:

• Create singularity container for warpkit: singularity build warpkit.simg docker://ghcr.io/vanandrew/warpkit:latest

Code (example for first 5 frames):

magnitude_data=( \
    "/local_raid/data/pbautin/data/pilot_dataset/sub-Pilot014/ses-01/func/sub-Pilot014_ses-01_task-rest_run-2_echo-1_bold.nii.gz" \
    "/local_raid/data/pbautin/data/pilot_dataset/sub-Pilot014/ses-01/func/sub-Pilot014_ses-01_task-rest_run-2_echo-2_bold.nii.gz" \
    "/local_raid/data/pbautin/data/pilot_dataset/sub-Pilot014/ses-01/func/sub-Pilot014_ses-01_task-rest_run-2_echo-3_bold.nii.gz" \
)

phase_data=( \
    "/local_raid/data/pbautin/data/pilot_dataset/sub-Pilot014/ses-01/func/sub-Pilot014_ses-01_task-rest_run-1_echo-1_part-phase_bold.nii.gz" \
    "/local_raid/data/pbautin/data/pilot_dataset/sub-Pilot014/ses-01/func/sub-Pilot014_ses-01_task-rest_run-1_echo-2_part-phase_bold.nii.gz" \
    "/local_raid/data/pbautin/data/pilot_dataset/sub-Pilot014/ses-01/func/sub-Pilot014_ses-01_task-rest_run-1_echo-3_part-phase_bold.nii.gz" \
)

metadata=( \
    "/local_raid/data/pbautin/data/pilot_dataset/sub-Pilot014/ses-01/func/sub-Pilot014_ses-01_task-rest_run-2_echo-1_bold.json" \
    "/local_raid/data/pbautin/data/pilot_dataset/sub-Pilot014/ses-01/func/sub-Pilot014_ses-01_task-rest_run-2_echo-2_bold.json" \
    "/local_raid/data/pbautin/data/pilot_dataset/sub-Pilot014/ses-01/func/sub-Pilot014_ses-01_task-rest_run-2_echo-3_bold.json" \
)

warpkit_img='/local_raid/data/pbautin/container/warpkit.simg'
out_prefix='/local_raid/data/pbautin/results/medic/medic_'


singularity exec --writable-tmpfs --containall  \
    -B /local_raid/data/pbautin/data/pilot_dataset/sub-Pilot014 \
    -B /local_raid/data/pbautin/results/medic \
    -B /home/bic/pbautin/.julia \
    $warpkit_img \
    --magnitude ${magnitude_data[0]} ${magnitude_data[1]} ${magnitude_data[2]} \
    --phase ${phase_data[0]} ${phase_data[1]} ${phase_data[2]} \
    --metadata ${metadata[0]} ${metadata[1]} ${metadata[2]} \
    --out_prefix $out_prefix -n 10

# Loop over the frames
for frame in {0..5}
do
  singularity exec --writable-tmpfs --containall  \
      -B /local_raid/data/pbautin/data/pilot_dataset/sub-Pilot014 \
      -B /local_raid/data/pbautin/results/medic \
      -B /home/bic/pbautin/.julia \
      $warpkit_img \
      extract_field_from_maps \
      ${out_prefix}_displacementmaps.nii \
      ${out_prefix}_displacementfield_${frame}.nii.gz \
      --frame_number ${frame} \
      --phase_encoding_axis j \
      --format fsl
  # Extract the specific frame using fslroi
  fslroi ${magnitude_data[0]} ${out_prefix}_${frame}_img.nii.gz 0 -1 0 -1 0 -1 ${frame} 1
  applywarp -i ${out_prefix}_${frame}_img.nii.gz -r ${out_prefix}_${frame}_img.nii.gz -w ${out_prefix}_displacementfield_${frame}.nii.gz -
done

# Concatenate all the warped frames into a 4D NIfTI image
fslmerge -t ${out_prefix}warped_img_4D.nii.gz ${out_prefix}_warped_{0..5}.nii.gz