ENH dwifslpreproc: Support phase encoding in IS direction

[Lestropie]

topup fails if the component of the phase encoding direction on the third spatial axis is non-zero. Ideally dwifslpreproc would internally permute spatial axes to shift the phase encoding direction onto one of the first two spatial axes, run topup / eddy, then revert the permutation for the output data. However all internal data generated by dwifslpreproc for FSL would need to reflect this appropriately.

[celstark]

Is there any reasonable work-around? The issue, as I understand it, is that MRTrix always has IS as the 3rd dimension and topup is assuming the axes are more-or-less as acquired. So, there would be no reason why the through-plane axis would be the phase encode (it must be one of the first two -- the in-plane ones).

Could this be handled right line ~800 in dwifslpreproc as you make the topup_in.nii, re-orienting the nifti file there and then re-orienting the output before applying it?

[Lestropie]

There's a combination of a few issues at play. MRtrix is happy for IS to be something other than the third axis in an input image; it internally rejigs its interpretation of the image at load time so that the third axis is whichever is closest to IS. But whenever within a script we export to NIfTI, including in preparation for topup, we explicitly specify the order of the strides such that eg. ijk = RAS or LAS in the output image.

So it's conceivable that in some circumstances an image that can be used with topup directly won't work with dwifslpreproc because it's the manipulations performed by the latter that results in the phase encoding direction being in the third axis. But I'm not 100% sure that this is what's happening from your description.

Modifying code to support this case without breaking any other use case is very difficult. There's previously been issues with the interaction between that internal axis rejigging and the external representation of phase encoding information (see #2308 if you want to make your eyes bleed), and perhaps full resolution of that since this issue was listed means that such explicit axis reordering is no longer necessary. I'm nevertheless paranoid that any attempt to support this use case will need extensive testing, and not just of this use case but every other as well. Also it would almost certainly not be sufficient to just modify the input to topup, because the field coefficients are calculated based on the orientation of that image, and that's then fed to eddy, which makes a stringent assumption of spatial correspondence between that pre-calculated field and its own input data. So everything downstream of that would need to be modified as well, which means that the outstanding issue with interpreting external diffusion gradient information in the presence of MRtrix internal transform realignment (#2416) would then rear its own ugly head.

Indeed even the statement that "the first two axes are in-plane and the third is through-plane" is not strictly true. It may be how the data come out of DICOM for a 2D sequence, but anything after that can be permuted to your heart's content. While this will be the same logic by which topup refuses to accept phase encoding along the third image axis, IMO there isn't actually any need for such an assumption, and so I was silently hoping that the FSL devs would simply remove that check, in which case this issue listing would resolve itself upstream.

If you have data that currently won't process because of this issue, the "reasonable" workaround is to manually do the permutation yourself prior to running the script (making sure that interpretation of any orientation-specific information like diffusion gradients and phase encoding are correct) and then reverting that change after script completion if you choose. At least that way if things go wrong it won't be my fault for once 😬

[celstark]

Thanks for the detailed note. On the through-plane, I was thinking that was FSL's assumption and therefore why they have that check. Looking at their code, it does seem this is quite hard-coded:

if (_pevec.Nrows() != 3) throw TopupException("TopupScan::TopupScan: pevec must have three elements");
  if (_pevec(3) != 0.0) throw TopupException("TopupScan::TopupScan: third element of pevec must be zero");
  _orig = boost::shared_ptr<NEWIMAGE::volume<float> >(new NEWIMAGE::volume<float>(scan));
  _orig->setextrapolationmethod(NEWIMAGE::periodic);
  if (_pevec(1) && !_pevec(2)) _orig->setextrapolationvalidity(true,false,false);
  else if (!_pevec(1) && _pevec(2)) _orig->setextrapolationvalidity(false,true,false);
  else _orig->setextrapolationvalidity(false,false,false);

They're pretty explicitly treating the first two dimensions one way and the third dimension differently. So, phase encoding must be in one of the first two. I'm a bit curious to see what happens if we take out their check and setextrapolationvalidity(false,false,true) there, but, as you note, there are a lot of potential downstream issues.

[celstark]

If you have a moment, can you peek over my shoulder here to check my thinking? I think I can just lie and leverage off of the fact that bvecs are image-axis based to make for a reasonably easy fix. By changing around the NIFTI header we change how the image will be interpreted, but the bvecs still map on just fine. Process it that way, swap it back, and then re-mif-ify the output. Is this basically what you were thinking?

Exploit limitation of bval/bvec

The bvec files are WRT image x,y,z, not real-world x,y,z. We can use this to our advantage here. In truth, our image is RIP (AFNI RIP - seems to be LSA in MRtrix [ -1 3 2 4 ]), but dwifslpreproc will dump it out as LAS (they use an explicit -strides -1,+2,+3,+4 when making the topup_in.nii - note AFNI calls this RPI). So, even when we create our NIFTI / bval / bvec version and run on that:

mrconvert -export_grad_fsl step0.bvec step0.bval step0.mif step0.nii.gz

We get something that MRtrix calls [ -1 3 2 4 ] (MRtrix LSA, AFNI RIP). But, MRtrix converts it back and again dumps it out as LAS. But, if we lie and recode the NIFTI file header first to something like LAS (MRtrix aka RPI in AFNI), we should be able to keep the bvecs (the image hasn't changed, just the header) it should work. The only trick is that we need to do this for the b0 images and we need to lie about which direction the phase encoding was on. It's not on IS anymore, but rather on "PA".

cp step0.nii.gz step1.nii.gz
3drefit -orient RPI step1.nii.gz

mrconvert  ${der_path}/../b0/sub-${sid}/sub-${sid}_dir-both_b0.mif step0_b0.nii.gz
cp step0_b0.nii.gz step1_b0.nii.gz
3drefit -orient RPI step1_b0.nii.gz

dwifslpreproc -config BZeroThreshold 20 \
    step1.nii.gz step1_edd.nii.gz \
    -fslgrad step0.bvec step0.bval \
    -rpe_pair -se_epi step1_b0.nii.gz \
    -pe_dir pa -nocleanup -scratch ${der_path}/preprocessed/sub-${sid}

[Lestropie]

Looking at their code, it does seem this is quite hard-coded:

OK, a couple of observations:

• I understand the logic of having through-plane data treated differently to within-plane. I've encountered that issue in the context of convolutional filters at image edges (#2278). My guess is that "extrapolation" here is referring to performing wrap-around when attempting to access the data beyond the FoV. What I don't understand is the necessity to constrain this capability to only the phase encoding axis; it's equally applicable to the readout axis, it's just unlikely to actually have an effect.
• The final line suggests that it's permissible to have phase encoding directions of [0, 0, 0] or [1, 1, 0], since otherwise that final line would never be executed.

By changing around the NIFTI header we change how the image will be interpreted, but the bvecs still map on just fine. Process it that way, swap it back, and then re-mif-ify the output. Is this basically what you were thinking?

Personally I think that any manual "tricks" performed here are more likely to be wrong than right. I prefer to lean on the API code that's been carefully tested. I don't have any experience with 3drefit so can't advise on its use, but any manipulation of NIfTI axes that does not concurrently modify bvecs makes me very nervous. The advantage of MRtrix3 internally having the representation of the gradient table as internal metadata of the image (whether actually stored within the image header on the filesystem, or explicitly loaded using eg. -fslgrad) is that whenever any modification is made to the image, any orientation-dependent metadata internally stored alongside the image is immediately modified accordingly.

I'd be doing: (assuming that the axes of the b=0 data are the same as the DWI series)

mrconvert step0.mif -axes 0,2,1,3 step1.mif
mrconvert ${der_path}/../b0/sub-${sid}/sub-${sid}_dir-both_b0.mif -axes 0,2,1,3 step1_bzero.mif
dwifslpreproc -config BZeroThreshold 20 \
    step1.mif step2.mif 
    -rpe_pair -se_epi step1_bzero.mif \
    -pe_dir pa -nocleanup -scratch ${der_path}/preprocessed/sub-${sid}
mrconvert step2.mif -axes 0,2,1,3 step3.mif

You would need to double-check whether the phase encoding direction is still appropriate once everything is imported into scratch. "j-" would be more faithful here since we're dealing more explicitly with image axes, even though the result should be the same. The tricky bit with that field is that they correspond to orientations / axes as they are interpreted by MRtrix3*, ie. following internal transform realignment, not* as they are stored on disk. Which can cause all sorts of problems.

[celstark]

LOL... a developer preferring to use the code he knows inside and out as he wrote it? Go figure ;) A user who has used a package for a quarter century preferring to lean on it rather than one he's just getting to know? Go figure ;) I've probably 3drefit thousands of datasets by now, so yea, it's one I know very well and I'm still just getting my head around strides vs axes and such.

Your version certainly is cleaner and isn't leaning on the fact that bvec is so image-axis based. The good news is that I ran both and they overlay perfectly. I've got a thread up on the FSL email list (their first thought was to ask why the heck phase would be encoded in the through-slice direction) and maybe some recoding there will come out of it. If not, this thread should help anyone else with the issue.

Thanks!

[Lestropie]

:grinning:

Couple of tangents:

1. In BIDS land, the axis corresponding to slice encoding direction gets explicitly specified in metadata, just as is the phase encoding direction. So in that scenario you know exactly what's what. Though the situation gets more complex for MRtrix3 then because that sidecar information could be either embedded in the image header or in a sidecar JSON file, and in either case if MRtrix3 does an internal transform rejigging then the same transformation needs to be applied to those metadata; hence the difficulty of #2308.

2. For eddy's slice-to-volume motion correction, there is an implicit assumption that the third image axis corresponds to slices. So the scope of this issue should technically be expanded a little. It's not just "permit phase encoding along the third axis", but more like "ensure that the third axis is slice encoding", but where the latter information is absent it must revert to at least "make sure phase encoding is not the third axis, since we know the same axis can't be used for both phase and slice encoding".

3. The issue of on-disk vs. interpreted strides / transform / metadata is listed as #2526. I think I need to prioritise that at some point as it will be a far better communication of the presence & consequences of the internal header transform realignment that MRtrix3 does by default.