New classes for fixel dataset handling

[Lestropie]

This PR is part of what the set of changes aggregated in #2644 have been working towards.
At time of initial posting the changes apply only to tcksift and tcksift2, but that list will increase over time.
Rather than receiving as input an FOD image and performing FOD segmentation at commencement of the command, they will instead take as input a fixel data file encoding fibre density, and automatically load all other information relevant to the fixel dataset.

---

It has taken quite some time since the fixel directory format superseded the initial sparse fixel file format for me to get around to addressing how it conflicts with both pre-existing handling of, and intentions for future developments involving, fixel data.

• In FBA, one needs full user flexibility over the number and names of files stored as well as how to interpret them, whereas in quantitative fixel-wise models there is a known set of parameters and corresponding interpretations.
• If a given command loops over the data for each fixel just once, the overhead of image indexing / invoking get-set functions from function pointers / data being spread across multiple images is fine; but for complex computations, one wants native access & RAM contiguity for all data in each fixel.

Previously, these preferable attributes were achieved using the DWI::Fixel_map<> class and derivatives.
This predated the fixel directory format, the legacy fixel file format, and even the term "fixel".
But it:

• Was populated using the results of FOD segmentation, such that one pipeline could involve segmenting one FOD image multiple times.
• Presented an unnecessary clash in terms of fixel data handling.
• Precluded certain desired future features.

The primary goal of this PR is to supersede the DWI::Fixel_map<> class with new MR::Fixel::Dataset.
This class is more consistent with the philosophy behind the fixel directory format.
However it also provides the relevant quantitative information per fixel in a native, memory-contiguous manner.
Further, it maintains utilisation of the dixel mask information in assigning streamlines to fixels; this was a key advantage of keeping the relevant commands utilising FOD segmentation rather than the reduced resolution of fixel mean direction information only.

---

PR is in a highly preliminary state; intend to branch other developments off of here onto a private repository.
Below is list of items that should largely be addressed before merging to base branch of #2644.

• afdconnectivity: Take fixel FD file as input
• fixelconnectivity: Make use of TrackMapperBase's new FixelMappingPlugin
• tckgen: Have dynamic seeding take fixel FD file as input
• tckmap: Integrate fixel dataset target destination
  • Correspondingly remove tck2fixel entirely
• Finish reworking of SIFT model weights
  • Rename "processing mask" to "model weights" throughout code / interface / documentation
  • Permit user to provide either voxel-wise or fixel-wise image data
  • No longer store voxel-wise image as member variable
  • Add output function(s) (fixel-wise and voxel-wise)
• Rename "total variation" variance, as it is in fact not "total variation" at all
• Change MR::Fixel helper functions to no longer be FORCE_INLINE
• Modify constructor for SIFT::ModelBase::FixelBase class and derivatives to grab ConstRowXprs for fixel direction and dixelmask
• Implement looping constructs for SIFT::Model class and derivatives
  • Just as MR::Fixel::IndexImage provides a looping construct that yields an iterable list of integer indices, this class should be able to do the same for an iterable list of the relevant Fixel class
  • In addition to looping over the set of fixels within a voxel, it should also be possible to loop over all fixels in the model
    (possibly automatically skipping over any fixels for which the model weight is zero)
  • Define separate mutable and immutable Fixel classes, enabling looping over fixels within member functions declared const
    (I spent ages trying to get this working without success...)
• Implement iterator construct for SIFT::TrackIndexRange
• Decide whether to essentially remove all manual utilisations of Image<index_type> as a fixel index image in favour of MR::Fixel::IndexImage
  In many circumstances it will not be desirable to construct a full MR::Fixel::Dataset instance.
  However rather than MR::Fixel::Loop in its current form, it may be preferable to instead modify that class to be initialised exclusively from an MR::Fixel::IndexImage.
• Move the functor responsible for resampling the ACT 5TT image out of the SIFT namespace and into the ACT namespace.
• Regarding the exclusion of fixels:
  • Change default behaviour to exclude untracked fixels rather than retain them, even for SIFT1 (this includes renaming the command-line option)
  • Add paragraph to Description explaining impact of fixel exclusion being done on a per-tractogram basis
  • For SIFT2, better explain (and disambiguate in code) distinction between fixel being excluded from model and being excluded from optimisation
• Speed up implementation mapping streamlines to fixels
  Currently involves boolean check for presence of dixel masks on a per-streamline-voxel-intersection basis; would prefer to reduce the number of unnecessary branch executions

Optionals:

• For afdconnectivity, tcksift and tcksift2: Allow user-specified input to be either an FOD image or a fixel FD file.
  This would require the ability to populate Fixel::Dataset from the contents of DWI::FMLS::FOD_lobes. This is conceptually not too difficult, and indeed the existing population of DWI::Fixel_map<> does something very similar. The catch is that the total number of fixels is not known a priori, and Eigen::Array<> is very expensive to do repeated .conservativeResize() calls (as we found previously with amp2response: amp2response speedup #1406). So what would be necessary is either:
  • A buffer class that would store the incoming fixel-wise data in something more friendly like std::vector<>'s and populate the Eigen::Array<>'s upon competion of segmentation; or:
  • A wrapper class that would do a more intelligent dynamic allocation of Fixel::Dataset members, with a .resize() to fit just the populated data at completion of segmentation.
  • In addition, would likely need to change the mechanism for construction of the IndexImage class. Currently the constructor receives a std::string, meaning the class is constructed from image only. If permitting creation of an image from scratch, would likely want to have static IndexImage::open() and IndexImage::scratch() functions (maybe an IndexImage::create() if wanting to homologate fod2fixel & others).
• Re-introduce capacity for "null lobes"
  This is a capability that exists within the current SIFT(2) implementation, is disabled by default, but I don't quite know what I want to do about it going forward.
  In its pre-existing form, this creates one additional "fixel" per voxel, that "claims" all directions not ascribed to any segmented fixel, and has an FD of zero. In this way, streamlines for which the tangent does not lie within the subset of the sphere corresponding to any fixel are assigned to this "extra" fixel, and always contribute to an excess in reconstructed density there. It may be possible to do something a bit more clever, for instance by comparing the total AFD of all retained fixels to the FOD l=0 term and attributing the difference between them to all directions not attributed to any fixel. But there's still uncertainty around how such data should be used in a downstream optimisation, rather than just assigning all streamlines to the nearest fixel (which is the current default behaviour).
  There's also the question of whether this should be handled at processing time, or whether fod2fixel could include such data in the fixel output directory, with the corresponding dixel mask but with a non-finite direction. All downstream applications would need to be tested for robustness against the latter case.
• Improve handling of non-WM tissue content
  A known issue with SIFT compared to COMMIT is the model fitting in areas of partial volume between GM and WM.
  With COMMIT, a separate forward model component can be provided for GM, and in partial volume voxels the optimiser can choose how much of the DWI signal to attribute to streamlines vs. GM.
  With SIFT, there is a kind of baked-in presumption of correspondence between the proportion of the voxel traversed by streamlines (which itself is likely influenced by ACT), and the proportion of the DWI signal attributed to WM vs. GM by MSMT (and that's not further accounting for which fixels may vs. may not have been successfully tracked).
  This can introduce large model errors at partial volume voxels, which can undesirably drive the weights of individual streamlines very high.
  While I have various tricks such as downregulation of such fixels in the model and regularisation of streamline weights, it remains a weakness.
  What I would prefer is something that better takes into account the presence of non-WM tissue, whether from MSMT or from a 5TT image. While I haven't figured out yet exactly how I would want this to work, I think it should be possible to do something to stabilise the fit.

[Arshiyasan]

Super excited about the changes to tckmap @Lestropie ! Mapping tracts and weights to a fixel dataset would be an amazing feature. Hoping to see it soon!