Problem with dwipreproc output

1

Dear MRtrix3 community,

I’m facing a problem with the output of the dwipreproc script.

By replacing the old “revpe_dwicombine” with “dwipreproc” into my pipeline I’ve noted that the streamlines produced were not pretty as before. The imaging output (result.mif) of both scripts is very similar but the bvecs are different, see screenshot below (the bvecs_OLD, from revpe_dwicombine, is more similar to the original table)

b-matrix_Comparison.png1128×211 31.2 KB

I’m running MRtrix under ubuntu 12.04 and python 2.7.3, FSL V 5.0.9.

Hoping this helps, below I list some issues I have tackled with next release; reading the forum, they were not generally encountered by other users:

1- The ‘‘compact’’ form of number of streamlines does not work (I’m still using 1000000 instead of 1M); “1M” is recognized as “1”.

2- within dwipreproc script I had to replace the occurences of str(rpe_b0_count[0]+1) with str(int(rpe_b0_count[0])+1)

3- dwipreproc: I’ve added ‘’+fsl_suffix+" everywhere a FSL-generated file is invoked (i.e. " runCommand(‘mrcat dwi1_pre_topup’ + fsl_suffix +’ dwi2_pre_topup’ + fsl_suffix +’ dwi_pre_eddy’ + fsl_suffix +’ -axis 3’)". Otherwise, the script was looking at files “.nii” instead of "nii.gz’’.

4- dwipreproc: In a case (other dwi protocol) the script got ‘‘division by 0’’ error at line 241 “factor = 1.0 / math.sqrt(norm2)”

5- eddy_cuda is not present in my FSL bin (I would be really pleased to get it)

Thanks!

Marco

2

Hi Marco,

Thanks for reporting these issues.

Are you sure your version is up-to-date? As far as I can remember, issues 1, 2 and 3 should be fixed in the current master branch.

About issue number 5: I suspect you are using neurodebian. As far as I can tell, FSL is no longer bundeling the source code of recent eddy binaries (eddy_openmp and eddy_cuda), causing them to be missing in the latest neurodebian package. The only way around this is to use the binary CentOS packages provided by fmrib, which does include them. I have been using the CentOS package on various other linux distributions without much problems.

Cheers,
Ben

3

Thank you Ben,

I’m using the updated version of the script (I’m working on the dwipreproc version released after “Fix FSL grad file export bug”).

I’m grateful for your advice on FSL distribution, I’ll try it.

Cheers,

Marco

4

I just realised you are using -rpe_all. Issues 2, 3 and 4 are specific to that case and are indeed not fixed.

Would you mind filing a bug report at https://github.com/MRtrix3/mrtrix3/issues, documenting issues 2, 3 and 4, together with your fixes? Then they will be addressed by the devs ASAP.

Issue 1 is already being tackled by @rsmith.

Cheers,
Ben

5

Perfect! I will do it ASAP. Do you have any idea on how to fix the main issue of wrong directions? It seems to be related to the sign of the components but I do not know why it occurs and how to fix it.

Thank you again,

Marco

6

I have a strong suspicion that the issue with the bvecs will be related to a bug in our handling of the bvecs format that we’ve recently been investigating (see this issue on GitHub for details). We’re about to push the fix out very shortly…

However, I wouldn’t expect it to make a huge amount of difference unless the subject moved considerably during the acquisition. Was this the case…?

7

It does appear to be another issue with the bvecs, as compared to what we’re fixing now: looking at the original post, there appears to be a sign flip in both x and y components in between the old and new bvecs… so I wouldn’t relate it to the amount of motion in this case.

8

Hi Marco,

Thanks for the detailed feedback!

  1. This should have been fixed here. I believe the issue was specific to tckgen, specifically because of how tractography properties were loaded from the command-line by that particular command; if you manage to find any other specific command options where this functionality still fails, please let me know.

  2. Nice spot. I wonder if this might have done an implicit conversion in Python3, for me to have not found it in testing. I’ll push the fix for this imminently.

  3. I’ve tried to use .nii wherever MRtrix is producing a file to feed to FSL, but use the fsl_suffix variable to predict what type of image FSL coimmands are going to provide as output. Even if FSLOUTPUTTYPE is set to NIFTI_GZ, I don’t see the point in performing image compression / decompression for temporary files if the script’s working directory is going to be deleted at script completion anyway. I’ve gone through and found some errors, so I’ll now test all possible scenarios (NIFTI / NIFTI_GZ, -rpe_none / -rpe_pair / -rpe_all).

  4. Your b=0 images must have a null vector direction. It’s more common for b=0 images to have a direction of e.g. [0 0 1] and a b-value of 0; but we have seen this case crop up occasionally. From meory some FSL commands have trouble interpreting this also. It can also cause problems with multi-b-value acquisitions, where the actual b-value needs to be calculated based on the b-value reported in the DICOM and the amplitude of the non-unit-vector direction (even though this conflicts with the DICOM standard). I’ll push a fix soon. Out of curiosity, what processing / conversion steps were performed prior to dwipreproc?

Regarding the main issue: As @Thijs noted, this appears to be different to the other bvecs issue we’re currently addressing, since it’s not just the x axis that’s flipped - and it looks like the flips may have occurred in scanner space rather than image space. Similar to point 4, what conversions and/or processing was performed to obtain the input data for dwipreproc?

Cheers
Rob

9

Hi Robert,
thank you for the support.

DWI data (from a 3T Siemens scanner) go quite directly to dwipreproc script:

I’m first converting DICOMs to nifti and then to .mif:

mrconvert $examFolder/HARDIAP_0* DWI_AP1.nii -export_grad_fsl bvec bval
mrconvert DWI_PA1.nii DWI_PA1.mif -fslgrad bvec bval

Old ‘‘combine’’ command: revpe_dwicombine ap DWI_AP1.mif DWI_PA1.mif DWIeddy.mif

Current ‘‘combine’’ command: dwipreproc AP DWI_AP1.mif DWIeddy.mif -rpe_all DWI_PA1.mif -export_grad_fsl bvecEDDY1 bvalEDDY1 -nocleanup -tempdir /tmp/-$sub-eddyPROC

Is it possible that the flipping of directions is due to a conversion problem such as the issue 2?

Thank you again

Marco

10

Any reason you’re converting to NIfTI first…? What’s wrong with just:

$ mrconvert $examFolder/*HARDI*AP_0*  DWI_PA1.mif

?

11

No particular reason! It is just a refuse from a early version.

Cheers

Marco

12

All good - just wondering whether there was some bug we needed to know about…

13

All good - just wondering whether there was some bug we needed to know about…

Well, I’d still be suspicious given the nature of the bvecs error reported.

@Marco_Aiello:

  1. Can you try converting directly from DICOM to .mif, then run dwipreproc? Though I wouldn’t trust comparing against the result of an older script to determine whether the result is right or wrong; it’s more useful to know this directly from the data. A couple of ways you can check to see if your gradient table is significantly corrupted are:
  • Visualise the raw diffusion data directly by loading your DWI series as a ‘dixel-based ODF’ in the mrview ODF tool.
  • Perform a quick run of deterministic tensor tractography, and make sure that in deep WM single-fibre structures, the streamlines are following the anatomy.
  • Perform CSD and visualise the FODs, make sure the fibre orientations follow continuous paths between voxels.
  1. Not sure if it was just a mis-quote in your post:
mrconvert $examFolder/*HARDI*AP_0*  DWI_AP1.nii -export_grad_fsl  bvec bval 
mrconvert DWI_PA1.nii DWI_PA1.mif  -fslgrad bvec bval

You should absolutely avoid exporting a bvecs/bvals file pair from one image, and importing it into another. Even if the two images have the same transform / voxel size / dimensions, if the data strides of the two images differ then their gradient tables in FSL format will be different.

Besides, there’s no reason to be using bvecs/bvals either at this point, or at the output of dwipreproc: in both cases, if you simply use the .mif format, the gradient table will be there in the header ready to go. You can subsequently export the final gradient table to bvecs/bvals using mrinfo -export_grad_fsl if you wish.

If taking these into account, your gradient table appears correct before dwipreproc but incorrect afterwards, this is definitely something we need to know about.

Cheers
Rob

14

Hi Robert,

during last days I’ve tried my pipeline following your advices after updating to the last release of dwipreproc. Now the output of dwipreproc is as expected (really nice!).

I had just to work around the following error that occurred after the ‘‘eddy’’ step of the script:

File “/usr/local/mrtrix3/scripts/dwipreproc”, line 370, in
if (abs(i-t) / (0.5*(i+t))) > 0.001:
ZeroDivisionError: float division by zero

In my case, it worked fine by replacing "if (abs(i-t) / (0.5*(i+t))) > 0.001 " with an always false condition such as " if 10< 0.001 " but I’m sure that is not a suitable solution for all cases.

Many thanks to the MRtrix team for your amazing work and support,

Marco

15

Ah yes, thanks for reporting that; that should definitely be an absolute threshold rather than a relative one, don’t know what I was thinking… I’ll push a fix so you don’t need that hack. Thanks!

16

Dear MRtrix developers,

I’m still facing the main issue of this topic (flipping of x and y directions of the output of dwipreproc -rpe_all script) . In the previous message I’ve erroneously written that it was solved (I was looking to data processed with revpe_dwicombine).

Checking for a solution, I’ve found that a possible mistake may arise from line 278 of dwipreproc:

if PE_design == ‘Pair’:
runCommand(‘mrconvert ’ + series_path + ’ - -stride -1,+2,+3,+4 | mrinfo - -export_grad_fsl bvecs bvals’)

There is the -stride options that, I guess, is strictly data-dependent and the conversion to FSL bvecs get wrong in my case where -stride from original data is [ -1, -2, 3, 4].

Is it feasible to omit the “-stride” option at this line?

Again, during early attempts I’ve reconstructed good streamlines from DWIeddy.mif (resulting from dwipreprocess) by manually flipping the gradient table in a quite contorted manner:

mrtransform -flip 0,1 DWIeddy.mif DWIeddyFlip.mif
mrinfo DWIeddyFlip.mif -export_grad_mrtrix bvecsF
rm DWIeddyFlip.mif
mrconvert DWIeddy.mif DWIeddyf.mif -grad bvecsF

Is there a way to directly flip the axes of a gradient table?

Thank you again,

Marco

17

The bvecs import/export should now be correct, as of 11 days ago. I’ve spent quite a bit of time getting to the bottom of it, in consultation with the FSL team, and tested things worked correctly within FSL. So I’m pretty sure there’s no actual mistake at this point, assuming you’re running a recent build.

However, these changes may have invalidated existing bvecs files, under certain conditions (namely, when the image associated with the bvecs is stored using neurological convention, with a positive determinant for the transform matrix). In this case, you might need to invert the values in the first row of the associated bvecs file - your procedure would appear to do just that, from what I can tell. But to be honest, it’s probably just as simple to import your data afresh, especially since any modification to the DW gradient table should be viewed as suspicious: the data should be imported correctly, especially if straight from DICOM. And in your case, flipping the directions after EDDY means that a recent version of EDDY would have performed the wrong reorientation of the gradients - including during its own internal processing. So even if you don’t use the reoriented gradients produced by a recent version of EDDY, there’s a chance the processing itself might have been subtly affected (although to be fair, that would probably be a very minor effect).

As to providing ways to directly modify the gradient table: I don’t think I’ll ever do this. There’s too many ways to get things wrong, with all the different conventions and coordinate systems involved. In many cases, a straight flip of one of the axes may be the right thing to do, in other cases it won’t be. Often, it may seem to do the right thing, but if the images are angulated, it may be wrong in subtle ways that won’t be immediately apparent. Sometimes flipping a component of the table in the bvecs format might be the right thing, sometimes you might need to flip a component in the MRtrix format encoding file - the difference being that the former is specified with respect to the image axes, the latter with respect to the scanner axes. If you know what you’re doing enough to understand these implications, I’d expect you’d have no trouble writing a simple script to do this (I find something as simple as awk is often sufficient for this).

18

Checking for a solution, I’ve found that a possible mistake may arise from line 278 of dwipreproc:

if PE_design == 'Pair':
runCommand('mrconvert ' + series_path + ' - -stride -1,+2,+3,+4 | mrinfo - -export_grad_fsl bvecs bvals')

There is the -stride options that, I guess, is strictly data-dependent and the conversion to FSL bvecs get wrong in my case where -stride from original data is [ -1, -2, 3, 4].
Is it feasible to omit the “-stride” option at this line?

Definitely can’t make that change. Those strides are explicitly set in that way in order to put the image in FSL’s preferred LAS format prior to exporting the bvecs/bvals. Moreover, line 279 won’t even execute if you’re using the -rpe_all option; only with the -rpe_pair option.

The theory behind the script is that any data being sent to FSL will be converted to LAS (-stride -1,+2,+3(,+4)) before being processed by FSL, and will be converted back to the strides of the original image only at creation of the final output image.

Thinking about it as I type, the latter may be the cause of the issue. The final output image will mimic the strides of the input image; but with your input image strides (-1,-2,+3,+4), your output image is neither LAS nor RAS but LPS (presumably because it came from a DICOM).

So a few things you can try:

  • Run dwipreproc directly on the DICOM directory. When you pre-convert to NIfTI and bvecs/bvals, literally the first thing the dwipreproc script is doing is converting to .mif and embedding the gradient table in the header.
  • Use a .mif as the output of dwipreproc rather than a NIfTI.
  • Having generated a .mif from dwipreproc, convert to NIfTI and bvecs/bvals after the fact, using the correct syntax to produce a LAS image for FSL:
    mrconvert dwi.mif dwi.nii.gz -stride -1,+2,+3,+4 -export_grad_fsl bvecs bvals

Once you can run the script with non-erroneous output, it’s then possible to pinpoint precisely which aspect is introducing the fault. I’ll start running some tests at my end as well.

@jdtournier (and other devs): Think I suggested in the pull request: Are we better off forcing strides to be either +1,+2,+3 or -1,+2,+3 on NIfTI export?

19

Sorry, I think I’m getting confused with all the various issues that have been discussed so far… Given the y flip, this can’t relate to the NIfTI bvecs import/export issue, and hence shouldn’t be affected by the strides. I’m worried the initial data import might be the problem. The simpler way to check is by running a short whole-brain tensor tracking run directly from whatever you feed into the script, without running deipreproc. I often run:

$ tckgen -algo tensor_det DWI.mif -seed_image mask.mif -num 10000 test.tck

And display the tracks in MRView to check for consistency with the anatomy in various planes.

Can you try this and report back as to whether the input data are OK? If they are, then there’s definitely something that we’ll need to investigate further in the script. If not, then we’ll need to figure out whether the DICOM import is the problem…

20

I’m pretty sure it’s an error in the dwipreproc script that only applies to the -rpe_all option; have duplicated the issue locally. Waiting for eddy to finish processing to confirm whether or not I’ve found the source of the problem.