Phase-encoding table for motion correction using the SHARD reconstruction approach

Dear all,

I was wondering if anyone could give an example of what the pe_table looks like in mrtrix format. I tried the -export_pe_table command but I get the following error: “no phase-encoding information found within image”. Similarly, the mrinfo command output doesn’t show any phase encoding information or total readout time (but that is probably a Siemens/Philips issue in the way the information is stored in the header I guess?). If we’re able to somehow retrieve this information from the raw data we can in theory create the able manually, right?
Specifically we would like to use the pe_table as input to the motion correction SHARD reconstruction software implemented in mrtrix ( So if anyone has had any experience with that, comments/suggestions are welcome.

Any help is much apppreciated!

Many thanks

Dear Maria,

Nice to hear that you are trying out our new motion correction module. It’s still a bit experimental and not at all well enough documented yet, so don’t hesitate to get in touch if you run into any issues.

About your question: the PE table lists, for every volume, the phase encoding axis and readout time in the same format as FSL Topup’s acqp.txt file (see their user guide for details). The difference with the FSL format, is that we do not use the indexing file, but repeat the corresponding entries directly. For example, an image with 6 volumes and alternating AP-PA encoding could have a PE table:

0 -1  0  0.058
0  1  0  0.058
0 -1  0  0.058
0  1  0  0.058
0 -1  0  0.058
0  1  0  0.058

whereas the FSL format would be

0 -1  0  0.058
0  1  0  0.058

1 2 1 2 1 2

Note that both formats are supported with the -import_pe_table (MRtrix) and -import_pe_eddy (FSL) options. If you import your images from DICOM, MRtrix will try very hard to find the PE information in the headers and store it in the converted .mif file. If you did not start with the DICOMs or if you use, e.g., .nii files, you will need to keep track of the PE table separately.

Note that SHARD motion correction requires the PE table to apply a B0 field map during reconstruction; it does not estimate the field itself. Therefore, you will likely need to run topup first on a subset of b=0 images with reverse phase encoding, and then pass the field to dwimotioncorrect. In this use case, the exact value of the readout time is not important, it just needs to be the same in both configurations.

I hope this is clear, but I would understand if it’s not… So get in touch if you need more help.


Hi Daan,

Thanks a lot for the tips. It’s very exciting to read and learn more about your new approach!
I had indeed already ran FSL’s topup to generate the B0 field map (using b0 images acquired in both AP and PA directions). I then created a PE table with 119 identical entries (to match the number of volumes in our dwi series) which seems to have worked fine as input for dwimotioncorrect. Does that sound about right?

A couple of additional questions that came to mind:

  1. We have multi-shell data from a large group of 5 year-old children (4 shells: 0, 700, 1000, 2000 in 7, 20, 32 and 60 directions respectively). So far I’ve tested it with all the defaults (only changed the mb factor to match our acquisition). Is it recommended in your opinion to change some of the parameters and if so, how can we make an informed choice? Are the appropriate lmax/rlmax values determined automatically based on the data (like in mrtrix)? And what about regularization values?

  2. I was wondering how the SHARD approach fits into the overall dwi pre-processing pipeline. SHARD takes care of the motion part, but what about the eddy/inhomogeneity distortion corrections that are usually embedded in the dwipreproc script? Would additional pre-processing steps be needed before we go on to estimate the response function and/or do tractography or any other analysis?

  3. To get an idea of subject motion, the script motionstats can be used, using as input the motion parameters and weights generated with dwimotioncorrect - right? Could you maybe provide some more information about the output of dwimotioncorrect? It’s not so straightforward for me at the moment what exactly these text files contain!

Many thanks in advance!

Hi Maria,

I’m very sorry about the late reply. Here it goes:

  1. For your data, set -lmax to 0,4,6,8, i.e., the maximum SH order per shell given the number of directions that were acquired. You can then play around with different values for -rlmax. The ‘r’ is for the ‘rank-reduced’ representation used internally in the registration. I typically use -rlmax 4,2,0 as a starting point. Lower values (e.g. -rlmax 2,2,0 or 2,0) can make the motion correction more robust, higher values can potentially improve the accuracy of the slice alignment.

  2. As it stands, the code does not currently handle eddy-current distortion. This is something I would like to add by extending the (currently rigid) registration. Susceptibility distortion is corrected for, given a B0 field map at the input, which you can get either from a calibration scan or with FSL Topup if you have reverse phase encoding data. The pre-processing pipeline we use is: denoising – degibbs – topup – dwimotioncorrect. The output at that stage is a 5-D image of multi-shell SH coefficients. You can then use mssh2amp to map these back into dMRI volumes that you can use with CSD etc.

  3. Yes, motionstats calculates a measure of the amount of motion during the scan. The script will output 3 numbers: a measure of the amount of translation, a measure of the amount of rotation, and a measure of the slice outlier prevalence. The motion metrics are proportional to the first derivative, i.e., to the change in subject position during the scan. The outlier metric is a simple ratio (between 0 and 1) of affected slices. Details are in the preprint.

I hope this helps,