Dwi2response dhollander with single-shell data

Hi,

I am trying to obtain 3-tissue responses using the dhollander algorithm on a clinical single-shell dataset (b=1000,64 directions, 2mm iso resolution), and will follow the ss3t-csd pipeline to continue my analysis.

Here, I used B0-T2 registration to correct the distortion.
Although I am able to obtain the 3-tissue responses for most of the cases, there are still a fraction of data that are not able to be processed (the errors below).
From the intermediate files, I noticed that the refined_gm.mif is mostly 0 in the image.

I am wondering if there are certain parameter that I can tune to successfully process these images? or it is really just a death sentence.

Thank you!

dwi2response: [ERROR] mrcalc refined_gm.mif safe_sdm.mif 3.14924455 -subtract -abs 1 -add 0 -if - | mrthreshold - - -bottom 0 -ignorezero | mrcalc refined_gm.mif - 0 -if - -dataty

pe bit | mrconvert - voxels_gm.mif -axes 0,1,2 (dhollander.py:208)

dwi2response: [ERROR] Information from failed command:
dwi2response:
              mrcalc: computing: (refined_gm.mif ? (|(safe_sdm.mif - 3.14924)| + 1) : 0)...  [=====================================================================================]
              mrthreshold: [ERROR] value supplied for option "bottom" is out of bounds (valid range: 1 to 9223372036854775807, value supplied: 0)
              mrcalc: [ERROR] Could not interpret string "-" as either an image path or a numerical value
              mrconvert: [ERROR] no filename supplied to standard input (broken pipe?)
              mrconvert: [ERROR] error opening image "-"

1 Like

Hi @xiaobird,

That should generally work out of the box; so something unusual must be going on. Let’s see if we can figure it out.

Good to know. Generally, I don’t see a problem per se with that, but you never know it might’ve done something strange to the data!

Yep, that shouldn’t be happening. But the problem likely originates before it gets to that point. But yeah, this point itself does reveal something has gone wrong in an earlier step.

Maybe, but I’m guessing this might have to do with another issue in the data. If this is a regular clinical dataset, then this should absolutely not be a death sentence, we just need to figure out what went wrong and why.

…so thanks for posting that output, but: from that I can only see that indeed at some point refined_gm.mif has dropped very close to 0 voxels, but I think not entirely to 0 at least, due to what I can see in that output. That’s weird for sure. Let’s follow up on 2 fronts to help you out. Can you post the entire output you get from dwi2response dhollander at the command line? That might help me get a bit of insight into the relevant numbers for your data and how the algorithm runs through them. It might give me the right hint to figure out what’s the issue even.

Apart from that, I’ll be in touch privately; maybe it’s useful to arrange a little data transfer of one dataset/subject where you encounter exactly this. I can then take a look at it myself (when I find some time; pretty busy with project meetings), see if I can reproduce the scenario and investigate.

Cheers,
Thijs

Hi,

I forgot to mention that I was running dwi2response on up-sampled dwi data (from 2mm to 1mm using mrresize). I was using the upsampled version for some other experiments, and thought dwi2response should work despite of the resolution difference.

I re-run dwi2response at 2mm resolution on one failed case, and it worked. Maybe the upsampling introduced some irregular values to the images. I will try these out on all failed subjects and see.

Thanks a lot!

Best,
Xiao

Ok, to report back and close this off: we’ve looked into it and solved it offline. The full output until before the earlier reported error was:

dwi2response: -------
dwi2response: 2 unique b-value(s) detected: 0,1000 with 1,64 volumes
dwi2response: -------
dwi2response: Preparation:
dwi2response: * Eroding brain mask by 3 pass(es)...
dwi2response:   [ mask: 1499226 -> 1280147 ]
dwi2response: * Computing signal decay metric (SDM):
dwi2response:  * b=0...
dwi2response:  * b=1000...
dwi2response: * Removing erroneous voxels from mask and correcting SDM...
dwi2response:   [ mask: 1280147 -> 1279835 ]
dwi2response: -------
dwi2response: Crude segmentation:
dwi2response: * Crude WM versus GM-CSF separation (at FA=0.2)...
dwi2response:   [ 1279835 -> 1279785 (WM) & 50 (GM-CSF) ]
dwi2response: * Crude GM versus CSF separation...
dwi2response:   [ 50 -> 17 (GM) & 33 (CSF) ]
dwi2response: -------
dwi2response: Refined segmentation:
dwi2response: * Refining WM...
dwi2response:   [ WM: 1279785 -> 1003999 ]
dwi2response: * Refining GM...
dwi2response:   [ GM: 17 -> 10 ]
dwi2response: * Refining CSF...
dwi2response:   [ CSF: 33 -> 3064 ]
dwi2response: -------
dwi2response: Final voxel selection and response function estimation:
dwi2response: * CSF:
dwi2response:  * Selecting final voxels (10.0% of refined CSF)...
dwi2response:    [ CSF: 3064 -> 306 ]
dwi2response:  * Estimating response function...
dwi2response: * GM:
dwi2response:  * Selecting final voxels (2.0% of refined GM)...

…which I could exactly reproduce on my end on the same upsampled data. The data also came with a nice pre-computed mask. No worries there.

The initial clue here was in:

dwi2response: * Crude WM versus GM-CSF separation (at FA=0.2)...
dwi2response:   [ 1279835 -> 1279785 (WM) & 50 (GM-CSF) ]

This is the step that separates WM and non-WM crudely initially (which is then later refined, etc…). In any case, the separation here looks very much out of balance, indicating that way too many voxels had an FA value above 0.2. However, looking into the upsampled data itself then, we realised a number of volumes at the end where all zeros. E.g. mrstats of the upsampled data within the brain mask gives:

      volume       mean     median        std        min        max      count
       [ 0 ]    494.738 388.090057    301.132   -117.253    3791.53    1499226
       [ 1 ]     151.67 160.248291    56.3384   -19.7745    863.292    1499226
       [ 2 ]    152.909 161.636276    57.8888   -28.6813    802.051    1499226
       [ 3 ]    154.691 163.452454    56.8032   -16.6526    912.083    1499226
       [ 4 ]    148.126 156.810303    56.2617    -25.557    811.781    1499226
       [ 5 ]    154.134 163.085571     56.868   -22.2393    865.215    1499226
       [ 6 ]    151.709 159.929688     55.836   -17.5558    851.449    1499226
       [ 7 ]    152.727 161.113647    56.5254   -18.3122     897.37    1499226
       [ 8 ]     152.41 161.774902    57.0294   -29.8998    868.062    1499226
       [ 9 ]    153.687 162.161835    56.4684     -13.26    884.924    1499226
      [ 10 ]    149.487 158.265808    56.3392   -17.7466    828.059    1499226
      [ 11 ]    149.845 158.797089    56.2865   -14.8034    847.888    1499226
      [ 12 ]    149.179 157.764954    56.2122   -15.9005    890.661    1499226
      [ 13 ]    152.913 161.815659    55.9334   -18.3439    832.262    1499226
      [ 14 ]    146.906 155.168762    55.5958   -29.5223    857.216    1499226
      [ 15 ]    149.081 157.652863    55.3762   -5.74525    933.559    1499226
      [ 16 ]    148.969 157.750107    55.4795   -23.2521    830.679    1499226
      [ 17 ]    150.163 158.642563    56.0775   -7.94573    876.392    1499226
      [ 18 ]    154.027 162.568237    56.2677   -11.8833    950.344    1499226
      [ 19 ]    150.185 159.171448    57.6265   -24.6021    1002.39    1499226
      [ 20 ]    151.181 159.403931    55.9195   -31.5592     881.49    1499226
      [ 21 ]    153.933 163.305298     56.899   -29.2348    883.925    1499226
      [ 22 ]    154.261 162.684387    56.3751   -30.3851    915.124    1499226
      [ 23 ]    149.102 157.824921    55.5241   -26.4253    819.713    1499226
      [ 24 ]    153.055 162.450211    56.2618   -23.6808    894.056    1499226
      [ 25 ]    156.925 166.234192    56.5952   -18.5176    977.849    1499226
      [ 26 ]    152.641 161.784088    57.4353   -24.1844    836.776    1499226
      [ 27 ]    150.062 158.706757     55.683   -29.5699    828.574    1499226
      [ 28 ]    152.358 161.353195    55.8603   -27.0974    905.059    1499226
      [ 29 ]    149.982 159.510071    56.4315   -14.0581    821.621    1499226
      [ 30 ]    149.213  158.12854    55.8289   -16.6214    892.409    1499226
      [ 31 ]    154.774 163.897568    55.9401   -19.6115    921.185    1499226
      [ 32 ]    155.856 164.580597    56.0664   -19.2309    933.737    1499226
      [ 33 ]     159.19 168.703522    57.0455   -26.0583    968.283    1499226
      [ 34 ]    155.719 165.070358      55.71   -18.1726    897.377    1499226
      [ 35 ]    153.395 162.080444    55.4853   -22.4217    868.023    1499226
      [ 36 ]    149.216 157.887573    54.7311    -19.353    860.021    1499226
      [ 37 ]     145.98 154.940552     55.168   -23.2812    830.195    1499226
      [ 38 ]    145.852 154.467316    54.4531   -17.9308    838.842    1499226
      [ 39 ]    151.661 160.197876    56.7448   -31.7362    929.559    1499226
      [ 40 ]    148.803 156.756592    55.2557   -19.2926    898.824    1499226
      [ 41 ]    148.465 157.243225    55.4671   -16.2449    830.806    1499226
      [ 42 ]    152.409  161.18512    57.4395   -20.1523    786.045    1499226
      [ 43 ]    145.729 153.784531    54.2876   -18.8578    904.257    1499226
      [ 44 ]    150.248 159.324127     56.789   -19.3947    907.398    1499226
      [ 45 ]    150.504 159.493057    56.1103   -26.5536    875.271    1499226
      [ 46 ]    155.168 164.796051    55.9632   -12.6542     879.15    1499226
      [ 47 ]    155.245  164.57254    56.7716    -22.087    824.504    1499226
      [ 48 ]    148.829 157.496735    55.4244   -27.8474    840.875    1499226
      [ 49 ]    71.2541          0    86.0476   -35.6985    866.412    1499226
      [ 50 ]          0          0          0          0          0    1499226
      [ 51 ]          0          0          0          0          0    1499226
      [ 52 ]          0          0          0          0          0    1499226
      [ 53 ]          0          0          0          0          0    1499226
      [ 54 ]          0          0          0          0          0    1499226
      [ 55 ]          0          0          0          0          0    1499226
      [ 56 ]          0          0          0          0          0    1499226
      [ 57 ]          0          0          0          0          0    1499226
      [ 58 ]          0          0          0          0          0    1499226
      [ 59 ]          0          0          0          0          0    1499226
      [ 60 ]          0          0          0          0          0    1499226
      [ 61 ]          0          0          0          0          0    1499226
      [ 62 ]          0          0          0          0          0    1499226
      [ 63 ]          0          0          0          0          0    1499226
      [ 64 ]          0          0          0          0          0    1499226

Volume 49 is already missing data, and volumes 50-64 are all zero. This naturally then caused the artificially increased FA values.

In the end dwi2response dhollander ran without issue on the non-upsampled data. On the upsampled data, it would’ve run alright too if those volumes weren’t zeroed. In general though, I advise people to run dwi2response dhollander on non-upsampled data: much faster, and slightly better designed for it.

So in the end, it wasn’t a problem with dwi2response dhollander, but with the upsampling step. I’ve seen this happen before for others even relatively recently, in particular on various cluster environments (which was the case here too), maybe with complicated network storage setups or file systems, where it ends up failing to write all the output, leaving zeros in the data.

1 Like

Hi Thijs,
I also encountered similar problems, when I computed the 3-tissue responses using the dhollander algorithm on a clinical single-shell dataset (b=1000,64 directions, 1×1×2mm )。


Could you please deal with the problems?
thank you so much
Liang

Hi Liang,

You appear to have a 4D image somewhere where a 3D image is expected; perhaps the mask image you provided is 4D? Try running mrinfo on your input images.

Cheers
Rob

Hi Rob,
Thank you for your reply and the problem was indeed in the brainmask. I have solved this problem.
Cheers,
Liang

Hello @ThijsDhollander

I have also the similar dataset ( one and only one b=0 volume and 64 volumes with b=1000 ) and have a suspicion regarding the use of the more suitable algorithm for estimating the response function using dwi2reponse (dhollander or tournier). As I read in the manual (Response function estimation — MRtrix 3.0 documentation), dhollander algorithm can be used for multi b-value (single-shell + b=0, or multi-shell) response functions but tournier is used for single b-value (single-shell).

Any suggestions there?

thanks in advance to your helps

bests,

-Milad

Hi Milad,

We may need to revise the documentation text regarding both response function estimation and CSD.

The text in the response function document refers specifically to single-tissue CSD, not to single-shell data. However even for single-shell data, we typically recommend use of the multi-tissue CSD algorithm with two tissues based on the two unique b-values. So unless one has a strong justification for strictly doing single-tissue CSD and using only one b-value, we would typically suggest using dwi2response dhollander across the board; and in the case of “single-shell” DWI data (b=0 data plus one b>0 shell), still using the dwi2fod mamt_csd algorithm.

The possible exception, as I had contemplated in a recent post, was the case where only one b=0 volume is available. My concern there was that any potential noise or artifact in that solitary b=0 volume would have considerable influence on the resulting deconvolution results (whereas if multiple b=0 volumes are available, solitary issues may get diluted across those volumes, depending on the nature of the artifact). In this particular scenario, there’s some chance that it may be preferable to not use the b=0 volume; but this is entirely speculative on my part. In such a scenario, using something like dwi2response tournier is more convenient, because the derived response function contains on the coefficients for the b>0 data, whereas if one uses dwi2response dhollander the WM-like response function also contains coefficients for the b=0 “shell”. But if you know what you’re doing, you could still in this scenario use dwi2response dhollander, strip out the b=0 coefficients, and then either use dwi2fod csd, or use dwi2fod msmt_csd providing only b>0 data as the input (so that the number of shells will match between the input data and the response function).

Cheers
Rob

Hi Rob

Thanks for your quick replay. As I find out from your comments, when two unique b-values ( b=0 pluse B>0) are available, two response functions for two tissues ( WM and CSF) can be obtained using dwi2response dhollander. I tested this strategy as follow:

dwi2response dhollander dwi_den_preproc_unbiased.mif -fslgrad  sub.bvec sub.bval wm.txt  csf.txt  -voxels voxels.mif -fa 0.2  

but I encountered with the following error massage:

Error: the following arguments are required: out_csf

how can i solve it? I must mention that with running the dwi2reponsd with 3 tissues in the same dwi data, I can get the response functions:

dwi2response dhollander dwi_den_preproc_unbiased.mif -fslgrad  sub.bvec sub-40451_session1_dwi1.bval wm_rf.txt gm_rf.txt csf_rf.txt  -voxels voxels_rf.mif -fa 0.2    
dwi2response: 
dwi2response: Note that this script makes use of commands / algorithms that have relevant articles for citation. Please consult the help page (-help option) for more information.
dwi2response: 
dwi2response: Generated scratch directory: /media/sn/5EB4CC35B4CC1207/run/sub-40451/dwi/session1/dwi2response-tmp-KLRWGK/
dwi2response: Importing DWI data (/media/sn/5EB4CC35B4CC1207/run/sub-40451/dwi/session1/dwi_den_preproc_unbiased.mif)...
dwi2response: Changing to scratch directory (/media/sn/5EB4CC35B4CC1207/run/sub-40451/dwi/session1/dwi2response-tmp-KLRWGK/)
dwi2response: Computing brain mask (dwi2mask)...
dwi2response: -------
dwi2response: 2 unique b-value(s) detected: 0,1000 with 1,64 volumes
dwi2response: -------
dwi2response: Preparation:
dwi2response: * Eroding brain mask by 3 pass(es)...
dwi2response:   [ mask: 171074 -> 129179 ]
dwi2response: * Computing signal decay metric (SDM):
dwi2response:  * b=0...
dwi2response:  * b=1000...
dwi2response: * Removing erroneous voxels from mask and correcting SDM...
dwi2response:   [ mask: 129179 -> 129043 ]
dwi2response: -------
dwi2response: Crude segmentation:
dwi2response: * Crude WM versus GM-CSF separation (at FA=0.2)...
dwi2response:   [ 129043 -> 72867 (WM) & 56176 (GM-CSF) ]
dwi2response: * Crude GM versus CSF separation...
dwi2response:   [ 56176 -> 39813 (GM) & 16363 (CSF) ]
dwi2response: -------
dwi2response: Refined segmentation:
dwi2response: * Refining WM...
dwi2response:   [ WM: 72867 -> 67772 ]
dwi2response: * Refining GM...
dwi2response:   [ GM: 39813 -> 23382 ]
dwi2response: * Refining CSF...
dwi2response:   [ CSF: 16363 -> 10374 ]
dwi2response: -------
dwi2response: Final voxel selection and response function estimation:
dwi2response: * CSF:
dwi2response:  * Selecting final voxels (10.0% of refined CSF)...
dwi2response:    [ CSF: 10374 -> 1037 ]
dwi2response:  * Estimating response function...
dwi2response: * GM:
dwi2response:  * Selecting final voxels (2.0% of refined GM)...
dwi2response:    [ GM: 23382 -> 468 ]
dwi2response:  * Estimating response function...
dwi2response: * Single-fibre WM:
dwi2response:  * Selecting final voxels (0.5% of refined WM)...
dwi2response:    Selecting WM single-fibre voxels using built-in (Dhollander et al., 2019) algorithm
dwi2response:    [ WM: 67772 -> 339 (single-fibre) ]
dwi2response:  * Estimating response function...
dwi2response: -------
dwi2response: Generating outputs...
dwi2response: -------
dwi2response: Changing back to original directory (/media/sn/5EB4CC35B4CC1207/run/sub-40451/dwi/session1)
dwi2response: Deleting scratch directory (/media/sn/5EB4CC35B4CC1207/run/sub-40451/dwi/session1/dwi2response-tmp-KLRWGK/)



regarding the dwi2fod msmt for single tissue, I use the following codes according to your comments(Fixel-based analysis results - #3 by Andrea_Morelli) for single shell data :

dwiextract -shell 1000 dwi_den_preproc_unbiased.mif dwi_ss.mif

tail -1 wm.txt > wm_ss.txt

dwi2fod msmt_csd dwi_ss.mif -mask mask.mif wm_ss.txt wmfod.mif

Am i doing this step correctly?

bests,

-Milad

when two unique b-values ( b=0 pluse B>0) are available, two response functions for two tissues ( WM and CSF) can be obtained using dwi2response dhollander.

dwi2response dhollander always estimates, and outputs, three tissue response functions, regardless of the number of unique b-values in the data. The number of unique b-values restricts how many response functions you can use at the deconvolution stage.

This is why you get the error message:

Error: the following arguments are required: out_csf

If we look at your usage:

dwi2response dhollander <options> dwi_den_preproc_unbiased.mif wm.txt csf.txt

compared to the command help page:

 dwi2response dhollander [ options ] input out_sfwm out_gm out_csf

The script demands that you specify three output paths, but you have specified only two.


dwiextract -shell 1000 dwi_den_preproc_unbiased.mif dwi_ss.mif
tail -1 wm.txt > wm_ss.txt
dwi2fod msmt_csd dwi_ss.mif -mask mask.mif wm_ss.txt wmfod.mif

Am i doing this step correctly?

If you wish to discard the b=0 data and utilise only the b=1000 data, in a single-tissue deconvolution, then yes, that should work. But as I said, I’m not strictly advocating for that approach; I’m merely pointing out that in the instance where one has only one b=0 image rather than several, it’s not strictly guaranteed that a two-shell two-tissue deconvolution will be preferable to a single-shell single-tissue deconvolution. Indeed making that assessment objectively would be quite difficult.

Cheers
Rob

Hi Rob

Thanks for your detailed response.

In this particular scenario, there’s some chance that it may be preferable to not use the b=0 volume; but this is entirely speculative on my part. In such a scenario, using something like dwi2response tournier is more convenient, because the derived response function contains on the coefficients for the b>0 data,
Sorry i forgot to read this paragraph accurately. Since my data have a solitary b0 volume ( one b0 volume and 64 volume with b=1000), I’ve decided to use the tournier algorithm for creating the response function in order to avoid the possible error in estimating the response function.

If you wish to discard the b=0 data and utilise only the b=1000 data, in a single-tissue deconvolution, then yes, that should work. But as I said, I’m not strictly advocating for that approach; I’m merely pointing out that in the instance where one has only one b=0 image rather than several, it’s not strictly guaranteed that a two-shell two-tissue deconvolution will be preferable to a single-shell single-tissue deconvolution.

Considering my dwi data, I use this approach for dwi2fod msmt to gain the benefit of this algorithm compare with csd as you mentioned previously (Fixel-based analysis results)