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.
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.
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.
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
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_amps.mif iter6_second_peaks.mif voxels.mif
The location of sf.mif is very unexpected to me. Here is the overlay of
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
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.