Dwi2response with msmt_5tt HCP tutorial

Just as an extra comment: while you can indeed already use this algorithm for 3-tissue response function estimation in single shell data, note that you can’t do much with these responses yet thereafter: the MSMT-CSD method still requires at least 3 b-values (including the b=0 data) if you want to separate WM, GM and CSF tissue classes. The best you could do on that front at the moment for single-shell data, would be to ask the MSMT-CSD algorithm just for, e.g., the WM and the CSF (and ditch the GM response). That would at least already do a (suboptimal) form of free-water elimination (and a-tiny-bit-of-GM-elimination, due to a slight bias).

That’s all about to change (hopefully not too long from now) though; I’ll provide an implementation of single-shell 3-tissue CSD at some point.

Yes, exactly! Sorry, re-reading my comment, was actually quite lacking of details :slight_smile:
You are right, I was referring to what I can see after applying dhollander or msmt_5tt AND the msmt_csd thereafter. And it seems exactly what you observed in your abstract, cool!

Thanks again,
cheers

Sara

Hooraay for reproducible research!! :sunglasses: :smile:

Thanks for confirming; this is really helpful! :thumbsup:

@Thijs, Unfortunately, I failed using the parameter 'dhollander '. When run the script:

dwi2response dhollander DWI.mif RF_WM.txt RF_GM.txt RF_CSF.txt

It produced more than 10GB temptory files before output errors:

Command: mrcalc _crudegmhighselect.mif 1 _crudegmlowselect.mif -if refined_gm.mif -datatype bit
Command: mrcalc _crudewmoutliers.mif safe_sdm.mif 0 -if 1.93704 -gt 1 crude_csf.mif -if _crudecsfextra.mif -datatype bit
Command: mrcalc _crudecsfextra.mif safe_sdm.mif 1.93704 -subtract 0 -if - | mrthreshold - - -mask _crudecsfextra.mif | mrcalc _crudecsfextra.mif - 0 -if refined_csf.mif -datatype bit
dwi2response: Running 'tournier' algorithm to select 1619 single-fibre WM voxels.
Command: dwi2response tournier dwi.mif _respsfwmss.txt -sf_voxels 1619 -iter_voxels 16190 -mask refined_wm.mif -voxels voxels_sfwm.mif -quiet -tempdir /tmp/dwi2response-tmp-P8M1DI/
dwi2response:
dwi2response: [ERROR] Command failed: dwi2response tournier dwi.mif _respsfwmss.txt -sf_voxels 1619 -iter_voxels 16190 -mask refined_wm.mif -voxels voxels_sfwm.mif -quiet -tempdir /tmp/dwi2response-tmp-P8M1DI/
dwi2response: Output of failed command:

 dwi2response: [ERROR] unknown option "-sf_voxels"**

PS: The multi-shell DWI data is HCP 100048
DWI.mif was produced by
mrconvert …/T1w/Diffusion/data.nii.gz DWI.mif -fslgrad …/T1w/Diffusion/bvecs …/T1w/Diffusion/bvals -datatype float32 -stride 0,0,0,1

Thanks

Hi @Haiyong_Wu,

It looks like the problem you’re running into is due to remnants of an older MRtrix version still being there. I suppose you’ve been an MRtrix user since before the big update we did in March?

The problem is that you still may have a command called dwi2response sitting around, and that is causing conflicts with our new dwi2response script.

If this is the case, try running the following at the command line:

./build clean
git pull
./configure
./build

…and see if that fixes things. This should clean all the old stuff up, get the latest version of everything, and set everything up fresh.

Cheers,
Thijs

Hey @Thijs!!

running “dhollander” algorithm for single shell I found this problem :confused:

dwi2response: [ERROR] Command failed: mrcalc _crudenonwm.mif safe_sdm.mif  -subtract 0 -if - | mrthreshold - - -mask _crudenonwm.mif | mrcalc _crudenonwm.mif - 0 -if crude_csf.mif -datatype bit
dwi2response: Output of failed command:
mrcalc: [ERROR] not enough operands in stack for operation "if"
mrthreshold: [ERROR] no filename supplied to standard input (broken pipe?)
mrthreshold: [ERROR] error opening image "-"
mrcalc: [ERROR] error converting string "-"
dwi2response: Changing back to original directory (/mnt/projects/LifeMabs/DWItest/mrtrix_2000/b2000_50dir)
dwi2response: Script failed while executing the command: mrcalc _crudenonwm.mif safe_sdm.mif  -subtract 0 -if - | mrthreshold - - -mask _crudenonwm.mif | mrcalc _crudenonwm.mif - 0 -if crude_csf.mif -datatype bit
dwi2response: For debugging, inspect contents of temporary directory: /tmp/dwi2response-tmp-46Y7GP/

:sweat_smile:

Hi @Nayo,

Looks to me like this might be a similar problem to what’s going on in this parallel discussion: One of the mask images derived during the response function estimation algorithm is empty, which is resulting in a call to mrstats not returning any result. In your particular case this subsequently results in a malformed call to mrcalc (See the double space between ‘safe_sdm.mif’ and ‘-subtract’? There’s supposed to be a number in the middle), but it’s the preceding empty mask that’s the fundamental problem. My guess is that there’s something particularly unusual about your data that’s breaking one of the assumptions within the dwi2response dhollander algorithm; so the best bet would be to make an example image available so that the problem can be reproduced at our end.

Cheers
Rob

Hey @Nayo,

Sorry, but it seems I had completely missed your post here! Rob is right in that this is a combination of a specific mask in the algorithm being empty, and then the way how this makes the script crash due to a particular technical caveat in a command that just surfaced for another user in another context too… I’ll eventually add in a useful error message for this particular event. You’d get that error message in that case any way, so you’d still not be able to proceed.

In your particular case, it seems the region within the initial brain mask that has FA values below 0.2 is empty. That would be very odd indeed… Maybe even the original brain masking went wrong. If you want, you could make the DWI image that you provides to the script available somewhere (via Dropbox of something), and let me know in private (you can use the messaging feature on the forum, if you’re ok with it); so I can take a look at it.

Anyway, it seems you’re after single-shell processing here; so for now (while you may share the dataset with me in between), you could still proceed using just the tournier algorithm in dwi2response, combined with dwi2fod csd afterwards, for simple single-shell single-tissue CSD.

Hope this helps a bit (for now). :slight_smile:

Cheers,
Thijs

Hi bjeurissen

I met the same problem using hcp data(the subject 163129,166438,169747), and I checked the 5tt.mif .I think the 5tt.mif is right.Are there other reasons causing the problem ?

Thanks

Hello,
I got a similar problem with my data, the mask is empty. Something that I know is NOT recommendable was performed in order to obtain and embed bval & bvec into dwi.mif file (not sure whether this is the cause of the problem). I got only raw data, converted to nii format by Bru2Nii, and extracted bvec and bval from $PVM_DwDir and $PVM_DwEffBval respectively from a file called “method” at the raw data, and arranged it in fsl format.

My data looks like this:
Carloss-MacBook-Pro:mrtrix_test carlosengutierrez$ mrinfo 14.mif
************************************************
Image: “14.mif”
************************************************
Dimensions: 128 x 30 x 128 x 130
Voxel size: 3.5 x 7 x 3.5 x 4
Data strides: [ 1 3 -2 4 ]
Format: MRtrix
Data type: signed 16 bit integer (little endian)
Intensity scaling: offset = 0, multiplier = 4.3775600000000001e-05
Transform: 1 -0.0005371 0.003784 -212.5
0.00053 1 0.001875 -86.65
-0.003785 -0.001873 1 -217.8
comments: 13_HR_DTI_128axis_b3000
dw_scheme: [ 130 entries ]
mrtrix_version: 005da18a

And the shells:
Carloss-MacBook-Pro:mrtrix_test carlosengutierrez$ mrinfo 14.mif -shells -shellcounts
0 3018.33 
2 128 

Carloss-MacBook-Pro:mrtrix_test carlosengutierrez$ mrinfo 14.mif -dwgrad
         -0           0           0           0
         -0           0           0           0
  -0.140642    0.989951   0.0147479     3041.35
  0.0827185    0.995958  -0.0350076     3042.13
  0.0306206    0.980789    0.192654     3040.53
  -0.375125    0.925683   0.0489144     3037.26
  -0.081265    0.973005   -0.216003     3034.51
   0.144104    0.952087   -0.269749     3031.95
    0.29863    0.950621  -0.0844939     3033.45
    0.24308    0.956728    0.159948     3032.08
 -0.0786406    0.924687    0.372518     3032.94

Then I tried the following:
a) dhollander:
Carloss-MacBook-Pro:mrtrix_test carlosengutierrez$ dwi2response dhollander 14.mif wm.txt gm.txt csf.txt

Command: mrcalc safe_mask.mif safe_fa.mif 0 -if 0.2 -gt crude_wm.mif -datatype bit
Command: mrcalc crude_wm.mif 0 safe_mask.mif -if _crudenonwm.mif -datatype bit
mrstats: [ERROR] Cannot output statistic of interest; no values read (empty mask?)
Command: mrcalc _crudenonwm.mif safe_sdm.mif -subtract 0 -if - | mrthreshold - - -mask _crudenonwm.mif | mrcalc _crudenonwm.mif - 0 -if crude_csf.mif -datatype bit
dwi2response:
dwi2response: [ERROR] Command failed: mrcalc _crudenonwm.mif safe_sdm.mif -subtract 0 -if - | mrthreshold - - -mask _crudenonwm.mif | mrcalc _crudenonwm.mif - 0 -if crude_csf.mif -datatype bit
dwi2response: Output of failed command:
mrcalc: [ERROR] not enough operands in stack for operation "if"
mrthreshold: [ERROR] no filename supplied to standard input (broken pipe?)
mrthreshold: [ERROR] error opening image "-"
mrcalc: [ERROR] error converting string "-"
dwi2response: Changing back to original directory (/Users/carlosengutierrez/datathon/mrtrix_test)
dwi2response: Script failed while executing the command: mrcalc _crudenonwm.mif safe_sdm.mif -subtract 0 -if - | mrthreshold - - -mask _crudenonwm.mif | mrcalc _crudenonwm.mif - 0 -if crude_csf.mif -datatype bit
dwi2response: For debugging, inspect contents of temporary directory: /tmp/dwi2response-tmp-ZWA1YX/

b) tournier

Carloss-MacBook-Pro:mrtrix_test carlosengutierrez$ dwi2response tournier 14.mif 14_tournier.mif -grad encode.b
....
Command: fod2fixel iter0_FOD.mif -peak iter0_peaks.msf -mask mask.mif -fmls_no_thresholds
dwi2response: Deleting file: iter0_FOD.mif
Command: fixel2voxel iter0_peaks.msf split_value iter0_amps.mif
dwi2response: 
dwi2response: [ERROR] Command failed: fixel2voxel iter0_peaks.msf split_value iter0_amps.mif
dwi2response: Output of failed command:
fixel2voxel: [ERROR] fixel image is empty
dwi2response: Changing back to original directory (/Users/carlosengutierrez/datathon/mrtrix_test)
dwi2response: Script failed while executing the command: fixel2voxel iter0_peaks.msf split_value iter0_amps.mif
dwi2response: For debugging, inspect contents of temporary directory: /tmp/dwi2response-tmp-EZILZ5/

Please advice,
Thanks, Carlos.

Hello,
About mask images getting empty, I found the problem for my case, the gradient data was not in MRTrix format.
I used commnad convert_bruker to obtain mih file with the correct gradient data format, then I converted mih to mif.

$ convert_bruker 14/pdata/1/2dseq 14.mih
$ mrconvert 14.mih 14_1.mif

You can see the before-after info of my data, dimensions, voxel size and Transform matrix changed


Image: “14.nii”


Dimensions: 128 x 30 x 128 x 130
Voxel size: 3.5 x 7 x 3.5 x 4
Data strides: [ 1 3 -2 4 ]
Format: NIfTI-1.1
Data type: signed 16 bit integer (little endian)
Intensity scaling: offset = 0, multiplier = 4.3775577069027349e-05
Transform: 1 -0.0005371 0.003784 -212.5
0.00053 1 0.001875 -86.65
-0.003785 -0.001873 1 -217.8
comments: 13_HR_DTI_128axis_b3000


Image: “14_1.mif”


Dimensions: 128 x 128 x 30 x 130
Voxel size: 0.35 x 0.35 x 0.7 x ?
Data strides: [ 1 2 3 4 ]
Format: MRtrix
Data type: signed 16 bit integer (little endian)
Intensity scaling: offset = 0, multiplier = 1
Transform: 1 0 0 -22.23
0 1 0 -22.23
0 0 1 -10.15
dw_scheme: [ 130 entries ]
mrtrix_version: 005da18a

Thanks, Carlos.

1 Like

Thanks for reporting back yourself already @Carlos_Gutierrez. I was already thinking it must’ve been something as fundamental as that; it definitely didn’t seem specific to (any algorithm for) response function selection.

Hi,

I have a similar issue with a 2D rat invivo dataset.

dwi2response with tournier works and dhollander gives the empty mask error!

However, when I convert the 2D dataset into a cube (matching the 3rd dimension with the other 2) dhollander works!

Karthik

Hi @Kar,

I suppose with “2D”, you’re referring to a single slice? The problem here will be that the mask is eroded by default, and that of course in 3D, which will make your “slice mask” vanish entirely. I suppose you’re also supplying your mask manually here via the -mask option? I can imagine dwi2mask could have issues with the situation as well. Let’s say you’ve manually defined a “2D” mask, the way to make the algorithm then work properly would be to specify -erode 0, so the mask doesn’t get eroded / obliterated. :wink:

However, your makeshift solution is just as good in practice, I reckon. If you’ve got good responses from artificially creating a 3D volume out of the 2D slice, don’t bother with the above: your responses should almost certainly come out just fine!

Cheers,
Thijs

Hi Thijs,

By ‘2D’ i meant 2D acquisition acquired at 0.125*0.125 mm to save tim :wink: The data has 5 slices in the 3rd dimension with a slice thickness of 0.5 mm.

Yes, I do supply a mask derived from dwi2mask which works quite well!

I tried using the -erode 0 with dhollander but it seems to take borders for the gm response

As suggested, I can use the dhollander on the cubed data! The only concern is that I would be interpolating the 3rd dimension. Besides. I also upsample the image just before obtaining the FOD’s so that would introduce further interpolation in the 3rd dimension.

Karthik

Just to add, I was wondering why is it working fine for isotropic data & not otherwise :thinking:

Ok, so not a slice, but definitely still a “slab”. :wink: The same “issue” applies though: the default erosion of 3 voxels that is applied to the mask will obliterate your mask (since it’s only 5 voxels wide, and 3 voxels get eroded at both sides). Not eroding, well, doesn’t erode of course; and you see indeed what that step is there for: to avoid those few voxels at the edge! :slightly_smiling_face: In principle, you would probably only want to erode in 2 dimensions (and not along the 3D dimension, i.e. along the thickness of your slab). But indeed, the easy solution here is to just replicate the slab a few times, and allow erosion as normal (i.e. don’t specify -erode explicitly, so the default of 3 voxels applies). It’s not “perfect” etc…, but you’ll get the good quality response functions you’ll need anyway (and that is the goal, after all).

No worries for the sake of response function selection: the algorithm never interpolates, only selects existing voxels. You can still up-sample afterwards, just before performing the CSD step.
This does bring up another possibility for your initial problem: up-sampling the data will also create extra “slices”. However, you’d then still end up only selecting voxels from the more middle slice of your slab, if erosion does what it does by default. Plus, the algorithm would run much slower on the up-sampled data. That’s why, typically, I don’t advise to up-sample the data (be)for(e) response function selection; it works perfectly fine on non-upsampled original resolution data.

Probably the reason it didn’t work on non-isotropic data is one of the above; so unrelated to the actual voxels being non-isotropic. I’ve actually run the dhollander algorithm with success here on data of a ridiculous low quality: adult human in-vivo brain, but highly anisotropic voxels, and low spatial resolution, and a b-value of less than 1000, and only 12 gradient directions (all of this in one dataset)… still selects perfectly fine response functions for all 3 tissue types; and the voxels it selects for it, also still make full sense. :smile:

Maybe one extra note though: the erosion step, and also a dilation step that happens within the single-fibre WM response selection (currently tournier algorithm), are performed in voxel space, using voxel units! So yes, severely anisotropic voxels result in that erosion also being “anisotropic” with respect to “real” (metric) space. Again, not an issue for a whole brain, but may become one if you’re working with a more 2D-ish slab. If you’re really working with some weird or exceptional data in terms of sizes and slab-by-ness, you could of course also first get a mask, and manually (using the ROI editor) edit it yourself to perform some custom erosion around the edges (or include/exclude other bits that the masking algorithm would fail to do so), and then provide it to dwi2response dhollander via the -mask option, and (if you’ve done hands-on erosion yourself already) explicitly state -erode 0. As long as whatever is in the area of the mask is still (and only) the actual tissue and (free) water, the dhollander algorithm will tackle a wide range of scenarios.

I’m actually working on improving the dhollander algorithm even further; but I found out that it’s not as trivial as I thought it would be… so depending on whether I can make the little trick I came up with work robustly (or not), that improvement may also eventually come along (or not, should I have to admit defeat at my attempt at some point). :grin:

Ok, that’s good to know that upsampling wouldn’t be an issue. I could even upsample the data prior to response function estimation since the file size is pretty nominal (compared to the high-resolution data sets I’m used to ;))

But, I tried manually eroding the mask and response function estimation stating -erode 0 explicitly. It still takes the very edge of the brain as gm (close to the crosshair in the pic).
Capture
So, I will just upsample the data and estimate the response functions as it seems to select the voxels neatly! I would still check if the fod’s look normal :slight_smile:

It’s nice to hear that there would be another version of the dhollander response function. In general, I’m a huge fan of MRtrix and an avid user too. Thanks to the MRtrix team for all the efforts!! Btw, Is the SS3T implementation also close :wink::stuck_out_tongue_winking_eye::sunglasses:

Cheers,
Karthik

1 Like

Yeah, so what I meant by “hands-on” erosion would be to manually already exclude some of those edges from the mask, but then in the “2D” sense: just shrinking the mask without making it thinner in the third dimension. But if upsampling similarly fixes the issue, that would be fine just as well!

Sounds perfect! As long as the majority of the selected voxels are in sensible locations, you should be good. One or two “wrong” voxels don’t really hurt the result as long as many more others are correct. :slightly_smiling_face:

…depending on how well I can make my little idea generalise. I’m keen on improving it, but of course not at the cost of broad applicability. It’s hard to test though; and I need to still get it robust in the first place.

I wish! :grimacing: I may do a beta-round at some point, among a few select collaborators where we’ve been running the prototype on their data before…

1 Like