MLEstimateComponentBasedNormalisation refactor #1499

robbietuk · 2024-08-29T20:54:55Z

Changes in this pull request

WIP

Refactor of ML_estimate_component_based_normalisation.

Created a class.
Simplified the system logic
Add the ability option to save the data or not

Remaining tasks

Add methods to get_efficiencies(), get_norm_geo_data() and get_norm_block_data()
Add methods to get and set the boolean values of do_geo, do_block, do_symmetry_per_block, do_KL, do_display, do_save_to_file
Improve documentation defining the efficiencies, norm_geo_data, norm_block_data variables (even after all this work I am still slightly confused).
Add unit tests

Testing performed

Added a new unit test for basic functionality.

Related issues

Fixes #1498

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 (if applicable)
[] The code builds and runs on my machine
[] documentation/release_XXX.md has been updated with any functionality change (if applicable)

…ions

… method arguments

robbietuk · 2024-08-29T20:58:22Z

Thus far, I have not (intentionally) changed any of the functionality of the find_ML_normfactors3D utility. I have only used refactoring to move the logic into a class, split functionality into component methods, created members.

The only new introduction is do_save_to_file.

robbietuk

A few comments. I am not 100% sure about updating to using shared_ptrs but the idea is that the data can be accessed without copying.

Additionally, I am not convinced there is any testing of this code.

src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h

src/recon_buildblock/ML_estimate_component_based_normalisation.cxx

src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h

KrisThielemans

some initial comments. I didn't review everything.

I'm not sure if we should provide access to the efficiencies_sptr etc. shared_ptr are nice, but to be avoided when we can. If we have a const auto get_efficiencies() member, there seems tobe no point in making efficiencies etc into a shared_ptr.

src/buildblock/ML_norm.cxx

src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h

src/recon_buildblock/ML_estimate_component_based_normalisation.cxx

src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h

KrisThielemans · 2024-08-30T08:18:02Z

Additionally, I am not convinced there is any testing of this code.

there are tests in the recon_test_pack I believe.

robbietuk · 2024-08-30T18:44:47Z

At this point, I want to add an interface with BinNormalisationPETFromComponents. Is the any specific reason for no direct setters of the eff, geo, block factors into BinNormalisationPETFromComponents?

STIR/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h

Lines 60 to 62 in 33a0106

    
             This class does not actually set the relevant factors. That is left to external 
        
             methods via the crystal_efficiencies(), geometric_factors() and block_factors() 
        
             members. If they are not set, the factor is not applied (i.e. assumed to be 1).

apply_normfactors3D uses streams to set these values

STIR/src/utilities/apply_normfactors3D.cxx

Lines 71 to 83 in 33a0106

    
           if (do_eff) 
        
             { 
        
               char* in_filename = new char[in_filename_prefix.size() + 30]; 
        
               sprintf(in_filename, "%s_%s_%d_%d.out", in_filename_prefix.c_str(), "eff", iter_num, eff_iter_num); 
        
               std::ifstream in(in_filename); 
        
               in >> norm.crystal_efficiencies(); 
        
               if (!in) 
        
                 { 
        
                   warning("Error reading %s, using all 1s instead\n", in_filename); 
        
                   do_eff = false; 
        
                 } 
        
               delete[] in_filename; 
        
             }

. I would prefer not to do this.

robbietuk · 2024-08-30T20:01:38Z

src/recon_buildblock/ML_estimate_component_based_normalisation.cxx

+    const ProjData& measured_projdata_v,
+    const ProjData& model_projdata_v,


Should there be a check that these ProjData are the same "shape". Are there any tests for that?

They should be. Not sure if it's checked anywhere.

Adds has_processed_data to MLEstimateComponentBasedNormalisation

KrisThielemans · 2024-08-30T20:22:10Z

Is the any specific reason for no direct setters of the eff, geo, block factors into BinNormalisationPETFromComponents?

I was just running out of time...

Note that BinNormalisationPETFromComponents was intended to form a parent class for other component-based norms, including Siemens and GE. Not 100% sure how feasible that'll be, but it'd remove a lot of obscure/copied code. I'm not saying this because you need to do it, but just to be aware of the longer term plan. Clearly, adding setters/getters is not a problem.

On another note, the "block" stuff really doesn't work and shouldn't be used. It was intended to cope with timing-alignment factors (for many scanners these are per block) but it introduces too much freedom (which I didn't know when I wrote that code). Timing alignment should be parametrised by 1 factor per block, not a norm-factor for every block-pair. Just saying that such that you don't spend too much time on that (i.e. maybe don't expose it)

robbietuk · 2024-08-30T20:24:40Z

src/recon_buildblock/BinNormalisationPETFromComponents.cxx

+void
+BinNormalisationPETFromComponents::allocate(shared_ptr<const ProjDataInfo> pdi_sptr, MLEstimateComponentBasedNormalisation& normalization_estimator)
+{
+  if (!normalization_estimator.has_processed_data())
+    {
+      error("BinNormalisationPETFromComponents: internal error: allocate called on MLEstimateComponentBasedNormalisation "
+            "without it having processed the data");
+    }
+  this->proj_data_info_sptr = pdi_sptr;
+  this->efficiencies = *normalization_estimator.get_efficiencies();
+  this->geo_data = *normalization_estimator.get_geo_data();
+  this->block_data = *normalization_estimator.get_block_data();
+  this->_already_allocated = true;
+}


I believe this is the cleanest way to implement all these changes. Here is a use case from my cpp program.

MLEstimateComponentBasedNormalisation ml_estimator(out_filename_tmp, *measured_projdata, *model_activity_projdata, num_eff_iterations, num_iterations, do_geo, do_block, do_symmetry_per_block, do_KL, do_display, do_save_to_file); ml_estimator.process(); // Apply the normalization factors BinNormalisationPETFromComponents bin_normalisation; bin_normalisation.allocate(measured_projdata->get_proj_data_info_sptr(), ml_estimator); bin_normalisation.set_up(measured_projdata->get_exam_info_sptr(), measured_projdata->get_proj_data_info_sptr()); // Compute the projdata ProjDataInMemory normalization_projdata(*measured_projdata); normalization_projdata.fill(1.F); bin_normalisation.undo(normalization_projdata);

I acknowledge #1498 (comment), however this does work and provides an alternative allocation of the data.

I'd rather do it the other way around I'm afraid. Your way is very different from anything else in STIR.

I very much prefer to let MLEstimateComponentBasedNormalisation have a shared_ptr<BinNormalisationPETFromComponents> normalisation_sptr member, which we can get via get_output() (returning a shared_ptr<const BinNorm...>). The MLEstimateComponentBasedNormalisation::get_efficiencies() would then just forward to normalisation_sptr->get_efficiencies() (note that I wouldn't let those return shared_ptr then, but references). Any reason not to do it this way?

robbietuk · 2024-08-31T00:01:19Z

Currently seeing the oddest error. The minimum tangential position is filled with inf.

TEST(TestNormalizationCalculator, TBD)
{
    auto scanner = PETScanner();
    auto measured_projdata = get_empty_projdata(scanner, 1);
    measured_projdata->fill(1.F);

    auto model_projdata = std::make_shared<ProjDataInMemory>(*measured_projdata);
    model_projdata->fill(1.F);


    int num_eff_iterations = 1;
    int num_iterations = 1;
    bool do_geo = false;
    bool do_block = false;
    bool do_symmetry_per_block = false;
    bool do_KL = false;
    bool do_display = false;
    bool do_save_to_file = false;
    stir::MLEstimateComponentBasedNormalisation ml_estimator("", *measured_projdata, *model_projdata,
                                                             num_eff_iterations, num_iterations, do_geo, do_block,
                                                             do_symmetry_per_block, do_KL, do_display, do_save_to_file);
    ml_estimator.process();

    // Apply the normalization factors
    BinNormalisationPETFromComponents bin_normalisation;
    bin_normalisation.allocate(measured_projdata->get_proj_data_info_sptr(), ml_estimator);
    bin_normalisation.set_up(measured_projdata->get_exam_info_sptr(), measured_projdata->get_proj_data_info_sptr());

    // Compute the projdata
    ProjDataInMemory normalization_projdata(*measured_projdata);
    normalization_projdata.fill(1.F);
    bin_normalisation.apply(normalization_projdata);
    std::cout << "Complete calculating normalization factors" << std::endl;
    std::cout << "Max: " << normalization_projdata.find_max() << std::endl;
    std::cout << "Min: " << normalization_projdata.find_min() << std::endl;
    auto filename = "test_norm";
    normalization_projdata.write_to_file(filename);

It outputs

Complete calculating normalization factors
Max: inf
Min: 1

Visualization

All voxels are 1 for every slice, except the tangential_position = 0.

I don't think this is a preexisting issue but I will check the following later

find_ML_normfactors3D.exe test_exe .\ones.hs .\twos.hs 1 1
apply_normfactors3D.exe test_apply_exe test_exe .\ones.hs 0 1 1

KrisThielemans · 2024-08-31T07:12:50Z

Infinities tend to come from 1/0. A candidate would be the TEMP FIX stuff.

KrisThielemans · 2024-08-31T07:13:29Z

Note that you'll have to rename class and file to MLEstimateComponentBasedNormalisation (no _ in class names according to STIR conventions)

robbietuk · 2024-09-05T17:56:05Z

Therefore, the first row tangential_pos is treated as a gap and set to 0. This 0 is later used in the denominator of a division that leads to the infinity values in the test.

yes, so the TEMP FIX could/should handle this (the model has to be zero in that bin). Or you could trim the proj-data first (with SSRB)

The TEMP FIX code sits in the MLEstimateComponentBasedNormalisation class. This class correctly populates the DetectorEfficiencies with uniform 0.707* values in the test. It also has the shape (24,384) = (num_rings, num_det_per_ring). I'm not sure that is the problem.

The issue is when the DetectorEfficiencies are applied to the invnorm_proj_data_sptr in BindNormalisationPETFromComponents::set_up. Because of the half_fan_size calculation, the min tangential position is not accounted for in the efficiencies.

I am not a huge fan of trimming the projdata with SSRB. I am trying to understand the code to figure out why the min_tangential_pos_num is not able to be populated

KrisThielemans · 2024-09-05T22:46:23Z

It's not the efficiencies themselves that are the problem, but the FanData

STIR/src/buildblock/ML_norm.cxx

Lines 745 to 746 in e168b07

    
           fan_indices[ra][a][rb] 
        
               = IndexRange<1>(a + num_detectors_per_ring / 2 - half_fan_size, a + num_detectors_per_ring / 2 + half_fan_size);

Moreover, there's plenty of loops doing

for (bin.tangential_pos_num() = -half_fan_size; bin.tangential_pos_num() <= half_fan_size; ++bin.tangential_pos_num())

This isn't going to be so easy to fix.

KrisThielemans · 2024-09-06T09:06:08Z

Note that #1430 updates the ML_norm files as well. We didn't manage to finish that one...

robbietuk · 2024-09-06T20:44:12Z

I am not up to date with the theory in the papers. Do you expect any issue if I fix the fan data by removing the half_fan_size and fan_size calculations and replace them with methods that better derive from get_num_tangential_poss(). I feel this is the problem now.

This allows the test to pass and ignores the to ignore even num_tangential_poss issue

robbietuk · 2024-09-12T20:45:51Z

I changed the scanner to "Discovery 690", with an odd number of tangential positions (381 I believe). This allows the test to pass.

The previous issue with 0 and +inf values in the resulting projection data is pushed to issue #1511

robbietuk

A few check before I believe this is done.

robbietuk · 2024-09-12T21:56:50Z

src/utilities/apply_normfactors3D.cxx

        sprintf(in_filename, "%s_%s_%d_%d.out", in_filename_prefix.c_str(), "eff", iter_num, eff_iter_num);
        std::ifstream in(in_filename);
-        in >> norm.crystal_efficiencies();
+        in >> norm.get_crystal_efficiencies();


Returns by reference. This allows backwards compatibility with this code snippet but in general get_crystal_efficiencies() method returning the efficiencies object by reference might be dangerous?

const & would do the trick

src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx

src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h

…ests function Also documentation

robbietuk · 2024-09-13T14:10:28Z

Pretty sure this PR is ready for review. I can't find any more outstanding issues. #1511 covers the other issue with the fan data and even numbers. Ill add to the release notes.

…based_normalisation_redesign

…redesign

robbietuk · 2024-12-22T23:24:34Z

This should be merged into v6.3. I have been using a branch derived from here for a while with no major issues.

KrisThielemans

Small changes required, but unfortunately quite a lot of them.

We cannot rename functions/files in 6.3, so I've flagged those up. Deprecations need to be listed in the release notes.

KrisThielemans · 2024-12-23T18:39:09Z

src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h

+    efficiencies = eff;
+  }
+
+  DetectorEfficiencies& get_crystal_efficiencies()


ok for changing to get/set pair, but then I think we should have

Suggested change

DetectorEfficiencies& get_crystal_efficiencies()

const DetectorEfficiencies& get_crystal_efficiencies() const

or would this break something?

However, cannot change backwards compatibility in a minor version update, so keep

#if STIR_VERSION < 070000 STIR_DEPRECATED DetectorEfficiencies& crystal_efficiencies() ... #endif Obviously the same for the others.

KrisThielemans · 2024-12-23T18:39:54Z

src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h

+    geo_data = geo;
+  }
+
+  GeoData3D& get_geometric_factors()


also const

KrisThielemans · 2024-12-23T18:40:04Z

src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h

+    block_data = block;
+  }
+
+  BlockData3D& get_block_factors()


also const

KrisThielemans · 2024-12-23T18:40:55Z

src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h

+/*!
+ \file
+ \ingroup recon_buildblock
+ \brief Declaration of stir::ML_estimate_component_based_normalisation


Suggested change

\brief Declaration of stir::ML_estimate_component_based_normalisation

\brief Declaration of stir::ML_estimate_component_based_normalisation and

stir::MLEstimateComponentBasedNormalisation

obviously, need to adjust when keeping ML_estimate_component_based_normalisation in its original file.

KrisThielemans · 2024-12-23T18:42:01Z

src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h

+    \param do_display_v Whether to display the progress of the algorithm.
+    \param do_save_to_file_v Whether to save the each iteration of the efficiencies, geo data and block data to file.
+  */
+  MLEstimateComponentBasedNormalisation(std::string out_filename_prefix_v,


generally STIR uses const std::string&

KrisThielemans · 2024-12-23T22:52:00Z

src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx

+  \ingroup recon_test
+  \brief Test program to ensure MLEstimateComponentBasedNormalisation works as expected and
+         produces the correct normalization factors from components.
+         This is a very basic test that ensure the basic functionality.


Suggested change

This is a very basic test that ensure the basic functionality.

This is a test on basic functionality only.

KrisThielemans · 2024-12-23T22:53:18Z

src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx

+    const auto scanner = std::make_shared<Scanner>(*Scanner::get_scanner_from_name("Discovery 690"));
+    const auto exam_info = std::make_shared<ExamInfo>(ImagingModality::PT);


rename to scanner_sptr and exam_info_sptr

KrisThielemans · 2024-12-23T22:53:45Z

src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx

+    const auto exam_info = std::make_shared<ExamInfo>(ImagingModality::PT);
+    exam_info->patient_position = PatientPosition(PatientPosition::HFS);
+
+    const shared_ptr<ProjDataInfo> projdata_info(ProjDataInfo::ProjDataInfoCTI(


use ProjDataInfo::construct_proj_data_info (same signature)

KrisThielemans · 2024-12-23T22:54:10Z

src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx

+        /*views*/ scanner->get_num_detectors_per_ring() / 2,
+        /*tang_pos*/ scanner->get_default_num_arccorrected_bins(),
+        /*arc_corrected*/ false,
+        /*tof_mash_factor*/ -1));


why -1? Set to 0 for non-TOF.

KrisThielemans · 2024-12-23T22:54:52Z

src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx

+    const auto model_projdata = std::make_shared<ProjDataInMemory>(*measured_projdata);
+    model_projdata->fill(2.F);
+
+    constexpr int num_eff_iterations = 6;


Can we make this smaller (to run faster)?

KrisThielemans · 2024-12-23T22:58:44Z

Also, you'll need to run pre-commit run --all (see the devel doc). Also check the Codacy report.

KrisThielemans · 2025-02-07T11:57:29Z

@robbietuk where are we with this? I think it might be best to postpone for after 6.3, as that needs to be completed soon.

robbietuk · 2025-02-07T17:39:53Z

Postpone for now unless you can handle the review items. Sorry, I don't have much free time for this

robbietuk added 5 commits August 29, 2024 11:08

Refactor ML_estimate_component_based_normalisation creating sub-funct…

3035db8

…ions

Create MLEstimateComponentBasedNormalisation class and start to strip…

bc58f39

… method arguments

Refactor to use class varables

0dca08e

Add do_save_to_file guards

9fed40f

Documentation and minor refactors

2579f81

robbietuk changed the title ~~Ml estimate component based normalisation redesign~~ MLEstimateComponentBasedNormalisation refactor Aug 29, 2024

robbietuk self-assigned this Aug 29, 2024

robbietuk added the enhancement label Aug 29, 2024

robbietuk added 2 commits August 29, 2024 14:18

Rename refactor

5d89071

Addition usage of shared_ptrs

86f6a82

robbietuk commented Aug 29, 2024

View reviewed changes

[ci skip] parameter rename, remove redundant logic, documentation

880e530

KrisThielemans requested changes Aug 30, 2024

View reviewed changes

robbietuk added 3 commits August 30, 2024 11:27

Remove model_data member and pass by reference into the constructor

dd711b1

Readd stir/steam.h to includes

2bcdc80

[ci skip] Update copywrite

0fa8455

robbietuk commented Aug 30, 2024

View reviewed changes

robbietuk added 2 commits August 30, 2024 13:03

[ci skip] Fix a bug with do_KL info print

90ea81b

Add allocate override using a MLEstimateComponentBasedNormalisation

f0c788b

Adds has_processed_data to MLEstimateComponentBasedNormalisation

robbietuk commented Aug 30, 2024

View reviewed changes

[ci skip] cleanup comments

e5e2cc7

robbietuk added 2 commits September 3, 2024 14:13

Rename files MLEstimateComponentBasedNormalisation

8adb620

Remove allocate with norm estimator

9ff77fe

robbietuk mentioned this pull request Sep 6, 2024

SAFIR updates #1502

Closed

3 tasks

robbietuk added 3 commits September 12, 2024 13:10

Change scanner to Discovery

7c18bc7

This allows the test to pass and ignores the to ignore even num_tangential_poss issue

Add additional test on efficiency values

42340f2

Ensure normalization_projdata max and min values are about 2.0

f424ca8

robbietuk mentioned this pull request Sep 12, 2024

ML FanData is limited to an odd number of detectors #1511

Open

robbietuk commented Sep 13, 2024

View reviewed changes

robbietuk added 4 commits September 13, 2024 06:43

Make test use check_if_equal

a5f4b91

Split test_normalization_calculation_with_efficiencies from the run_t…

1eefeff

…ests function Also documentation

Improve usage of data_is_processed()

dcee251

Add getters and setters

63f8560

robbietuk requested a review from KrisThielemans September 13, 2024 14:07

robbietuk added 5 commits September 13, 2024 09:26

[ci skip] Add const to args

3d7e9ff

Merge remote-tracking branch 'UCL/master' into ML_estimate_component_…

55e5f1a

…based_normalisation_redesign

[ci skip] Update release 6.3 notes

55897ee

[ci skip] Fix error message

cbe81ec

Merge branch 'master' into ML_estimate_component_based_normalisation_…

fd194fc

…redesign

robbietuk added this to the v6.3 milestone Dec 22, 2024

KrisThielemans requested changes Dec 23, 2024

View reviewed changes

KrisThielemans modified the milestones: v6.3, v6.4 Feb 7, 2025

		const ProjData& measured_projdata_v,
		const ProjData& model_projdata_v,

	DetectorEfficiencies& get_crystal_efficiencies()
	const DetectorEfficiencies& get_crystal_efficiencies() const

	\brief Declaration of stir::ML_estimate_component_based_normalisation
	\brief Declaration of stir::ML_estimate_component_based_normalisation and
	stir::MLEstimateComponentBasedNormalisation

	This is a very basic test that ensure the basic functionality.
	This is a test on basic functionality only.

		const auto scanner = std::make_shared<Scanner>(*Scanner::get_scanner_from_name("Discovery 690"));
		const auto exam_info = std::make_shared<ExamInfo>(ImagingModality::PT);

MLEstimateComponentBasedNormalisation refactor #1499

Are you sure you want to change the base?

MLEstimateComponentBasedNormalisation refactor #1499

Uh oh!

Conversation

robbietuk commented Aug 29, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Changes in this pull request

Remaining tasks

Testing performed

Related issues

Checklist before requesting a review

Uh oh!

robbietuk commented Aug 29, 2024

Uh oh!

robbietuk left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

KrisThielemans left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

KrisThielemans commented Aug 30, 2024

Uh oh!

robbietuk commented Aug 30, 2024

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KrisThielemans commented Aug 30, 2024

Uh oh!

Choose a reason for hiding this comment

Uh oh!

robbietuk Aug 30, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

robbietuk commented Aug 31, 2024

Uh oh!

KrisThielemans commented Aug 31, 2024

Uh oh!

KrisThielemans commented Aug 31, 2024

Uh oh!

robbietuk commented Sep 5, 2024

Uh oh!

KrisThielemans commented Sep 5, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

KrisThielemans commented Sep 6, 2024

Uh oh!

robbietuk commented Sep 6, 2024

Uh oh!

robbietuk commented Sep 12, 2024

Uh oh!

robbietuk left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

robbietuk commented Aug 29, 2024 •

edited

Loading

robbietuk left a comment •

edited

Loading

robbietuk Aug 30, 2024 •

edited

Loading

KrisThielemans commented Sep 5, 2024 •

edited

Loading