Issue with initial preprocessing in BATMAN pipeline on Philips 3T DWI data (reverse phase-encoding)

aritban · December 2, 2024, 4:22pm

Hi MRtrix3 Community,

After successfully completing the BATMAN tutorial dataset and getting it to run smoothly, I am now working on processing my own DWI data collected on a Philips 3T scanner. I have already converted it to BIDS format and organized properly. My plan is to eventually run the pipeline on multiple subjects, but for now, I am starting with just one subject to work through any issues. My main challenge is at the initial preprocessing stage, specifically in handling reverse phase-encoding (RPE) corrections.

Data Overview

DWI Directory:

The dwi folder contains two acquisitions collected in the AP direction:

sub-22_dir-AP_acq-1_dwi.bval  
sub-22_dir-AP_acq-1_dwi.bvec  
sub-22_dir-AP_acq-1_dwi.nii.gz  
sub-22_dir-AP_acq-1_dwi.json  
sub-22_dir-AP_acq-2_dwi.bval  
sub-22_dir-AP_acq-2_dwi.bvec  
sub-22_dir-AP_acq-2_dwi.nii.gz  
sub-22_dir-AP_acq-2_dwi.json

Both acquisitions have the following dimensions:

mrinfo sub-22_dir-AP_acq-1_dwi.nii.gz
Dimensions:        144 x 144 x 81 x 51
Voxel size:        1.66667 x 1.66667 x 1.7 x 5.3

Fieldmap Directory:

The fmap folder contains AP and PA b=0 images:

sub-22_dir-AP_epi.nii.gz  
sub-22_dir-AP_epi.json  
sub-22_dir-PA_epi.nii.gz  
sub-22_dir-PA_epi.json

Both fieldmaps have the following dimensions:

mrinfo sub-22_dir-AP_epi.nii.gz
Dimensions:        144 x 144 x 81
Voxel size:        1.66667 x 1.66667 x 1.7

Preprocessing Steps

I have followed a systematic approach to prepare my data for dwifslpreproc:

1. DWI Directory:

Concatenated Gradient Tables:

paste -d' ' sub-22_dir-AP_acq-1_dwi.bval sub-22_dir-AP_acq-2_dwi.bval > concatenated.bval
paste -d' ' sub-22_dir-AP_acq-1_dwi.bvec sub-22_dir-AP_acq-2_dwi.bvec > concatenated.bvec

Attached Gradients to the DWI File:

mrconvert dwi_den_unr_even.mif dwi_den_unr_even_with_gradients.mif \
          -fslgrad concatenated.bvec concatenated.bval

Extracted b=0 Volumes and Computed Mean:

dwiextract dwi_den_unr_even_with_gradients.mif --bzero b0_AP_even.mif
mrmath b0_AP_even.mif mean mean_b0_AP_even.mif -axis 3

2. Fieldmap Directory:

Prepared AP-PA Fieldmaps:
- Converted NIfTI to MIF:

mrconvert sub-22_dir-PA_epi.nii.gz b0_PA_even.mif
mrconvert sub-22_dir-AP_epi.nii.gz b0_AP_even.mif

Computed Mean Images:

mrmath b0_PA_even.mif mean mean_b0_PA_even.mif -axis 3
mrmath b0_AP_even.mif mean mean_b0_AP_even.mif -axis 3

Concatenated AP and PA:

mrcat mean_b0_AP_even.mif mean_b0_PA_even.mif -axis 3 b0_pair.mif

Processing with `dwifslpreproc`

I attempted to process my data using dwifslpreproc:

singularity exec ~/bin/mrtrix3_latest.sif dwifslpreproc \
  dwi_den_even_resized_with_gradients.mif \
  sub-02_den_preproc.mif \
  -nocleanup \
  -pe_dir AP \
  -rpe_pair \
  -se_epi b0_pair.mif \
  -eddy_options " --slm=linear --data_is_shelled"

The Problem

While the command runs without crashing or giving me any errors, the warping correction doesn’t seem to work as expected. Instead of producing a properly corrected AP-PA combined image, the output still resembles the original AP image and it doesn’t seem like the PA image contributes to the distortion correction.

Additional Information

I have run the same pipeline successfully using the BATMAN tutorial dataset.
I have manually run topup on similar data previously, and it worked.
The b0_pair.mif file appears valid, and the dimensions match.

Question

Is there something specific about Philips DWI data (e.g., odd slice dimensions, intensity scaling, etc.) that could affect the performance of dwifslpreproc or topup during preprocessing?

I have been trying to ensure that the RPE correction works as expected, specifically that the PA information is correctly applied to correct the AP data. However, I am running into issues where the resulting warp does not seem to properly resolve the distortions in the AP images. Instead, the output looks more like the original AP image, as if the PA information wasn’t correctly incorporated.

Here’s what I’ve tried so far:

Tested both odd and even slice dimensions, ensuring the images align properly.
Processed the files in two ways: by concatenating acquisitions into one file and also by processing them individually.
Confirmed the phase-encoding directions in the metadata (AP and PA).
Resized images to ensure matching dimensions for AP and PA images before concatenation for b0_pair.mif.

Despite these steps, the correction doesn’t seem to work as expected. Are there any known issues with Philips data or preprocessing steps (e.g., intensity scaling, qform/sform discrepancies, or metadata handling) that I might be missing? How can I ensure that the reverse phase-encoding (RPE) correction is correctly applied?

Any insights, suggestions, or things to check would be greatly appreciated. Thank you so much!

jdtournier · December 4, 2024, 7:58pm

Hi @aritban,

I have to admit, I’m not sure what might be going on here – I can see you’ve already spent quite a bit of time and effort trying to debug this… I just have a few suggestions, though I don’t think they’re necessarily the issue here:

For the initial data import, I’d recommend using dwicat to merge the two runs, which requires running mrconvert -fslgrad on each dataset first:

for f in sub-22_dir-AP_acq-{1,2}_dwi; do 
  mrconvert $f.nii.gz -json_import $f.json -fslgrad $f.bvec $f.bval $f.mif; 
done
dwicat sub-22_dir-AP_acq-{1,2}_dwi.mif dwi.mif

I’m using a bit of bash substitution to ease the handling of all these long filenames, but hopefully you get the gist. This allows you to import both the JSON and DW gradient table information into each DWI run before merging them, and dwicat is designed to ensure appropriate intensity matching across the two runs.

For your fieldmaps, I’m a bit confused as to which one you’re using for the AP run. You provide instructions to extract them from the main DWI series (using dwiextract), but also by using the sub-22_dir-AP_epi.nii.gz file directly. In general, it’s preferable to ensure the first b=0 in your fieldmaps is the same image as the first image in your DWI series, since topup computes the field with respect to that image, and eddy then relies on the assumption that the field provided is aligned to the first volume in the series.

I know @rsmith has included a lot of smarts in dwifslpreproc to meet these requirements, but I think this is only possible when using the -rpe_header option. Since it’s not uncommon for Philips data to have the b=0 image as the last volume in the DWI series, it may be a good idea to use the -align_seepi option to dwifslpreproc – have a look at your data and check if this applies to you.

Note that this could be problematic when there is significant motion – I’ve no reason to think this is what is causing the issue here.

There should be no need to pad out the data to ensure an even number of slices, dwifslpreproc should take care of that for you. I recommend keeping things as simple and close to the raw data as possible, so we can rule out as many issues as possible outside of dwifslpreproc.

I’m not sure where your dwi_den_even_resized_with_gradients.mif input to dwifslpreproc is coming from, but hopefully it’s all as expected. I assume from what you’re saying you’ve already tried quite a few different inputs there… I’m also not sure where the singularity container you’re using comes from, since we don’t provide such containers – but hopefully it’ll be running the latest version of MRtrix and invoking a recent (hopefully latest) version of FSL.

Assuming this all checks out, it would be good to inspect the contents of the temporary folder created by dwifslpreproc (which should be reported on the terminal, and preserved since you’re using the -nocleanup option). You should be able to step through all the commands the script executed (recorded in one of the logs), and inspect all the temporary outputs generated, including the field map produced by topup. If you’re able to inspect that map, and ideally apply it to the b=0 AP image to check it is as expected, that might give you a clue as to where the problem is occurring. At the end of the day, dwifslpreproc should be running a similar sequence of commands to what you might do if you were running the process manually, so if you already have experience with that, it should be relatively straightforward to step through and identify the problematic step (at least I hope!).

Feel free to post screenshots or any other information from the temporary folder if none of this helps, maybe someone will spot an issue we haven’t thought of…

All the best,
Donald.