ENH dwifslpreproc: Volume combination options

[Lestropie]

Combination of long-standing desire of mine and a recent feature request.

Once the geometric corrections have been performed, there are currently two code branches:

1. Do nothing; all volumes generated by eddy will form a part of the output series.
2. If -rpe_all is used (or -rpe_header is used and the input data conform to the corresponding requirements), perform explicit combination of volume pairs that have identical diffusion sensitisation but opposing phase encoding, as per Skare and Bammer 2010.

There are however multiple other prospects that are not currently possible:

1. Tentatively called "-rpe_split": A scheme where DWIs have been split between two phase encoding directions. The DWIs acquired within each phase encoding direction individually are each reasonably well distributed, but they are also designed such that when they are concatenated together the coverage becomes more dense, ie. each fills the gaps of the other. Where data are acquired in this way, I would like for there to be an additional heuristic step applied to the output reconstructed volumes. Where the data have been severely compressed by EPI distortions, and therefore there will be spatial smoothing in the reconstructed image, I would like for the output image to be a linear combination of the empirical data and the predicted intensity based on other volumes that are expanded rather than compressed by those same distortions.
2. If data have been acquired in a way intended for -rpe_all, but there is substantial motion in the data such that explicit combination of volumes with matched diffusion gradients but opposing phase encoding directions is no longer appropriate (as the direction of diffusion sensitisation relative to anatomy has changed), it should be possible to use the -rpe_all flag to allow for inference of the phase encoding table, but to disable the explicit recombination of volumes using another command-line flag.
3. It should be possible to activate eddy's Least-Squares Reconstruction (LSR). This had been disabled in a very early version of eddy, hence I had implemented something different; I believe it's now available, though it can produce some strange artifacts sometimes. (*Edit: This is already listed as #1095, but makes sense as part of this overall request)

[Lestropie]

On reflection there's an additional complexity here:

When we do the "-rpe_split" approach, we don't acquire exactly the same number of volumes with each of the two phase encoding directions. We deliberately get one more volume with one phase encoding direction than the other. This protects against the error of loading the same set of diffusion gradient directions for both phase encoding directions, as the # volumes & acquisition time will be different for the two series. But it means a simple boolean command-line flag is not sufficient: where the input series has an odd number of volumes, the script would need to know whether the first phase encoding direction has one more volume than the second, or the other way around.

[Lestropie]

RE "When we do the "-rpe_split" approach" above: I think that I would implement this specific command-line functionality in a similar way to how -rpe_all and -rpe_pair work, in that they necessitate assuming that the first half of all volumes have one phase encoding direction, and the second half have the opposite. In our use case, where there is not a perfectly equivalent number of volumes per phase encoding direction, this would necessitate the use of -rpe_header; the allocation of volumes to phase encoding groups would then be known from the metadata.

[Lestropie]

The -rpe_split approach could potentially become quite complex. Particularly in the -rpe_header approach where you are not necessarily just utilising data from two phase encoding groups with opposing polarity: you would have some number of phase encoding groups, and for every DWI intensity I would want to generate a new intensity based on not only that sample, but also the intensities predicted from all other phase encoding groups, and the Jacobians of all phase encoding groups. This would also become very clumsy if trying to manage within a Python script, potentially even more so in a workflow management engine environment (pinging @tclose).

So I am considering whether this should instead be a C++ binary.

In other MRtrix-related projects there's been a need for a C++ command that "reconstructs" diffusion data given a set of input DWI data. I think that what I'm requiring in the context of this Issue falls into the same category, and so it might make sense as a single command that has different operating modes or flags. So for instance:

• Leave-one-out: Dictates whether the output intensity of a given DWI volume should be done in a way that ignores the empirical intensity of that volume. (So in one use case would do something akin to dipy's Patch2Self; but could also potentially be used in pre-processing scripts)

• Reconstruction with varying phase encoding: For each phase encoding group, generate a predicted DWI signal; for each DWI volume, derive some heuristic for reconstructing the output intensity based on the empirical intensity as well as the estimated intensities from other phase encoding groups. This is something that I'm yet to derive, but have thought about a bit. I want to do something simple like the Skare 2010 ISMRM abstract that is used for volume recombination, ie. a linear combination of signals based on Jacobians, but the nature of the expression should be different here. Basically, if a sample has been spread out in space by susceptibility distortions, it should be Jacobian modulated but otherwise left as-is; if it is instead compressed, and therefore there is a risk of loss of spatial contrast, it should progressively get down-weighted in favour of predicted intensities from other phase encoding groups that have not experienced such compression. Yet to propose an analytic expression though.

• Volume recombination: Implements what currently happens within dwifslpreproc, where volume pairs are identified with equivalent diffusion sensitisation encoding but opposed phase encoding, computes per voxel the relative weights based on the Jacobian of the susceptibility distortion field, and optimally combines them into a single volume.

• Imputing basis: There should be a natural augmentation between using shell-wise SH fit and a SHARD fit.

• Imputing weighting: Typically in this context people have done a simple amp2SH -> SH2amp with lower lmax to generate predicted intensities. I'm also interested in evaluating whether using a weighted SH fit, with greater weights ascribed to more proximal samples, could be of benefit here.

I've a suspicion that using @dchristiaens's dwirecon as the starting point, and working backwards to implement more simplistic features, might be on the harder side. @jdtournier Did you have source code somewhere for simple DWI prediction that would be useful as a starting point?