Dwifslpreproc eddy slice-to-volume correction error

1

Dear All,

I try using dwifslpreproc with eddy-cuda and slice-to-volume correction, like


dwifslpreproc predwi_denoised.mif dwi_preprocessed.mif -rpe_header -eddy_options " --mporder=16"

As input, predwi_denoised.mif is generated using mrcat, with successive concatenations of 4 .mif vols generated from DWI sets via mrconvert.

I have 4 dwi vols sets/subject:
set 1: 10Xb0, 60 dirXb700 (70 vols total), pe dir: j-
set 2: 10Xb0, 60 dirXb700 (70 vols total), pe dir: j
set 3: 10Xb0, 60 dirXb3000 (70 vols total), pe dir: j-
set 4: 10Xb0, 60 dirXb3000 (70 vols total), pe dir: j

Below you have listed the slice timing info for individual sets:

mrinfo set1.mif:
SliceTiming: 4.0250,0.0000,4.1500,0.1250,4.2775,0.2500,4.4025,0.3775,4.5275,0.5025,4.6550,0.6275,4.7800,0.7550,4.9050,0.8800,5.0325,1.0050,5.1575,1.1325,5.2825,1.2575,5.4075,1.3825,5.5350,1.5100,5.6600,1.6350,5.7850,1.7600,5.9125,1.8850,6.0375,2.0125,6.1625,2.1375,6.2900,2.2625,6.4150,2.3900,6.5400,2.5150,6.6675,2.6400,6.7925,2.7675,6.9175,2.8925,7.0450,3.0175,7.1700,3.1450,7.2950,3.2700,7.4225,3.3950,7.5475,3.5225,7.6725,3.6475,7.8000,3.7725,7.9250,3.9000
mrinfo set2.mif:
SliceTiming: 4.0275,0.0000,4.1525,0.1275,4.2775,0.2525,4.4025,0.3775,4.5300,0.5050,4.6550,0.6300,4.7800,0.7550,4.9075,0.8800,5.0325,1.0075,5.1575,1.1325,5.2850,1.2575,5.4100,1.3850,5.5350,1.5100,5.6625,1.6350,5.7875,1.7625,5.9125,1.8875,6.0400,2.0125,6.1650,2.1400,6.2900,2.2650,6.4175,2.3900,6.5425,2.5175,6.6675,2.6425,6.7950,2.7675,6.9200,2.8950,7.0450,3.0200,7.1725,3.1450,7.2975,3.2725,7.4225,3.3975,7.5475,3.5225,7.6750,3.6500,7.8000,3.7750,7.9250,3.9000
mrinfo set3.mif:
SliceTiming: 4.8425,0.0000,4.9950,0.1500,5.1450,0.3025,5.2975,0.4525,5.4475,0.6050,5.6000,0.7550,5.7500,0.9075,5.9025,1.0575,6.0550,1.2100,6.2050,1.3600,6.3575,1.5125,6.5075,1.6650,6.6600,1.8150,6.8100,1.9675,6.9625,2.1175,7.1125,2.2700,7.2650,2.4200,7.4150,2.5725,7.5675,2.7225,7.7200,2.8750,7.8700,3.0275,8.0225,3.1775,8.1725,3.3300,8.3250,3.4800,8.4750,3.6325,8.6275,3.7825,8.7775,3.9350,8.9300,4.0850,9.0825,4.2375,9.2325,4.3875,9.3850,4.5400,9.5350,4.6925
mrinfo set4.mif:
SliceTiming: 4.8425,0.0000,4.9950,0.1500,5.1475,0.3025,5.2975,0.4550,5.4500,0.6050,5.6000,0.7575,5.7525,0.9075,5.9025,1.0600,6.0550,1.2100,6.2050,1.3625,6.3575,1.5125,6.5100,1.6650,6.6600,1.8150,6.8125,1.9675,6.9625,2.1200,7.1150,2.2700,7.2650,2.4225,7.4175,2.5725,7.5675,2.7250,7.7200,2.8750,7.8700,3.0275,8.0225,3.1775,8.1750,3.3300,8.3250,3.4825,8.4775,3.6325,8.6275,3.7850,8.7800,3.9350,8.9300,4.0875,9.0825,4.2375,9.2325,4.3900,9.3850,4.5400,9.5375,4.6925

From the above, set1.mif and set2.mif have almost similar (but not identical) slice timing data), but very different from set3.mif and set4.mif.
Interestingly,

mrcat set1.mif set2.mif set12.mif

creates set12.mif having listed the slice timing info of set2.mif (similarly set34.mif lists set4.mif slice timing if using mrinfo), and the overall concat (mrcat set12.mif set34.mif predwi_denoised.mif) specifies ‘SliceTiming: variable’, hence the error.

Please help, thank you, how can I combine the 4 sets so that slspec data is read properly by dwifslpreproc; if this is not possible, how should I create the slspec file based on such concatenation resulting in dwi_denoised. mif (below I include part of mrinfo output for multiband info etc)
Octavian

mrinfo predwi_denoised.mif -all

Image name: “predwi_denoised”.mif

Dimensions: 128 x 128 x 64 x 280
Voxel size: 2 x 2 x 2 x ?
Data strides: [ -2 -3 4 1 ]
Format: MRtrix
Data type: 32 bit float (little endian)
Intensity scaling: offset = 0, multiplier = 1
Transform: 1 0 -0 -130
0 1 -0 -106.4
-0 -0 1 -44.28
EchoTime: variable
FlipAngle: 90
MultibandAccelerationFactor: 1
PixelBandwidth: 1395
RepetitionTime: variable
SliceEncodingDirection: k
SliceTiming: variable

2

Hi Octavian,

From the above, set1.mif and set2.mif have almost similar (but not identical) slice timing data)
Interestingly, mrcat set1.mif set2.mif set12.mif creates set12.mif having listed the slice timing info of set2.mif

Yes, this comes from my own dealings with similar data; see GitHub Issue. Siemens reports slice times to the nearest 2.5us, but for acquisitions with opposing phase encoding sometimes the rounding goes the other way (:man_shrugging:), resulting in slightly different header content; previously to the change linked above, even concatenating set1.mif and set2.mif resulted in SliceTiming being set to “variable” (which arises from very generic code looking for differences in header content between merged image headers). So I added a bit of tolerance to specifically SliceTiming. The effect doesn’t arise if you adequately take this factor into account in selection of your TR.

If I’m looking at your data correctly, the order of the slice acquisition is identical between the two acquisitions, as is the number of slices; the only difference is that the time in between each successive slice is extended by some factor, presumably because you have minimised your TR independently for the b=700 and b=3000 acquisitions (which is problematic in itself; I’ll get to that later).

In this scenario, the format in which eddy receives the slice timing information actually works out in your benefit. It’s not currently possible to accurately embed this slice timing in the image header, as the timing is different between different volumes. But eddy doesn’t actually operate on slice timing information: it expects as input an “slspec” file, which indicates slice group membership and order. Your image headers contain slice timing, as this is what is extracted from DICOM headers. dwifslpreproc, when appropriate, converts the slice timing information from the header into the slspec format expected by eddy, so that the majority of MRtrix3 users need not learn the intricacies of such things. The reason this is fortuitous for you is that the slice group membership and order is in fact equivalent between these two slice timing vectors, which means that the slspec file derived from either case would actually be identical.

(It does mean that the relationship between the discrete cosine basis coefficients in which eddy represents within-volume movement and absolute time differs slightly; but I digress)

So I see two solutions for you here:

  1. Instead of relying on slice timing in the image header, generate the slspec file in the format that eddy expects, and provide it via the -eddy_slspec command-line option in dwifslprepreoc.

  2. After concatenating, simply re-insert any one of the slice timing vectors back into the image header (doesn’t matter which, since the resulting slspec file will be the same):

    mrcat set1.mif set2.mif set3.mif set4.mif -axis 3 - | mrconvert - -set_property SliceTiming $(mrinfo set1.mif -property SliceTiming) set1234.mif

The fact that your slice timing vectors between b=700 and b=3000 are different by some multiplicative factor means that you’ve used a different TR for each of your b-values. I also note that your echo time is reported as variable. So presumably you have minimised the TE and TR separately for each b-value in order to maximise signal and minimise acquisition time. This limits the types of analyses that can be done with these data. Anything that assumes a parametric relationship between b-value and signal intensity (e.g. the tensor model) cannot be applied to these data, as the b-value is not the only thing that varies between those volumes. In the case of spherical deconvolution, this isn’t a problem, since each b-value shell is treated independently and no relationship is assumed between them. Where this may cause a problem is in your b=0 volumes, which will not contain identical contrast between sets 1&2 and sets 3&4. This is quite likely to mess up estimation of the susceptibility field by topup; I don’t know what effect it may or may not have on eddy. I would suggest throwing away the b=0 volumes for one of those two pairs of sets prior to running dwifslpreproc.

Rob

3

Thank you, Rob, for a comprehensive, clear, and prompt help.
I sent a message to fsl re topup/eddy use on data with multiple TE/TRs, I paste it here:

"topup: I would not mix b0s across different TE/TR. My suggestion would be to use the one good b0 each from data sets 1 and 2. Ideally, provided there is no problem with it, the first b0 in the input to topup should be the first b0 from data set 1.

eddy: There are advantages to running eddy on all data simultaneously, so I would try that first. There is no problem with the different TE/TR for the dwis. The Gaussian process will deal with that. The only potential problem is with the b0s, but given that we are only using a rigid body model for them I think it will be fine. So my suggestion is to just concatenate all four data sets, and maybe look a little extra carefully at the b0 alignment in the results.
Jesper
"

  1. So, I will try to discard the b0s from set34.mif all together and feed the rest to dwifslpreproc, OR, if I understand correctly Jesper’s answer, I would do so only for the topup, and feed the whole data set to eddy, I assume for the latter the following will do:
dwiextract set12.mif b0_12.mif -bzero
dwifslpreproc set1234.mif set1234_proc.mif -rpe_header -se_epi b0_12.mif -eddy_options " --mporder=16"
  1. Would you know if different TE/TRs data would affect other processing steps, specifically dwidenoise, mrdegibbs, and dwibiascorrect?
    Octavian
4

That would indeed work.

Though while Jesper’s instructions give you the ability to preserve b=0 volumes with both TE/TRs, you would want to be sure that you actually want to preserve such. Unless you’re doing custom DWI processing, chances are that whatever you are subsequently feeding those data to are either going to be completely unaware of the TE/TR difference, or not be able to do anything useful with that difference. So the only consequence of keeping all b=0’s is that your “b=0 data” is going to have way more variance between volumes than what would be expected given the noise level, because the software would be interpreting volumes of different contrast as though they have the same contrast.

  • dwidenoise is based on the noise level, so is insensitive to the image contrast; however if you were to concatenate your series prior to dwidenoise, then it would be sensitive to any recalibrations that may have occurred between the two acquisitions.

  • mrdegibbs is a purely per-volume process;

  • dwibiascorrect is based on the mean b=0 image, so would be affected somehow. Exactly how depends on the response of the N4 algorithm (in the case of the recommended ants algorithm) to the averaging of those volumes, which is difficult to predict, but I wouldn’t expect it to introduce a systematic bias in any specific direction.

Rob

5

Thank you, Rob, as always. It becomes clear that I should throw out the b0s for set 3 and 4.
Another way… I wonder if there is a way to run dwidenoise (mrdegibbs), and dwifslpreproc separately on sets_12.mif, sets_34.mif, and then preserve from dwifslpreproc/eddy just the registration of sets_34.mif vols to sets_12.mif, I am not sure registration can be disjointed from topup output though, thank you. Octavian

6

Theoretically, what you would want to do is:

  • Determine the rigid-body transformation from the distortion-corrected Set12 to the distortion-corrected Set34;

  • Compose this transformation with all other transformations applied by eddy, which includes the susceptibility distortion correction using the field as estimated by topup;

  • Apply this set of transformations in a single resampling step.

This would give you DWI data from all sets resampled on a common voxel grid and mutually aligned. But actually doing so is very much a nontrivial exercise, unless somebody somewhere has already created a tool for doing so. The reason for the complexity is that while you could simply register the two corrected sets and transform one to the other, this would result in two transformation & regridding steps rather than one, which will lead to undesirably smooth data. But this wouldn’t be “preserving” the registration between Set12 and Set34 if you were to run pre-processing entirely independently between the two, since that’s precisely the information that yoy wouldn’t have in that case.