rotate bvecs in revpe_distcorr

[bjeurissen]

Currently, the revpe_distcorr script does not rotate the gradient directions.

This can be implemented fairly straightforwardly as explained in http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/EDDY/Faq#Will_eddy_rotate_my_bevcs_for_me.3F

It doesn't seem to be a lot of work (apart from getting the coordinate system conventions right) and might be a useful feature.

What do you guys think?

[Lestropie]

I didn't do it at the time because there was still ambiguity about which parameter was which. I also wasn't sure whether or not an FSL update would then include this step. If we can get the rotation right without ambiguity then it should be done.

From memory the script already re-imports the created bvecs / bvals files (in case of FSL updating them itself). Probably safest if the rotation is done on the bvecs data, then re-imported into the header; least likely to result in errors due to axis conventions etc..

[tjroine]

Hi,

Here's a script which can be used to rotate bvecs based on the results from FSL's eddy. Unfortunately I don't have a c++ version yet, but maybe this will be of help in implementing this step into the revpe_distcorr script. In any case this script can be used individually (in Matlab) before performing CSD. I personally use it in a bash script with an extra exit command included in the end of the function, to be able to continue after the rotation to the next step in the script.

Cheers, Timo

Script updated on 9.3.2015 (corrected rotation matrices, and works on *.bvec instead of grad.txt - includes also instructions on exporting and importing bvecs into mrtrix3): https://www.dropbox.com/s/co1drajyuv97ona/rotate_bvecs.m?dl=0

[bjeurissen]

The posted script does not correctly rotate around the x- and y-axis. Also, as Rob suggested, it should not rotate the mrtrix grad file, but the fsl bvec file and then re-import into the header to avoid issues with axis conventions.

@Lestropie, I could easily include this in the revpe_distrcorr python script (and/or the libraries it depends on), the only issue is that the matrix operations would result in a numpy dependency, which might be undesirable. What do you think?

[tjroine]

Hi,

I wasn't aware of this difference between mrtrix grads and fsl (or other software) bvecs coordinate system. If I understood correctly, although a same gradient table in the scanner is used, the resulting bvecs will depend on the possible rotation of the FoV? Does this mean that a separate bvecs file should always be saved for each subject? Can this be done if the scanner doesn't write gradient information into the DICOM header?

I think this difference is an important distinction, which will affect many researchers and could be emphasized more in the wiki.

Anyway, the new version of the function will require the intermediation through bvecs and should work correctly irrespective to the possible differences in the coordinate systems.

https://www.dropbox.com/s/co1drajyuv97ona/rotate_bvecs.m?dl=0

Cheers, Timo

P.S. If you want to implement this in mrtrix3 or in the revpe_distcorr script without Numpy, these c-functions may help: http://www.cs.rochester.edu/~brown/Crypto/assts/projects/adj.html

[jdtournier]

I wasn't aware of this difference between mrtrix grads and fsl (or other software) bvecs coordinate system. If I understood correctly, although a same gradient table in the scanner is used, the resulting bvecs will depend on the possible rotation of the FoV? Does this mean that a separate bvecs file should always be saved for each subject? Can this be done if the scanner doesn't write gradient information into the DICOM header?

Well, really the difference is between DICOM and FSL's bvec format. MRtrix follows the DICOM patient centred coordinate system model, whereas FSL chose a format that is relative to the image axis, presumably because it simplifies downstream processing applications. Unfortunately, it also makes these gradient tables extremely brittle, often in unexpected ways. The first problem is that images converted from DICOM are rarely stored using the NIfTI 'default' ordering, since the row & column vectors are inverted between DICOM & NIfTI (assuming the data were acquired in the transverse plane) - things get more complicated still for coronal or sagittal acquisitions. Anything that modifies the image axes by reformatting the data will also require a corresponding change to the bvecs. So reformatting an image so that its row & column vectors are as per the NIfTI standard will mean the bvecs x & y vectors will also need to be inverted.

Where the bvecs format is helpful is in handling rotations - it's insensitive to that since the gradients will rotate along with the image given the bvecs definition. This is the reason why some of the FSL applications can afford to drop the transform information altogether, without completely compromising downstream processing applications. But this will only hold true if no actual reformatting of the data takes place - i.e. only the transform matrix is modified, not the order of the data on file (at least for NIfTI).

To answer your question about whether this means each subject would need its own bvecs, that depends on the specifics of your scanner & acquisition protocol, but in general, yes. On most scanners, the gradients are played out in scanner (or at least patient-centred) coordinates, irrespective of the image acquisition plane. That said, Philips will play them out in image coordinates unless overplus is used (since with that strategy you couldn't play out two axes' worth of gradient along one of the axes alone). On Siemens, you can supply your own encoding and specify the "prs" coordinate system, in which case the gradients will be played out in image coordinates - in any other case, including the built-in schemes, you'll get them in world/scanner/pcs coordintates. I have a feeling with GE, the gradients are played out relative to the image axes (at least that's how they were stored until the DICOM standard was updated to standardize all this). Assuming your data will be acquired with arbitrary tilts (which you should always assume, irrespective of what you think the radiographers should be doing...), it is almost never safe to assume that the gradients you specified for the scanner can be translated as-is to the bvecs file.

If you use recent versions of mrconvert & mrinfo to convert from DICOM to NIfTI, this should all be taken care of for you. But the minute you produce a NIfTI image with its bvecs/bvals file, you should always assume that these vectors apply to that image and that image alone.

I think this difference is an important distinction, which will affect many researchers and could be emphasized more in the wiki.

Good suggestion. I'll try to write up a full description of all this when I have a minute...

[dchristiaens]

Since this issue came up on Google+ recently, I thought I'd send it back to the top of the list. I came across this recently too and I think a proper fix can and should be implemented.

According to @bjeurissen, eddy should handle the b-vector rotation directly since FSL 5.0.9. While this is obviously the best solution, Kerstin commented that corrected gradients are written to a file <basename>.eddy_rotated_bvecs. The current revpe_distcorr script does not yet deal with this, and re-imports the original bvec file instead. This should be a small fix, and I'm happy to take this on. Similarly, I would propose to at least output a warning when used with an older FSL version. It would then be up to the user to either ignore the issue, upgrade FSL, or rotate the b-vectors manually.

Additionally, I was thinking it may be useful to have a script that calls eddy for data without reverse phase encoding. I could just simplify Rob's script and call it distcorr, or if we really want one single script the interface could be adapted to handle such data. Any thoughts?

[jdtournier]

Seems sensible to me... Any objections from anyone?

[draffelt]

All sounds good to me On 3 Dec 2015 7:31 am, "J-Donald Tournier" notifications@github.com wrote:

Seems sensible to me... Any objections from anyone?

— Reply to this email directly or view it on GitHub https://github.com/MRtrix3/mrtrix3/issues/141#issuecomment-161424647.

[bjeurissen]

A potential issue might be that the FSL manual recommends against eddy without topup on gradient sets that don't sample the full sphere. Not sure how serious this is. All sounds good to me On 3 Dec 2015 7:31 am, "J-Donald Tournier" notifications@github.com wrote:

Seems sensible to me... Any objections from anyone?

— Reply to this email directly or view it on GitHub https://github.com/MRtrix3/mrtrix3/issues/141#issuecomment-161424647.

— Reply to this email directly or view it on GitHub https://github.com/MRtrix3/mrtrix3/issues/141#issuecomment-161428750.

[tjroine]

Sounds like a good plan to me. I have used eddy without topup for 64 directions uniformly sampled but without reverse phase-encoding. The results look good and they should be better than with alternative approaches. The rotated bvecs from eddy are almost equal to those provided by my script (most likely just a precision difference from the calculations in matlab vs eddy).

[dchristiaens]

I tested revpe_distcorr in both FSL 5.0.2 (produces a warning), and 5.0.9 (rotates the gradient directions). The new script distcorr is identical but stripped of all references to topup.