Hi,
I am doing FBA analysis on multiple subjects and I need to use average response function to estimate the FOD. If I use ‘tournier’ algorithm as per the read document then I get error in ‘dwi2fod msmt_csd’ during FOD estimation step as below:
dwi2fod: [ERROR] number of rows in response function must match number of b-value shells
Thanks
Sarbani
The tournier approach (feels weird writing that…) provides a single-shell response, and you’re trying to run a multi-shell CSD with that single-shell response. You’d need to use the msmt_5tt or the shiny new dhollander algorithm instead (although the latter is so new, I haven’t even seen it yet - it’s being presented this week at the ISMRM diffusion workshop…). Have a look at the Response function estimation and Multi-tissue constrained spherical deconvolution pages for further details.
The tournier approach (feels weird writing that…)
Yes, I never had the intention of response function estimation algorithms being named after authors; I just couldn’t figure out any other way of logically separating your approach from Chantal’s given their similarity 
dwi2fod: [ERROR] number of rows in response function must match number of b-value shells
Because dwi2fod msmt_csd is a multi-shell algorithm, it will automatically try to use data from all b-value shells (including b=0 volumes, which we do consider to be a ‘shell’). However the tournier algorithm is a purely single-shell algorithm, and so only provides response function coefficients for one shell. So even if your acquired data is ‘single-shell’, dwi2fod msmt_csd sees two ‘shells’ in your image data, but only one in your response function.
I suspect what you’re looking for is:
dwi2fod msmt_csd ... -shell bvalue
(replace ‘bvalue’ with the b-value of your outer shell)
This will instruct dwi2fod to essentially discard all image volumes that don’t belong to your b-value of interest, and then perform deconvolution of the image data in that one shell with your single-shell response function. This step is not required with dwi2fod csd because it is intrinsically a single-shell algorithm, and therefore will automatically select only the outer shell.
Hi Sarbani,
Step 3 in the FBA instructions suggests using the msmt_csd algorithm for even single shell data because we believe the hard non-negativity constraint works better than the ‘soft’ constraint in the original CSD paper.
The step 3 instructions have currently been written as if you only have single shell data (this involves extracting the DW volumes only with dwiextract, so that the msmt algorithm does not complain when there are no b=0 response function values supplied). However, if you have multi-shell data I would highly recommend performing a proper multi-tissue CSD to improve the FOD estimates in voxels with partial volume. Cleaning up these FODs is also likely to benefit downstream registration. I will update the documentation to have different steps for single-shell and multi-shell data ASAP.
If you do have multi-shell data, I would definitely recommend giving the new dwi2response dhollander method a go (abstract here) to get the response functions. This is quite a new option so you may have to update MRtrix. You will then need to average each WM, GM, CSF response across groups. The average_response script does not currently handle multi-shell response functions, however I’ll update that tomorrow.
Cheers,
Dave
Thanks Dave, Robert and Jdtournier for your valuable inputs!
Initially I thought to use the msmt_5tt algorithm for dwi2response, but now as you said I will try with dwi2response dhollander method and see what’s coming out.
Many thanks!
Sarbani
Hi @Sarbani ,
As @Dave mentioned, it really depends on what data you’re working with. I can see how some confusion may have happened due to the tutorial being on single-shell data, but still using the msmt_csd algorithm; as Dave already mentioned, we did that at the time for the purpose of using the hard constraints in that algorithm (which are not (yet) offered as part of the classical single-shell single-tissue CSD algorithm option).
If you do have multi-shell data though, definitely give the dhollander algorithm a go; and in that case, I’d love to hear about your experiences with it! 
If you’re sitting on single-shell data, go with the tutorial for now; but within some time, the SS3T-CSD (single-shell 3-tissue CSD) wil become available in MRtrix as well; which should allow you to get multi-tissue results from single-shell data too (as long as your b-value is reasonably high).
Cheers,
Thijs
Hi,
I generate the response function with dhollander algorithm and then I run dwi2fod msmt_csd, but I got same error:
dwi2fod: [ERROR] number of rows in response function must match number of b-value shells
Then I play safe with response function and try msmt_5tt algorithm in dwi2response command, but still the error is there while executing dwi2fod.
dwi2fod: [ERROR] number of rows in response function must match number of b-value shells
The value in my average response function is as below:
avg_out_csf.txt ==9771.80513
avg_out_gm.txt ==5144.05102
avg_out_sfwm ==4066.94191 0.00000 0.00000 0.00000 0.00000
My commands are as below:
(Here the algorithm is msmt_5tt as this one is the latest one, but before that I tried dhollander, but got same error in dwi2fod)
command: dwi2response msmt_5tt /home/koushik/HCP_Data/2016-09-11/output_normalised_dwi/156637_DWI.mif /home/koushik/HCP_Data/2016-09-11/input_dwi/156637_input/5TT.mif /home/koushik/HCP_Data/2016-09-11/response_function_dir/156637_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/156637_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/156637_out_csf.txt -voxels /home/koushik/HCP_Data/2016-09-11/response_function_dir/156637_voxel.mif
Command: average_response /home/koushik/HCP_Data/2016-09-11/response_function_dir/397760_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/149337_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/499566_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/161731_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/414229_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/151627_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/201111_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/192540_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/214423_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/135932_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/133928_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/221319_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/212318_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/136833_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/156637_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/avg_out_sfwm.txt
Command: average_response /home/koushik/HCP_Data/2016-09-11/response_function_dir/397760_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/149337_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/499566_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/161731_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/414229_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/151627_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/201111_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/192540_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/214423_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/135932_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/133928_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/221319_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/212318_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/136833_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/156637_out_gm.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/avg_out_gm.txt
Command: average_response /home/koushik/HCP_Data/2016-09-11/response_function_dir/397760_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/149337_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/499566_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/161731_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/414229_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/151627_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/201111_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/192540_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/214423_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/135932_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/133928_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/221319_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/212318_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/136833_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/156637_out_csf.txt /home/koushik/HCP_Data/2016-09-11/response_function_dir/avg_out_csf.txt
Command: dwi2fod msmt_csd /home/koushik/HCP_Data/2016-09-11/output_normalised_dwi/397760_DWI.mif /home/koushik/HCP_Data/2016-09-11/response_function_dir/avg_out_sfwm.txt /home/koushik/HCP_Data/2016-09-11/output_fod_dir/397760_WM_fod.mif /home/koushik/HCP_Data/2016-09-11/response_function_dir/avg_out_gm.txt /home/koushik/HCP_Data/2016-09-11/output_fod_dir/397760_GM_fod.mif /home/koushik/HCP_Data/2016-09-11/response_function_dir/avg_out_csf.txt /home/koushik/HCP_Data/2016-09-11/output_fod_dir/397760_CSF_fod.mif -mask /home/koushik/HCP_Data/2016-09-11/input_brain_mask/397760nodif_brain_mask.nii.gz
Many Thanks!
Sarbani
I think the problem is that you are using the “average_response” script with multi-shell responses (that should have one row per shell) but this script calculates and saves the average of the first rows so you are “losing” the other ones.
Hi @Thijs, @diagiraldo
I found the average_response script is using the highest b-value of each individual response function which is eventually the first row value of response function and calculating single shell average response function, wheres the value is asked to use for multi shell multi tissue csd calculation in FBA instructions. I am bit confused here which is right way to do it.
Many thanks!
Sarbani
Hi Sarbina,
@diagiraldo is right. I’ve now updated the average_response script to handle multi-shell responses (i.e. multiple lines). Sorry I didn’t get to it earlier as promised. Use git pull && ./build to update your MRtrix install.
Cheers,
Dave
Just for reference: it was actually only using the first line, but that happened to (typically) be the b=0 “shell”. You could clearly see that from your average SF WM response: all the higher order terms were 0.00000, and only the first DC term was non-zero. That’s basically the b=0 “shell” of the SF WM response.
Seeing that you’re working with the HCP data, which definitely is multi-shell, means that you can happily use either msmt_5tt and/or the dhollander algorithm. I’d encourage you to use the dhollander algorithm though (but myself being that “Dhollander”, I may be a bit biased); I’ve seen cases where --even with the HCP data-- it performs better than the msmt_5tt algorithm. Not to mention the fact that it involves less hassle (i.e. no need to bring in the 5tt file).
Thanks @Dave, for updating the script for multi shell data. Its working fine now with dwi2fod msmt_csd.
Cheers!
Sarbani
Yes @Thijs, that’s correct, to me too dhollander algorithm looks hassle free, and thanks for the update that it gives better result with HCP data.
But now I am trying to figure out while generating population_template for fod’s, should I generate separate population_template(s) for CSF-fod images, WM-fod images and GM-fod images and same for the later on steps.
Thanks
Sarbani
Unless you want to analyse the tissue types separately then you need to align them to the same space. I would probably combine the tissue types before registration. Alternatively, you could run population_template for instance only on WM FODs and apply the resulting transformation to the GM FODs to bring them to the same space. This would assume that the WM only registration is also valid for the GM which I am not sure about. A third way would be to align one tissue type (WM), transform the others based on that, combine, and register again.
Cheers,
Max
If you’re only after fixel-based analysis (i.e. of apparent fibre density, fibre cross-section, etc.), then it would currently be sufficient to only use the WM-FOD images to build the population template, and simply ignore the CSF and GM “FOD” images. Note that most value of doing multi-tissue CSD has not been lost because of that: the WM-FOD image itself is still very much cleaned up by removing the contributions of GM and CSF. I suppose this is what you were probably after…? This is what we are currently doing in our lab for all our fixel-based analyses. The first one of those, using my new single-shell 3-tissue CSD (SS3T-CSD), was recently presented at ECTRIMS. Here’s the poster: https://www.researchgate.net/publication/307889845_Apparent_Fiber_Density_A_novel_method_to_detect_axonal_degeneration_in_patients_with_MS . See the right column for images of the template; note the nice depiction of the WM boundary in that template: this is due to the GM and CSF contributions being taken out. This by itself already greatly improves registration and template building outcomes, due to a much better contrast for those algorithms to work with!
But if you’re after looking at the GM and/or CSF itself though, then it comes down to whether you want the analysis of WM, GM and CSF to be in the same template space (to make joint conclusions); if not, separate templates would be sufficient. But if you want a single template space (which probably makes the most sense in this case), then a combined template building strategy of sorts would be what you’re after. There’s ways to modify the script to achieve such things (but not perse unique “trivially best” ways). @Dave and I are actually looking into this, as we understand that people would want to have a combined template (and benefit of the combined WM, GM and CSF information to build such a template).
Dear Max,
I have the outputs from the dhollander method - the WM’s, GM’s and CSF’s FODs. How can I combine the tissue types now?
Thanks
Inês
Hello Inês,
It very much depends on what you want to do. For most analysis of white matter properties, it is reasonable using the WM compartment to co-register the images and using the resulting warps to transform the other tissue compartments if needed. Currently, mrtrix does not offer simultaneous multi-tissue registration but we are working on it.
Cheers,
Max