Moving SyneRBI-Challenge nema-data utilities to SIRF

evgueni-ovtchinnikov commented 2 months ago

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
[x] 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:

[x] 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 commented 2 months ago

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 commented 2 months ago

@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 commented 2 months ago

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 commented 2 months ago

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

evgueni-ovtchinnikov commented 1 month ago

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 commented 1 month ago

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

evgueni-ovtchinnikov commented 1 month ago

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 commented 1 month ago

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

KrisThielemans commented 1 month ago

@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 commented 1 month ago

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

SyneRBI / SIRF