Moving SyneRBI-Challenge nema-data utilities to SIRF

[evgueni-ovtchinnikov]

Changes in this pull request

Testing performed

Related issues

Checklist before requesting a review

• I have performed a self-review of my code
• I have added docstrings/doxygen in line with the guidance in the developer guide
• I have implemented unit tests that cover any new or modified functionality
• The code builds and runs on my machine
• CHANGES.md has been updated with any functionality change

Contribution Notes

Please read and adhere to the contribution guidelines.

Please tick the following:

• The content of this Pull Request (the Contribution) is intentionally submitted for inclusion in SIRF (the Work) under the terms and conditions of the Apache-2.0 License.

[KrisThielemans]

I'm not 100% sure if this should be in C++ TBH. I think most people will never check the C++ code, while they might do for the Python side. However, of course we can have it in C++ and an equivalent in Python as illustration.

Regarding the "incompatible projdata", if you post an error, I can try to help.

[evgueni-ovtchinnikov]

@KrisThielemans
The following line in cSTIR/tests/test7.cpp

		asm_.normalise(*background_sptr);

gives me

WARNING: BinNormalisationFromProjData: incompatible projection data:

Norm projdata info:
ProjDataInfoCylindricalNoArcCorr := 
Scanner parameters:=
  Scanner type := Siemens mMR
  Number of rings                          := 64
  Number of detectors per ring             := 504
  Inner ring diameter (cm)                 := 65.6
  Average depth of interaction (cm)        := 0.7
  Distance between rings (cm)              := 0.40625
  Default bin size (cm)                    := 0.208626
  View offset (degrees)                    := 0
  Maximum number of non-arc-corrected bins := 344
  Default number of arc-corrected bins     := 344
  Energy resolution         := 0.145
  Reference energy (in keV) := 511
  Number of blocks per bucket in transaxial direction         := 1
  Number of blocks per bucket in axial direction              := 2
  Number of crystals per block in axial direction             := 8
  Number of crystals per block in transaxial direction        := 9
  Number of detector layers                                   := 1
  Number of crystals per singles unit in axial direction      := 16
  Number of crystals per singles unit in transaxial direction := 9
  Scanner geometry (BlocksOnCylindrical/Cylindrical/Generic)  := Cylindrical
End scanner parameters:=

start vertical bed position (mm) := 0
start horizontal bed position (mm) := 0

TOF mashing factor in data:      0
Number of TOF positions in data: 1

Segment_num range:           (-1, 1)
Number of Views:                126
Number of axial positions per seg: {115 127 115 }
Number of tangential positions: 344
Azimuthal angle increment (deg):   1.42857
Azimuthal angle extent (deg):      180
ring differences per segment: 
(-16,-6)(-5,5)(6,16)
End :=

Emission projdata info (made non-TOF if norm is non-TOF):

ProjDataInfoCylindricalNoArcCorr := 
Scanner parameters:=
  Scanner type := Siemens mMR
  Number of rings                          := 64
  Number of detectors per ring             := 504
  Inner ring diameter (cm)                 := 65.6
  Average depth of interaction (cm)        := 0.7
  Distance between rings (cm)              := 0.40625
  Default bin size (cm)                    := 0.208626
  View offset (degrees)                    := 0
  Maximum number of non-arc-corrected bins := 344
  Default number of arc-corrected bins     := 344
  Energy resolution         := 0.145
  Reference energy (in keV) := 511
  Number of blocks per bucket in transaxial direction         := 1
  Number of blocks per bucket in axial direction              := 2
  Number of crystals per block in axial direction             := 8
  Number of crystals per block in transaxial direction        := 9
  Number of detector layers                                   := 1
  Number of crystals per singles unit in axial direction      := 16
  Number of crystals per singles unit in transaxial direction := 9
  Scanner geometry (BlocksOnCylindrical/Cylindrical/Generic)  := Cylindrical
End scanner parameters:=

start vertical bed position (mm) := 0
start horizontal bed position (mm) := -19

TOF mashing factor in data:      0
Number of TOF positions in data: 1

Segment_num range:           (-1, 1)
Number of Views:                126
Number of axial positions per seg: {115 127 115 }
Number of tangential positions: 344
Azimuthal angle increment (deg):   1.42857
Azimuthal angle extent (deg):      180
ring differences per segment: 
(-16,-6)(-5,5)(6,16)
End :=

--- (end of incompatible projection data info)---

(note the difference in start horizontal bed position (mm)), and the rest of the code is not executed.

I tried different declarations of background_sptr to no avail.

[KrisThielemans]

Cutting out some stuff

std::string f_template = append_path(path, "mMR_template_span11_small.hs");
CREATE_OBJECT(STIRAcquisitionData, STIRAcquisitionDataInFile,
		acq_data_template, templ_sptr, f_template.c_str());
...
converter.sinograms_and_randoms_from_listmode
		(lm_data, 0, 10, acq_data_template, sinograms_sptr, randoms_sptr);
...
PETAttenuationModel::compute_ac_factors(templ_sptr, mu_map_sptr, acf_sptr, iacf_sptr);

The last line seems incorrect. the template should no longer be used, as that doesn't know about bed positions, time frames etc actually used. Instead, use sinograms_sptr

[evgueni-ovtchinnikov]

@KrisThielemans thanks a lot - seems to work ok now!

[KrisThielemans]

not sure about the direction this is going I'm afraid. Some more detail in the comments.

[KrisThielemans]

let's rename sino to prompts. I'd prefer to rename rand to randoms_estimate (or at least, say that in a comment). They are not the randoms of course.

Also rename the function to get_prompts_and_randoms(...). However, I'm not a great fan of this I'm afraid. We have set_input etc, and here a (presumably static) function that takes an lm_data and time-frame info. I don't see a good reason to have both interfaces. So, I suppose I suggest to have no arguments for this function (i.e. bunch of sets, a set_up(), and then call this). It raises the question why we have it, aside from having a one-liner.

We're mixing 2 kinds of interfaces here, which creates confusion.

If you don't like the set functions. we could have a constructor function for Listmode2Sinograms

def __init__(self, list_file, interval, template)

Same comments will apply on the C++ side of course.

[evgueni-ovtchinnikov]

This method is just a one-line alternative to already available usage

lm2sino = pet.ListmodeToSinograms()
lm2sino.set_input(lm_data)
lm2sino.set_output_prefix("prompts")
...

But this PR is about adding Challenge nema-data utilities to SIRF, and all those utilities are one-liners, as can be seen in try_data_preparation.py:

sinograms_and_randoms_from_listmode(listmode_filename, (0, 600), template_filename, sinograms_filename, randoms_filename)

ac_factors = acquisition_sensitivity_from_attenuation(siemens_attn_header, sinograms_filename, stir_attn_header, ac_sens_filename)

estimate_scatter(sinograms_filename, randoms_filename, stir_norm_header, stir_attn_header, ac_factors, scatter_filename)

mult_factors = create_multfactors(stir_norm_header, ac_factors, mult_factors_filename)

create_background_and_additive(randoms_filename, scatter_filename,  mult_factors, background_filename, additive_term_filename)

So I tried to make respective SIRF one-liners, just replacing files with objects in memory.

[KrisThielemans]

I am not a fan of poluting class interfaces with multiple options to do the same thing (although we have a fair amount of that in SIRF already), but I suppose others might prefer oneliners. Nevertheless, I certainly don't want to have all the above one-liners in class interfaces (it wouldn't make sense anyway).

I still think this particular function has a few problem:

• it is PET specific (randoms don't exist in SPECT). This was already slightly problematic in ListModeToSinograms but at least people would have to explicitly call estimate_randoms, which nobody in SPECT would do.
• for dynamic PET, people might need to call get_time_at_which_num_prompts_exceeds_threshold first. To do that, they need to call set_input first. It's contour-intuititive (and suboptimal) that this function then needs the listmode data again, so you'd need to remove the lm_data argument.
• it calls set_time_interval etc, so will modify the underlying ListModeToSinograms object, which could be surprising, at least without documentation.
• (as stated above, the name is wrong)

See general comment in the review.

[KrisThielemans]

This PR introduces convenience functions for PET data preparation. In particular, it adds sirf.STIR.AcquisitionSensitivityModel.compute_attenuation_factors and sirf.STIR.ListmodeToSinograms.prompts_and_randoms_from_listmode (and corresponding C/C++ functions).

For the challenge, we could have just added them in the Challenge repo, but it's obviously better to have them in SIRF itself.

However, I think the current PR has several problems.

• it adds "one-liner" to class interfaces, introducing duplication in functionality
• most of it is PET specific but they are at "top-level" sirf.STIR.
  • sirf.STIR.AcquisitionSensitivityModel.compute_attenuation_factors does not make sense for SPECT, where they attenuation image has to be entered explicitly in the AcquisitionModel. So, that function (if we keep it) needs to be moved to PETAttenuationModel (where it already is in C++).
  • randoms do not exist in SPECT

My preferred solution is to have sirf.STIR.PET_data_preparation Python module (it could be on the C++ side as well of course, but more work) where these kind of one-liners are located (as functions, not as class members). I'm not sure if that's easy to do quickly though.

Possibly best things to do in this PR is to:

• move sirf.STIR.AcquisitionSensitivityModel.compute_attenuation_factors to a PET specific model. If that doesn't exist on the Python side, I'm afraid I would not add this to the class.
• drop the list-mode one-liner

Aside from this, there are huge diffs in stir_x.h, probably just due to moving things around. I haven't reviewed any of that.

[KrisThielemans]

I am not a fan of poluting class interfaces with multiple options to do the same thing (although we have a fair amount of that in SIRF already), but I suppose others might prefer oneliners. Nevertheless, I certainly don't want to have all the above one-liners in class interfaces (it wouldn't make sense anyway).

I still think this particular function has a few problem:

• it is PET specific (randoms don't exist in SPECT). This was already slightly problematic in ListModeToSinograms but at least people would have to explicitly call estimate_randoms, which nobody in SPECT would do.
• for dynamic PET, people might need to call get_time_at_which_num_prompts_exceeds_threshold first. To do that, they need to call set_input first. It's contour-intuititive (and suboptimal) that this function then needs the listmode data again, so you'd need to remove the lm_data argument.
• it calls set_time_interval etc, so will modify the underlying ListModeToSinograms object, which could be surprising, at least without documentation.
• (as stated above, the name is wrong)

See general comment in the review.

[evgueni-ovtchinnikov]

randoms don't exist in SPECT

So, current ListmodeToSinograms will not work for SPECT data? Does STIR now has a converter that works for SPECT?

[KrisThielemans]

Listmode 2 sinograms wil work for SPECT as long as you don't ask for the randoms.

[evgueni-ovtchinnikov]

Listmode 2 sinograms wil work for SPECT as long as you don't ask for the randoms.

So far as I can see, ListmodeToSinograms::set_up assigns randoms.sptr to an empty STIRAcquisitionData object. Should we then add checking the size of randoms data in estimate_randoms(), save_randoms() and get_randoms_sptr() and issuing a warning if it is zero? In SIRF 4.0 we could do a more proper fix via inheritance from a common base class.

[KrisThielemans]

Can be an error I think, but yes, good plan

[KrisThielemans]

@evgueni-ovtchinnikov let's try to finish this off and merge. Given above problems with Listmode2Sinograms::sinograms_and_randoms_from_listmode I suggest to cut it out, and keep relevant lines in the python code. you could make it a function there obviously.

Then add something to CHANGES.md and go ahead.

[KrisThielemans]

I'm sorry, but although this last commit is fine, I don't like it. I find introducing multiple ways to the class interface to do something in the end less productive. How are we going to handle the randoms now? Add some other member to attempt a one liner? Where do we stop?

We need to move on, and we keep circling. so please abandon this PR and just add equivalent functions in the Challenge PR. Sorry.

We could pick it up later, but it'd need design decisions.

• compute_ac_factors
  • C++: it is a static function of PETAttenuationModel, that takes a PETAttenuationModel as argument, that's pretty weird.
  • Python: move sirf.STIR.AcquisitionSensitivityModel.compute_attenuation_factors to a sirf.STIR.PETAttenuationModel. However, that doesn't exist on the Python side yet, so I don't know how much work that'd be. This would be worth doing though, as really the current AcquisitionSensitivityModel constructor taking an ImageData as argument should sit in PETAttenuationModel. As that breaks backwards compatibility, it's for 4.0.
• remove prompts_*_from_listmode completely.
• remove test_data_preparation.py
• change CHANGES.md to say you've added a convenience function compute_af_factors
• make sure that the attenuation model is NOT returned from compute_af_factors
• rename test7.cxx to something that says what it's testing.

[KrisThielemans]

this won't make sense to users.

[KrisThielemans]

Right, I see this was merged. Please revert the merge. Sorry