Subset Handling Classes

[robbietuk]

I have been looking into the use of subsets in STIR, with a focus on sinogram reconstruction.

Background

STIR subsets only relate to the "ordered" subsets methodology. The two key variables are num_subsets and subset_num, which are defined at IterativeReconstruction and in GeneralisedObjectiveFunction classes. The number of subsets is used in the projectors (symmetries etc.).

At each iteration of an iterative algorithm, subset_num = this->get_subset_num(). subset_num is parsed through various methods to distributable_computation where both variables are used in find_basic_vs_nums_in_subset (also takes symmetries). find_basic_vs_nums_in_subset returns a vector of ViewSegmentNumbers, a class composed of 3 variables (view, segment, and tof), which corresponds to all view, segment (and tof) values associated with the particular subset specified.

Note: find_basic_vs_nums_in_subset also requires proj_data_info and symmetries to be parsed to iterate from min to max values in each view, segment, and tof dimension

Proposed

It may be of interest to collect these methods into a common class. I therefore propose a new class hierachy (temporarily denoted SubsetNumbersClass) that would contain all information related to subsets, including find_basic_vs_nums_in_subset and num_subsets.

Secondly, I propose this class be parsed seperately to OSMAPOSL and OSSPSwith its own add_parsing_key and that num_subsets and subset_num be renamed num_view_subsets and subset_view_num respectively. This would allow for furture expansion into subset in additional dimentions, e.g. TOF and segments.

Furthermore, additional classes can be introduced to handle get_subset_num() to allow for an expansion into additional subset sampling methods (e.g. stochastic).

Points of discussion

• Class hierarchy and abstraction of classes to allow for projection and list mode data from a common abstract class
• Impact on STIR-python/matlab and general backwards compatability.

[AnderBiguri]

Extra info: the (still in progress) GPU-projector code in PR #817 can do arbitrary subsets (the GPU code itself, not the PR). Currently, the subset are composed by projections that share some symmetry (I am I am not wrong) so they can not be fully arbitrary, otherwise the CPU code will be much slower (because not exploiting symmetries means much much slower code). However, GPU code does not have this limitation, you can compute any arbitrary projection/bin on the same subset without any loss of generality or speed. Allowing to do have arbitrary subsets can be of benefit here.

So, adding this idea would have the benefit of also unlocking extra capabilities of the soon to come GPU code.

[robbietuk]

This is correct and the GPU projector will infuture be able to handle subsets in the same way STIR can handle them. Currently it can only handle 1 subset.

[KrisThielemans]

Note that for SIRF/CIL, with @paskino we were discussing another way of handling subsets: you encode the subset info into the acquisition data (e.g. if on disk and read-only, by sharing the file, if in memory, presumably not). You can then write the objective function as a sum of objective functions (in STIR using SumOfGeneralisedObjectiveFunctions, and generic subset algos for any sum. This is quite interesting, but also somewhat different. Not so easy to write OSEM et al. that way (as it has to check balancing etc, and it's very unclear how to do that generically).

Of course, both ways (subset of data, or subset of objective function) need a generic Subset class.

[robbietuk]

It is definitely an interesting way to approach the problem, and probably moreso correct based upon optimisation literature than PET EM like algorithms, although the concept still holds. The biggest trouble is that this would probably require A LOT more work as it probably involves rewriting a lot of STIR algorithms, which probably needs doing to if we were to add more.

We definitely need a discussion with @paskino and CIL people on this.