I’m trying to replicate the method of Mirzaalian et al., 2016, in which spherical harmonic deconvolution is used to harmonize DMRI data from different scanners (with similar acquisition schemes). Basically they calculated scaling factors on the SH coefficient images using RISH values from their test subjects at each site. I can see that the dwi2fod command in MRtrix gets me part of the way there to implementing this – I now have coefficient images. But is there a command to go in the opposite direction? That is, to recreate the DWI data from the coefficients after I modify them? The idea would be that, using this method, any DMRI processing software should be able to handle the harmonized images.
I have a feeling you’re after the amp2sh and sh2amp commands for this. dwi2fod does a fair bit more than a straight SH fit, which I think is what you’re after here?
I have two questions regarding amp2sh and sh2amp to project back and forth between DWI and SH:
When I apply amp2sh and try to reconstruct then DWI from the resulting SH, it seems to me I do not get the original DWI volumes. Is this discrepancy expected, or am I doing something wrong? I do not normalise to b0 volumes, and I use the max number of l (i.e. for that data 6)
When applying sh2amp, I do either get no b0 volumes at all (that is the case when I feed a gradient file with [x,y,z,b] to sh2amp), or the b0 volumes are completely empty (when I feed a file with [x,y,z], I assume a zero vector just means no signal). How can I recover b0 volumes from SH, or is that not possible at all?
The amplitude to SH transform generally entails rank reduction: the number of SH coefficients will typically be fewer than the number of DW signals you’re fitting to – unless the number of DW volumes matches exactly (i.e. 28 volumes for lmax=6). This means that there’s always some degree of ‘smoothing’ or information loss involved, and you won’t in general be able to reproduce the input signals exactly (unless, as I mentioned above, the number of SH coefficients matches the number of DW signals exactly).
The SH fit is done on a per-shell basis. So by default, the fit is done to the highest shell available, all other shells are ignored. To get around this would require extending the SH basis to also apply along the radial / b-value domain – which is what we tried to do with this recent paper. Watch this space…
Okay, that is what I expected, but I am happy you kindly confirmed.
You are referring to DWI -> SH on multi shell data if I am not mistaken, but I was wondering if I had single shell data, applying amp2sh I get my SH coeff. Then I want to retrieve my DWI signal from SH again (sh2amp), but then the b0 volumes cannot be reconstructed from SH coeff? Or am I misunderstanding…
Regarding point 2: the fit is done per shell, and a shell means: the set of volumes acquired with the same b-value (and different directions). So by that definition, the b=0 volumes constitute a separate shell – albeit a very small one – so they won’t be part of the outer shell’s SH fit, and can’t be regenerated from its coefficients.
If you think about it, that also means that a single-shell HARDI acquisition is already actually multi-shell… Which is why we can do 2-tissue multi-shell multi-tissue (MSMT) CSD with regular single-shell acquisitions.
How can I recover b0 volumes from SH, or is that not possible at all?
To answer this question in the most literal sense, independently of the capabilities and interfaces of the relevant MRtrix3 commands:
Extract the b=0 volumes from your DWI dataset.
Compute the mean image intensity across volumes.
Congratulations! You just did the equivalent of performing a lmax=0 least-squares SH fit to your b=0 “shell” (which is the maximum spherical harmonic degree achievable for such data), and mapped the result back to the predicted “DWI” (b=0) image intensities.
A further question regarding sh2amp. When i use am2sh to construct SH and then reverse the operation with sh2amp, it seem that the diffusion sensitisation is flipped for the x-axis. I use a gradient file for the directions ([x,y,z]), when I adjust the directions to [-x,y,z] it seems I get nearly the same signal as in the original DWI volume. Am I doing something wrong, or is this maybe an inconsistency with FSL and MRtrix?
When i use amp2sh to construct SH and then reverse the operation with sh2amp, it seem that the diffusion sensitisation is flipped for the x-axis. I use a gradient file for the directions ([x,y,z]), when I adjust the directions to [-x,y,z] it seems I get nearly the same signal as in the original DWI volume. Am I doing something wrong, or is this maybe an inconsistency with FSL and MRtrix?
To know for sure, we would need to know exactly:
Where the diffusion gradients are obtained for the initial amp2sh call;
How you derive the external gradient file that you use for the sh2amp call;
Which version of MRtrix3 you are using (particularly whether it is newer than this fix).
How you derive the external gradient file that you use for the sh2amp call
This is derived simply by taking the b = 1000 b-vectors from the bvecs file and transposing them so that the resulting b1000_bvecs file will have 3 columns and 90 rows.
then using the command
sh2amp sh_coeffs.nii.gz b1000_bvecs out.nii.gz
Exactly the same as @stefan, the x axis attenuation seems flipped, and simply by modifying the b1000_bvecs from [x, y, z] to [-x, y, z], the resultant image from sh2amp is almost exactly perfect.
OK, the issue here will relate to the conventions assumed in how the DW directions are stored in MRtrix and FSL formats. The sh2amp call expects the directions to be provided in MRtrix format, which does not make the same assumptions as the FSL bvecs/bvals format. You’ll find all the details in our docs, but the main point is that directions in MRtrix format are assumed to be provided with respect to the scanner/word coordinates, while FSL assumes them to be provided with respect to the image frame – with an additional (and unexpected!) potential flip of the x-axis depending on the handedness of the image coordinate system (i.e. whether the image is stored in ‘neurological’ or ‘radiological’ convention, or more precisely, the sign of the determinant of the transform matrix). That last gotcha has caught us out in the past…
If you’re using MRtrix 3.0.2, then you should be able to provide the original image as the source for the directions, along with the shell you’re interested in, i.e. something like this:
Thanks @jdtournier for the quick and in-depth response! This is certainly not the first time that I’ve been caught out by the x-axis flip
The command you suggested did work as expected. Although sh2amp did complain that -shells was an unknown option; not too much of a problem though as you can easily export a single shell using mrconvert.
Which is why we can do 2-tissue multi-shell multi-tissue (MSMT) CSD with regular single-shell acquisitions. 
