New command: dwirecon

[Lestropie]

Creating a draft PR to facilitate external tracking.

This work was done in collaboration with @tclose and @arkiev for creation of a DWI pre-processing pipeline in Pydra. We've ported much of dwifslpreproc over to that environment. But one of the final components, being the explicit weighted recombination of DWI volume pairs with identical diffusion sensitisation but opposing phase encoding direction, would have been exceptionally difficult to port to that environment.

Moreover, I am commonly utilising an acquisition strategy that is currently not explicitly handled by dwifslpreproc. I utilise a subset of the capabilities of the scattered-slice acquisition work, where I take the full gradient table, split it into two individually-well-distributed subsets, and acquire one subset with one phase encoding direction and then the other with the opposite phase encoding direction, such that after pre-processing and concatenation the aggregated data is well distributed. However in regions of strong susceptibility distortion, half of the images will be blurred. Ideally, this would be handled either by something like what the SHARD recon does, or with a weighted diffusion model fit. What I've done here instead, as an incremental improvement that's consistent with the current workflow, is to introduce an additional capability into the pre-processing, that will appear in the same place as the explicit volume pair combination. Here, where a volume is expected to be erroneously smooth due to EPI distortion, the reconstructed image intensity will be heuristically weighted with the intensity predicted by other volume groups that have different effects of distortions. This allows slightly better results of DWI model fitting in regions of strong susceptibility distortions. Unlike the SHARD recon this is done on a per-shell basis. This would have been very difficult to manage as Python code within dwifslpreproc.

Therefore here I'm proposing a new C++ command dwirecon, which can perform different operations involving reconstruction of DWI data depending on what is required for any given processing pipeline.

The implemented combine_pairs operation duplicates exactly what is currently done in dwifslpreproc. The combine_predicted operation seems to be doing a sensible job in terms of utilising predicted data in regions of strong susceptibility distortions but leaving data elsewhere unmodified; I'll present some examples here when I get the chance to come back to it.

---

• [ ] Integrate use of dwirecon into dwifslpreproc. dwirecon using the "combine_pairs" option supersedes a reasonably large block of code in dwifslpreproc, which can be erased and replaced with a single command execution.

• [ ] Implement -rpe_split option to dwifslpreproc. This will allow manual user specification that within the input DWI series, the first half of the volumes have the phase encoding direction as specified at the command-line, and the second half have the opposite phase encoding direction; but unlike -rpe_all, where this results in explicit volume recombination, here it would instead trigger the "combine_predicted" operation of dwirecon.

• [ ] Implement -rpe_header detection that the combine_predicted operation should be performed. This will require some form of heuristic that looks at the different volume groups, their gradient tables, and the gradient table of the full concatenated set. I will require this functionality to be in place myself given that:

  • We are intentionally acquiring two phase encoding directions, with the number of volumes in the two sets differing by one, such that the natural user interface of -rpe_split will not apply
  • I want to be able to process these data using either my existing, or upcoming, pre-processing pipeline, which are based on automatically determining the appropriate pre-processing steps to apply given the input data.

• [ ] Implement -volume_pairs option for combine_pairs operation. Given that it is expected that this command will be executed within dwifslpreproc after eddy, it will be taking the motion-corrected gradient table. This introduces ambiguity into the determination of matching pairs of volumes with equivalent diffusion sensitisation but opposing phase encoding directions (assuming you want to do explicit matching rather than assuming equal order). Instead, dwifslpreproc or some other processing pipeline might determine upstream that this is the type of acquisition that has occurred, store that volume index matching in a file, and provide that file as input to dwirecon.

• [ ] Clean up excessive code comments from development process.

• [ ] (optional) Implement "leave_one_out" operation. For each DWI volume, generate an A->SH->A transform using all other volumes in that shell, and write the predictions to the output image. This will do something akin to the "Patch2Self" algorithm, just exclusively on a per-voxel basis rather than a sliding window.

• [ ] (optional) Implement proximity weighting for A->SH transforms When generating a predicted intensity for a given direction, something more similar to eddy's GP predictor could be achieved by doing a weighted least squares fit, where the weights are determined by eg. cosine distance to the direction to be predicted. This means that there would be a separate transform for every direction to be predicted, so would be more computationally expensive; but I'm nevertheless interested to see what kind of impact this might have.

• [ ] Discuss command name. While dwirecon I think makes the most sense given the underlying operations, it conflicts with an existing command in external project https://github.com/dchristiaens/shard-recon. Longer term, I could see the operation of the SHARD dwirecon command being incorporated into this dwirecon command as one operation among many. However in the absence of such a merge, the command name conflict exists. Open to thoughts on how to best resolve this.

[dchristiaens]

The term "reconstruction" is so widely used medical imaging that I suppose a potential naming conflict for dwirecon was only a matter of time. That said, I would have preferred to merge SHARD-recon into the main repo before other commands with this name are added.

Moreover, the overarching purpose you describe could, I think, be achieved directly with dwirecon (SHARD edition). You would simply need to provide the field map and the PE table, and set the appropriate slice thickness for the split acquisition you describe. It is perfectly possible to run this single-shell, if you want. From what I understand, the only part that is not supported in dwirecon (SHARD edition) is the leave-one-out approach you describe, but this is not yet implemented in this PR either. Did you compare both approaches?

Longer term, I could see the operation of the SHARD dwirecon command being incorporated into this dwirecon command as one operation among many.

May I suggest to do it the other way around? The additional functionality to deal with subject motion, the slice profile, outlier weights, etc., bring about a lot more complexity.

[Lestropie]

dwi2dwi? :upside_down_face:

the overarching purpose you describe could, I think, be achieved directly with dwirecon (SHARD edition).

Hmmm, potentially. The SHARD recon behaviour in the -rpe_all context would maybe be closer in behaviour to eddy's least-squares reconstruction, rather than the very simple weighted averaging currently in place. For the split PE design, it would certainly be more principled than what I'm doing here. But how the outcomes may differ, I don't know, I haven't tried it. For the Pydra pipeline we'd like to integrate SHARD pre-processing as an alternative workflow option in the future (was going to email you about that separately once we had a minimal working product). Here I'm currently only trying to facilitate utilisation of the current pipeline in that environment, and am conscious of potential road blocks given external time pressures.

What would be the prospect of integration of a subset of the SHARD capabilities into 3.1: enough to facilitate the functionalities required here, but not the entire pre-processing pipeline?

One of the tasks on my to-do list is to acquire single-session data using all of the different phase encoding strategies that can be handled by dwifslpreproc, to use as exemplar workshop data. Perhaps in the process of generating that data I can finally try getting my hands dirty with SHARD, not only as a potential substitute for dwifslpreproc but potentially also just for estimation of the final DWI data post-eddy.

May I suggest to do it the other way around?

My thinking was that if these were separate "operations" within the same command, akin to something like mrfilter, then the order of implementation / distribution wouldn't matter; there may just be some operations that have many command-line options only applicable to that operation. Whereas if I'm following your train of thought, there would be no need for such distinct "operations"; they would all fall under the banner of "given these input DWI data, and other input parameters, give me output DWI data".