Dwi2response tax does not work and the result of dwi2response tournier is strange

Dear all,

I would like to your expert opinion on what’s can be wrong in the following case.
One sample in our dataset gave me strange result in dwi2response tournier (similar result created by dwi2response fa).

Meanwhile, it failed to run dwi2response tax.

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: Generated temporary directory: /tmp/dwi2response-tmp-FYYRBC/
Command: mrconvert /media/meng/Data/DBS/Tremor_bk/Result/Gaibler_Elisabeth/mrtrix/dwi_upsample_preproc.mif - -stride 0,0,0,1 | dwiextract - /tmp/dwi2response-tmp-FYYRBC/dwi.mif
dwi2response: Changing to temporary directory (/tmp/dwi2response-tmp-FYYRBC/)
Command: dwi2mask dwi.mif mask.mif
Command: amp2sh dwi.mif dwiSH.mif
Command: dwiextract dwi.mif shell.mif
Command: dwi2fod csd dwi.mif init_RF.txt iter0_FOD.mif -mask mask.mif
Command: fod2fixel iter0_FOD.mif -peak iter0_peaks.msf -mask mask.mif -fmls_no_thresholds
Command: fixel2voxel iter0_peaks.msf split_value iter0_amps.mif
Command: mrconvert iter0_amps.mif iter0_first_peaks.mif -coord 3 0 -axes 0,1,2
Command: mrconvert iter0_amps.mif iter0_second_peaks.mif -coord 3 1 -axes 0,1,2
Command: fixel2voxel iter0_peaks.msf split_dir iter0_all_dirs.mif
Command: mrconvert iter0_all_dirs.mif iter0_first_dir.mif -coord 3 0:2
Command: mrcalc iter0_second_peaks.mif iter0_first_peaks.mif -div iter0_peak_ratio.mif
Command: mrcalc iter0_peak_ratio.mif 0.1 -lt mask.mif -mult iter0_SF.mif
Command: sh2response dwiSH.mif iter0_SF.mif iter0_first_dir.mif iter0_RF.txt
Command: dwi2fod csd dwi.mif iter0_RF.txt iter1_FOD.mif -mask iter0_SF.mif
Command: fod2fixel iter1_FOD.mif -peak iter1_peaks.msf -mask iter0_SF.mif -fmls_no_thresholds
Command: fixel2voxel iter1_peaks.msf split_value iter1_amps.mif
Command: mrconvert iter1_amps.mif iter1_first_peaks.mif -coord 3 0 -axes 0,1,2
Command: mrconvert iter1_amps.mif iter1_second_peaks.mif -coord 3 1 -axes 0,1,2
Command: fixel2voxel iter1_peaks.msf split_dir iter1_all_dirs.mif
Command: mrconvert iter1_all_dirs.mif iter1_first_dir.mif -coord 3 0:2
Command: mrcalc iter1_second_peaks.mif iter1_first_peaks.mif -div iter1_peak_ratio.mif
Command: mrcalc iter1_peak_ratio.mif 0.1 -lt iter0_SF.mif -mult iter1_SF.mif
Traceback (most recent call last):
File “/media/meng/Data/Soft/mrtrix3/scripts/dwi2response”, line 120, in
File “/media/meng/Data/Soft/mrtrix3/scripts/src/dwi2response/tax.py”, line 88, in execute
errorMessage(‘Aborting: All voxels have been excluded from single-fibre selection’)
NameError: global name ‘errorMessage’ is not defined

Any suggestion what’s wrong with this case?

Thank you in advance.
Best, Meng

Not sure what’s going on, but using the -sf_voxels option to output the final single-fibre mask should at least give us an idea as to whether the voxels used to estimate the response make anatomical sense (they should consist of voxels containing coherent white matter).

Not too surprised about the Tax algorithm, it does fail quite regularly, which is why we’ve been moving away from it recently.

Also, given your other posts about getting the DW information, there may be issue with the gradient directions causing this issue. You’d need to make sure all this is fixed up properly to rule that out as the source of the problem.

Given the results here, as well as the presence of your other post, this is almost certainly an issue with the diffusion gradient orientations. If your image axes are permuted in any way, the gradient vectors need to be correspondingly modified.

  • The tournier algorithm gives the result it does because it relies on direct averaging of the signal in the top 300 single-fibre voxels, after rotating those signals such that the fibre direction is along the z-axis. So if the gradient directions are incorrect, that averaging process goes awry. But there’s a good chance that the selected voxels themselves (as shown by the -sf_voxels option) will in fact be appropriate.

  • The tax algorithm relies on iteratively rejecting voxels from the single-fibre selection, based on the response function estimated during the previous iteration. But if calculation of said response function from single-voxel averaging (see point above) is going awry, you end up with a bad response function, concomitantly bad FODs, and all voxels rejected from the single-fibre mask. The message isn’t quite right (it should just print ‘Aborting: All voxels have been excluded from single-fibre selection’, but I’ve failed to import the library function to print said message; I’ll fix that), but the point stands nonetheless.


@jdtournier, @rsmith, thank you all for your replies.

First of all, thi gradient table for this sample is correct. This sample comes from a routine clinic scan. And the encodings is the same as the others in the same project (I have double checked that). The data in my another post (on orientation issues) comes from a postmortem pilot study study. Sorry for the confusions.

I got some sf images using this code, dwi2response tournier -sf_voxels 100 dwi_upsample_preproc.mif dwi_response_tournier.txt -nocleanup

meng@Meng:/tmp/dwi2response-tmp-4OWD2H$ ls
command.txt             iter3_CF.mif            iter6_SF_diff.mif
cwd.txt                 iter3_first_dir.mif     iter6_SF_dilated.mif
dwi.mif                 iter3_first_peaks.mif   iter6_SF.mif
dwiSH.mif               iter3_FOD.mif           iter7_all_dirs.mif
encoding.b              iter3_peaks.msf         iter7_amps.mif
init_RF.txt             iter3_RF.txt            iter7_CF.mif
iter0_all_dirs.mif      iter3_second_peaks.mif  iter7_first_dir.mif
iter0_amps.mif          iter3_SF_diff.mif       iter7_first_peaks.mif
iter0_CF.mif            iter3_SF_dilated.mif    iter7_FOD.mif
iter0_first_dir.mif     iter3_SF.mif            iter7_peaks.msf
iter0_first_peaks.mif   iter4_all_dirs.mif      iter7_RF.txt
iter0_FOD.mif           iter4_amps.mif          iter7_second_peaks.mif
iter0_peaks.msf         iter4_CF.mif            iter7_SF_diff.mif
iter0_RF.txt            iter4_first_dir.mif     iter7_SF_dilated.mif
iter0_second_peaks.mif  iter4_first_peaks.mif   iter7_SF.mif
iter0_SF_dilated.mif    iter4_FOD.mif           iter8_all_dirs.mif
iter0_SF.mif            iter4_peaks.msf         iter8_amps.mif
iter1_all_dirs.mif      iter4_RF.txt            iter8_CF.mif
iter1_amps.mif          iter4_second_peaks.mif  iter8_first_dir.mif
iter1_CF.mif            iter4_SF_diff.mif       iter8_first_peaks.mif
iter1_first_dir.mif     iter4_SF_dilated.mif    iter8_FOD.mif
iter1_first_peaks.mif   iter4_SF.mif            iter8_peaks.msf
iter1_FOD.mif           iter5_all_dirs.mif      iter8_RF.txt
iter1_peaks.msf         iter5_amps.mif          iter8_second_peaks.mif
iter1_RF.txt            iter5_CF.mif            iter8_SF_diff.mif
iter1_second_peaks.mif  iter5_first_dir.mif     iter8_SF_dilated.mif
iter1_SF_diff.mif       iter5_first_peaks.mif   iter8_SF.mif
iter1_SF_dilated.mif    iter5_FOD.mif           iter9_all_dirs.mif
iter1_SF.mif            iter5_peaks.msf         iter9_amps.mif
iter2_all_dirs.mif      iter5_RF.txt            iter9_CF.mif
iter2_amps.mif          iter5_second_peaks.mif  iter9_first_dir.mif
iter2_CF.mif            iter5_SF_diff.mif       iter9_first_peaks.mif
iter2_first_dir.mif     iter5_SF_dilated.mif    iter9_FOD.mif
iter2_first_peaks.mif   iter5_SF.mif            iter9_peaks.msf
iter2_FOD.mif           iter6_all_dirs.mif      iter9_RF.txt
iter2_peaks.msf         iter6_amps.mif          iter9_second_peaks.mif
iter2_RF.txt            iter6_CF.mif            iter9_SF_diff.mif
iter2_second_peaks.mif  iter6_first_dir.mif     iter9_SF_dilated.mif
iter2_SF_diff.mif       iter6_first_peaks.mif   log.txt
iter2_SF_dilated.mif    iter6_FOD.mif           mask.mif
iter2_SF.mif            iter6_peaks.msf         response.txt
iter3_all_dirs.mif      iter6_RF.txt            
iter3_amps.mif          iter6_second_peaks.mif  voxels.mif

The location of sf.mif is very unexpected to me. Here is the overlay of iter0_SF.mif on dwiSH.mif. It seems the script choose an unreliable region to estimate FOD response.

Sorry, I gave you the wrong instructions! The -sf_voxels option specifies how many voxels to include in the estimation, and by default this is 300. I don’t recommend you change that - if anything, you’ll want to increase that number if these data were up-sampled. What I wanted was the -voxels option, which outputs the final single-fibre mask.

What you’re showing here is the single-fibre mask used in the first iteration - what we want is the one used in the final iteration (presumably voxels.mif in your listing). Otherwise this would indeed be far from ideal…

By the way, are your data phase-encoded along the inferior-superior axis? The distortions seem very unusual to me…

Thank you. I rerun the dwi2response with the followings: dwi2response tournier -voxel voxels.mif dwi_upsample_preproc.mif dwi_response_tournier.txt -nocleanup
Here is the overlay of voxels.mif on dwiSH.mif.

I don’t think we have phase-encoded along the inferior-superior axis. But there is stripes noises in the original images. I probably need to drop this sample out.

OK, the voxels in the region with high susceptibility-induced distortions may throw off the estimation - but I wouldn’t have expected it to throw it off that much, provided most of the voxels are in sensible locations (which seems to be the case, by and large). What you could try if you want is to use a more conservative initial mask than what dwi2response computes by default (it uses a simple dwi2mask call), by eroding what would have been the default mask more aggressively and passing that through to dwi2response:

dwi2mask dwi_upsample_preproc.mif - | maskfilter - erode -npass 5 eroded_mask.mif
dwi2response tournier -mask eroded_mask.mif -voxels voxels.mif dwi_upsample_preproc.mif dwi_response_tournier.txt

This erodes 5 voxels off the edge of the mask, hopefully that’ll get rid of anything near the edge of the brain. You may want to experiment with higher values if those problematic voxels remain included in the final selection…

If that doesn’t fix it, then there’s a good chance your stripe artefacts are too severe, you would need to focus on removing those before proceeding. You can use something like:

mrconvert -coord 3 1:10,12:end dwi.mif dwi_reduced.mif

to remove volume 11, for instance. mrconvert will take care of updating the gradient table to match.

Hi @jdtournier, @rsmith, the problem is solved. You are right, it’s data issue – the last two volumes were missing.
Thank you for your patience! I do learn a lot on MRTrix3 these days.