1dDW_Grad_O_Mat++ and missing -proc_dset flag

Hi Afni Experts,
I have started using FATCAT tools for preprocessing of diffusion weighted data. One of the things I wanted to do was the following:

  1. Create an average B0 image by linearly registering multiple B0’s to create an optimized mean image as opposed to just creating a mean (since this does not account for head motion and may make the mean look worse).
  2. Put the mean image as the first volume for eddy current reference

Using an older afni package, there was just 1dDW_Grad_O_Mat which would create a mean B0 image (although not linearly registered i’m guessing) when the -proc_dset flag is used and would output the new dataset with the -pref_dset flag. When updating and using 1dDW_Grad_O_Mat++ these flags are no longer available, and so I wanted to see if this is handled elsewhere or if this functionality is not being replaced? The issue is that even though the bvec/bval files are updated by removing other B0 values, I do not see any flags to output a corrected image. It looks like bref_mean_top will handle the text files but I do not see how to output a corrected image (with the mean B0 in place which I can replace with a linearly registered version).

Any help would be much appreciated.

Thanks,
Ajay

Hi, Ajay-

1) Create an average B0 image by linearly registering multiple B0’s to create an optimized mean image as opposed to just creating a mean (since this does not account for head motion and may make the mean look worse).
→ The choice is yours, but mine own thinking on this kind of point has changed over time. And that is one reason for the change from 1dDW_Grad_o_Mat → 1dDW_Grad_o_Mat++. Basically, for the quality of fit, I believe you are better off not averaging the b=0 volumes together. Those volumes have the highest SNR of the acquired DWI set, and averaging them together effectively reduces their weight/influence on the tensor fit, in the manner used by 3dDWItoDT (if using other software with different fitting, then it might make sense to average together, I s’pose). This point basically propagates to answer other points you raise in this posting.

In terms of performing motion correction across b=0s, I would think that you would actually want to motion correct across all the acquired volumes, b=0s and DWIs alike. And to do that, you are likely best off putting your data into TORTOISE’s DIFFPREP tool for both motion and eddy current distortion, and possibly into their DR-BUDDI tool if you have a dual phase encoded acquisition (e.g., AP-PA acquired data). That is a more rigorous way to handle motion correction, and you benefit by getting other distortion reductions (which exist in all acquired DWI data). NB: TORTOISE 3.0 is now command line processable (THANK YOU, OKAN!), and can be downloaded for free from their website:
https://tortoisedti.nichd.nih.gov/stbb/login.html

2) Put the mean image as the first volume for eddy current reference
If you are using TORTOISE, this wouldn’t be necessary to do.

Using an older afni package, there was just 1dDW_Grad_O_Mat which would create a mean B0 image (although not linearly registered i’m guessing) when the -proc_dset flag is used and would output the new dataset with the -pref_dset flag. When updating and using 1dDW_Grad_O_Mat++ these flags are no longer available, and so I wanted to see if this is handled elsewhere or if this functionality is not being replaced?
As per noted in the first part, the reason this no longer exists is because I believe it would be a non-preferred step on the way to fitting tensors. Note that 1dDW_Grad_o_Mat does still exist in its native form, for backward-compatability reasons, but I don’t think it should be used, really-- the basic functionality of the 1dDW*++ function is easier to use and more to my liking. (Of course, that is just a personal opinion.)

The issue is that even though the bvec/bval files are updated by removing other B0 values, I do not see any flags to output a corrected image. It looks like bref_mean_top will handle the text files but I do not see how to output a corrected image (with the mean B0 in place which I can replace with a linearly registered version).
Well, if you’re going to be doing motion correction on DWIs, you’d want to update the gradient/matrix table appropriately; for b=0 volumes, this isn’t an issue, but, again, I think motion correction on the full dset really makes sense-- TORTOISE handles that part of the processing, as well. NB: I’m not paid by the TORTOISE folks (really), but I would strongly recommend using that for these kinds of processing steps.

Happy to hear your thoughts on that. The above are only recommendations.

–pt

Thanks for the quick reply!

  1. I can see what you mean. The reason I was looking at the average is that there was better definition of the pons as there is more signal on the average (and less distortion), however the shape of the best SNR shape may match the individual directions better for an improved overlap/fit. I will look into this further.

  2. For tortoise 3.0 do you know if we need IDL license in order to use the command line tools? Since I am going to be running the scripts on a grid computer system, I do not think that would necessarily be feasible. In any case, I have been getting very good results for eddy current and warping using AFNI tools (ie 3dAllineate/3dQwarp etc).

Thanks,
Ajay

Hi, Ajay-

1) I can see what you mean. The reason I was looking at the average is that there was better definition of the pons as there is more signal on the average (and less distortion), however the shape of the best SNR shape may match the individual directions better for an improved overlap/fit. I will look into this further.

OK. Note that TORTOISE (primarily+preferably) uses registration to reference anatomical (a T2w one), which would have higher SNR, likely.


2) For tortoise 3.0 do you know if we need IDL license in order to use the command line tools? Since I am going to be running the scripts on a grid computer system, I do not think that would necessarily be feasible. In any case, I have been getting very good results for eddy current and warping using AFNI tools (ie 3dAllineate/3dQwarp etc).

No, TORTOISE DIFFPREP and DR-BUDDI are written in C++ and you don’t need an IDL license. Hallelujah.
I’m glad that the AFNI tools have been serving you well, but note that TORTOISE is designed to deal with the differing contrasts of the DWIs. I’d be curious to know what steps are you using to process in AFNI.

–pt

Hi,
Sorry for the delay, I just noticed your question regarding steps at the end…

  1. 3dAllineate to perform eddy current correction. I created an average B0 as a reference frame. I also aligned my T2 to my B0 reference scan rigidly.
    For subjects with lots of motion you can use -twopass…its robust but takes a while to run.

  2. 3dQwarp for B0 to pseudo-T2 distortion correction but with a slightly lower penalty (1 instead of 5). I also limited number of iterations so that it doesn’t get weird results from over warping since they are close to begin with.

It seemed to give me reasonable results for the corrections. Plus I could use wsinc5 to write out files for improved interpolation.

QUESTIONS:

  1. The part that is a bit more difficult is how to rotate the bvecs. I used MRtrix3 to import the nifti and output the bval/bvec in “FSL format” which I assume is just the row vector and has nothing to do with sign of values. I want to apply my 3dAllineate matrix but I only want to apply motion related 6dof from the 12dof matrix for the rotation of bvecs, which I DO NOT THINK the following afni blogs perform. It seems that I would want to keep rotation and translation and throw out scaling and shear to accomplish this goal, but the scripts below throw out translation (not sure why).

  2. Another question is regarding afni implementations I found:
    http://blog.cogneurostats.com/?p=302
    Vecwarp -matvec tmp.transform.final.${aRow} -input tmp.bvecs.${aRow} -forward >> bvec_rotated.1D

When comparing the post above to one of Daniel’s old post, there was a conflict in the following line:
https://afni.nimh.nih.gov/afni/community/board/read.php?1,52909,54478#msg-54478
Vecwarp -matvec tempmatvec0.1D -backward -input tempvector.1D -force -output temprotgrad.1D > & /dev/null

The difference between these two afni suggestions could be due to the way 3dAllineate and WarpDrive function differently , or there is a mistake in one of these two versions. Is the -forward option correct for applying 3dAllineate matrix?

  1. For Afni’s DTI Fit is the reweight/nonlinear option the same as Restore/iRestore algorithm or different?

Any suggestions would be much appreciated.

Thanks,
Ajay

Just as a followup regarding MRTrix3 for the FSL format (http://mrtrix.readthedocs.io/en/latest/concepts/dw_scheme.html)

“It is important to note that in this format, the gradient vectors are provided with respect to the image axes, not in real or scanner coordinates (actually, it’s a little bit more complicated than that, refer to the FSL wiki for details). This is a rich source of confusion, since seemingly innocuous changes to the image can introduce inconsistencies in the b-vectors. For example, simply reformatting the image from sagittal to axial will effectively rotate the b-vectors, since this operation changes the image axes. It is also important to remember that a particular bvals/bvecs pair is only valid for the particular image that it corresponds to.”

I can output the files in " real formats".
“It is important to note that in this format, the direction vectors are assumed to be provided with respect to real or scanner coordinates. This is the same convention as is used in the DICOM format.”

Does the “FSL” or real format match better to 3dAllineate to apply rotation of bvecs?

Thanks,
Ajay

Hi, Ajay-

I reply to various points below, but in general it sounds like a lot of these things could be accomplished in a single stroke using TORTOISE, which is publicly available and command line-processable now; no IDL license needed. Is there a reason you aren’t using that for the various preprocessing steps? It does a real form of both eddy and EPI distortion correction, with consistent motion correction; it also rotates the DW gradients as it goes.

In general, I haven’t seen a very big difference between 3dDWItoDT and DIFF_CALC for most tensor fitting in the WM tissue. I use 3dDWItoDT (and sure, using -reweight is fine). Nonlinear fitting should always be done for the main tensor fit results.


1) 3dAllineate to perform eddy current correction. I created an average B0 as a reference frame. I also aligned my T2 to my B0 reference scan rigidly.
For subjects with lots of motion you can use -twopass...its robust but takes a while to run.

I don’t see how affine registration could adequately account for eddy current distortions. Affine registration comes from a global fit with a fixed number of levers to pull at for matching, and it doesn’t have the ability to deal with local/twisty/pulling distortions. I know that some other eddy correction tools do just use affine registration, but I don’t see how that kind of alignment could provide enough appropriate correction.


2) 3dQwarp for B0 to pseudo-T2 distortion correction but with a slightly lower penalty (1 instead of 5). I also limited number of iterations so that it doesn't get weird results from over warping since they are close to begin with.
It seemed to give me reasonable results for the corrections. Plus I could use wsinc5 to write out files for improved interpolation.

I’m not sure what pseudo-T2 distortion correction is-- is that EPI distortion correction?

Sidenote: you can use 3dedge3 to make “edge” volumes to help see how alignments really line up across the brain.


QUESTIONS:
1) The part that is a bit more difficult is how to rotate the bvecs. I used MRtrix3 to import the nifti and output the bval/bvec in "FSL format" which I assume is just the row vector and has nothing to do with sign of values. I want to apply my 3dAllineate matrix but I only want to apply motion related 6dof from the 12dof matrix for the rotation of bvecs, which I DO NOT THINK the following afni blogs perform. It seems that I would want to keep rotation and translation and throw out scaling and shear to accomplish this goal, but the scripts below throw out translation (not sure why).

I’d have to read more about the scripts; this would be another case where using TORTOISE would handle all of this consistently and internally.


2) Another question is regarding afni implementations I found:
[blog.cogneurostats.com]
Vecwarp -matvec tmp.transform.final.${aRow} -input tmp.bvecs.${aRow} -forward >> bvec_rotated.1D
When comparing the post above to one of Daniel's old post, there was a conflict in the following line:
[afni.nimh.nih.gov]
Vecwarp -matvec tempmatvec0.1D -backward -input tempvector.1D -force -output temprotgrad.1D > & /dev/null
The difference between these two afni suggestions could be due to the way 3dAllineate and WarpDrive function differently , or there is a mistake in one of these two versions. Is the -forward option correct for applying 3dAllineate matrix?

The different functions do work differently in many respects; one difference between forward/backward might be whether A is aligned to B in one case, and B aligned to A in the other. I’m not certain offhand; I’d have to look at the scripts closely, and probably even test data. But the fact that the scripts use very different functions would lead me to not be surprised that subsequent steps could be different.


3) For Afni's DTI Fit is the reweight/nonlinear option the same as Restore/iRestore algorithm or different?

I think they have similar aspects of weighting and refitting (and can be repeated a given/large number of times); I would have to look at the full Restore paper, though, to see what other tricks might be involved. Basically, ‘-reweight’ just repeats weighting and refitting; my guess is that details differ, such as cost function used and things like that, at a minimum.

So, in summary, my main advice would be to use TORTOISE to address pretty much all of the above issues…

–pt

Hi, Ajay–


Just as a followup regarding MRTrix3 for the FSL format (http://mrtrix.readthedocs.io/en/latest/concepts/dw_scheme.html)
"It is important to note that in this format, the gradient vectors are provided with respect to the image axes, not in real or scanner coordinates (actually, it’s a little bit more complicated than that, refer to the FSL wiki for details). This is a rich source of confusion, since seemingly innocuous changes to the image can introduce inconsistencies in the b-vectors. For example, simply reformatting the image from sagittal to axial will effectively rotate the b-vectors, since this operation changes the image axes. It is also important to remember that a particular bvals/bvecs pair is only valid for the particular image that it corresponds to."

I can output the files in " real formats".
"It is important to note that in this format, the direction vectors are assumed to be provided with respect to real or scanner coordinates. This is the same convention as is used in the DICOM format."

Does the "FSL" or real format match better to 3dAllineate to apply rotation of bvecs?

Offhand, I don’t know. Changing formats across softwares can be really tough, indeed, just because different assumptions are made-- not even with one being right/wrong, just different (and those seem to be the hardest to catch). One has to look at data very closely and know what to expect in some test cases when doing this.

I mainly use dcm2nii(x) for converting DICOM to NIFTI+text files. I can see where changing an image from RAI → AIR or something would have a big effect, because the order of coordinates is effectively changed. I would only want to deal with things where everything is in (x,y,z)-type order, I think, for ease of envisioning what is happening, and I would imagine that that would be the easiest for dealing with 3dAllineate with the grads, as well. Note that the “expected” order for transformation values is usually specified in a program, though; I think, for example, that 3dvolreg uses a slightly different order for its size parameters than those six would appear in standard 3dAllineate format. However, I don’t know that I have had to do this for anything.

–pt

Thanks for the info regarding these differences. I wanted to follow up with you as I looked into it and tested the conversion using dcm2nii (I believe fatcat uses this as well for importing dicoms), and saw that the magnitudes are the same for all three rows as compared with the mrtrix3 output (fsl format only). The first two rows have the opposite sign of the dcm2nii bvec output so this probably has to deal with FSL RAI (dicom order) from MRtrix3 vs LPI (neurological order) from dcm2nii. Using 1dDW_O_grad_Mat -flip_x -flip_y ended up fixing the output to match that of dcm2nii. The reason why I used mrtrix3 for importing is that we have some Philips data using the new dicom format which imported correctly using this tool and was not able to be converted using dcm2nii. I think that if I flip the signs to match dcm2nii (as well as verifying the fiber directions are in the correct orientation) this may be the fix.

Best,
Ajay

Hi, Ajay–

Re. flipping gradients: different softwares and scanner systems just have different ways of referring to positive/negative x,y,z directions. I don’t know why it’s not standard, but many things in MRI like this aren’t (very unfortunately!). I have usually found dcm2nii/dcm2niix to be reliable and low maintenance: I always check, though, to be sure, esp. when I get data from a new scanner, protocol, center, etc.

Basically, +/- differences don’t affect magnitude values of DTI fitting, but they do affect the relative orientations in space (which itself is pretty important). I use tractography to help me feel confident that the system is consistent, and the tool in AFNI-FATCAT for this is @GradFlipTest-- I’d recommend using that. It basically tests out WB tractographic results for each flip (including no flip), and counts tracts to guess a winner: usually, then difference is quite large and clear, and you also get prompted to look at each result to verify.

Note that moving to/from software tools will change the flipping, likely: for example, TORTOISE internally does use a different system than AFNI (such is life), but one can always fix this pretty straightforwardly with @GradFlipTest and 1dDW_Grad_o_Mat++.

–pt

It’s been a while since I wrote that post (http://blog.cogneurostats.com/?p=302) on rotating b-vectors. 3dAllineate and 3dWarpDrive both output a base-to-input matrix for the transform. However, while align_epi_anat.py uses 3dAllineate in my command, it will output two matrices, and I choose the one that seemed to make more sense to me - the matrix aligning the DWI to the anatomical dataset. One could use the other matrix and just invert the transform.

I am in full agreement with Paul: use TORTOISE. It has the ability to 1) deal with EPI distortions via a quadratic warp test (above and beyond what 3dAllineate, Eddy/eddy_correct); 2) handle all the registration to anatomical datasets with rotations of b-vectors. It does other things too, but those two are probably the two that apply most to your questions.

Hi Peter and Paul,
Thank you once again for the clarification and help.

In my case where the B0ref is the base and the dwi scan is the input, would you think the forward transform from 3dAllineate is appropriate or backwards? When I took my T2 (aligned to B0ref) as my input and applied the inverse of the 3dAllineate matrix I was able to get it aligned to my original DWIs. This seems similar to your situation, in which case the forward may make more sense.

Thanks,
Ajay

Hi,
Based on the fact that the documentation (as well as you) pointed out the 3dAllineate matrix goes from base to source instead of source to base I tried using rotate_bvec.sh “backwards” but my values seemed a bit large (i.e original vector [1.17166 0.03995 -.98435] and affine matrix [.999809 .0144843 -.00778 .325539 -.01837 .9851 -.01444 .09296 0.007527 0.01142 .99987 1.33466]. I replaced the first column with 1 0 0, second column with 0 1 0 and translation with 0 0 0 so that the rotation is the only component accounted for. (originally I tried replacing all 3 columns with 0 0 0 but got -nan error when using backwards mode)

Forward yielded new vector as -1.1413 1.0013 -.021987
Backwards yielded new vector as 133.62 -0.468 -133.833

My B0ref is my base and my dti are my inputs. Based on this setup does backwards make the most sense to use or forward?

Thanks

Ajay

I haven’t double checked the direction of the warp. If your tests look correct then you could go with that.

Importantly, I think the correct thing (in terms of matrices) is to leave the identity matrix on the first two columns (or just leave them as they are). Placing actual zeros doesn’t make much sense. Aligning an object with itself would yield the following matrix:


1 0 0 0
0 1 0 0
0 0 1 0

By zeroing out numbers on the diagonal, you’re probably doing something a bit odd that my brain isn’t going to ponder too deeply on a Friday afternoon. In case you’re curious, I wrote some instructions recently for doing preprocessing of DTI data in TORTOISE on the blog. You might check your results against those run through TORTOISE, particularly with tensors fit.

Hi Peter,
Thanks again for the help before. In looking through Tortoise as well as the other diffusion afni blogs I had a general question. In tortoise they will upsample files to 1.5 mm resolution, then calculate eddy parameters and then apply the transforms for eddy current correction to the bvecs file.

If I performed this same operation within AFNI, and use 3dAllineate for eddy current correction of an upsampled scan (i.e. 1.5mm iso from the original 2mm iso), I would then extract the rotational components (invert to make it source to base instead of base to source) using Vecwarp -backward, and apply it to the original bvecs.

Question: Since in theory rotational components should be independent of voxel size, I assume I do not have to take the different grid size into account when multiplying the rotational eddy component with the original bvec file. I wanted to make sure I am correct in this assumption and not overlooking anything in regards to how 3dAllineate,Vecwarp etc work.

Also, does going from a non-isotropic voxel dimension to iso affect this value?

If there is an issue, could it be handled by importing a rotational matrix (derived on 1.5mm iso data), putting a master dataset (orig dti space) and then output recalculated values (if this is needed based on answers above). Is there an AFNI tool poised to do such as process if needed?

Thanks,
Ajay