Hi MRtrix3 experts,
I want to evaluate the performance of ‘dwi2response dhollander’ script on single-shell data with low b-values. I used the classic Fiber cup datasets, which include 3 different b-value, e.g. 650, 1500, 2000. My command looks like:
dwi2response dhollander …/dwi/acq-averaged_b-650.nii.gz RF_WM.txt RF_GM.txt RF_CSF.txt -grad gradients_650.txt
it outputs:
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 temporary directory: /Users/hywu/data/fibercup/result/dwi2response-tmp-BHX0GW/
Command: mrconvert /Users/hywu/data/fibercup/dwi/acq-averaged_b-650.nii.gz /Users/hywu/data/fibercup/result/dwi2response-tmp-BHX0GW/dwi.mif -stride 0,0,0,1 -grad /Users/hywu/data/fibercup/result/gradients_650.txt
dwi2response: Changing to temporary directory (/Users/hywu/data/fibercup/result/dwi2response-tmp-BHX0GW/)
Command: dwi2mask dwi.mif mask.mif
dwi2response: 2 unique b-value(s) detected: 0,650 with 1,64 volumes.
Command: maskfilter mask.mif erode eroded_mask.mif -npass 3
Command: dwiextract dwi.mif -shell 0 - | mrmath - mean mean_b0.mif -axis 3
Command: mrcalc mean_b0.mif -finite mean_b0.mif 0 -if 0 -le err_b0.mif -datatype bit
Command: dwiextract dwi.mif -shell 650 - | mrmath - mean mean_b650.mif -axis 3
Command: mrcalc mean_b650.mif -finite mean_b650.mif 0 -if 0 -le err_b650.mif -datatype bit
Command: mrcalc mean_b0.mif mean_b650.mif -divide -log sdm_b650.mif
Command: mrcalc sdm_b650.mif 64 -mult 64 -divide full_sdm.mif
Command: mrcalc full_sdm.mif -finite full_sdm.mif 0 -if 0 -le err_sdm.mif -datatype bit
Command: mrcalc err_b0.mif err_b650.mif -add err_sdm.mif -add 0 eroded_mask.mif -if safe_mask.mif -datatype bit
Command: mrcalc safe_mask.mif full_sdm.mif 0 -if 10 -min safe_sdm.mif
Command: dwi2tensor dwi.mif - -mask safe_mask.mif | tensor2metric - -fa safe_fa.mif -vector safe_vecs.mif -modulate none -mask safe_mask.mif
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 (dhollander.py:105)
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 “-”
Any suggestions will be deeply appreciated.
Best,
Haiyong
I resolved this issue.
The reason may lie in that the Fiber Cup dataset is ‘2D’ data, so the script cannot produce a mask for it. I duplicated the data at z-axis to achieve a ‘3D’ cube, then it worked well.
Many thanks to the MRtrix3 experts for providing such a powerful DWI tool, and always explaining kinds of naive problems patiently.
Hi @Haiyong_Wu,
Thanks for giving the dhollander algorithm a go. Even though you technically resolved the issue that was causing the error, I should note that the dhollander algorithm is not really designed for something like that particular fiber cup phantom. It’s really designed for more “realistic” datasets that have 3 actual tissue types which relate in certain specific ways to each other (such as white matter, grey matter and CSF do). It does generalise surprisingly well in those cases, across acquisition settings (b-values, etc…), pathologies (even gross ones), ex vivo under certain conditions, and various animal models (again under certain conditions). For the latter two challenging cases (ex vivo; animal), there may be a need to tune one particular parameter: the -fa option to the dhollander algorithm. In those cases, the anisotropy in the dataset can sometimes be quite low compared to what typically goes in an in vivo human. The purpose of the -fa parameter is to perform a very rough separation of WM versus GM/CSF. It doesn’t have to be perfect (and it never will, since some extreme crossings, especially 3-way crossings, have very low FA), but if the result of that step roughly indicates the WM-GM boundary, the rest of the algorithm takes cares of the details automatically.
In case of the fiber cup phantom, you’re really only after a single “WM” (fiber cup synthetic fibers, that is) response, there’s noting else in the phantom. Evaluating the performance of the dhollander algorithm in this setting is a bit like evaluating the performance of a small aircraft to drive on a (car) highway… if that makes sense (well, the point is: it doesn’t really). 
If you set the -fa parameter to something decent (the default may be ok for that phantom, but I’m not sure), you will probably end up with a decent multi-shell “WM” response though; but the GM and CSF ones are almost certainly meaningless with respect to the dataset.
I hope that all makes some sense. 
Cheers,
Thijs
Maybe to add a bit to your specific purpose: you could just as well evaluate the tournier algorithm in this setting, it’ll do pretty much the same thing (given you had a decent -fa setting to dhollander before). In my experience, tournier works decently on the fiber cup phantom, even at the lowest b-values.
If you want to evaluate the dhollander algorithm at low b-values, it is probably best done on an in vivo human dataset. I will hopefully soon get to publishing this algorithm decently (apart from the original abstract), where I will definitely include a comprehensive analysis at low b-values as well. No worries by the way: I’ve tested it in some pretty extreme conditions (things like b-values below 1000, only 12 gradient directions, and wildly anisotropic voxel sizes; all combined in one dataset). My ultimate focus for these basic steps is always robustness; to provide users with something they can (“safely”) rely upon. If we ever aim to bring these things into a clinical reality, robustness is certainly key. I wouldn’t want to see the thing randomly mess up if real patient(s) are involved.
Hi @ThijsDhollander,
I am working on the DWI data from the infants’ brain, with low b-value and low SNR. Moreover, their WM, GM, and CSF are not distinguishable easily.
I really appreciate your detail explanation, and look forward to reading your new paper.
Best,
Haiyong
In the mean time, the abstract may help give some insights into how the algorithm works (and what it targets and aims to find in a dataset): https://www.researchgate.net/publication/307863133_Unsupervised_3-tissue_response_function_estimation_from_single-shell_or_multi-shell_diffusion_MR_data_without_a_co-registered_T1_image
The majority of the text is quite boring (mostly listing of all individual steps for reference/documentation), but the figures on the right give some qualitative insight in what the algorithm filters the data down to, to extract responses from. If you want to see which final voxels get selected, give the -voxels option a go: the output from that one can be opened (or overlaid) directly in mrview, and will use the same colours (red = CSF, green = GM, blue = single-fibre WM).
I have learned much from this abstract and your script at Github. I prefer the form of the abstract since it is not boring but able to concisely illustrate a new idea. 
Thank you @ThijsDhollander,
1 Like
