Dwifslpreproc eddy_mask.nii looks faulty

[HelenaV]

Dear MRtrix experts,

I am currently analysing multishell data (b=0, b=700, b=1200 and b=2800) and I just ran dwifslpreproc (3.0_RC4) on our HPC cluster (the people from HPC installed this version, I don’t really get how, but apparently, this is the newest of newest?). I wanted to check data quality of the output in the eddy output folder and the first thing I opened was the eddy_mask.nii that is created during the process. It looks like a whole part of the brain is cut off, but I have no idea why. As this mask is being used during eddy, I guess this is not good news for the analysis?

This is what the mask looks like (and it is like this for every subject):

For the analysis, I used the following command

dwifslpreproc ${input}/dwi_degibbsed.mif ${input}/geomcorr.mif -pe_dir AP -rpe_pair -se_epi ${input}/b0_revpe.mif -eddy_options " --repol " -eddyqc_all ${input}/eddy/

The information that was saved, is the following:

dwifslpreproc: 
dwifslpreproc: Note that this script makes use of commands / algorithms that have relevant articles for citation; INCLUDING FROM EXTERNAL SOFTWARE PACKAGES. Please consult the help page (-help option) for more information.
dwifslpreproc: 
dwifslpreproc: Generated scratch directory: /vscmnt/gent_kyukon_scratch/_kyukon_scratch_gent/vo/000/gvo00047/vsc41214/LAT_study/dwifslpreproc-tmp-ZNUOJV/
Command:  mrconvert /user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/dwi_degibbsed.mif /vscmnt/gent_kyukon_scratch/_kyukon_scratch_gent/vo/000/gvo00047/vsc41214/LAT_study/dwifslpreproc-tmp-ZNUOJV/dwi.mif -json_export /vscmnt/gent_kyukon_scratch/_kyukon_scratch_gent/vo/000/gvo00047/vsc41214/LAT_study/dwifslpreproc-tmp-ZNUOJV/dwi.json
Command:  mrconvert /user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/b0_revpe.mif /vscmnt/gent_kyukon_scratch/_kyukon_scratch_gent/vo/000/gvo00047/vsc41214/LAT_study/dwifslpreproc-tmp-ZNUOJV/se_epi.mif
dwifslpreproc: Changing to scratch directory (/vscmnt/gent_kyukon_scratch/_kyukon_scratch_gent/vo/000/gvo00047/vsc41214/LAT_study/dwifslpreproc-tmp-ZNUOJV/)
Command:  dirstat dwi.mif -output asym
Command:  mrinfo dwi.mif -export_grad_mrtrix grad.b
dwifslpreproc: 1 spatial axis of DWIs has non-even size; this will be automatically padded for compatibility with topup, and the extra slice erased afterwards
Command:  mrconvert se_epi.mif -coord 2 74 - | mrcat se_epi.mif - se_epi_pad2.mif -axis 2
Command:  mrconvert dwi.mif -coord 2 74 - | mrcat dwi.mif - dwi_pad2.mif -axis 2
Command:  mrconvert se_epi_pad2.mif topup_in.nii -strides -1,+2,+3,+4 -export_pe_table topup_datain.txt
Command:  topup --imain=topup_in.nii --datain=topup_datain.txt --out=field --fout=field_map.nii.gz --config=/apps/gent/CO7/skylake-ib/software/FSL/6.0.3-foss-2019b-Python-3.7.4/fsl/etc/flirtsch/b02b0.cnf --verbose
Command:  mrinfo dwi_pad2.mif -export_pe_eddy applytopup_config.txt applytopup_indices.txt
Command:  dwiextract dwi_pad2.mif -pe 0.0,-1.0,0.0,0.0336 - | mrconvert - dwi_pad2_pe_1.nii -json_export dwi_pad2_pe_1.json
Command:  applytopup --imain=dwi_pad2_pe_1.nii --datain=applytopup_config.txt --inindex=1 --topup=field --out=dwi_pad2_pe_1_applytopup.nii --method=jac
Command:  mrconvert dwi_pad2_pe_1_applytopup.nii.gz dwi_pad2_pe_1_applytopup.mif -json_import dwi_pad2_pe_1.json
Command:  dwi2mask dwi_pad2_pe_1_applytopup.mif - | maskfilter - dilate - | mrconvert - eddy_mask.nii -datatype float32 -strides -1,+2,+3
Command:  mrconvert dwi_pad2.mif eddy_in.nii -strides -1,+2,+3,+4 -export_grad_fsl bvecs bvals -export_pe_eddy eddy_config.txt eddy_indices.txt
Command:  eddy_openmp --imain=eddy_in.nii --mask=eddy_mask.nii --acqp=eddy_config.txt --index=eddy_indices.txt --bvecs=bvecs --bvals=bvals --topup=field --repol --out=dwi_post_eddy --verbose
dwifslpreproc: Command 'eddy_quad' not found in PATH; not running automated eddy qc tool QUAD
Command:  mrconvert dwi_post_eddy.nii.gz result.mif -coord 2 0:74 -strides -2,-3,4,1 -fslgrad dwi_post_eddy.eddy_rotated_bvecs bvals
Function: os.makedirs('/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy')
Function: shutil.copy('dwi_post_eddy.eddy_parameters', '/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy/eddy_parameters')
Function: shutil.copy('dwi_post_eddy.eddy_movement_rms', '/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy/eddy_movement_rms')
Function: shutil.copy('dwi_post_eddy.eddy_restricted_movement_rms', '/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy/eddy_restricted_movement_rms')
Function: shutil.copy('dwi_post_eddy.eddy_post_eddy_shell_alignment_parameters', '/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy/eddy_post_eddy_shell_alignment_parameters')
Function: shutil.copy('dwi_post_eddy.eddy_post_eddy_shell_PE_translation_parameters', '/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy/eddy_post_eddy_shell_PE_translation_parameters')
Function: shutil.copy('dwi_post_eddy.eddy_outlier_report', '/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy/eddy_outlier_report')
Function: shutil.copy('dwi_post_eddy.eddy_outlier_map', '/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy/eddy_outlier_map')
Function: shutil.copy('dwi_post_eddy.eddy_outlier_n_stdev_map', '/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy/eddy_outlier_n_stdev_map')
Function: shutil.copy('dwi_post_eddy.eddy_outlier_n_sqr_stdev_map', '/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy/eddy_outlier_n_sqr_stdev_map')
Function: shutil.copy('dwi_post_eddy.eddy_outlier_free_data.nii.gz', '/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy/eddy_outlier_free_data.nii.gz')
Function: shutil.copy('eddy_mask.nii', '/user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/eddy/eddy_mask.nii')
Command:  mrconvert result.mif /user/gent/412/vsc41214/scratch_vo_user/LAT_study/subjects/LAT_072/geomcorr.mif
dwifslpreproc: Changing back to original directory (/vscmnt/gent_kyukon_scratch/_kyukon_scratch_gent/vo/000/gvo00047/vsc41214/LAT_study)
dwifslpreproc: Deleting scratch directory (/vscmnt/gent_kyukon_scratch/_kyukon_scratch_gent/vo/000/gvo00047/vsc41214/LAT_study/dwifslpreproc-tmp-ZNUOJV/)

To me, the output did not look weird.

This is before dwifslpreproc:

And after dwifslpreproc:

Did something go wrong? And if so: what went wrong and how can I resolve this?

Thank you in advance for checking my post!

Best,
Helena

[jdtournier]

I’ve no idea what might be causing this specific issue (sorry), but on the following topic:

@HelenaV

I just ran dwifslpreproc (3.0_RC4) on our HPC cluster (the people from HPC installed this version, I don’t really get how, but apparently, this is the newest of newest?)

This is currently being discussed on GitHub – and it’s all my fault, sorry…

There’s nothing wrong with using this tag, but I hadn’t expected people to use it for deployment (I’d created it for internal testing purposes, I should have been more careful). What you get is a pretty up to date version of our development dev branch – it’s perfectly stable for use, but not an official release as such, and there’s no corresponding documentation (although the documentation for the dev branch would be close enough). It’s essentially the version that will be released officially in the hopefully very near future…

[ThijsDhollander]

@jdtournier

There’s nothing wrong with using this tag, but I hadn’t expected people to use it for deployment

It does have the change in FOD amplitude evaluation for tractography, but not the “corresponding” change in default thresholds and parameters, nor matching updated documentation or advice for how relevant thresholds would likely have to be altered. Depending on what pipeline users follow, this is not unlikely to have a substantial impact on results. To each their own opinion, but I would personally not risk it. You cannot know whether it mattered until you re-run with e.g. RC3 / master / …

@jdtournier

(although the documentation for the dev branch would be close enough)

Is also not yet updated in some places with regard to the above; maybe other things too (?). Also, documentation for dev can change quite dynamically over time…? So this depends on when the documentation was or is looked at by users.

Again, I wouldn’t take any risks. I personally know from experience e.g. the tractography cutoff matters substantially for both connectomics as well as FBA pipelines, and results can be altered non-trivially and even non-monotonically.

[jdtournier]

All good points, and I fully agree you shouldn’t install 3.0_RC4 if you can help it.

Nonetheless, there’s nothing inherently wrong with the code in this tag – no overt bugs that you won’t also find on previous releases, at least that we know of. It’s completely true that it won’t be equivalent to other releases, but then that’s true of any release: no two releases should ever be considered equivalent in terms of the precise results they will yield.

But the documentation for that tag is indeed unlikely to be completely up to date, and the documentation for the dev branch will progressively diverge over time. So yes, best to steer clear.

I’ll make a separate announcement on the forum about this now.

[HelenaV]

Thank you for your answers, I get what happened with the installation of the MRtrix 3.0_RC4 on our HPC cluster and how it is not ideal to use it right now.

I am just wondering…

1. Do you suggest to continue analyses (I am at the tckgen step) using the 3.0_RC3_latest, given the changes in tractography in the 3.0_RC4 release mentioned by @ThijsDhollander? So switching between the two version?

2. Is there any way to find out what actually happened during eddy? Because if the code of dwifslpreproc is the same as for dwipreproc, I am wondering what is causing this eddy_mask.nii to look like this. I checked in the eddy output folder of some older analyses to see if it happened there as well, but it appears this mask isn’t saved when I use the 3.0_RC3 version.

Thank you for your suggestions!

[ThijsDhollander]

@HelenaV

Do you suggest to continue analyses (I am at the tckgen step) using the 3.0_RC3_latest, given the changes in tractography in the 3.0_RC4 release mentioned by @ThijsDhollander? So switching between the two version?

Opinions might again differ… personally I would in most scenarios re-run, unless maybe this is always supposed to be a “one-off” pre-processing. The typical problems tend to arise when at a later stage, you need to re-visit this preprocessing; e.g. adding subjects, changing a step, having to address some reviewer’s particular requests, … I can imagine you’d remove the RC4 version from those systems so others don’t run into this, but your analysis does depend partially on it. Not that there’s no ways to deal with it, but it might get messy. Even if your initial preprocessing took some time to set up and get it all together, I suppose once you had the design of it figured out, re-running it might be more simple now; i.e. requires far less “human effort” and just waiting for computations to finish? If that’s the case, it might be a no-brainer to simply fire up the preprocessing steps again and mostly just wait for them to finish again? In turn, you get some peace of mind and it potentially makes life easier further down the track. This mostly from real-life experience with scenarios where this had gone wrong, and basically results couldn’t be reproduced anymore; only to result in re-running the whole lot anyway. If you’re still early in the pipeline, you’ve got the opportunity to not allow this to become an issue down the track in any case.

@HelenaV

Is there any way to find out what actually happened during eddy? Because if the code of dwifslpreproc is the same as for dwipreproc, I am wondering what is causing this eddy_mask.nii to look like this. I checked in the eddy output folder of some older analyses to see if it happened there as well, but it appears this mask isn’t saved when I use the 3.0_RC3 version.

I’m also not sure about this one. One way to at least test whether this could’ve had a substantial impact, is to indeed re-run preprocessing with another version (maybe even first just on a single subject, only for the sake of this test). If things differ substantially, it’s certainly worthwhile taking a much closer look…

[jdtournier]

Generally agree with the advice given by @ThijsDhollander here. There’s other reasons you might want to stick to more ‘official’ releases than whether the code is any less buggy.

Ultimately though, you may need to re-run the preprocessing anyway if the issue with the mask that you originally reported does indeed turn out to be a bug in 3.0_RC4 (it’s certainly unexpected). Any chance you could share the input data to one of your dwifslpreproc calls so we can try to reproduce and investigate?

[HelenaV]

Thank you @jdtournier and @ThijsDhollander for your responses. Indeed, I think rerunning everything is the proper thing to do. I just wonder how I can check if things change “substantially”. Just visually? I don’t think the 3.0_RC3 automatically saves this eddy_mask.nii, so is there a way to do this anyway, in order to check if the mask looks better using this version?

@jdtournier, I don’t mind sharing the input data. How can I share this with you?

[jdtournier]

Sorry about the delay - got sidetracked and forgot to press send on my reply…

@HelenaV

I just wonder how I can check if things change “substantially”

The simplest is to generate the sum of squared difference image, using something like:

mrcalc dwi1.mif dwi2.mif -sub 2 - | mrmath - rms -axis 3 rms_diff.mif

You can then open this in mrview and check for how big these differences are. If the image is strictly zero, then you have nothing to worry about. Anything else, and you might be better re-running again.

@HelenaV

How can I share this with you?

Simplest is dropbox or Google Drive – you can send the invite to my account: jdtournier@gmail.com.

[rsmith]

I am wondering what is causing this eddy_mask.nii to look like this.

Me too Could you please send the data to me as well?

I don’t think the 3.0_RC3 automatically saves this eddy_mask.nii

You would only get access to it by using the -nocleanup option.

[HelenaV]

Hi @rsmith, @jdtournier,

I uploaded 1 dataset a week ago on wetransfer. I just saw that it will expire in 19 hours I will send it to your email addresses. Unfortunately I am not able to re-upload or send it via another way at the moment as I am traveling right now. Will be back behind a computer on the 28th.

@rsmith, drop me your email address if you would like it as well.

Thank you for checking this out!

Helena

[jdtournier]

OK, problem confirmed on my system too. I’ll investigate further…

[jdtournier]

Right, problem is due to mismatched strides between the dwi_pad2_pe_1.nii image and the field map in the applytoptup call. I find it surprising that it works as well as it does, to be honest. The output seems to be corrected pretty well, despite this mask being completely wrong. Maybe I just don’t understand what eddy's --mask option is supposed to do…

@rsmith: any reason we don’t set the strides to -1,+2,+3,+4 immediately in the initial mrconvert calls? I reckon that would fix it.

@HelenaV: quick fix in the meantime is to set the strides for all input images to -1,2,3,4 (the default Analyze orientation), and process those:

mrconvert -stride -1,2,3,4 dwi_degibbsed.mif dwi_degibbsed_mod.mif
mrconvert -stride -1,2,3,4 b0_revpe.mif b0_revpe_mod.mif
dwifslpreproc dwi_degibbsed_mod.mif geomcorr_mod.mif -pe_dir AP -rpe_pair -se_epi b0_revpe.mif -eddy_options " --repol " -eddyqc_all eddy_mod/ 

I’m running this as we speak to check whether this makes any difference to the eddy_mask.nii – and/or to the final output…

[jdtournier]

Just to feed back on this idea:

@jdtournier

quick fix in the meantime is to set the strides for all input images to -1,2,3,4

This doesn’t seem to fix it… I get a different issue with the mask: now the front of the brain is missing:

I’m still running the processing with (what we think is) the correct fix applied.

---

We’ve established that the bug dates back to version 3.0_RC1, so this is a long-standing issue. I’m baffled that it would have gone unnoticed for so long… In our defence, this might be because:

Thankfully, as far as I can tell, the difference it makes to the output seems negligible (I’m seeing differences of ≪1% in your data) – though there is one volume where it looks like the outlier rejection was triggered in one case but not the other… So this does make a difference, but it looks like it remains within the expected range of performance. I expect the biggest difference in the output will be confined to borderline outlier slices. We’ll need to look into this further to verify…

[jdtournier]

@jdtournier

I’m still running the processing with (what we think is) the correct fix applied

Just to finish this off: the fix works (the mask is now correct), the difference it makes to the output is equally negligible. As far as I can tell, its impact is limited to minor changes in the outlier detection for borderline slices (i.e. slices that are just on the edge of suspicious). Here’s an example of an affected volume from your dataset (most volumes are unaffected):

Before the fix (current 3.0_RC4 code):

With the fix (current tip of dev branch):

Difference at the affected slice (scaled up for display):

You’d be hard-pressed to spot the slice that was affected in the original data – and what’s more, in this instance, the pre-fix version looks more consistent than the post-fix version (the signal is lower in the corrected, post-fix version, and less consistent with its neighbours, if you look closely). So it seems the pre-fix version decided to outlier reject and replace this slice, while the post-fix version did not – presumably because it is teetering around the threshold for rejection (I don’t think you could pick it as an outlier by visual inspection).

So all in all, while we are undoubtedly fairly unimpressed with ourselves at having failed to spot this for so long, it looks like the impact of this bug are fairly minimal, if not altogether negligible. Obviously, you may want to look into the difference in your own results to satisfy yourselves that this does not warrant reprocessing your data…

[daveyin]

Hello,

I upgraded mrtrix to the latest version, but now eddy mask looks completely off when I overlap it on the AP_PA output image from Dwifslpreproc.

[jdtournier]

Hi, what are you using to display the overlay? What does it look like in mrview? Can you also post the output of mrinfo on both the mask and underlying data images (eddy_mask.nii & AP_PA_out.nii, from what you’re showing)? There’s a chance this is simply a mismatch in the strides, which may or may not be an issue here…

[daveyin]

Yup, using mrview showed proper alignment, so thanks for pointing that out! I think I will provide my own mask to remove the eye. Thanks agin.