Eddy's estimate move by susceptibility option causing crash on denoised data?

Hi there! Disclaimer: Completely new to mrtrix, working on my preprocessing pipeline for first time.

I’ve noticed that when I run dwipreproc on denoised data and include eddy’s --estimate_move_by_susceptibility=true option, eddy will fail. However it will not fail if I skip the denoising step.

Here is the message I get upon failure:
Basisfield:: msg=hadamard: Image dimension mismatch
EDDY::: MoveBySuscCF.cpp::: NEWMAT::ReturnMatrix EDDY::MoveBySuscCFImpl::grad(const ColumnVector&):Exception thrown

My dataset is nothing out of the ordinary: a 3-shell 150 direction AP dataset, with an additional scan of PA bzeros for topup. I’m running the cuda implementation of Eddy. And again, no issues if I leave out that option from Eddy (or if I use non-denoised data with that option).

Thanks!
Mark

Hi Mark,

Congratulations, your data might be too clean to be processed by eddy :fireworks: Eddy’s Gaussian process expects some degree of noise in the data and if denoising works well, eddy might not – (discussion here). As far as I can tell, the error message and dependency on estimate_move_by_susceptibility are different from those reported so far though so there might be a problem with the denoised data or dwipreproc.

Might be worth going into the temporary directories created by dwipreproc -nocleanup and checking the input data to eddy and the exact commands. In particular, are there any differences in image dimensions and voxel sizes (mrinfo)? You could also try to verify if eddy works on denoised data to which you added some Gaussian noise (mrcalc denoised.mif randn $(mrmath 1.mif max -axis 3 - | mrstats - -output mean) -mult 0.5 -mult -add - | mrview - denoised.mif or similar).

Hi Max!

So mrinfo shows all my images dimensions and voxels sizes are matching, so no issue there. Just to be clear, I can run eddy on denoised data with no issues if I exclude that one option from the eddy command.

I can next try to follow your suggestion to add back some noise to the denoised data to see if that will make a difference. Can you walk me through that mrcalc command you provided? What is 1.mif in your example?

Understood, just wanted to make our option parsing does not change the input to eddy.

This might be helpful in understanding our command line fu. 1.mif should have been denoised.mif. The command breaks down to first estimating the intensity level and then applying scaled Gaussian noise to each voxel. This command heuristic here is quite arbitrary which is why I pipe the image to mrview to tweak the noise level visually. Alternatively, you can use the noise map from dwidenoise to calibrate the noise level quantitatively.

mrmath denoised.mif max -axis 3 - | mrstats - -output mean estimates the maximum intensity in denoised.mif across volumes and averages that within the resulting 3D image and prints the result to stdout to get a rough estimate of the image intensity scale. The command mrcalc denoised.mif randn $scale -mult 0.5 -mult -add - | mrview - denoised.mif adds Gaussian noise to each voxel: denoised.mif + ((randn() * $scale) * 0.5)) and pipes the resulting image to mrview which also opens the original for comparison.

1 Like

Got it! So I added noise back into the denoised.mif data and tried again. Same result: eddy fails if I include the --estimate_move_by_susceptibility=true option, and succeeds if I remove it.

I’ve not actually experimented with that capability yet.

Ultimately as with other FSL issue reports, I would recommend trying to process the data independently of MRtrix3. If you can prepare the requisite data and images for topup and eddy without actually invoking dwipreproc, and eddy still fails when that option is included, then the issue is entirely independent of MRtrix3 and needs to be reported to FSL to investigate.

Can you double-check that that capability works fine on non-denoised data but crashes on denoised data, and there are definitely no other variables that differ between those two images? Need to be absolutely sure that that detail is not a red herring.

Hi @MrRibbits

Just a wild random guess based on nothing more …than guessing; anyways: could it be that eddy with that particular option doesn’t play well with negative intensities? That’s the one thing of the top of my head I can think of that the denoising generally introduces in the data, and that some subsequent methods or processing steps might not be (prepared to) dealing with.

One way to test this in the most direct way possible would be to run denoising, and then cap negative intensities to zero. Unless it also doesn’t deal with zeros… if even that would be a sensible explanation, you could try and cap to a very small (but strictly positive) intensity.

Cheers,
Thijs

Hi @ThijsDhollander !
I think you’re onto something. I just thresholded the input to remove negative voxels and it seems to have worked. I’m going go try it once more on another dataset to be certain and report back once complete.

@rsmith: Indeed I can run topup/eddy independent of MRtrix3 successfully without a failure on that same data.

Best,
Mark

1 Like

Makes good sense in some ways: negative signals aren’t physically sensible. Applications relying on such properties would indeed run into issues at times.

Good idea! :+1:

I can easily threshold my eddy_in.nii to remove negative voxels using a fslmaths -thr command, but any chance someone give me the equivalent mrcalc command so I can apply the thresholding to the .mif file prior to invoking the dwipreproc script? Still in the learning phase here, sorry!

The key to working with mrcalc in general, is to become used to working with “postfix” notation, also known as “Reverse Polish Notation” (RPN). The Wikipedia article on this is actually pretty decent, with a nice example; even though it mostly focuses on binary operators. Basically, rather than putting your operator in between your operands (“infix notation”), as you’re probably more used to, you put it after your operands instead. So rather than writing “5 + 3”, you write “5 3 +”. See the Wikipedia article to convince yourself if needed, but essentially one of the benefits is that you suddenly don’t need brackets anymore to write down any expression you can dream of. But instead, the position of your operators in the whole expression will implicitly define “where the brackets sit” for an equivalent infix notation. Also, postfix notation allows in a more generalised way to have operators with different numbers of operands, e.g. unary, binary, ternary, etc… (whereas infix notation strictly only applies to binary operators, because the operator needs to sit between the first and the second operand). For your command, you might for instance be thinking of using the -if operator: it’s basically an if-then-else operator that can check for something (1st operand), and then return the 2nd operand if the checked thing is true, or return the 3rd operand if the check thing is false.

The most explicit way to write your command I can think of (there’s many ways you could achieve the same result; but I personally prefer those that are as explicitly readable as possible), would be:

mrcalc dwi_denoised.mif 0 -lt 0 dwi_denoised.mif -if dwi_denoised_capped.mif

The first little bit up until and including the -lt (“less than”) operator there compares the input image with 0, and returns true if the input image value is less than 0. The result of that (true or false) then is the 1st operand to -if. If the test is returns true, the -if will return the 2nd operand (0), and else it’ll return the 3rd operand (the input image value itself). That is then in the end stored in the output image, which always sits at the end of the expression. So basically, for a postfix-trained reader, this mrcalc expression reads as “if dwi_denoised.mif is less than zero, return zero, else return dwi_denoised.mif”. Or zero-capping indeed. :slightly_smiling_face:

Another, far shorter, way of achieving exactly the same thing though:

mrcalc dwi_denoised.mif 0 -max dwi_denoised_capped.mif

In some peoples’ minds, this immediately yells “capping”, but in others’, it doesn’t. Regardless, you might prefer it because it’s shorter. Up to personal taste, I’d say. :wink:

1 Like

Hi again!

I spent some more time with this and have more info to provide. So the denoising was indeed a red herring. I got lost somewhere in all the file shuffles, but long story short: when I ran eddy (successfully) on the original non-denoised data, the input file I used was created with dcm2niix and not mrconvert (my mistake for not catching that difference sooner and leading everyone astray).

So this is where I’m at:
If I create a .mif from dicom with mrconvert (then mrconvert to .nii for eddy), eddy fails when I include that --move… option.
If I convert dicom directly to .nii using mrconvert, eddy still fails when I include that option.
And if convert dicom to .nii using dcm2niix, eddy succeeds when I include that option.

Mrinfo on both nii files (from dcm2nii and from mrconvert) look identical.

That’s a bit worrying… What’s the output of fslhd for both images?

Pasting below, I spot some differences…

mrconvert dicom to nii

filename mrconvert-10-HARDI-AP.nii.gz

size of header 348

data_type UINT16

dim0 4

dim1 96

dim2 96

dim3 64

dim4 151

dim5 1

dim6 1

dim7 1

vox_units mm

time_units s

datatype 512

nbyper 2

bitpix 16

pixdim0 1.000000

pixdim1 2.500000

pixdim2 2.500000

pixdim3 2.500000

pixdim4 nan

pixdim5 0.000000

pixdim6 0.000000

pixdim7 0.000000

vox_offset 352

cal_max 0.000000

cal_min 0.000000

scl_slope 1.000000

scl_inter 0.000000

phase_dim 0

freq_dim 0

slice_dim 0

slice_name Unknown

slice_code 0

slice_start 0

slice_end 0

slice_duration 0.000000

toffset 0.000000

intent Unknown

intent_code 0

intent_name

intent_p1 0.000000

intent_p2 0.000000

intent_p3 0.000000

qform_name Scanner Anat

qform_code 1

qto_xyz:1 -2.490381 0.122238 0.181823 106.471519

qto_xyz:2 -0.138247 -2.486280 -0.222036 117.008377

qto_xyz:3 0.169969 -0.231236 2.483474 -53.856674

qto_xyz:4 0.000000 0.000000 0.000000 1.000000

qform_xorient Right-to-Left

qform_yorient Anterior-to-Posterior

qform_zorient Inferior-to-Superior

sform_name Scanner Anat

sform_code 1

sto_xyz:1 -2.490381 0.122236 0.181823 106.471519

sto_xyz:2 -0.138245 -2.486280 -0.222036 117.008377

sto_xyz:3 0.169969 -0.231236 2.483473 -53.856674

sto_xyz:4 0.000000 0.000000 0.000000 1.000000

sform_xorient Right-to-Left

sform_yorient Anterior-to-Posterior

sform_zorient Inferior-to-Superior

file_type NIFTI-1+

file_code 1

descrip MRtrix version: 3.0_RC3-204-g266a5b4f

aux_file

dcm2niix dicom to nii

filename 10_HARDI-AP.nii.gz

size of header 348

data_type INT16

dim0 4

dim1 96

dim2 96

dim3 64

dim4 151

dim5 0

dim6 0

dim7 0

vox_units mm

time_units s

datatype 4

nbyper 2

bitpix 16

pixdim0 -1.000000

pixdim1 2.500000

pixdim2 2.500000

pixdim3 2.500000

pixdim4 5.000000

pixdim5 0.000000

pixdim6 0.000000

pixdim7 0.000000

vox_offset 352

cal_max 0.000000

cal_min 0.000000

scl_slope 1.000000

scl_inter 0.000000

phase_dim 2

freq_dim 1

slice_dim 3

slice_name Unknown

slice_code 0

slice_start 0

slice_end 0

slice_duration 0.000000

toffset 0.000000

intent Unknown

intent_code 0

intent_name

intent_p1 0.000000

intent_p2 0.000000

intent_p3 0.000000

qform_name Scanner Anat

qform_code 1

qto_xyz:1 -2.490381 -0.122236 0.181824 118.083931

qto_xyz:2 -0.138245 2.486280 -0.222036 -119.188232

qto_xyz:3 0.169970 0.231236 2.483473 -75.824120

qto_xyz:4 0.000000 0.000000 0.000000 1.000000

qform_xorient Right-to-Left

qform_yorient Posterior-to-Anterior

qform_zorient Inferior-to-Superior

sform_name Scanner Anat

sform_code 1

sto_xyz:1 -2.490381 -0.122236 0.181823 118.083931

sto_xyz:2 -0.138245 2.486280 -0.222036 -119.188232

sto_xyz:3 0.169969 0.231236 2.483474 -75.824120

sto_xyz:4 0.000000 0.000000 0.000000 1.000000

sform_xorient Right-to-Left

sform_yorient Posterior-to-Anterior

sform_zorient Inferior-to-Superior

file_type NIFTI-1+

file_code 1

descrip TE=1.1e+02;Time=134216.158;phase=1;dwell=0.600

aux_file imgComments

OK, looking at the differences between them, there’s a few that stand out:

  • datatype is INT16 for dcm2niix, UINT16 for mrconvert. The latter will stick with what’s specified in the DICOM images, which is UINT16. Add -datatype int16 to match dcm2niix – but I doubt that’s the issue.

  • The data provided by dcm2niix has been regridded to LAS (i.e. left->right, posterior->anterior, inferior->superior), whereas mrconvert leaves it as LPS, as it was in the DICOM. Add -stride -1,2,3 to match dcm2niix – but again, that’s unlikely to be the issue here.

  • There’s information in the decript field provided by dcm2iix (TE=1.1e+02;Time=134216.158;phase=1;dwell=0.600), which is not there with mrconvert. Not sure whether this information is necessary for that option to work.

  • There’s an external aux file specified with dcm2niix (called ImgComments, but no such file with mrconvert. Again, not sure whether this contains information necessary for that option to work.

  • The phase_dim, freq_dim and slice_dim entries are filled in using dcm2niix, but zero-filled with mrconvert. Given the option you’re trying to use, I reckon this is the most likely culprit – although still not all that likely since eddy should have all the information it needs (the per-volume PE direction) from the user-provided acquisition parameters file.

If it is due to that last point, it’s also going to be the hardest to fix, since it’s difficult / impossible to reliably extract this information from the DICOM headers across all manufacturers. Furthermore, it’s a subset of the PE direction handling backend that @rsmith has put in place, and would require quite a bit of work to ensure some form of compatibility so that simple cases like this (where the PE direction is unique across the entire dataset) can be handled seamlessly. On top of which, we don’t currently store the slice direction in our headers (haven’t needed to so far), so we’d need to implement that too. Not so say that we shouldn’t do all of this, but it won’t be an overnight fix… No doubt @rsmith will be able to comment on this, he’s closest to that side of things at the moment.

In the meantime, I suspect getting MRtrix3-generated data to work with that option will require some form of custom handling. Probably the ‘simplest’ might be to let it fail, then go into the temporary folder that the script generated, modifying the header using the FSL tools (e.g. fsledithd) to add the missing information, then running the script again with the -continue option to allow it to carry on from where it left off. Not the most user-friendly approach, I have to admit, but it might at least allow you (& us) to figure out what the cause of the problem is exactly…

I just thresholded the input to remove negative voxels and it seems to have worked.

Given the confusion about the sources of other images being tested, can you re-confirm whether or not this was indeed the case? This is something I could get dwipreproc to do in cases where that particular eddy option is provided in order to prevent others from encountering the same issue, whether or not FMRIB resolve at their end.

There’s information in the decript field provided by dcm2iix

This will almost certainly be ignored; both the fields and format would be specific to dcm2niix, and I highly doubt they would write code that is tied to a specific conversion tool.

we don’t currently store the slice direction in our headers (haven’t needed to so far),

Fake news.

Okay I reconfirmed with a 2nd dataset. Also it’s not the INT16 vs UINT16 difference or the -stride difference.

I can post a download link to one dataset, along w/ the commands I’m using in a script if that’s helpful at all:

I edited the phase_dim , freq_dim and slice_dim entries, and tried again and it ran successfully. Will try again today to be sure that’s what just happened.

So I’ve replicated this again, so I’m pretty sure this is the cause.

  1. I let dwipreproc crash at the eddy step.
  2. I used AFNI’s nifti_tool to edit the eddy_in.nii header.
  3. I started dwipreproc again with -continue, specifying eddy_in.nii.

The only thing I noticed is the 2nd dwipreproc run will crash here:
dwipreproc: [ERROR] Number of b-values reported by mrinfo (4) does not match number of outputs provided by dirstat (0)

I commented out lines 235-251 to avoid that error for now. It’s the section that queries the quality of the diffusion acq scheme. Not sure why it’s complaining or even running that section of code that happens much earlier than eddy.

Hope this helps!