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.

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

there are tests in the recon_test_pack I believe.

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?

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)

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

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

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

I rebased onto master but GitHub doesn't handle that well it seems. I have reverted the rebase.

This is failing because of a divide by 0. This because the projdata is 192 tangential positions but the fan size is calculated as 191 by

   const int half_fan_size = min(proj_data_info.get_max_tangential_pos_num(), -proj_data_info.get_min_tangential_pos_num()); 
   fan_size = 2 * half_fan_size + 1; 

in here. There seems to be some requirement for the fan size to be odd?

Regardless, during BindNormalisationPETFromComponents::set_up is called it constructs the projdata, which eventually reaches

The issue is in the last for loop and the bin.tangential_pos_num() being initialized at -half_fan_size = -95 whereas the proj_data.get_min_tangential_pos_num() returns -96. 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.

There seems to be some requirement for the fan size to be odd?

yes, that is the case at the moment. Code could be enhanced, but I was too lazy (having 1% extra data won't make or break the efficiency estimation)

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)

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

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

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.

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

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.

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

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.