The data was converted from Nifti (dcm2nii) to MRtrix format including the FSL style b-vectors and b-values. I checked the intermediate results from tensor2metric and amp2sh - they look reasonable.
Something that irritates me is the b-values.
mrinfo -shells gives me some odd values:
Firstly, the msdwi2response script is now deprecated; it should disappear on update. The equivalent functionality can now be accessed using dwi2response manual. You may find that some aspects change when using the newer script: for instance, I believe I’ve tweaked the heuristics that trigger a warning regarding the selection of shells, and I would expect it to no longer give a warning in this instance.
Furthermore, if you are able to perform EPI distortion correction and T1 tissue segmentation, the dwi2response msmt_5tt script is much more convenient for deriving multi-shell, multi-tissue response functions if you’re happy with a more automated approach.
Assuming that you want to continue with a manual approach of selecting exemplar voxels for response function estimation, I’d definitely test with the most up-to-date code first. The old msdwi2response script was written in bash and didn’t include sanity checks on the input, so the newer script may flag something. If you’ve still got the issue with the most up-to-date code, we might need to see the data; my only hypothesis was that there was no gradient table embedded in the image header, but you’ve verified that’s not the case.
Regarding the slightly-off b-values, these are calculated directly from the gradient table - but by default, will include scaling of the b-value if the gradient direction is not of unit amplitude. So even if the table quotes exactly 1000 or 2700 for every volume, it’s possible that some or all of your directions don’t have unit norm. Try running mrinfo -dwgrad with and without the -raw_dwgrad option and see if there’s a difference. Also, if this is a problem, you can try bypassing dcm2nii and use MRtrix3 to convert straight from DICOM to .mif.
That’s all for now; let us know if you have the same issue with dwi2response manual
Cheers
Rob
Updating to the latest version indeed removed the msdwi2response script. So I used dwi2response manual but still get-nans in the response file but for the first value.
I used gdb to debug sh2response and found that it is indeed a problem with my data: one of the voxels shows an invalid direction.
Tweaking sh2response to ignore voxels with invalid direction solved the problem:
diff --git a/cmd/sh2response.cpp b/cmd/sh2response.cpp
index 8e53897..024ed56 100644
--- a/cmd/sh2response.cpp
+++ b/cmd/sh2response.cpp
@@ -90,6 +90,11 @@ void run ()
continue;
Eigen::Vector3d d = dir.row(3);
+ // skip voxels with invalid direction
+ if (std::isnan(d(0)) || std::isnan(d(1)) || std::isnan(d(2))) {
+ WARN("invalid direction in \"" + std::string (argument[2]) + "\" - skipping voxel");
+ continue;
+ }
d.normalize();
Math::SH::delta (delta, d, lmax);
I noticed that dwi2response manual produces a different output compared to msdwi2response. First, before I had to specify shells/ b-values and lmax including 0: 0,1000,2700 and 0,6,8. With dwi2response manual it looks like I must not include 0. That might be because in the old version the number of coefficients were padded with 0s to the maximum number of coefficients while dwi2response manual seems to truncate to the least number of coefficients. Or should I compute the response for every shell separately?
How should the final response file look like?
Btw: Is it possible to use the parameters -force and -verbose? I always got the error Error: unrecognized arguments: -verbose/ -force…
Unfortunately, we do not have the data for EPI distortion correction yet so I currently have to go with the simple versions.
Actually, my gradient directions are slightly off unit amplitude. But when I average the product of b-values and b-vector-amplitudes I get 997 for b1000 and 2694 for b27000 while mrinfo -shells returns 994 and 2689. Apart from that, does MRtrix by default use the b-values scaled by b-vector or the raw b-values? Do you know if there is a meaning behind the b-vector amplitude on a Siemens TIM Trio scanner?
I think I came across the same issue just before the ISMRM. I think there’s a bug in there, we’ll need to look into it.
Actually, the scaling applied is by the square of the gradient norm: see here for the exact code. That probably explains the discrepancy.
By default, MRtrix3 will scale the b-values - it’s usually the right thing to do. You can disable it in commands that expect a DW encoding using the -bvalue_scaling option.
There usually is… This used to be the only way to run a multi-shell acquisition in a single scan using a custom gradient file. I think this may have been relaxed now, but I’m not sure. In any case, where the gradient amplitudes differ significantly from unity, you have to assume there is some meaning attached to that, otherwise they would have just stored unit vectors…
Was going through my unread forum posts and realised that I didn’t alert you to the changes I made to the script, so we may not have fully covered this one off.
Tweaking sh2response to ignore voxels with invalid direction solved the problem
Good idea. I’ll add that change to the main code.
I noticed that dwi2response manual produces a different output compared to msdwi2response. First, before I had to specify shells/ b-values and lmax including 0: 0,1000,2700 and 0,6,8. With dwi2response manual it looks like I must not include 0. That might be because in the old version the number of coefficients were padded with 0s to the maximum number of coefficients while dwi2response manual seems to truncate to the least number of coefficients. Or should I compute the response for every shell separately? How should the final response file look like?
dwi2response manualshould require inclusion of 0 in the b-value shell listing in order for the b=0 images to be included as a ‘shell’, and this is certainly the behaviour I’m getting right now. Let me know if this is not the case for you.
It should also be padding any shells with a lower lmax with zeros so that every row contains the same number of columns; this is required in order for msdwi2fod to be able to import the files using standardized matrix file import code. Again I’m testing it right now and getting that result, so please let me know if it’s any different at your end with up-to-date code.
Btw: Is it possible to use the parameters -force and -verbose? I always got the error Error: unrecognized arguments: -verbose/ -force…
This is a limitation of my current adoption of Python’s argparse functionality for command-line parsing. In those scripts that require selection of an algorithm as the first command-line argument, standardized options (those that are available in all scripts) must be specified before the name of the algorithm; e.g. dwi2response -verbose manual ...
Hopefully I’ll get around to finding a fix for this at some point.
just pulled the latest version (0.3.14-619-gbb24a238) and it works like charm.
I used the tournier algorithm with the -voxels option on a single shell subset of my data to generate a mask of voxels with parallel fibres. This mask was then used with the manual algorithm to compute the response for all shells. The coefficients of the response are nicely padded with zeros.
Is it correct to specify only a single value for -lmax in msdwi2fod when using only a single tissue type but multiple shells?
Thank you very much for the great support!
Cheers,
Jan
Is it correct to specify only a single value for -lmax in msdwi2fod when using only a single tissue type but multiple shells?
Correct. In the case of msdwi2fod, the -lmax option applies to the output FODs, not the input responses. So you need a comma-separated list with an lmax value for each tissue; or just one value if you’re fitting only one tissue.
This is indeed the correct specification if you want to fit only a single tissue using data from multiple shells.
But… I would highly recommend against actually doing that. The “WM” FODs that you’ll get when including low b-value and b=0 shells in a single-tissue fit, will be large spheres of random peaks in GM, and enormous spheres of random peaks in CSF. The principle here is that, the more data (i.e., low b-values and b=0 in this case) you try to fit, the more you should take into account to actually represent that data correctly in the fit. A single-tissue model is not suited to fit multiple b-values over the whole brain… you need to account for GM and CSF by including responses for them as well (even if you don’t need the GM and CSF parts of the fit in the end).
