From 3035db86beb58ad048eebc71821c30081289a0e1 Mon Sep 17 00:00:00 2001 From: rts Date: Thu, 29 Aug 2024 11:08:02 -0700 Subject: [PATCH 01/34] Refactor ML_estimate_component_based_normalisation creating sub-functions --- ...estimate_component_based_normalisation.cxx | 461 +++++++++++------- 1 file changed, 274 insertions(+), 187 deletions(-) diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index a22ed6b0b..298e19b43 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -35,6 +35,196 @@ START_NAMESPACE_STIR +//! Helper function to write efficiencies to a file +void +write_efficiencies_to_file(const std::string& out_filename_prefix, + int iter_num, + int eff_iter_num, + const DetectorEfficiencies& efficiencies) +{ + char* out_filename = new char[out_filename_prefix.size() + 30]; + sprintf(out_filename, "%s_%s_%d_%d.out", out_filename_prefix.c_str(), "eff", iter_num, eff_iter_num); + std::ofstream out(out_filename); + out << efficiencies; +} + +//! Helper function to write geo data to a file +void +write_geo_data_to_file(const std::string& out_filename_prefix, int iter_num, const GeoData3D& norm_geo_data) +{ + char* out_filename = new char[out_filename_prefix.size() + 30]; + sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "geo", iter_num); + std::ofstream out(out_filename); + out << norm_geo_data; +} + +//! Helper function to write block data to a file +void +write_block_data_to_file(const std::string& out_filename_prefix, int iter_num, const BlockData3D& norm_block_data) +{ + char* out_filename = new char[out_filename_prefix.size() + 30]; + sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "block", iter_num); + std::ofstream out(out_filename); + out << norm_block_data; +} + +// Function to compute factors dependent on the data +void +compute_initial_data_dependent_factors(const ProjData& measured_data, + const FanProjData& model_fan_data, + FanProjData& measured_fan_data, + DetectorEfficiencies& data_fan_sums, + GeoData3D& measured_geo_data, + BlockData3D& measured_block_data, + bool do_display, + float& threshold_for_KL) +{ + make_fan_data_remove_gaps(measured_fan_data, measured_data); + + /* TEMP FIX */ + for (int ra = model_fan_data.get_min_ra(); ra <= model_fan_data.get_max_ra(); ++ra) + { + for (int a = model_fan_data.get_min_a(); a <= model_fan_data.get_max_a(); ++a) + { + for (int rb = std::max(ra, model_fan_data.get_min_rb(ra)); rb <= model_fan_data.get_max_rb(ra); ++rb) + { + for (int b = model_fan_data.get_min_b(a); b <= model_fan_data.get_max_b(a); ++b) + if (model_fan_data(ra, a, rb, b) == 0) + measured_fan_data(ra, a, rb, b) = 0; + } + } + } + + threshold_for_KL = measured_fan_data.find_max() / 100000.F; + // display(measured_fan_data, "measured data"); + + make_fan_sum_data(data_fan_sums, measured_fan_data); + make_geo_data(measured_geo_data, measured_fan_data); + make_block_data(measured_block_data, measured_fan_data); + if (do_display) + display(measured_block_data, "raw block data from measurements"); +} + +// Function to handle efficiency iteration +void +efficiency_iteration(FanProjData& fan_data, + const FanProjData& model_fan_data, + const GeoData3D& norm_geo_data, + const BlockData3D& norm_block_data, + DetectorEfficiencies& efficiencies, + const DetectorEfficiencies& data_fan_sums, + const std::string& out_filename_prefix, + int iter_num, + int eff_iter_num, + bool do_KL, + bool do_display, + const FanProjData& measured_fan_data, + float threshold_for_KL) +{ + iterate_efficiencies(efficiencies, data_fan_sums, fan_data); + write_efficiencies_to_file(out_filename_prefix, iter_num, eff_iter_num, efficiencies); + if (do_KL) + { + apply_efficiencies(fan_data, efficiencies); + std::cerr << "measured*norm min " << measured_fan_data.find_min() << " ,max " << measured_fan_data.find_max() << std::endl; + std::cerr << "model*norm min " << fan_data.find_min() << " ,max " << fan_data.find_max() << std::endl; + if (do_display) + display(fan_data, "model_times_norm"); + info(boost::format("KL %1%") % KL(measured_fan_data, fan_data, threshold_for_KL)); + // now restore for further iterations + fan_data = model_fan_data; + apply_geo_norm(fan_data, norm_geo_data); + apply_block_norm(fan_data, norm_block_data); + } + if (do_display) + { + fan_data.fill(1); + apply_efficiencies(fan_data, efficiencies); + display(fan_data, "eff norm"); + // now restore for further iterations + fan_data = model_fan_data; + apply_geo_norm(fan_data, norm_geo_data); + apply_block_norm(fan_data, norm_block_data); + } +} + +// Function to handle geo normalization +void +geo_normalization_iteration(FanProjData& fan_data, + const FanProjData& model_fan_data, + const DetectorEfficiencies& efficiencies, + const BlockData3D& norm_block_data, + GeoData3D& norm_geo_data, + const GeoData3D& measured_geo_data, + const std::string& out_filename_prefix, + int iter_num, + bool do_geo, + bool do_KL, + bool do_display, + const FanProjData& measured_fan_data, + float threshold_for_KL) +{ + fan_data = model_fan_data; + apply_efficiencies(fan_data, efficiencies); + apply_block_norm(fan_data, norm_block_data); + + if (do_geo) + { + iterate_geo_norm(norm_geo_data, measured_geo_data, fan_data); + } + + write_geo_data_to_file(out_filename_prefix, iter_num, norm_geo_data); + if (do_KL) + { + apply_geo_norm(fan_data, norm_geo_data); + info(boost::format("KL %1%") % KL(measured_fan_data, fan_data, threshold_for_KL)); + } + if (do_display) + { + fan_data.fill(1); + apply_geo_norm(fan_data, norm_geo_data); + display(fan_data, "geo norm"); + } +} + +// Function to handle block normalization +void +block_normalization_iteration(FanProjData& fan_data, + const FanProjData& model_fan_data, + const DetectorEfficiencies& efficiencies, + const GeoData3D& norm_geo_data, + BlockData3D& norm_block_data, + const BlockData3D& measured_block_data, + const std::string& out_filename_prefix, + int iter_num, + bool do_block, + bool do_KL, + bool do_display, + const FanProjData& measured_fan_data, + float threshold_for_KL) +{ + fan_data = model_fan_data; + apply_efficiencies(fan_data, efficiencies); + apply_geo_norm(fan_data, norm_geo_data); + if (do_block) + { + iterate_block_norm(norm_block_data, measured_block_data, fan_data); + } + write_block_data_to_file(out_filename_prefix, iter_num, norm_block_data); + if (do_KL) + { + apply_block_norm(fan_data, norm_block_data); + info(boost::format("KL %1%") % KL(measured_fan_data, fan_data, threshold_for_KL)); + } + if (do_display) + { + fan_data.fill(1); + apply_block_norm(fan_data, norm_block_data); + display(norm_block_data, "raw block norm"); + display(fan_data, "block norm"); + } +} + void ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, const ProjData& measured_data, @@ -102,207 +292,104 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix BlockData3D norm_block_data(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); make_fan_data_remove_gaps(model_fan_data, model_data); - { - // next could be local if KL is not computed below - FanProjData measured_fan_data; - float threshold_for_KL; - // compute factors dependent on the data - { - make_fan_data_remove_gaps(measured_fan_data, measured_data); - /* TEMP FIX */ - for (int ra = model_fan_data.get_min_ra(); ra <= model_fan_data.get_max_ra(); ++ra) - { - for (int a = model_fan_data.get_min_a(); a <= model_fan_data.get_max_a(); ++a) - { - for (int rb = std::max(ra, model_fan_data.get_min_rb(ra)); rb <= model_fan_data.get_max_rb(ra); ++rb) - { - for (int b = model_fan_data.get_min_b(a); b <= model_fan_data.get_max_b(a); ++b) - if (model_fan_data(ra, a, rb, b) == 0) - measured_fan_data(ra, a, rb, b) = 0; - } - } - } + // next could be local if KL is not computed below + FanProjData measured_fan_data; + float threshold_for_KL; + compute_initial_data_dependent_factors(measured_data, + model_fan_data, + measured_fan_data, + data_fan_sums, + measured_geo_data, + measured_block_data, + do_display, + threshold_for_KL); - threshold_for_KL = measured_fan_data.find_max() / 100000.F; - // display(measured_fan_data, "measured data"); + // std::cerr << "model min " << model_fan_data.find_min() << " ,max " << model_fan_data.find_max() << std::endl; + if (do_display) + display(model_fan_data, "model"); - make_fan_sum_data(data_fan_sums, measured_fan_data); - make_geo_data(measured_geo_data, measured_fan_data); - make_block_data(measured_block_data, measured_fan_data); + for (int iter_num = 1; iter_num <= std::max(num_iterations, 1); ++iter_num) + { + if (iter_num == 1) + { + efficiencies.fill(sqrt(data_fan_sums.sum() / model_fan_data.sum())); + norm_geo_data.fill(1); + norm_block_data.fill(1); + } + fan_data = model_fan_data; + apply_geo_norm(fan_data, norm_geo_data); + apply_block_norm(fan_data, norm_block_data); if (do_display) - display(measured_block_data, "raw block data from measurements"); - - /* { - char *out_filename = new char[20]; - sprintf(out_filename, "%s_%d.out", - "fan", ax_pos_num); - std::ofstream out(out_filename); - out << data_fan_sums; - delete[] out_filename; + { + display(fan_data, "model*geo*block"); } - */ - } - // std::cerr << "model min " << model_fan_data.find_min() << " ,max " << model_fan_data.find_max() << std::endl; - if (do_display) - display(model_fan_data, "model"); -#if 0 + // Efficiency iterations + for (int eff_iter_num = 1; eff_iter_num <= num_eff_iterations; ++eff_iter_num) { - shared_ptr out_proj_data_ptr = - new ProjDataInterfile(model_data.get_proj_data_info_sptr()->clone, - output_file_name); - - set_fan_data(*out_proj_data_ptr, model_fan_data); + efficiency_iteration(fan_data, + model_fan_data, + norm_geo_data, + norm_block_data, + efficiencies, + data_fan_sums, + out_filename_prefix, + iter_num, + eff_iter_num, + do_KL, + do_display, + measured_fan_data, + threshold_for_KL); } -#endif - - for (int iter_num = 1; iter_num <= std::max(num_iterations, 1); ++iter_num) - { - if (iter_num == 1) - { - efficiencies.fill(sqrt(data_fan_sums.sum() / model_fan_data.sum())); - norm_geo_data.fill(1); - norm_block_data.fill(1); - } - // efficiencies - { - fan_data = model_fan_data; - apply_geo_norm(fan_data, norm_geo_data); - apply_block_norm(fan_data, norm_block_data); - if (do_display) - display(fan_data, "model*geo*block"); - for (int eff_iter_num = 1; eff_iter_num <= num_eff_iterations; ++eff_iter_num) - { - iterate_efficiencies(efficiencies, data_fan_sums, fan_data); - { - char* out_filename = new char[out_filename_prefix.size() + 30]; - sprintf(out_filename, "%s_%s_%d_%d.out", out_filename_prefix.c_str(), "eff", iter_num, eff_iter_num); - std::ofstream out(out_filename); - out << efficiencies; - delete[] out_filename; - } - if (do_KL) - { - apply_efficiencies(fan_data, efficiencies); - std::cerr << "measured*norm min " << measured_fan_data.find_min() << " ,max " << measured_fan_data.find_max() - << std::endl; - std::cerr << "model*norm min " << fan_data.find_min() << " ,max " << fan_data.find_max() << std::endl; - if (do_display) - display(fan_data, "model_times_norm"); - info(boost::format("KL %1%") % KL(measured_fan_data, fan_data, threshold_for_KL)); - // now restore for further iterations - fan_data = model_fan_data; - apply_geo_norm(fan_data, norm_geo_data); - apply_block_norm(fan_data, norm_block_data); - } - if (do_display) - { - fan_data.fill(1); - apply_efficiencies(fan_data, efficiencies); - display(fan_data, "eff norm"); - // now restore for further iterations - fan_data = model_fan_data; - apply_geo_norm(fan_data, norm_geo_data); - apply_block_norm(fan_data, norm_block_data); - } - } - } // end efficiencies - - // geo norm - fan_data = model_fan_data; - apply_efficiencies(fan_data, efficiencies); - apply_block_norm(fan_data, norm_block_data); + // geo norm + geo_normalization_iteration(fan_data, + model_fan_data, + efficiencies, + norm_block_data, + norm_geo_data, + measured_geo_data, + out_filename_prefix, + iter_num, + do_geo, + do_KL, + do_display, + measured_fan_data, + threshold_for_KL); - if (do_geo) - iterate_geo_norm(norm_geo_data, measured_geo_data, fan_data); + // block norm + block_normalization_iteration(fan_data, + model_fan_data, + efficiencies, + norm_geo_data, + norm_block_data, + measured_block_data, + out_filename_prefix, + iter_num, + do_block, + do_KL, + do_display, + measured_fan_data, + threshold_for_KL); -#if 0 + //// print KL for fansums + if (do_KL) + { + DetectorEfficiencies fan_sums(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); + GeoData3D geo_data(num_physical_axial_crystals_per_basic_unit, + num_physical_transaxial_crystals_per_basic_unit / 2, + num_physical_rings, + num_physical_detectors_per_ring); // inputes have to be modified + BlockData3D block_data(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); - { // check - for (int a=0; a Date: Thu, 29 Aug 2024 11:26:58 -0700 Subject: [PATCH 02/34] Create MLEstimateComponentBasedNormalisation class and start to strip method arguments --- ...L_estimate_component_based_normalisation.h | 78 ++++++- ...estimate_component_based_normalisation.cxx | 193 +++++++++--------- 2 files changed, 169 insertions(+), 102 deletions(-) diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h index 9fcbc0650..ae1e56f4c 100644 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h @@ -14,16 +14,18 @@ */ #include "stir/common.h" #include +#include "stir/ProjData.h" +#include "stir/ML_norm.h" START_NAMESPACE_STIR -class ProjData; - /*! \brief Find normalisation factors using a maximum likelihood approach \ingroup recon_buildblock */ + +// The function that calls the class MLEstimateComponentBasedNormalisation and runs the normalisation estimation void ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, const ProjData& measured_data, const ProjData& model_data, @@ -35,4 +37,76 @@ void ML_estimate_component_based_normalisation(const std::string& out_filename_p bool do_KL, bool do_display); +class MLEstimateComponentBasedNormalisation +{ +public: + MLEstimateComponentBasedNormalisation(const std::string& out_filename_prefix, + const ProjData& measured_data, + const ProjData& model_data, + int num_eff_iterations, + int num_iterations, + bool do_geo, + bool do_block, + bool do_symmetry_per_block, + bool do_KL, + bool do_display); + + void run(); + +private: + void write_efficiencies_to_file(int iter_num, int eff_iter_num, const DetectorEfficiencies& efficiencies); + + void write_geo_data_to_file(int iter_num, const GeoData3D& norm_geo_data); + + void write_block_data_to_file(int iter_num, const BlockData3D& norm_block_data); + + void compute_initial_data_dependent_factors(const FanProjData& model_fan_data, + FanProjData& measured_fan_data, + DetectorEfficiencies& data_fan_sums, + GeoData3D& measured_geo_data, + BlockData3D& measured_block_data); + + void efficiency_iteration(FanProjData& fan_data, + const FanProjData& model_fan_data, + const GeoData3D& norm_geo_data, + const BlockData3D& norm_block_data, + DetectorEfficiencies& efficiencies, + const DetectorEfficiencies& data_fan_sums, + int iter_num, + int eff_iter_num, + const FanProjData& measured_fan_data); + + void geo_normalization_iteration(FanProjData& fan_data, + const FanProjData& model_fan_data, + const DetectorEfficiencies& efficiencies, + const BlockData3D& norm_block_data, + GeoData3D& norm_geo_data, + const GeoData3D& measured_geo_data, + int iter_num, + const FanProjData& measured_fan_data); + + void block_normalization_iteration(FanProjData& fan_data, + const FanProjData& model_fan_data, + const DetectorEfficiencies& efficiencies, + const GeoData3D& norm_geo_data, + BlockData3D& norm_block_data, + const BlockData3D& measured_block_data, + int iter_num, + const FanProjData& measured_fan_data); + // Arguments + std::string out_filename_prefix; + const ProjData& measured_data; + const ProjData& model_data; + int num_eff_iterations; + int num_iterations; + bool do_geo; + bool do_block; + bool do_symmetry_per_block; + bool do_KL; + bool do_display; + + // Calculated values + float threshold_for_KL; +}; + END_NAMESPACE_STIR diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index 298e19b43..9069227f7 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -35,12 +35,58 @@ START_NAMESPACE_STIR +void +ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, + const ProjData& measured_data, + const ProjData& model_data, + int num_eff_iterations, + int num_iterations, + bool do_geo, + bool do_block, + bool do_symmetry_per_block, + bool do_KL, + bool do_display) +{ + MLEstimateComponentBasedNormalisation estimator(out_filename_prefix, + measured_data, + model_data, + num_eff_iterations, + num_iterations, + do_geo, + do_block, + do_symmetry_per_block, + do_KL, + do_display); + estimator.run(); +} + +MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(const std::string& out_filename_prefix, + const ProjData& measured_data, + const ProjData& model_data, + int num_eff_iterations, + int num_iterations, + bool do_geo, + bool do_block, + bool do_symmetry_per_block, + bool do_KL, + bool do_display) + : measured_data(measured_data), + model_data(model_data), + out_filename_prefix(out_filename_prefix), + num_eff_iterations(num_eff_iterations), + num_iterations(num_iterations), + do_geo(do_geo), + do_block(do_block), + do_symmetry_per_block(do_symmetry_per_block), + do_KL(do_KL), + do_display(do_display) +{} + //! Helper function to write efficiencies to a file void -write_efficiencies_to_file(const std::string& out_filename_prefix, - int iter_num, - int eff_iter_num, - const DetectorEfficiencies& efficiencies) +MLEstimateComponentBasedNormalisation::write_efficiencies_to_file(int iter_num, + int eff_iter_num, + const DetectorEfficiencies& efficiencies) { char* out_filename = new char[out_filename_prefix.size() + 30]; sprintf(out_filename, "%s_%s_%d_%d.out", out_filename_prefix.c_str(), "eff", iter_num, eff_iter_num); @@ -50,7 +96,7 @@ write_efficiencies_to_file(const std::string& out_filename_prefix, //! Helper function to write geo data to a file void -write_geo_data_to_file(const std::string& out_filename_prefix, int iter_num, const GeoData3D& norm_geo_data) +MLEstimateComponentBasedNormalisation::write_geo_data_to_file(int iter_num, const GeoData3D& norm_geo_data) { char* out_filename = new char[out_filename_prefix.size() + 30]; sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "geo", iter_num); @@ -60,7 +106,7 @@ write_geo_data_to_file(const std::string& out_filename_prefix, int iter_num, con //! Helper function to write block data to a file void -write_block_data_to_file(const std::string& out_filename_prefix, int iter_num, const BlockData3D& norm_block_data) +MLEstimateComponentBasedNormalisation::write_block_data_to_file(int iter_num, const BlockData3D& norm_block_data) { char* out_filename = new char[out_filename_prefix.size() + 30]; sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "block", iter_num); @@ -70,14 +116,11 @@ write_block_data_to_file(const std::string& out_filename_prefix, int iter_num, c // Function to compute factors dependent on the data void -compute_initial_data_dependent_factors(const ProjData& measured_data, - const FanProjData& model_fan_data, - FanProjData& measured_fan_data, - DetectorEfficiencies& data_fan_sums, - GeoData3D& measured_geo_data, - BlockData3D& measured_block_data, - bool do_display, - float& threshold_for_KL) +MLEstimateComponentBasedNormalisation::compute_initial_data_dependent_factors(const FanProjData& model_fan_data, + FanProjData& measured_fan_data, + DetectorEfficiencies& data_fan_sums, + GeoData3D& measured_geo_data, + BlockData3D& measured_block_data) { make_fan_data_remove_gaps(measured_fan_data, measured_data); @@ -107,22 +150,18 @@ compute_initial_data_dependent_factors(const ProjData& measured_data, // Function to handle efficiency iteration void -efficiency_iteration(FanProjData& fan_data, - const FanProjData& model_fan_data, - const GeoData3D& norm_geo_data, - const BlockData3D& norm_block_data, - DetectorEfficiencies& efficiencies, - const DetectorEfficiencies& data_fan_sums, - const std::string& out_filename_prefix, - int iter_num, - int eff_iter_num, - bool do_KL, - bool do_display, - const FanProjData& measured_fan_data, - float threshold_for_KL) +MLEstimateComponentBasedNormalisation::efficiency_iteration(FanProjData& fan_data, + const FanProjData& model_fan_data, + const GeoData3D& norm_geo_data, + const BlockData3D& norm_block_data, + DetectorEfficiencies& efficiencies, + const DetectorEfficiencies& data_fan_sums, + int iter_num, + int eff_iter_num, + const FanProjData& measured_fan_data) { iterate_efficiencies(efficiencies, data_fan_sums, fan_data); - write_efficiencies_to_file(out_filename_prefix, iter_num, eff_iter_num, efficiencies); + write_efficiencies_to_file(iter_num, eff_iter_num, efficiencies); if (do_KL) { apply_efficiencies(fan_data, efficiencies); @@ -150,19 +189,14 @@ efficiency_iteration(FanProjData& fan_data, // Function to handle geo normalization void -geo_normalization_iteration(FanProjData& fan_data, - const FanProjData& model_fan_data, - const DetectorEfficiencies& efficiencies, - const BlockData3D& norm_block_data, - GeoData3D& norm_geo_data, - const GeoData3D& measured_geo_data, - const std::string& out_filename_prefix, - int iter_num, - bool do_geo, - bool do_KL, - bool do_display, - const FanProjData& measured_fan_data, - float threshold_for_KL) +MLEstimateComponentBasedNormalisation::geo_normalization_iteration(FanProjData& fan_data, + const FanProjData& model_fan_data, + const DetectorEfficiencies& efficiencies, + const BlockData3D& norm_block_data, + GeoData3D& norm_geo_data, + const GeoData3D& measured_geo_data, + int iter_num, + const FanProjData& measured_fan_data) { fan_data = model_fan_data; apply_efficiencies(fan_data, efficiencies); @@ -173,7 +207,7 @@ geo_normalization_iteration(FanProjData& fan_data, iterate_geo_norm(norm_geo_data, measured_geo_data, fan_data); } - write_geo_data_to_file(out_filename_prefix, iter_num, norm_geo_data); + write_geo_data_to_file(iter_num, norm_geo_data); if (do_KL) { apply_geo_norm(fan_data, norm_geo_data); @@ -189,19 +223,14 @@ geo_normalization_iteration(FanProjData& fan_data, // Function to handle block normalization void -block_normalization_iteration(FanProjData& fan_data, - const FanProjData& model_fan_data, - const DetectorEfficiencies& efficiencies, - const GeoData3D& norm_geo_data, - BlockData3D& norm_block_data, - const BlockData3D& measured_block_data, - const std::string& out_filename_prefix, - int iter_num, - bool do_block, - bool do_KL, - bool do_display, - const FanProjData& measured_fan_data, - float threshold_for_KL) +MLEstimateComponentBasedNormalisation::block_normalization_iteration(FanProjData& fan_data, + const FanProjData& model_fan_data, + const DetectorEfficiencies& efficiencies, + const GeoData3D& norm_geo_data, + BlockData3D& norm_block_data, + const BlockData3D& measured_block_data, + int iter_num, + const FanProjData& measured_fan_data) { fan_data = model_fan_data; apply_efficiencies(fan_data, efficiencies); @@ -210,7 +239,7 @@ block_normalization_iteration(FanProjData& fan_data, { iterate_block_norm(norm_block_data, measured_block_data, fan_data); } - write_block_data_to_file(out_filename_prefix, iter_num, norm_block_data); + write_block_data_to_file(iter_num, norm_block_data); if (do_KL) { apply_block_norm(fan_data, norm_block_data); @@ -226,18 +255,8 @@ block_normalization_iteration(FanProjData& fan_data, } void -ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, - const ProjData& measured_data, - const ProjData& model_data, - int num_eff_iterations, - int num_iterations, - bool do_geo, - bool do_block, - bool do_symmetry_per_block, - bool do_KL, - bool do_display) +MLEstimateComponentBasedNormalisation::run() { - const int num_transaxial_blocks = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); const int num_axial_blocks = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); const int virtual_axial_crystals @@ -295,15 +314,9 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix // next could be local if KL is not computed below FanProjData measured_fan_data; - float threshold_for_KL; - compute_initial_data_dependent_factors(measured_data, - model_fan_data, - measured_fan_data, - data_fan_sums, - measured_geo_data, - measured_block_data, - do_display, - threshold_for_KL); + + compute_initial_data_dependent_factors( + model_fan_data, measured_fan_data, data_fan_sums, measured_geo_data, measured_block_data); // std::cerr << "model min " << model_fan_data.find_min() << " ,max " << model_fan_data.find_max() << std::endl; if (do_display) @@ -334,29 +347,14 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix norm_block_data, efficiencies, data_fan_sums, - out_filename_prefix, iter_num, eff_iter_num, - do_KL, - do_display, - measured_fan_data, - threshold_for_KL); + measured_fan_data); } // geo norm - geo_normalization_iteration(fan_data, - model_fan_data, - efficiencies, - norm_block_data, - norm_geo_data, - measured_geo_data, - out_filename_prefix, - iter_num, - do_geo, - do_KL, - do_display, - measured_fan_data, - threshold_for_KL); + geo_normalization_iteration( + fan_data, model_fan_data, efficiencies, norm_block_data, norm_geo_data, measured_geo_data, iter_num, measured_fan_data); // block norm block_normalization_iteration(fan_data, @@ -365,13 +363,8 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix norm_geo_data, norm_block_data, measured_block_data, - out_filename_prefix, iter_num, - do_block, - do_KL, - do_display, - measured_fan_data, - threshold_for_KL); + measured_fan_data); //// print KL for fansums if (do_KL) From 0dca08ef8cff1781a29d03fd66f387853eee013e Mon Sep 17 00:00:00 2001 From: rts Date: Thu, 29 Aug 2024 11:46:05 -0700 Subject: [PATCH 03/34] Refactor to use class varables --- ...L_estimate_component_based_normalisation.h | 56 +++--- ...estimate_component_based_normalisation.cxx | 187 +++++++----------- 2 files changed, 92 insertions(+), 151 deletions(-) diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h index ae1e56f4c..1d1874a10 100644 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h @@ -60,39 +60,13 @@ class MLEstimateComponentBasedNormalisation void write_block_data_to_file(int iter_num, const BlockData3D& norm_block_data); - void compute_initial_data_dependent_factors(const FanProjData& model_fan_data, - FanProjData& measured_fan_data, - DetectorEfficiencies& data_fan_sums, - GeoData3D& measured_geo_data, - BlockData3D& measured_block_data); - - void efficiency_iteration(FanProjData& fan_data, - const FanProjData& model_fan_data, - const GeoData3D& norm_geo_data, - const BlockData3D& norm_block_data, - DetectorEfficiencies& efficiencies, - const DetectorEfficiencies& data_fan_sums, - int iter_num, - int eff_iter_num, - const FanProjData& measured_fan_data); - - void geo_normalization_iteration(FanProjData& fan_data, - const FanProjData& model_fan_data, - const DetectorEfficiencies& efficiencies, - const BlockData3D& norm_block_data, - GeoData3D& norm_geo_data, - const GeoData3D& measured_geo_data, - int iter_num, - const FanProjData& measured_fan_data); - - void block_normalization_iteration(FanProjData& fan_data, - const FanProjData& model_fan_data, - const DetectorEfficiencies& efficiencies, - const GeoData3D& norm_geo_data, - BlockData3D& norm_block_data, - const BlockData3D& measured_block_data, - int iter_num, - const FanProjData& measured_fan_data); + void compute_initial_data_dependent_factors(); + + void efficiency_iteration(int iter_num, int eff_iter_num); + + void geo_normalization_iteration(int iter_num); + + void block_normalization_iteration(int iter_num); // Arguments std::string out_filename_prefix; const ProjData& measured_data; @@ -107,6 +81,22 @@ class MLEstimateComponentBasedNormalisation // Calculated values float threshold_for_KL; + + // Calculated data + FanProjData model_fan_data; + FanProjData fan_data; + DetectorEfficiencies data_fan_sums; + DetectorEfficiencies efficiencies; + BlockData3D norm_block_data; + BlockData3D measured_block_data; + GeoData3D norm_geo_data; + GeoData3D measured_geo_data; + + // do_KL specific varaibles + FanProjData measured_fan_data; + DetectorEfficiencies fan_sums; + GeoData3D geo_data; + BlockData3D block_data; }; END_NAMESPACE_STIR diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index 9069227f7..eb16a3557 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -80,7 +80,68 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(con do_symmetry_per_block(do_symmetry_per_block), do_KL(do_KL), do_display(do_display) -{} +{ + const int num_transaxial_blocks = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); + const int num_axial_blocks = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); + const int virtual_axial_crystals + = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_axial_crystals_per_block(); + const int virtual_transaxial_crystals + = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_transaxial_crystals_per_block(); + const int num_physical_rings = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_rings() + - (num_axial_blocks - 1) * virtual_axial_crystals; + const int num_physical_detectors_per_ring + = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_detectors_per_ring() + - num_transaxial_blocks * virtual_transaxial_crystals; + const int num_transaxial_buckets = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_buckets(); + const int num_axial_buckets = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_buckets(); + const int num_transaxial_blocks_per_bucket + = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks_per_bucket(); + const int num_axial_blocks_per_bucket + = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks_per_bucket(); + + int num_physical_transaxial_crystals_per_basic_unit + = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_crystals_per_block() + - virtual_transaxial_crystals; + int num_physical_axial_crystals_per_basic_unit + = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_crystals_per_block() - virtual_axial_crystals; + // If there are multiple buckets, we increase the symmetry size to a bucket. Otherwise, we use a block. + if (do_symmetry_per_block == false) + { + if (num_transaxial_buckets > 1) + { + num_physical_transaxial_crystals_per_basic_unit *= num_transaxial_blocks_per_bucket; + } + if (num_axial_buckets > 1) + { + num_physical_axial_crystals_per_basic_unit *= num_axial_blocks_per_bucket; + } + } + + data_fan_sums = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); + efficiencies = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); + + measured_geo_data = GeoData3D(num_physical_axial_crystals_per_basic_unit, + num_physical_transaxial_crystals_per_basic_unit / 2, + num_physical_rings, + num_physical_detectors_per_ring); // inputes have to be modified + norm_geo_data = GeoData3D(num_physical_axial_crystals_per_basic_unit, + num_physical_transaxial_crystals_per_basic_unit / 2, + num_physical_rings, + num_physical_detectors_per_ring); // inputes have to be modified + + measured_block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); + norm_block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); + + make_fan_data_remove_gaps(model_fan_data, model_data); + compute_initial_data_dependent_factors(); + + fan_sums = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); + geo_data = GeoData3D(num_physical_axial_crystals_per_basic_unit, + num_physical_transaxial_crystals_per_basic_unit / 2, + num_physical_rings, + num_physical_detectors_per_ring); // inputes have to be modified + block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); +} //! Helper function to write efficiencies to a file void @@ -116,11 +177,7 @@ MLEstimateComponentBasedNormalisation::write_block_data_to_file(int iter_num, co // Function to compute factors dependent on the data void -MLEstimateComponentBasedNormalisation::compute_initial_data_dependent_factors(const FanProjData& model_fan_data, - FanProjData& measured_fan_data, - DetectorEfficiencies& data_fan_sums, - GeoData3D& measured_geo_data, - BlockData3D& measured_block_data) +MLEstimateComponentBasedNormalisation::compute_initial_data_dependent_factors() { make_fan_data_remove_gaps(measured_fan_data, measured_data); @@ -150,15 +207,7 @@ MLEstimateComponentBasedNormalisation::compute_initial_data_dependent_factors(co // Function to handle efficiency iteration void -MLEstimateComponentBasedNormalisation::efficiency_iteration(FanProjData& fan_data, - const FanProjData& model_fan_data, - const GeoData3D& norm_geo_data, - const BlockData3D& norm_block_data, - DetectorEfficiencies& efficiencies, - const DetectorEfficiencies& data_fan_sums, - int iter_num, - int eff_iter_num, - const FanProjData& measured_fan_data) +MLEstimateComponentBasedNormalisation::efficiency_iteration(const int iter_num, const int eff_iter_num) { iterate_efficiencies(efficiencies, data_fan_sums, fan_data); write_efficiencies_to_file(iter_num, eff_iter_num, efficiencies); @@ -189,14 +238,7 @@ MLEstimateComponentBasedNormalisation::efficiency_iteration(FanProjData& fan_dat // Function to handle geo normalization void -MLEstimateComponentBasedNormalisation::geo_normalization_iteration(FanProjData& fan_data, - const FanProjData& model_fan_data, - const DetectorEfficiencies& efficiencies, - const BlockData3D& norm_block_data, - GeoData3D& norm_geo_data, - const GeoData3D& measured_geo_data, - int iter_num, - const FanProjData& measured_fan_data) +MLEstimateComponentBasedNormalisation::geo_normalization_iteration(int iter_num) { fan_data = model_fan_data; apply_efficiencies(fan_data, efficiencies); @@ -223,14 +265,7 @@ MLEstimateComponentBasedNormalisation::geo_normalization_iteration(FanProjData& // Function to handle block normalization void -MLEstimateComponentBasedNormalisation::block_normalization_iteration(FanProjData& fan_data, - const FanProjData& model_fan_data, - const DetectorEfficiencies& efficiencies, - const GeoData3D& norm_geo_data, - BlockData3D& norm_block_data, - const BlockData3D& measured_block_data, - int iter_num, - const FanProjData& measured_fan_data) +MLEstimateComponentBasedNormalisation::block_normalization_iteration(const int iter_num) { fan_data = model_fan_data; apply_efficiencies(fan_data, efficiencies); @@ -257,67 +292,6 @@ MLEstimateComponentBasedNormalisation::block_normalization_iteration(FanProjData void MLEstimateComponentBasedNormalisation::run() { - const int num_transaxial_blocks = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); - const int num_axial_blocks = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); - const int virtual_axial_crystals - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_axial_crystals_per_block(); - const int virtual_transaxial_crystals - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_transaxial_crystals_per_block(); - const int num_physical_rings = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_rings() - - (num_axial_blocks - 1) * virtual_axial_crystals; - const int num_physical_detectors_per_ring - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_detectors_per_ring() - - num_transaxial_blocks * virtual_transaxial_crystals; - const int num_transaxial_buckets = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_buckets(); - const int num_axial_buckets = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_buckets(); - const int num_transaxial_blocks_per_bucket - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks_per_bucket(); - const int num_axial_blocks_per_bucket - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks_per_bucket(); - - int num_physical_transaxial_crystals_per_basic_unit - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_crystals_per_block() - - virtual_transaxial_crystals; - int num_physical_axial_crystals_per_basic_unit - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_crystals_per_block() - virtual_axial_crystals; - // If there are multiple buckets, we increase the symmetry size to a bucket. Otherwise, we use a block. - if (do_symmetry_per_block == false) - { - if (num_transaxial_buckets > 1) - { - num_physical_transaxial_crystals_per_basic_unit *= num_transaxial_blocks_per_bucket; - } - if (num_axial_buckets > 1) - { - num_physical_axial_crystals_per_basic_unit *= num_axial_blocks_per_bucket; - } - } - - FanProjData model_fan_data; - FanProjData fan_data; - DetectorEfficiencies data_fan_sums(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); - DetectorEfficiencies efficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); - - GeoData3D measured_geo_data(num_physical_axial_crystals_per_basic_unit, - num_physical_transaxial_crystals_per_basic_unit / 2, - num_physical_rings, - num_physical_detectors_per_ring); // inputes have to be modified - GeoData3D norm_geo_data(num_physical_axial_crystals_per_basic_unit, - num_physical_transaxial_crystals_per_basic_unit / 2, - num_physical_rings, - num_physical_detectors_per_ring); // inputes have to be modified - - BlockData3D measured_block_data(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); - BlockData3D norm_block_data(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); - - make_fan_data_remove_gaps(model_fan_data, model_data); - - // next could be local if KL is not computed below - FanProjData measured_fan_data; - - compute_initial_data_dependent_factors( - model_fan_data, measured_fan_data, data_fan_sums, measured_geo_data, measured_block_data); - // std::cerr << "model min " << model_fan_data.find_min() << " ,max " << model_fan_data.find_max() << std::endl; if (do_display) display(model_fan_data, "model"); @@ -341,41 +315,18 @@ MLEstimateComponentBasedNormalisation::run() // Efficiency iterations for (int eff_iter_num = 1; eff_iter_num <= num_eff_iterations; ++eff_iter_num) { - efficiency_iteration(fan_data, - model_fan_data, - norm_geo_data, - norm_block_data, - efficiencies, - data_fan_sums, - iter_num, - eff_iter_num, - measured_fan_data); + efficiency_iteration(iter_num, eff_iter_num); } // geo norm - geo_normalization_iteration( - fan_data, model_fan_data, efficiencies, norm_block_data, norm_geo_data, measured_geo_data, iter_num, measured_fan_data); + geo_normalization_iteration(iter_num); // block norm - block_normalization_iteration(fan_data, - model_fan_data, - efficiencies, - norm_geo_data, - norm_block_data, - measured_block_data, - iter_num, - measured_fan_data); + block_normalization_iteration(iter_num); //// print KL for fansums if (do_KL) { - DetectorEfficiencies fan_sums(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); - GeoData3D geo_data(num_physical_axial_crystals_per_basic_unit, - num_physical_transaxial_crystals_per_basic_unit / 2, - num_physical_rings, - num_physical_detectors_per_ring); // inputes have to be modified - BlockData3D block_data(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); - make_fan_sum_data(fan_sums, fan_data); make_geo_data(geo_data, fan_data); make_block_data(block_data, measured_fan_data); From 9fed40f642e96558fab54900a94b085508b88085 Mon Sep 17 00:00:00 2001 From: rts Date: Thu, 29 Aug 2024 12:59:00 -0700 Subject: [PATCH 04/34] Add do_save_to_file guards --- ...L_estimate_component_based_normalisation.h | 7 +++-- ...estimate_component_based_normalisation.cxx | 29 +++++++++++++------ src/utilities/find_ML_normfactors3D.cxx | 3 +- 3 files changed, 27 insertions(+), 12 deletions(-) diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h index 1d1874a10..4a4f06977 100644 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h @@ -35,7 +35,8 @@ void ML_estimate_component_based_normalisation(const std::string& out_filename_p bool do_block, bool do_symmetry_per_block, bool do_KL, - bool do_display); + bool do_display, + bool do_save_to_file); class MLEstimateComponentBasedNormalisation { @@ -49,7 +50,8 @@ class MLEstimateComponentBasedNormalisation bool do_block, bool do_symmetry_per_block, bool do_KL, - bool do_display); + bool do_display, + bool do_save_to_file); void run(); @@ -78,6 +80,7 @@ class MLEstimateComponentBasedNormalisation bool do_symmetry_per_block; bool do_KL; bool do_display; + bool do_save_to_file; // Calculated values float threshold_for_KL; diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index eb16a3557..da29a627c 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -45,7 +45,8 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix bool do_block, bool do_symmetry_per_block, bool do_KL, - bool do_display) + bool do_display, + bool do_save_to_file) { MLEstimateComponentBasedNormalisation estimator(out_filename_prefix, measured_data, @@ -56,7 +57,8 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix do_block, do_symmetry_per_block, do_KL, - do_display); + do_display, + do_save_to_file); estimator.run(); } @@ -69,7 +71,8 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(con bool do_block, bool do_symmetry_per_block, bool do_KL, - bool do_display) + bool do_display, + bool do_save_to_file) : measured_data(measured_data), model_data(model_data), out_filename_prefix(out_filename_prefix), @@ -79,7 +82,8 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(con do_block(do_block), do_symmetry_per_block(do_symmetry_per_block), do_KL(do_KL), - do_display(do_display) + do_display(do_display), + do_save_to_file(do_save_to_file) { const int num_transaxial_blocks = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); const int num_axial_blocks = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); @@ -210,7 +214,10 @@ void MLEstimateComponentBasedNormalisation::efficiency_iteration(const int iter_num, const int eff_iter_num) { iterate_efficiencies(efficiencies, data_fan_sums, fan_data); - write_efficiencies_to_file(iter_num, eff_iter_num, efficiencies); + if (do_save_to_file) + { + write_efficiencies_to_file(iter_num, eff_iter_num, efficiencies); + } if (do_KL) { apply_efficiencies(fan_data, efficiencies); @@ -248,8 +255,10 @@ MLEstimateComponentBasedNormalisation::geo_normalization_iteration(int iter_num) { iterate_geo_norm(norm_geo_data, measured_geo_data, fan_data); } - - write_geo_data_to_file(iter_num, norm_geo_data); + if (do_save_to_file) + { + write_geo_data_to_file(iter_num, norm_geo_data); + } if (do_KL) { apply_geo_norm(fan_data, norm_geo_data); @@ -274,7 +283,10 @@ MLEstimateComponentBasedNormalisation::block_normalization_iteration(const int i { iterate_block_norm(norm_block_data, measured_block_data, fan_data); } - write_block_data_to_file(iter_num, norm_block_data); + if (do_save_to_file) + { + write_block_data_to_file(iter_num, norm_block_data); + } if (do_KL) { apply_block_norm(fan_data, norm_block_data); @@ -330,7 +342,6 @@ MLEstimateComponentBasedNormalisation::run() make_fan_sum_data(fan_sums, fan_data); make_geo_data(geo_data, fan_data); make_block_data(block_data, measured_fan_data); - info(boost::format("KL on fans: %1%, %2") % KL(measured_fan_data, fan_data, 0) % KL(measured_geo_data, geo_data, 0)); } } diff --git a/src/utilities/find_ML_normfactors3D.cxx b/src/utilities/find_ML_normfactors3D.cxx index 2e78c0cbf..0dc8abf52 100644 --- a/src/utilities/find_ML_normfactors3D.cxx +++ b/src/utilities/find_ML_normfactors3D.cxx @@ -114,7 +114,8 @@ main(int argc, char** argv) do_block, do_symmetry_per_block, do_KL, - do_display); + do_display, + /*do_save_to_file*/ true); timer.stop(); info(boost::format("CPU time %1% secs") % timer.value()); From 2579f81110533b5a2342e2a9554327d5763179d9 Mon Sep 17 00:00:00 2001 From: rts Date: Thu, 29 Aug 2024 13:36:24 -0700 Subject: [PATCH 05/34] Documentation and minor refactors --- ...L_estimate_component_based_normalisation.h | 112 ++++++-- ...estimate_component_based_normalisation.cxx | 262 +++++++++--------- 2 files changed, 229 insertions(+), 145 deletions(-) diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h index 4a4f06977..7e82bf798 100644 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h @@ -1,5 +1,6 @@ /* Copyright (C) 2022, University College London + Copyright (C) 2024 This file is part of STIR. SPDX-License-Identifier: Apache-2.0 @@ -20,12 +21,20 @@ START_NAMESPACE_STIR /*! - \brief Find normalisation factors using a maximum likelihood approach - + \brief Find normalisation factors using a maximum likelihood approach \ingroup recon_buildblock + \param out_filename_prefix The prefix for the output files + \param measured_data The measured projection data + \param model_data The model projection data + \param num_eff_iterations The number of (sub-)efficiency iterations to perform per iteration of the algorithm + \param num_iterations The number of algorithm iterations to perform + \param do_geo Whether to perform geo normalization calculations + \param do_block Whether to perform block normalization calculations + \param do_symmetry_per_block Whether to perform block normalization calculations with symmetry + \param do_KL Whether to perform Kullback-Leibler divergence calculations and display the KL value. This can be slow. + \param do_display Whether to display the progress of the algorithm. + \param do_save_to_file Whether to save the each iteration of the efficiencies, geo data and block data to file. */ - -// The function that calls the class MLEstimateComponentBasedNormalisation and runs the normalisation estimation void ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, const ProjData& measured_data, const ProjData& model_data, @@ -36,12 +45,30 @@ void ML_estimate_component_based_normalisation(const std::string& out_filename_p bool do_symmetry_per_block, bool do_KL, bool do_display, - bool do_save_to_file); + bool do_save_to_file = true); +/*! + \brief Find normalisation factors using a maximum likelihood approach + \ingroup recon_buildblock +*/ class MLEstimateComponentBasedNormalisation { public: - MLEstimateComponentBasedNormalisation(const std::string& out_filename_prefix, + /*! + \brief Constructor + \param out_filename_prefix The prefix for the output files + \param measured_data The measured projection data + \param model_data The model projection data + \param num_eff_iterations The number of (sub-)efficiency iterations to perform per iteration of the algorithm + \param num_iterations The number of algorithm iterations to perform + \param do_geo Whether to perform geo normalization calculations + \param do_block Whether to perform block normalization calculations + \param do_symmetry_per_block Whether to perform block normalization calculations with symmetry + \param do_KL Whether to perform Kullback-Leibler divergence calculations and display the KL value. This can be slow. + \param do_display Whether to display the progress of the algorithm. + \param do_save_to_file Whether to save the each iteration of the efficiencies, geo data and block data to file. + */ + MLEstimateComponentBasedNormalisation(std::string out_filename_prefix, const ProjData& measured_data, const ProjData& model_data, int num_eff_iterations, @@ -53,39 +80,86 @@ class MLEstimateComponentBasedNormalisation bool do_display, bool do_save_to_file); - void run(); + /*! + \brief Run the normalisation estimation algorithm using the parameters provided in the constructor. + */ + void process(); private: - void write_efficiencies_to_file(int iter_num, int eff_iter_num, const DetectorEfficiencies& efficiencies); - - void write_geo_data_to_file(int iter_num, const GeoData3D& norm_geo_data); - - void write_block_data_to_file(int iter_num, const BlockData3D& norm_block_data); - - void compute_initial_data_dependent_factors(); - + /*! + \brief Performs an efficiency iteration to update the efficiancies from the data_fan_sums and model. + Additionally, handles the saving of the efficiencies iteration to file, KL calculation and display. + \param[in] iter_num The iteration number + \param[in] eff_iter_num The efficiency iteration number + */ void efficiency_iteration(int iter_num, int eff_iter_num); + /*! + \brief Performs a geo normalization iteration to update the geo data from the measured_geo_data and model_data. + Additionally, handles the saving of the geo data iteration to file, KL calculation and display. + \param[in] iter_num The iteration number + */ void geo_normalization_iteration(int iter_num); + /*! + \brief Performs a block normalization iteration to update the block data from the measured_block_data and model_data. + Additionally, handles the saving of the block data iteration to file, KL calculation and display. + * @param iter_num The iteration number + */ void block_normalization_iteration(int iter_num); - // Arguments + + /*! + \brief Write the efficiencies to a file (regardless of the value of do_save_to_file) + \param[in] iter_num The iteration number + \param[in] eff_iter_num The efficiency iteration number + */ + void write_efficiencies_to_file(int iter_num, int eff_iter_num) const; + + /*! + \brief Write the efficiencies to a file (regardless of the value of do_save_to_file) + \param[in] iter_num The iteration number + */ + void write_geo_data_to_file(int iter_num) const; + + /*! + \brief Write the efficiencies to a file (regardless of the value of do_save_to_file) + \param[in] iter_num The iteration number + */ + void write_block_data_to_file(int iter_num) const; + + /*! + \brief Computes the threshold for the Kullback-Leibler divergence calculation. This is a purely heuristic value. + \return The threshold value + */ + float compute_threshold_for_KL(); + + // Constructor parameters + //! The prefix for the output files std::string out_filename_prefix; + //! The measured projection data const ProjData& measured_data; + //! The model projection data const ProjData& model_data; + //! The number of (sub-)efficiency iterations to perform per iteration of the algorithm int num_eff_iterations; + //! The number of algorithm iterations to perform int num_iterations; + //! Whether to perform geo normalization calculations bool do_geo; + //! Whether to perform block normalization calculations bool do_block; + //! Whether to perform block normalization calculations with symmetry bool do_symmetry_per_block; + //! Whether to perform Kullback-Leibler divergence calculations and display the KL value. This can be slow. bool do_KL; + //! Whether to display the progress of the algorithm. bool do_display; + //! Whether to save the each iteration of the efficiencies, geo data and block data to file. bool do_save_to_file; - // Calculated values + // Calculated variables + //! The threshold for the Kullback-Leibler divergence calculation float threshold_for_KL; - - // Calculated data FanProjData model_fan_data; FanProjData fan_data; DetectorEfficiencies data_fan_sums; diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index da29a627c..ff31a93e3 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -3,6 +3,7 @@ Copyright (C) 2019-2020, 2022, University College London Copyright (C) 2016-2017, PETsys Electronics Copyright (C) 2021, Gefei Chen + Copyright (C) 2024 This file is part of STIR. SPDX-License-Identifier: Apache-2.0 @@ -18,6 +19,7 @@ \author Kris Thielemans \author Tahereh Niknejad \author Gefei Chen + \author Robert Twyman Skelly */ #include "stir/recon_buildblock/ML_estimate_component_based_normalisation.h" @@ -32,6 +34,7 @@ #include #include #include +#include START_NAMESPACE_STIR @@ -39,14 +42,14 @@ void ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, const ProjData& measured_data, const ProjData& model_data, - int num_eff_iterations, - int num_iterations, - bool do_geo, - bool do_block, - bool do_symmetry_per_block, - bool do_KL, - bool do_display, - bool do_save_to_file) + const int num_eff_iterations, + const int num_iterations, + const bool do_geo, + const bool do_block, + const bool do_symmetry_per_block, + const bool do_KL, + const bool do_display, + const bool do_save_to_file) { MLEstimateComponentBasedNormalisation estimator(out_filename_prefix, measured_data, @@ -59,23 +62,23 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix do_KL, do_display, do_save_to_file); - estimator.run(); + estimator.process(); } -MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(const std::string& out_filename_prefix, +MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std::string out_filename_prefix, const ProjData& measured_data, const ProjData& model_data, - int num_eff_iterations, - int num_iterations, - bool do_geo, - bool do_block, - bool do_symmetry_per_block, - bool do_KL, - bool do_display, - bool do_save_to_file) - : measured_data(measured_data), + const int num_eff_iterations, + const int num_iterations, + const bool do_geo, + const bool do_block, + const bool do_symmetry_per_block, + const bool do_KL, + const bool do_display, + const bool do_save_to_file) + : out_filename_prefix(std::move(out_filename_prefix)), + measured_data(measured_data), model_data(model_data), - out_filename_prefix(out_filename_prefix), num_eff_iterations(num_eff_iterations), num_iterations(num_iterations), do_geo(do_geo), @@ -121,6 +124,7 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(con } } + // Setup the data structures given the PET scanner geometry data_fan_sums = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); efficiencies = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); @@ -137,8 +141,17 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(con norm_block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); make_fan_data_remove_gaps(model_fan_data, model_data); - compute_initial_data_dependent_factors(); + make_fan_data_remove_gaps(measured_fan_data, measured_data); + + threshold_for_KL = compute_threshold_for_KL(); + + make_fan_sum_data(data_fan_sums, measured_fan_data); + make_geo_data(measured_geo_data, measured_fan_data); + make_block_data(measured_block_data, measured_fan_data); + if (do_display) + display(measured_block_data, "raw block data from measurements"); + // Compute the do_KL specific varaibles from the measured data fan_sums = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); geo_data = GeoData3D(num_physical_axial_crystals_per_basic_unit, num_physical_transaxial_crystals_per_basic_unit / 2, @@ -147,76 +160,60 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(con block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); } -//! Helper function to write efficiencies to a file -void -MLEstimateComponentBasedNormalisation::write_efficiencies_to_file(int iter_num, - int eff_iter_num, - const DetectorEfficiencies& efficiencies) -{ - char* out_filename = new char[out_filename_prefix.size() + 30]; - sprintf(out_filename, "%s_%s_%d_%d.out", out_filename_prefix.c_str(), "eff", iter_num, eff_iter_num); - std::ofstream out(out_filename); - out << efficiencies; -} - -//! Helper function to write geo data to a file void -MLEstimateComponentBasedNormalisation::write_geo_data_to_file(int iter_num, const GeoData3D& norm_geo_data) +MLEstimateComponentBasedNormalisation::process() { - char* out_filename = new char[out_filename_prefix.size() + 30]; - sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "geo", iter_num); - std::ofstream out(out_filename); - out << norm_geo_data; -} + if (do_display) + display(model_fan_data, "model"); -//! Helper function to write block data to a file -void -MLEstimateComponentBasedNormalisation::write_block_data_to_file(int iter_num, const BlockData3D& norm_block_data) -{ - char* out_filename = new char[out_filename_prefix.size() + 30]; - sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "block", iter_num); - std::ofstream out(out_filename); - out << norm_block_data; -} + efficiencies.fill(sqrt(data_fan_sums.sum() / model_fan_data.sum())); + norm_geo_data.fill(1); + norm_block_data.fill(1); -// Function to compute factors dependent on the data -void -MLEstimateComponentBasedNormalisation::compute_initial_data_dependent_factors() -{ - make_fan_data_remove_gaps(measured_fan_data, measured_data); - - /* TEMP FIX */ - for (int ra = model_fan_data.get_min_ra(); ra <= model_fan_data.get_max_ra(); ++ra) + for (int iter_num = 1; iter_num <= std::max(num_iterations, 1); ++iter_num) { - for (int a = model_fan_data.get_min_a(); a <= model_fan_data.get_max_a(); ++a) + fan_data = model_fan_data; + apply_geo_norm(fan_data, norm_geo_data); + apply_block_norm(fan_data, norm_block_data); + if (do_display) { - for (int rb = std::max(ra, model_fan_data.get_min_rb(ra)); rb <= model_fan_data.get_max_rb(ra); ++rb) - { - for (int b = model_fan_data.get_min_b(a); b <= model_fan_data.get_max_b(a); ++b) - if (model_fan_data(ra, a, rb, b) == 0) - measured_fan_data(ra, a, rb, b) = 0; - } + display(fan_data, "model*geo*block"); } - } - threshold_for_KL = measured_fan_data.find_max() / 100000.F; - // display(measured_fan_data, "measured data"); + // Internal Efficiency iterations loop + for (int eff_iter_num = 1; eff_iter_num <= num_eff_iterations; ++eff_iter_num) + { + efficiency_iteration(iter_num, eff_iter_num); + } - make_fan_sum_data(data_fan_sums, measured_fan_data); - make_geo_data(measured_geo_data, measured_fan_data); - make_block_data(measured_block_data, measured_fan_data); - if (do_display) - display(measured_block_data, "raw block data from measurements"); + if (do_geo) + { + // geo norm iteration + geo_normalization_iteration(iter_num); // Calculate geo norm iteration + } + if (do_block) + { + // block norm iteration + block_normalization_iteration(iter_num); // Calculate block norm iteration + } + //// print KL for fansums + if (do_KL) + { + make_fan_sum_data(fan_sums, fan_data); + make_geo_data(geo_data, fan_data); + make_block_data(block_data, measured_fan_data); + info(boost::format("KL on fans: %1%, %2") % KL(measured_fan_data, fan_data, 0) % KL(measured_geo_data, geo_data, 0)); + } + } } -// Function to handle efficiency iteration void MLEstimateComponentBasedNormalisation::efficiency_iteration(const int iter_num, const int eff_iter_num) { iterate_efficiencies(efficiencies, data_fan_sums, fan_data); if (do_save_to_file) { - write_efficiencies_to_file(iter_num, eff_iter_num, efficiencies); + write_efficiencies_to_file(iter_num, eff_iter_num); } if (do_KL) { @@ -243,21 +240,18 @@ MLEstimateComponentBasedNormalisation::efficiency_iteration(const int iter_num, } } -// Function to handle geo normalization void MLEstimateComponentBasedNormalisation::geo_normalization_iteration(int iter_num) { - fan_data = model_fan_data; - apply_efficiencies(fan_data, efficiencies); - apply_block_norm(fan_data, norm_block_data); - if (do_geo) - { - iterate_geo_norm(norm_geo_data, measured_geo_data, fan_data); - } + fan_data = model_fan_data; // Reset fan_data to model_data + apply_efficiencies(fan_data, efficiencies); // Apply efficiencies + apply_block_norm(fan_data, norm_block_data); // Apply block norm + iterate_geo_norm(norm_geo_data, measured_geo_data, fan_data); + if (do_save_to_file) { - write_geo_data_to_file(iter_num, norm_geo_data); + write_geo_data_to_file(iter_num); } if (do_KL) { @@ -272,20 +266,18 @@ MLEstimateComponentBasedNormalisation::geo_normalization_iteration(int iter_num) } } -// Function to handle block normalization void MLEstimateComponentBasedNormalisation::block_normalization_iteration(const int iter_num) { - fan_data = model_fan_data; - apply_efficiencies(fan_data, efficiencies); - apply_geo_norm(fan_data, norm_geo_data); - if (do_block) - { - iterate_block_norm(norm_block_data, measured_block_data, fan_data); - } + + fan_data = model_fan_data; // Reset fan_data to model_data + apply_efficiencies(fan_data, efficiencies); // Apply efficiencies + apply_geo_norm(fan_data, norm_geo_data); // Apply geo norm + iterate_block_norm(norm_block_data, measured_block_data, fan_data); // Iterate block norm calculation + if (do_save_to_file) { - write_block_data_to_file(iter_num, norm_block_data); + write_block_data_to_file(iter_num); } if (do_KL) { @@ -302,49 +294,67 @@ MLEstimateComponentBasedNormalisation::block_normalization_iteration(const int i } void -MLEstimateComponentBasedNormalisation::run() +MLEstimateComponentBasedNormalisation::write_efficiencies_to_file(const int iter_num, const int eff_iter_num) const { - // std::cerr << "model min " << model_fan_data.find_min() << " ,max " << model_fan_data.find_max() << std::endl; - if (do_display) - display(model_fan_data, "model"); + char* out_filename = new char[out_filename_prefix.size() + 30]; + sprintf(out_filename, "%s_%s_%d_%d.out", out_filename_prefix.c_str(), "eff", iter_num, eff_iter_num); + std::ofstream out(out_filename); + out << efficiencies; +} - for (int iter_num = 1; iter_num <= std::max(num_iterations, 1); ++iter_num) - { - if (iter_num == 1) - { - efficiencies.fill(sqrt(data_fan_sums.sum() / model_fan_data.sum())); - norm_geo_data.fill(1); - norm_block_data.fill(1); - } - fan_data = model_fan_data; - apply_geo_norm(fan_data, norm_geo_data); - apply_block_norm(fan_data, norm_block_data); - if (do_display) - { - display(fan_data, "model*geo*block"); - } +void +MLEstimateComponentBasedNormalisation::write_geo_data_to_file(const int iter_num) const +{ + char* out_filename = new char[out_filename_prefix.size() + 30]; + sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "geo", iter_num); + std::ofstream out(out_filename); + out << norm_geo_data; +} - // Efficiency iterations - for (int eff_iter_num = 1; eff_iter_num <= num_eff_iterations; ++eff_iter_num) - { - efficiency_iteration(iter_num, eff_iter_num); - } +void +MLEstimateComponentBasedNormalisation::write_block_data_to_file(const int iter_num) const +{ + char* out_filename = new char[out_filename_prefix.size() + 30]; + sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "block", iter_num); + std::ofstream out(out_filename); + out << norm_block_data; +} - // geo norm - geo_normalization_iteration(iter_num); +float +MLEstimateComponentBasedNormalisation::compute_threshold_for_KL() +{ + // Set the max found value to -inf + float max_elem = -std::numeric_limits::infinity(); - // block norm - block_normalization_iteration(iter_num); + const auto min_ra = model_fan_data.get_min_ra(); + const auto max_ra = model_fan_data.get_max_ra(); + const auto min_a = model_fan_data.get_min_a(); + const auto max_a = model_fan_data.get_max_a(); - //// print KL for fansums - if (do_KL) + for (auto ra = min_ra; ra <= max_ra; ++ra) + { + const auto min_rb = std::max(ra, model_fan_data.get_min_rb(ra)); + const auto max_rb = model_fan_data.get_max_rb(ra); + + for (auto a = min_a; a <= max_a; ++a) { - make_fan_sum_data(fan_sums, fan_data); - make_geo_data(geo_data, fan_data); - make_block_data(block_data, measured_fan_data); - info(boost::format("KL on fans: %1%, %2") % KL(measured_fan_data, fan_data, 0) % KL(measured_geo_data, geo_data, 0)); + const auto min_b = model_fan_data.get_min_b(a); + const auto max_b = model_fan_data.get_max_b(a); + + for (auto rb = min_rb; rb <= max_rb; ++rb) + { + for (auto b = min_b; b <= max_b; ++b) + { + if (model_fan_data(ra, a, rb, b) == 0) + { + max_elem = std::max(max_elem, measured_fan_data(ra, a, rb, b)); + } + } + } } } + + return max_elem / 100000.F; } END_NAMESPACE_STIR From 5d89071b1d37da93338b5ab24d9f06ee50bfd121 Mon Sep 17 00:00:00 2001 From: rts Date: Thu, 29 Aug 2024 14:18:08 -0700 Subject: [PATCH 06/34] Rename refactor --- ...L_estimate_component_based_normalisation.h | 4 +- ...estimate_component_based_normalisation.cxx | 41 ++++++++++--------- 2 files changed, 24 insertions(+), 21 deletions(-) diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h index 7e82bf798..13b20560b 100644 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h @@ -137,9 +137,9 @@ class MLEstimateComponentBasedNormalisation //! The prefix for the output files std::string out_filename_prefix; //! The measured projection data - const ProjData& measured_data; + const ProjData& measured_projdata; //! The model projection data - const ProjData& model_data; + const ProjData& model_projdata; //! The number of (sub-)efficiency iterations to perform per iteration of the algorithm int num_eff_iterations; //! The number of algorithm iterations to perform diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index ff31a93e3..1db8fb94d 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -66,8 +66,8 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix } MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std::string out_filename_prefix, - const ProjData& measured_data, - const ProjData& model_data, + const ProjData& measured_projdata, + const ProjData& model_projdata, const int num_eff_iterations, const int num_iterations, const bool do_geo, @@ -77,8 +77,8 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std const bool do_display, const bool do_save_to_file) : out_filename_prefix(std::move(out_filename_prefix)), - measured_data(measured_data), - model_data(model_data), + measured_projdata(measured_projdata), + model_projdata(model_projdata), num_eff_iterations(num_eff_iterations), num_iterations(num_iterations), do_geo(do_geo), @@ -88,29 +88,29 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std do_display(do_display), do_save_to_file(do_save_to_file) { - const int num_transaxial_blocks = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); - const int num_axial_blocks = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); + const int num_transaxial_blocks = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); + const int num_axial_blocks = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); const int virtual_axial_crystals - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_axial_crystals_per_block(); + = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_axial_crystals_per_block(); const int virtual_transaxial_crystals - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_transaxial_crystals_per_block(); - const int num_physical_rings = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_rings() + = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_transaxial_crystals_per_block(); + const int num_physical_rings = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_rings() - (num_axial_blocks - 1) * virtual_axial_crystals; const int num_physical_detectors_per_ring - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_detectors_per_ring() + = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_detectors_per_ring() - num_transaxial_blocks * virtual_transaxial_crystals; - const int num_transaxial_buckets = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_buckets(); - const int num_axial_buckets = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_buckets(); + const int num_transaxial_buckets = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_buckets(); + const int num_axial_buckets = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_buckets(); const int num_transaxial_blocks_per_bucket - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks_per_bucket(); + = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks_per_bucket(); const int num_axial_blocks_per_bucket - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks_per_bucket(); + = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks_per_bucket(); int num_physical_transaxial_crystals_per_basic_unit - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_crystals_per_block() + = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_crystals_per_block() - virtual_transaxial_crystals; int num_physical_axial_crystals_per_basic_unit - = measured_data.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_crystals_per_block() - virtual_axial_crystals; + = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_crystals_per_block() - virtual_axial_crystals; // If there are multiple buckets, we increase the symmetry size to a bucket. Otherwise, we use a block. if (do_symmetry_per_block == false) { @@ -140,8 +140,8 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std measured_block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); norm_block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); - make_fan_data_remove_gaps(model_fan_data, model_data); - make_fan_data_remove_gaps(measured_fan_data, measured_data); + make_fan_data_remove_gaps(model_fan_data, model_projdata); + make_fan_data_remove_gaps(measured_fan_data, measured_projdata); threshold_for_KL = compute_threshold_for_KL(); @@ -164,8 +164,11 @@ void MLEstimateComponentBasedNormalisation::process() { if (do_display) - display(model_fan_data, "model"); + { + display(model_fan_data, "model"); + } + // Initialize the efficiencies, geo data and block data to 1 efficiencies.fill(sqrt(data_fan_sums.sum() / model_fan_data.sum())); norm_geo_data.fill(1); norm_block_data.fill(1); From 86f6a82183477b14219a445ebc168cb3919a1aed Mon Sep 17 00:00:00 2001 From: rts Date: Thu, 29 Aug 2024 15:12:41 -0700 Subject: [PATCH 07/34] Addition usage of shared_ptrs --- src/buildblock/ML_norm.cxx | 12 +- ...L_estimate_component_based_normalisation.h | 24 +-- ...estimate_component_based_normalisation.cxx | 147 +++++++++++------- src/utilities/find_ML_normfactors3D.cxx | 4 +- 4 files changed, 115 insertions(+), 72 deletions(-) diff --git a/src/buildblock/ML_norm.cxx b/src/buildblock/ML_norm.cxx index bc8c10910..986ed3b7a 100644 --- a/src/buildblock/ML_norm.cxx +++ b/src/buildblock/ML_norm.cxx @@ -987,6 +987,10 @@ get_fan_info(int& num_rings, int& num_detectors_per_ring, int& max_ring_diff, in { error("Can only process data without axial compression (i.e. span=1)\n"); } + if (proj_data_info_ptr->is_tof_data()) + { + error("make_fan_data: Incompatible with TOF data. Abort."); + } num_rings = proj_data_info.get_scanner_ptr()->get_num_rings(); num_detectors_per_ring = proj_data_info.get_scanner_ptr()->get_num_detectors_per_ring(); const int half_fan_size = min(proj_data_info.get_max_tangential_pos_num(), -proj_data_info.get_min_tangential_pos_num()); @@ -1141,23 +1145,19 @@ make_fan_data_remove_gaps(FanProjData& fan_data, const ProjData& proj_data) int num_detectors_per_ring; int fan_size; int max_delta; - - if (proj_data.get_proj_data_info_sptr()->is_tof_data()) - error("make_fan_data: Incompatible with TOF data. Abort."); - const ProjDataInfo& proj_data_info = *proj_data.get_proj_data_info_sptr(); get_fan_info(num_rings, num_detectors_per_ring, max_delta, fan_size, proj_data_info); if (proj_data.get_proj_data_info_sptr()->get_scanner_ptr()->get_scanner_geometry() == "Cylindrical") { - auto proj_data_info_ptr = dynamic_cast(&proj_data_info); + const auto proj_data_info_ptr = dynamic_cast(&proj_data_info); make_fan_data_remove_gaps_help( fan_data, num_rings, num_detectors_per_ring, max_delta, fan_size, *proj_data_info_ptr, proj_data); } else { - auto proj_data_info_ptr = dynamic_cast(&proj_data_info); + const auto proj_data_info_ptr = dynamic_cast(&proj_data_info); make_fan_data_remove_gaps_help( fan_data, num_rings, num_detectors_per_ring, max_delta, fan_size, *proj_data_info_ptr, proj_data); diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h index 13b20560b..b6af36f73 100644 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h @@ -36,8 +36,8 @@ START_NAMESPACE_STIR \param do_save_to_file Whether to save the each iteration of the efficiencies, geo data and block data to file. */ void ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, - const ProjData& measured_data, - const ProjData& model_data, + const std::shared_ptr& measured_data, + const std::shared_ptr& model_data, int num_eff_iterations, int num_iterations, bool do_geo, @@ -69,8 +69,8 @@ class MLEstimateComponentBasedNormalisation \param do_save_to_file Whether to save the each iteration of the efficiencies, geo data and block data to file. */ MLEstimateComponentBasedNormalisation(std::string out_filename_prefix, - const ProjData& measured_data, - const ProjData& model_data, + const std::shared_ptr& measured_data, + const std::shared_ptr& model_data, int num_eff_iterations, int num_iterations, bool do_geo, @@ -85,6 +85,10 @@ class MLEstimateComponentBasedNormalisation */ void process(); + std::shared_ptr get_efficiencies() const; + std::shared_ptr get_geo_data() const; + std::shared_ptr get_block_data() const; + private: /*! \brief Performs an efficiency iteration to update the efficiancies from the data_fan_sums and model. @@ -137,9 +141,9 @@ class MLEstimateComponentBasedNormalisation //! The prefix for the output files std::string out_filename_prefix; //! The measured projection data - const ProjData& measured_projdata; + const std::shared_ptr& measured_projdata; //! The model projection data - const ProjData& model_projdata; + const std::shared_ptr& model_projdata; //! The number of (sub-)efficiency iterations to perform per iteration of the algorithm int num_eff_iterations; //! The number of algorithm iterations to perform @@ -158,17 +162,19 @@ class MLEstimateComponentBasedNormalisation bool do_save_to_file; // Calculated variables + std::shared_ptr efficiencies_ptr; + std::shared_ptr norm_block_data_ptr; + std::shared_ptr norm_geo_data_ptr; //! The threshold for the Kullback-Leibler divergence calculation float threshold_for_KL; FanProjData model_fan_data; FanProjData fan_data; DetectorEfficiencies data_fan_sums; - DetectorEfficiencies efficiencies; - BlockData3D norm_block_data; BlockData3D measured_block_data; - GeoData3D norm_geo_data; GeoData3D measured_geo_data; + bool data_processed = false; + // do_KL specific varaibles FanProjData measured_fan_data; DetectorEfficiencies fan_sums; diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index 1db8fb94d..6aa4c73dd 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -40,8 +40,8 @@ START_NAMESPACE_STIR void ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, - const ProjData& measured_data, - const ProjData& model_data, + const std::shared_ptr& measured_data, + const std::shared_ptr& model_data, const int num_eff_iterations, const int num_iterations, const bool do_geo, @@ -66,8 +66,8 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix } MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std::string out_filename_prefix, - const ProjData& measured_projdata, - const ProjData& model_projdata, + const std::shared_ptr& measured_projdata_sptr, + const std::shared_ptr& model_projdata_sptr, const int num_eff_iterations, const int num_iterations, const bool do_geo, @@ -77,8 +77,8 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std const bool do_display, const bool do_save_to_file) : out_filename_prefix(std::move(out_filename_prefix)), - measured_projdata(measured_projdata), - model_projdata(model_projdata), + measured_projdata(measured_projdata_sptr), + model_projdata(model_projdata_sptr), num_eff_iterations(num_eff_iterations), num_iterations(num_iterations), do_geo(do_geo), @@ -88,29 +88,32 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std do_display(do_display), do_save_to_file(do_save_to_file) { - const int num_transaxial_blocks = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); - const int num_axial_blocks = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); + const int num_transaxial_blocks = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); + const int num_axial_blocks = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); const int virtual_axial_crystals - = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_axial_crystals_per_block(); + = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_axial_crystals_per_block(); const int virtual_transaxial_crystals - = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_transaxial_crystals_per_block(); - const int num_physical_rings = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_rings() + = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_transaxial_crystals_per_block(); + const int num_physical_rings = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_rings() - (num_axial_blocks - 1) * virtual_axial_crystals; const int num_physical_detectors_per_ring - = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_detectors_per_ring() + = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_detectors_per_ring() - num_transaxial_blocks * virtual_transaxial_crystals; - const int num_transaxial_buckets = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_buckets(); - const int num_axial_buckets = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_buckets(); + const int num_transaxial_buckets + = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_buckets(); + const int num_axial_buckets = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_buckets(); const int num_transaxial_blocks_per_bucket - = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks_per_bucket(); + = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks_per_bucket(); const int num_axial_blocks_per_bucket - = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks_per_bucket(); + = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks_per_bucket(); int num_physical_transaxial_crystals_per_basic_unit - = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_crystals_per_block() + = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_crystals_per_block() - virtual_transaxial_crystals; int num_physical_axial_crystals_per_basic_unit - = measured_projdata.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_crystals_per_block() - virtual_axial_crystals; + = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_crystals_per_block() + - virtual_axial_crystals; + // If there are multiple buckets, we increase the symmetry size to a bucket. Otherwise, we use a block. if (do_symmetry_per_block == false) { @@ -126,22 +129,23 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std // Setup the data structures given the PET scanner geometry data_fan_sums = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); - efficiencies = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); + efficiencies_ptr = std::make_shared(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); measured_geo_data = GeoData3D(num_physical_axial_crystals_per_basic_unit, num_physical_transaxial_crystals_per_basic_unit / 2, num_physical_rings, num_physical_detectors_per_ring); // inputes have to be modified - norm_geo_data = GeoData3D(num_physical_axial_crystals_per_basic_unit, - num_physical_transaxial_crystals_per_basic_unit / 2, - num_physical_rings, - num_physical_detectors_per_ring); // inputes have to be modified + norm_geo_data_ptr = std::make_shared(num_physical_axial_crystals_per_basic_unit, + num_physical_transaxial_crystals_per_basic_unit / 2, + num_physical_rings, + num_physical_detectors_per_ring); // inputes have to be modified measured_block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); - norm_block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); + norm_block_data_ptr + = std::make_shared(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); - make_fan_data_remove_gaps(model_fan_data, model_projdata); - make_fan_data_remove_gaps(measured_fan_data, measured_projdata); + make_fan_data_remove_gaps(model_fan_data, *model_projdata_sptr); + make_fan_data_remove_gaps(measured_fan_data, *measured_projdata_sptr); threshold_for_KL = compute_threshold_for_KL(); @@ -149,7 +153,9 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std make_geo_data(measured_geo_data, measured_fan_data); make_block_data(measured_block_data, measured_fan_data); if (do_display) - display(measured_block_data, "raw block data from measurements"); + { + display(measured_block_data, "raw block data from measurements"); + } // Compute the do_KL specific varaibles from the measured data fan_sums = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); @@ -163,21 +169,22 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std void MLEstimateComponentBasedNormalisation::process() { + data_processed = true; if (do_display) { display(model_fan_data, "model"); } // Initialize the efficiencies, geo data and block data to 1 - efficiencies.fill(sqrt(data_fan_sums.sum() / model_fan_data.sum())); - norm_geo_data.fill(1); - norm_block_data.fill(1); + efficiencies_ptr->fill(sqrt(data_fan_sums.sum() / model_fan_data.sum())); + norm_geo_data_ptr->fill(1); + norm_block_data_ptr->fill(1); for (int iter_num = 1; iter_num <= std::max(num_iterations, 1); ++iter_num) { fan_data = model_fan_data; - apply_geo_norm(fan_data, norm_geo_data); - apply_block_norm(fan_data, norm_block_data); + apply_geo_norm(fan_data, *norm_geo_data_ptr); + apply_block_norm(fan_data, *norm_block_data_ptr); if (do_display) { display(fan_data, "model*geo*block"); @@ -210,17 +217,47 @@ MLEstimateComponentBasedNormalisation::process() } } +std::shared_ptr +MLEstimateComponentBasedNormalisation::get_efficiencies() const +{ + if (!data_processed) + { + return nullptr; + } + return efficiencies_ptr; +} + +std::shared_ptr +MLEstimateComponentBasedNormalisation::get_geo_data() const +{ + if (!data_processed) + { + return nullptr; + } + return norm_geo_data_ptr; +} + +std::shared_ptr +MLEstimateComponentBasedNormalisation::get_block_data() const +{ + if (!data_processed) + { + return nullptr; + } + return norm_block_data_ptr; +} + void MLEstimateComponentBasedNormalisation::efficiency_iteration(const int iter_num, const int eff_iter_num) { - iterate_efficiencies(efficiencies, data_fan_sums, fan_data); + iterate_efficiencies(*efficiencies_ptr, data_fan_sums, fan_data); if (do_save_to_file) { write_efficiencies_to_file(iter_num, eff_iter_num); } if (do_KL) { - apply_efficiencies(fan_data, efficiencies); + apply_efficiencies(fan_data, *efficiencies_ptr); std::cerr << "measured*norm min " << measured_fan_data.find_min() << " ,max " << measured_fan_data.find_max() << std::endl; std::cerr << "model*norm min " << fan_data.find_min() << " ,max " << fan_data.find_max() << std::endl; if (do_display) @@ -228,18 +265,18 @@ MLEstimateComponentBasedNormalisation::efficiency_iteration(const int iter_num, info(boost::format("KL %1%") % KL(measured_fan_data, fan_data, threshold_for_KL)); // now restore for further iterations fan_data = model_fan_data; - apply_geo_norm(fan_data, norm_geo_data); - apply_block_norm(fan_data, norm_block_data); + apply_geo_norm(fan_data, *norm_geo_data_ptr); + apply_block_norm(fan_data, *norm_block_data_ptr); } if (do_display) { fan_data.fill(1); - apply_efficiencies(fan_data, efficiencies); + apply_efficiencies(fan_data, *efficiencies_ptr); display(fan_data, "eff norm"); // now restore for further iterations fan_data = model_fan_data; - apply_geo_norm(fan_data, norm_geo_data); - apply_block_norm(fan_data, norm_block_data); + apply_geo_norm(fan_data, *norm_geo_data_ptr); + apply_block_norm(fan_data, *norm_block_data_ptr); } } @@ -247,10 +284,10 @@ void MLEstimateComponentBasedNormalisation::geo_normalization_iteration(int iter_num) { - fan_data = model_fan_data; // Reset fan_data to model_data - apply_efficiencies(fan_data, efficiencies); // Apply efficiencies - apply_block_norm(fan_data, norm_block_data); // Apply block norm - iterate_geo_norm(norm_geo_data, measured_geo_data, fan_data); + fan_data = model_fan_data; // Reset fan_data to model_data + apply_efficiencies(fan_data, *efficiencies_ptr); // Apply efficiencies + apply_block_norm(fan_data, *norm_block_data_ptr); // Apply block norm + iterate_geo_norm(*norm_geo_data_ptr, measured_geo_data, fan_data); if (do_save_to_file) { @@ -258,13 +295,13 @@ MLEstimateComponentBasedNormalisation::geo_normalization_iteration(int iter_num) } if (do_KL) { - apply_geo_norm(fan_data, norm_geo_data); + apply_geo_norm(fan_data, *norm_geo_data_ptr); info(boost::format("KL %1%") % KL(measured_fan_data, fan_data, threshold_for_KL)); } if (do_display) { fan_data.fill(1); - apply_geo_norm(fan_data, norm_geo_data); + apply_geo_norm(fan_data, *norm_geo_data_ptr); display(fan_data, "geo norm"); } } @@ -273,10 +310,10 @@ void MLEstimateComponentBasedNormalisation::block_normalization_iteration(const int iter_num) { - fan_data = model_fan_data; // Reset fan_data to model_data - apply_efficiencies(fan_data, efficiencies); // Apply efficiencies - apply_geo_norm(fan_data, norm_geo_data); // Apply geo norm - iterate_block_norm(norm_block_data, measured_block_data, fan_data); // Iterate block norm calculation + fan_data = model_fan_data; // Reset fan_data to model_data + apply_efficiencies(fan_data, *efficiencies_ptr); // Apply efficiencies + apply_geo_norm(fan_data, *norm_geo_data_ptr); // Apply geo norm + iterate_block_norm(*norm_block_data_ptr, measured_block_data, fan_data); // Iterate block norm calculation if (do_save_to_file) { @@ -284,14 +321,14 @@ MLEstimateComponentBasedNormalisation::block_normalization_iteration(const int i } if (do_KL) { - apply_block_norm(fan_data, norm_block_data); + apply_block_norm(fan_data, *norm_block_data_ptr); info(boost::format("KL %1%") % KL(measured_fan_data, fan_data, threshold_for_KL)); } if (do_display) { fan_data.fill(1); - apply_block_norm(fan_data, norm_block_data); - display(norm_block_data, "raw block norm"); + apply_block_norm(fan_data, *norm_block_data_ptr); + display(*norm_block_data_ptr, "raw block norm"); display(fan_data, "block norm"); } } @@ -302,7 +339,7 @@ MLEstimateComponentBasedNormalisation::write_efficiencies_to_file(const int iter char* out_filename = new char[out_filename_prefix.size() + 30]; sprintf(out_filename, "%s_%s_%d_%d.out", out_filename_prefix.c_str(), "eff", iter_num, eff_iter_num); std::ofstream out(out_filename); - out << efficiencies; + out << *efficiencies_ptr; } void @@ -311,7 +348,7 @@ MLEstimateComponentBasedNormalisation::write_geo_data_to_file(const int iter_num char* out_filename = new char[out_filename_prefix.size() + 30]; sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "geo", iter_num); std::ofstream out(out_filename); - out << norm_geo_data; + out << *norm_geo_data_ptr; } void @@ -320,7 +357,7 @@ MLEstimateComponentBasedNormalisation::write_block_data_to_file(const int iter_n char* out_filename = new char[out_filename_prefix.size() + 30]; sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "block", iter_num); std::ofstream out(out_filename); - out << norm_block_data; + out << *norm_block_data_ptr; } float diff --git a/src/utilities/find_ML_normfactors3D.cxx b/src/utilities/find_ML_normfactors3D.cxx index 0dc8abf52..8293bb245 100644 --- a/src/utilities/find_ML_normfactors3D.cxx +++ b/src/utilities/find_ML_normfactors3D.cxx @@ -106,8 +106,8 @@ main(int argc, char** argv) timer.start(); ML_estimate_component_based_normalisation(out_filename_prefix, - *measured_data, - *model_data, + measured_data, + model_data, num_eff_iterations, num_iterations, do_geo, From 880e5303684c7cb0fe52c4216527799bc3912c7b Mon Sep 17 00:00:00 2001 From: rts Date: Thu, 29 Aug 2024 15:50:21 -0700 Subject: [PATCH 08/34] [ci skip] parameter rename, remove redundant logic, documentation --- ...L_estimate_component_based_normalisation.h | 36 ++++---- ...estimate_component_based_normalisation.cxx | 91 +++++++++---------- 2 files changed, 61 insertions(+), 66 deletions(-) diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h index b6af36f73..b7e469f29 100644 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h @@ -12,6 +12,7 @@ \ingroup recon_buildblock \brief Declaration of stir::ML_estimate_component_based_normalisation \author Kris Thielemans + \author Robert Twyman Skelly */ #include "stir/common.h" #include @@ -59,34 +60,37 @@ class MLEstimateComponentBasedNormalisation \param out_filename_prefix The prefix for the output files \param measured_data The measured projection data \param model_data The model projection data - \param num_eff_iterations The number of (sub-)efficiency iterations to perform per iteration of the algorithm - \param num_iterations The number of algorithm iterations to perform - \param do_geo Whether to perform geo normalization calculations - \param do_block Whether to perform block normalization calculations - \param do_symmetry_per_block Whether to perform block normalization calculations with symmetry - \param do_KL Whether to perform Kullback-Leibler divergence calculations and display the KL value. This can be slow. - \param do_display Whether to display the progress of the algorithm. - \param do_save_to_file Whether to save the each iteration of the efficiencies, geo data and block data to file. + \param num_eff_iterations_v The number of (sub-)efficiency iterations to perform per iteration of the algorithm + \param num_iterations_v The number of algorithm iterations to perform + \param do_geo_v Whether to perform geo normalization calculations + \param do_block_v Whether to perform block normalization calculations + \param do_symmetry_per_block_v Whether to perform block normalization calculations with symmetry + \param do_KL_v Whether to perform Kullback-Leibler divergence calculations and display the KL value. This can be slow. + \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, const std::shared_ptr& measured_data, const std::shared_ptr& model_data, - int num_eff_iterations, - int num_iterations, - bool do_geo, - bool do_block, - bool do_symmetry_per_block, - bool do_KL, - bool do_display, - bool do_save_to_file); + int num_eff_iterations_v, + int num_iterations_v, + bool do_geo_v, + bool do_block_v, + bool do_symmetry_per_block_v, + bool do_KL_v, + bool do_display_v, + bool do_save_to_file_v); /*! \brief Run the normalisation estimation algorithm using the parameters provided in the constructor. */ void process(); + //! Get the efficiencies, nullptr if not calculated std::shared_ptr get_efficiencies() const; + //! Get the geo data, nullptr if not calculated std::shared_ptr get_geo_data() const; + //! Get the block data, nullptr if not calculated std::shared_ptr get_block_data() const; private: diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index 6aa4c73dd..37c49a47b 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -22,18 +22,13 @@ \author Robert Twyman Skelly */ #include "stir/recon_buildblock/ML_estimate_component_based_normalisation.h" - #include "stir/ML_norm.h" #include "stir/Scanner.h" -#include "stir/stream.h" #include "stir/display.h" #include "stir/info.h" -#include "stir/warning.h" #include "stir/ProjData.h" -#include -#include + #include -#include #include START_NAMESPACE_STIR @@ -65,66 +60,62 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix estimator.process(); } -MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std::string out_filename_prefix, - const std::shared_ptr& measured_projdata_sptr, - const std::shared_ptr& model_projdata_sptr, - const int num_eff_iterations, - const int num_iterations, - const bool do_geo, - const bool do_block, - const bool do_symmetry_per_block, - const bool do_KL, - const bool do_display, - const bool do_save_to_file) - : out_filename_prefix(std::move(out_filename_prefix)), - measured_projdata(measured_projdata_sptr), - model_projdata(model_projdata_sptr), - num_eff_iterations(num_eff_iterations), - num_iterations(num_iterations), - do_geo(do_geo), - do_block(do_block), - do_symmetry_per_block(do_symmetry_per_block), - do_KL(do_KL), - do_display(do_display), - do_save_to_file(do_save_to_file) +MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation( + std::string out_filename_prefix_v, + const std::shared_ptr& measured_projdata_sptr_v, + const std::shared_ptr& model_projdata_sptr_v, + const int num_eff_iterations_v, + const int num_iterations_v, + const bool do_geo_v, + const bool do_block_v, + const bool do_symmetry_per_block_v, + const bool do_KL_v, + const bool do_display_v, + const bool do_save_to_file_v) + : out_filename_prefix(std::move(out_filename_prefix_v)), + measured_projdata(measured_projdata_sptr_v), + model_projdata(model_projdata_sptr_v), + num_eff_iterations(num_eff_iterations_v), + num_iterations(num_iterations_v), + do_geo(do_geo_v), + do_block(do_block_v), + do_symmetry_per_block(do_symmetry_per_block_v), + do_KL(do_KL_v), + do_display(do_display_v), + do_save_to_file(do_save_to_file_v) { - const int num_transaxial_blocks = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); - const int num_axial_blocks = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); + const int num_transaxial_blocks + = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); + const int num_axial_blocks = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); const int virtual_axial_crystals - = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_axial_crystals_per_block(); + = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_axial_crystals_per_block(); const int virtual_transaxial_crystals - = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_transaxial_crystals_per_block(); - const int num_physical_rings = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_rings() + = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_transaxial_crystals_per_block(); + const int num_physical_rings = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_rings() - (num_axial_blocks - 1) * virtual_axial_crystals; const int num_physical_detectors_per_ring - = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_detectors_per_ring() + = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_detectors_per_ring() - num_transaxial_blocks * virtual_transaxial_crystals; const int num_transaxial_buckets - = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_buckets(); - const int num_axial_buckets = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_buckets(); + = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_buckets(); + const int num_axial_buckets = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_buckets(); const int num_transaxial_blocks_per_bucket - = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks_per_bucket(); + = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks_per_bucket(); const int num_axial_blocks_per_bucket - = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks_per_bucket(); + = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks_per_bucket(); int num_physical_transaxial_crystals_per_basic_unit - = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_crystals_per_block() + = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_crystals_per_block() - virtual_transaxial_crystals; int num_physical_axial_crystals_per_basic_unit - = measured_projdata_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_crystals_per_block() + = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_crystals_per_block() - virtual_axial_crystals; // If there are multiple buckets, we increase the symmetry size to a bucket. Otherwise, we use a block. if (do_symmetry_per_block == false) { - if (num_transaxial_buckets > 1) - { - num_physical_transaxial_crystals_per_basic_unit *= num_transaxial_blocks_per_bucket; - } - if (num_axial_buckets > 1) - { - num_physical_axial_crystals_per_basic_unit *= num_axial_blocks_per_bucket; - } + num_physical_transaxial_crystals_per_basic_unit *= num_transaxial_blocks_per_bucket; + num_physical_axial_crystals_per_basic_unit *= num_axial_blocks_per_bucket; } // Setup the data structures given the PET scanner geometry @@ -144,8 +135,8 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation(std norm_block_data_ptr = std::make_shared(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); - make_fan_data_remove_gaps(model_fan_data, *model_projdata_sptr); - make_fan_data_remove_gaps(measured_fan_data, *measured_projdata_sptr); + make_fan_data_remove_gaps(model_fan_data, *model_projdata_sptr_v); + make_fan_data_remove_gaps(measured_fan_data, *measured_projdata_sptr_v); threshold_for_KL = compute_threshold_for_KL(); From dd711b1f1dec2f284d0a8d355995b04c5f4de06b Mon Sep 17 00:00:00 2001 From: rts Date: Fri, 30 Aug 2024 11:27:19 -0700 Subject: [PATCH 09/34] Remove model_data member and pass by reference into the constructor --- ...L_estimate_component_based_normalisation.h | 12 ++---- ...estimate_component_based_normalisation.cxx | 42 +++++++++---------- src/utilities/find_ML_normfactors3D.cxx | 4 +- 3 files changed, 26 insertions(+), 32 deletions(-) diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h index b7e469f29..8c36fd3a3 100644 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h @@ -37,8 +37,8 @@ START_NAMESPACE_STIR \param do_save_to_file Whether to save the each iteration of the efficiencies, geo data and block data to file. */ void ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, - const std::shared_ptr& measured_data, - const std::shared_ptr& model_data, + const ProjData& measured_projdata, + const ProjData& model_projdata, int num_eff_iterations, int num_iterations, bool do_geo, @@ -70,8 +70,8 @@ class MLEstimateComponentBasedNormalisation \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, - const std::shared_ptr& measured_data, - const std::shared_ptr& model_data, + const ProjData& measured_data_v, + const ProjData& model_data_v, int num_eff_iterations_v, int num_iterations_v, bool do_geo_v, @@ -144,10 +144,6 @@ class MLEstimateComponentBasedNormalisation // Constructor parameters //! The prefix for the output files std::string out_filename_prefix; - //! The measured projection data - const std::shared_ptr& measured_projdata; - //! The model projection data - const std::shared_ptr& model_projdata; //! The number of (sub-)efficiency iterations to perform per iteration of the algorithm int num_eff_iterations; //! The number of algorithm iterations to perform diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index 37c49a47b..dd7b5a265 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -35,8 +35,8 @@ START_NAMESPACE_STIR void ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, - const std::shared_ptr& measured_data, - const std::shared_ptr& model_data, + const ProjData& measured_projdata, + const ProjData& model_projdata, const int num_eff_iterations, const int num_iterations, const bool do_geo, @@ -47,8 +47,8 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix const bool do_save_to_file) { MLEstimateComponentBasedNormalisation estimator(out_filename_prefix, - measured_data, - model_data, + measured_projdata, + model_projdata, num_eff_iterations, num_iterations, do_geo, @@ -62,8 +62,8 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation( std::string out_filename_prefix_v, - const std::shared_ptr& measured_projdata_sptr_v, - const std::shared_ptr& model_projdata_sptr_v, + const ProjData& measured_projdata_v, + const ProjData& model_projdata_v, const int num_eff_iterations_v, const int num_iterations_v, const bool do_geo_v, @@ -73,8 +73,6 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation( const bool do_display_v, const bool do_save_to_file_v) : out_filename_prefix(std::move(out_filename_prefix_v)), - measured_projdata(measured_projdata_sptr_v), - model_projdata(model_projdata_sptr_v), num_eff_iterations(num_eff_iterations_v), num_iterations(num_iterations_v), do_geo(do_geo_v), @@ -85,30 +83,30 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation( do_save_to_file(do_save_to_file_v) { const int num_transaxial_blocks - = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); - const int num_axial_blocks = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); + = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); + const int num_axial_blocks = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); const int virtual_axial_crystals - = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_axial_crystals_per_block(); + = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_axial_crystals_per_block(); const int virtual_transaxial_crystals - = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_transaxial_crystals_per_block(); - const int num_physical_rings = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_rings() + = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_transaxial_crystals_per_block(); + const int num_physical_rings = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_rings() - (num_axial_blocks - 1) * virtual_axial_crystals; const int num_physical_detectors_per_ring - = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_detectors_per_ring() + = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_detectors_per_ring() - num_transaxial_blocks * virtual_transaxial_crystals; const int num_transaxial_buckets - = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_buckets(); - const int num_axial_buckets = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_buckets(); + = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_buckets(); + const int num_axial_buckets = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_buckets(); const int num_transaxial_blocks_per_bucket - = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks_per_bucket(); + = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks_per_bucket(); const int num_axial_blocks_per_bucket - = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks_per_bucket(); + = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks_per_bucket(); int num_physical_transaxial_crystals_per_basic_unit - = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_crystals_per_block() + = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_crystals_per_block() - virtual_transaxial_crystals; int num_physical_axial_crystals_per_basic_unit - = measured_projdata_sptr_v->get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_crystals_per_block() + = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_crystals_per_block() - virtual_axial_crystals; // If there are multiple buckets, we increase the symmetry size to a bucket. Otherwise, we use a block. @@ -135,8 +133,8 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation( norm_block_data_ptr = std::make_shared(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); - make_fan_data_remove_gaps(model_fan_data, *model_projdata_sptr_v); - make_fan_data_remove_gaps(measured_fan_data, *measured_projdata_sptr_v); + make_fan_data_remove_gaps(model_fan_data, model_projdata_v); + make_fan_data_remove_gaps(measured_fan_data, measured_projdata_v); threshold_for_KL = compute_threshold_for_KL(); diff --git a/src/utilities/find_ML_normfactors3D.cxx b/src/utilities/find_ML_normfactors3D.cxx index 8293bb245..0dc8abf52 100644 --- a/src/utilities/find_ML_normfactors3D.cxx +++ b/src/utilities/find_ML_normfactors3D.cxx @@ -106,8 +106,8 @@ main(int argc, char** argv) timer.start(); ML_estimate_component_based_normalisation(out_filename_prefix, - measured_data, - model_data, + *measured_data, + *model_data, num_eff_iterations, num_iterations, do_geo, From 2bcdc80e86d2ee9f13ccc76b87bc1cf5321928a8 Mon Sep 17 00:00:00 2001 From: rts Date: Fri, 30 Aug 2024 11:32:34 -0700 Subject: [PATCH 10/34] Readd stir/steam.h to includes --- .../ML_estimate_component_based_normalisation.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index dd7b5a265..02bbda81a 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -27,7 +27,7 @@ #include "stir/display.h" #include "stir/info.h" #include "stir/ProjData.h" - +#include "stir/stream.h" #include #include From 0fa8455f44e46afa56b2165a9c090ca562ddfbce Mon Sep 17 00:00:00 2001 From: rts Date: Fri, 30 Aug 2024 11:32:47 -0700 Subject: [PATCH 11/34] [ci skip] Update copywrite --- .../ML_estimate_component_based_normalisation.h | 2 +- .../ML_estimate_component_based_normalisation.cxx | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h index 8c36fd3a3..79bf94065 100644 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h @@ -1,6 +1,6 @@ /* Copyright (C) 2022, University College London - Copyright (C) 2024 + Copyright (C) 2024, Robert Twyman Skelly This file is part of STIR. SPDX-License-Identifier: Apache-2.0 diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index 02bbda81a..a1e8ac8c0 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -3,7 +3,7 @@ Copyright (C) 2019-2020, 2022, University College London Copyright (C) 2016-2017, PETsys Electronics Copyright (C) 2021, Gefei Chen - Copyright (C) 2024 + Copyright (C) 2024, Robert Twyman Skelly This file is part of STIR. SPDX-License-Identifier: Apache-2.0 From 90ea81bb9e89219e03c5f3757e9bd0236f0a6da9 Mon Sep 17 00:00:00 2001 From: rts Date: Fri, 30 Aug 2024 13:03:44 -0700 Subject: [PATCH 12/34] [ci skip] Fix a bug with do_KL info print --- .../ML_estimate_component_based_normalisation.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index a1e8ac8c0..bbb4736bc 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -201,7 +201,7 @@ MLEstimateComponentBasedNormalisation::process() make_fan_sum_data(fan_sums, fan_data); make_geo_data(geo_data, fan_data); make_block_data(block_data, measured_fan_data); - info(boost::format("KL on fans: %1%, %2") % KL(measured_fan_data, fan_data, 0) % KL(measured_geo_data, geo_data, 0)); + info(boost::format("KL on fans: %1%, %2%") % KL(measured_fan_data, fan_data, 0) % KL(measured_geo_data, geo_data, 0)); } } } From f0c788b554f4c77a367452379aae44d9d73a1426 Mon Sep 17 00:00:00 2001 From: rts Date: Fri, 30 Aug 2024 13:19:41 -0700 Subject: [PATCH 13/34] Add allocate override using a MLEstimateComponentBasedNormalisation Adds has_processed_data to MLEstimateComponentBasedNormalisation --- .../BinNormalisationPETFromComponents.h | 3 +++ .../ML_estimate_component_based_normalisation.h | 3 +++ .../BinNormalisationPETFromComponents.cxx | 14 ++++++++++++++ .../ML_estimate_component_based_normalisation.cxx | 7 ++++++- 4 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h b/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h index 714ffca9b..355c712ed 100644 --- a/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h +++ b/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h @@ -20,6 +20,7 @@ #ifndef __stir_recon_buildblock_BinNormalisationPETFromComponents_H__ #define __stir_recon_buildblock_BinNormalisationPETFromComponents_H__ +#include "ML_estimate_component_based_normalisation.h" #include "stir/recon_buildblock/BinNormalisation.h" #include "stir/ML_norm.h" #include "stir/ProjData.h" @@ -126,6 +127,8 @@ class BinNormalisationPETFromComponents : public BinNormalisation void allocate(shared_ptr, bool do_eff, bool do_geo, bool do_block = false, bool do_symmetry_per_block = false); + void allocate(shared_ptr pdi_sptr, MLEstimateComponentBasedNormalisation &normalization_estimator); + DetectorEfficiencies& crystal_efficiencies() { return efficiencies; diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h index 79bf94065..35417a8d7 100644 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h @@ -86,6 +86,9 @@ class MLEstimateComponentBasedNormalisation */ void process(); + //! Check if the data has been processed + bool has_processed_data() const; + //! Get the efficiencies, nullptr if not calculated std::shared_ptr get_efficiencies() const; //! Get the geo data, nullptr if not calculated diff --git a/src/recon_buildblock/BinNormalisationPETFromComponents.cxx b/src/recon_buildblock/BinNormalisationPETFromComponents.cxx index f07d50552..bcfbbd1d7 100644 --- a/src/recon_buildblock/BinNormalisationPETFromComponents.cxx +++ b/src/recon_buildblock/BinNormalisationPETFromComponents.cxx @@ -173,6 +173,20 @@ BinNormalisationPETFromComponents::allocate( this->_already_allocated = true; } +void +BinNormalisationPETFromComponents::allocate(shared_ptr 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; +} void BinNormalisationPETFromComponents::create_proj_data() diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index bbb4736bc..8fb9546f6 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -158,7 +158,6 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation( void MLEstimateComponentBasedNormalisation::process() { - data_processed = true; if (do_display) { display(model_fan_data, "model"); @@ -204,6 +203,12 @@ MLEstimateComponentBasedNormalisation::process() info(boost::format("KL on fans: %1%, %2%") % KL(measured_fan_data, fan_data, 0) % KL(measured_geo_data, geo_data, 0)); } } + data_processed = true; +} +bool +MLEstimateComponentBasedNormalisation::has_processed_data() const +{ + return data_processed; } std::shared_ptr From e5e2cc7be10cd1fbd1099730384ca990f765c35a Mon Sep 17 00:00:00 2001 From: rts Date: Fri, 30 Aug 2024 16:06:53 -0700 Subject: [PATCH 14/34] [ci skip] cleanup comments --- .../ML_estimate_component_based_normalisation.cxx | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx index 8fb9546f6..35b688d37 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx @@ -123,11 +123,11 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation( measured_geo_data = GeoData3D(num_physical_axial_crystals_per_basic_unit, num_physical_transaxial_crystals_per_basic_unit / 2, num_physical_rings, - num_physical_detectors_per_ring); // inputes have to be modified + num_physical_detectors_per_ring); norm_geo_data_ptr = std::make_shared(num_physical_axial_crystals_per_basic_unit, num_physical_transaxial_crystals_per_basic_unit / 2, num_physical_rings, - num_physical_detectors_per_ring); // inputes have to be modified + num_physical_detectors_per_ring); measured_block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); norm_block_data_ptr @@ -151,7 +151,7 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation( geo_data = GeoData3D(num_physical_axial_crystals_per_basic_unit, num_physical_transaxial_crystals_per_basic_unit / 2, num_physical_rings, - num_physical_detectors_per_ring); // inputes have to be modified + num_physical_detectors_per_ring); block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); } @@ -186,17 +186,15 @@ MLEstimateComponentBasedNormalisation::process() if (do_geo) { - // geo norm iteration geo_normalization_iteration(iter_num); // Calculate geo norm iteration } if (do_block) { - // block norm iteration block_normalization_iteration(iter_num); // Calculate block norm iteration } - //// print KL for fansums if (do_KL) { + // print KL for fansums make_fan_sum_data(fan_sums, fan_data); make_geo_data(geo_data, fan_data); make_block_data(block_data, measured_fan_data); From 8adb620bcda8c3e6606fa288d7e68556f3011aaf Mon Sep 17 00:00:00 2001 From: rts Date: Tue, 3 Sep 2024 14:04:46 -0700 Subject: [PATCH 15/34] Rename files MLEstimateComponentBasedNormalisation --- .../stir/recon_buildblock/BinNormalisationPETFromComponents.h | 4 ++-- ...ormalisation.h => MLEstimateComponentBasedNormalisation.h} | 0 src/recon_buildblock/CMakeLists.txt | 2 +- ...lisation.cxx => MLEstimateComponentBasedNormalisation.cxx} | 2 +- src/utilities/find_ML_normfactors3D.cxx | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) rename src/include/stir/recon_buildblock/{ML_estimate_component_based_normalisation.h => MLEstimateComponentBasedNormalisation.h} (100%) rename src/recon_buildblock/{ML_estimate_component_based_normalisation.cxx => MLEstimateComponentBasedNormalisation.cxx} (99%) diff --git a/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h b/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h index 355c712ed..0970f0e45 100644 --- a/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h +++ b/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h @@ -20,7 +20,7 @@ #ifndef __stir_recon_buildblock_BinNormalisationPETFromComponents_H__ #define __stir_recon_buildblock_BinNormalisationPETFromComponents_H__ -#include "ML_estimate_component_based_normalisation.h" +#include "MLEstimateComponentBasedNormalisation.h" #include "stir/recon_buildblock/BinNormalisation.h" #include "stir/ML_norm.h" #include "stir/ProjData.h" @@ -35,7 +35,7 @@ START_NAMESPACE_STIR Components currently supported are crystal efficiencies, geometric factors (constrained by symmetry) and block data. The latter were introduced to cope with timing alignment issues between blocks, but are generally - not recommended in the current estimation process (by ML_estimate_component_based_normalisation) + not recommended in the current estimation process (by MLEstimateComponentBasedNormalisation) as the model allows for too much freedom. The detection efficiency of a crystal pair is modelled as diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h similarity index 100% rename from src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h rename to src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index 61286275f..42948aec6 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -61,7 +61,7 @@ set(${dir_LIB_SOURCES} BinNormalisationFromAttenuationImage.cxx BinNormalisationSPECT.cxx BinNormalisationPETFromComponents.cxx - ML_estimate_component_based_normalisation.cxx + MLEstimateComponentBasedNormalisation.cxx GeneralisedPrior.cxx ProjDataRebinning.cxx FourierRebinning.cxx diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx similarity index 99% rename from src/recon_buildblock/ML_estimate_component_based_normalisation.cxx rename to src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx index 35b688d37..62073add7 100644 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx @@ -21,7 +21,7 @@ \author Gefei Chen \author Robert Twyman Skelly */ -#include "stir/recon_buildblock/ML_estimate_component_based_normalisation.h" +#include "stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h" #include "stir/ML_norm.h" #include "stir/Scanner.h" #include "stir/display.h" diff --git a/src/utilities/find_ML_normfactors3D.cxx b/src/utilities/find_ML_normfactors3D.cxx index 0dc8abf52..f7bdcd998 100644 --- a/src/utilities/find_ML_normfactors3D.cxx +++ b/src/utilities/find_ML_normfactors3D.cxx @@ -17,7 +17,7 @@ Just a wrapper around ML_estimate_component_based_normalisation \author Kris Thielemans */ -#include "stir/recon_buildblock/ML_estimate_component_based_normalisation.h" +#include "stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h" #include "stir/CPUTimer.h" #include "stir/info.h" #include "stir/ProjData.h" From 9ff77fedd03924cef4d53db4524108509524fd27 Mon Sep 17 00:00:00 2001 From: rts Date: Tue, 3 Sep 2024 14:14:12 -0700 Subject: [PATCH 16/34] Remove allocate with norm estimator --- .../BinNormalisationPETFromComponents.h | 3 --- .../BinNormalisationPETFromComponents.cxx | 14 -------------- 2 files changed, 17 deletions(-) diff --git a/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h b/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h index 0970f0e45..a2e2f11dc 100644 --- a/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h +++ b/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h @@ -20,7 +20,6 @@ #ifndef __stir_recon_buildblock_BinNormalisationPETFromComponents_H__ #define __stir_recon_buildblock_BinNormalisationPETFromComponents_H__ -#include "MLEstimateComponentBasedNormalisation.h" #include "stir/recon_buildblock/BinNormalisation.h" #include "stir/ML_norm.h" #include "stir/ProjData.h" @@ -127,8 +126,6 @@ class BinNormalisationPETFromComponents : public BinNormalisation void allocate(shared_ptr, bool do_eff, bool do_geo, bool do_block = false, bool do_symmetry_per_block = false); - void allocate(shared_ptr pdi_sptr, MLEstimateComponentBasedNormalisation &normalization_estimator); - DetectorEfficiencies& crystal_efficiencies() { return efficiencies; diff --git a/src/recon_buildblock/BinNormalisationPETFromComponents.cxx b/src/recon_buildblock/BinNormalisationPETFromComponents.cxx index bcfbbd1d7..f07d50552 100644 --- a/src/recon_buildblock/BinNormalisationPETFromComponents.cxx +++ b/src/recon_buildblock/BinNormalisationPETFromComponents.cxx @@ -173,20 +173,6 @@ BinNormalisationPETFromComponents::allocate( this->_already_allocated = true; } -void -BinNormalisationPETFromComponents::allocate(shared_ptr 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; -} void BinNormalisationPETFromComponents::create_proj_data() From 5c1a0f75acc201d6da0c28b148250dedc13d55e6 Mon Sep 17 00:00:00 2001 From: rts Date: Tue, 3 Sep 2024 14:34:40 -0700 Subject: [PATCH 17/34] Add construct_bin_normfactors_from_components to MLEstimateComponentBasedNormalisation --- .../BinNormalisationPETFromComponents.h | 24 ++++++++++++++++--- .../MLEstimateComponentBasedNormalisation.h | 7 ++++++ .../MLEstimateComponentBasedNormalisation.cxx | 22 +++++++++++++++++ src/utilities/apply_normfactors3D.cxx | 6 ++--- 4 files changed, 53 insertions(+), 6 deletions(-) diff --git a/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h b/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h index a2e2f11dc..67f6e0a18 100644 --- a/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h +++ b/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h @@ -126,15 +126,33 @@ class BinNormalisationPETFromComponents : public BinNormalisation void allocate(shared_ptr, bool do_eff, bool do_geo, bool do_block = false, bool do_symmetry_per_block = false); - DetectorEfficiencies& crystal_efficiencies() + + void set_crystal_efficiencies(const DetectorEfficiencies& eff) + { + efficiencies = eff; + } + + DetectorEfficiencies& get_crystal_efficiencies() { return efficiencies; } - GeoData3D& geometric_factors() + + void set_geometric_factors(const GeoData3D& geo) + { + geo_data = geo; + } + + GeoData3D& get_geometric_factors() { return geo_data; } - BlockData3D& block_factors() + + void set_block_factors(const BlockData3D& block) + { + block_data = block; + } + + BlockData3D& get_block_factors() { return block_data; } diff --git a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h index 35417a8d7..2c08f169b 100644 --- a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h +++ b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h @@ -14,10 +14,12 @@ \author Kris Thielemans \author Robert Twyman Skelly */ +#include "BinNormalisationPETFromComponents.h" #include "stir/common.h" #include #include "stir/ProjData.h" #include "stir/ML_norm.h" +#include "stir/ProjDataInMemory.h" START_NAMESPACE_STIR @@ -96,6 +98,8 @@ class MLEstimateComponentBasedNormalisation //! Get the block data, nullptr if not calculated std::shared_ptr get_block_data() const; + BinNormalisationPETFromComponents construct_bin_normfactors_from_components() const; + private: /*! \brief Performs an efficiency iteration to update the efficiancies from the data_fan_sums and model. @@ -164,6 +168,9 @@ class MLEstimateComponentBasedNormalisation //! Whether to save the each iteration of the efficiencies, geo data and block data to file. bool do_save_to_file; + // The projdata info of the measured data + std::shared_ptr projdata_info; + // Calculated variables std::shared_ptr efficiencies_ptr; std::shared_ptr norm_block_data_ptr; diff --git a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx index 62073add7..baed4921d 100644 --- a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx @@ -82,6 +82,7 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation( do_display(do_display_v), do_save_to_file(do_save_to_file_v) { + projdata_info = measured_projdata_v.get_proj_data_info_sptr(); const int num_transaxial_blocks = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); const int num_axial_blocks = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_axial_blocks(); @@ -239,6 +240,27 @@ MLEstimateComponentBasedNormalisation::get_block_data() const return norm_block_data_ptr; } +BinNormalisationPETFromComponents +MLEstimateComponentBasedNormalisation::construct_bin_normfactors_from_components() const +{ + if (!data_processed) + { + error("MLEstimateComponentBasedNormalisation::construct_bin_normfactors_from_components: data has not been processed yet"); + } + auto bin_norm = BinNormalisationPETFromComponents(); + bin_norm.allocate(projdata_info, true, do_geo, do_block, do_symmetry_per_block); + bin_norm.set_crystal_efficiencies(*efficiencies_ptr); + if (do_geo) + { + bin_norm.set_geometric_factors(*norm_geo_data_ptr); + } + if (do_block) + { + bin_norm.set_block_factors(*norm_block_data_ptr); + } + return bin_norm; +} + void MLEstimateComponentBasedNormalisation::efficiency_iteration(const int iter_num, const int eff_iter_num) { diff --git a/src/utilities/apply_normfactors3D.cxx b/src/utilities/apply_normfactors3D.cxx index a235aea2b..0be120d05 100644 --- a/src/utilities/apply_normfactors3D.cxx +++ b/src/utilities/apply_normfactors3D.cxx @@ -73,7 +73,7 @@ main(int argc, char** argv) 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(); + in >> norm.get_crystal_efficiencies(); if (!in) { warning("Error reading %s, using all 1s instead\n", in_filename); @@ -89,7 +89,7 @@ main(int argc, char** argv) char* in_filename = new char[in_filename_prefix.size() + 30]; sprintf(in_filename, "%s_%s_%d.out", in_filename_prefix.c_str(), "geo", iter_num); std::ifstream in(in_filename); - in >> norm.geometric_factors(); + in >> norm.get_geometric_factors(); if (!in) { warning("Error reading %s, using all 1s instead\n", in_filename); @@ -106,7 +106,7 @@ main(int argc, char** argv) char* in_filename = new char[in_filename_prefix.size() + 30]; sprintf(in_filename, "%s_%s_%d.out", in_filename_prefix.c_str(), "block", iter_num); std::ifstream in(in_filename); - in >> norm.block_factors(); + in >> norm.get_block_factors(); if (!in) { warning("Error reading %s, using all 1s instead\n", in_filename); From a114d14cb289bab5e767159668abfc16828816d3 Mon Sep 17 00:00:00 2001 From: rts Date: Tue, 3 Sep 2024 14:42:40 -0700 Subject: [PATCH 18/34] Change efficiencies, norm_block_data, norm_geo_data values from sptr --- .../MLEstimateComponentBasedNormalisation.h | 18 +-- .../MLEstimateComponentBasedNormalisation.cxx | 125 +++++++++--------- 2 files changed, 75 insertions(+), 68 deletions(-) diff --git a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h index 2c08f169b..53748b82b 100644 --- a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h +++ b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h @@ -91,12 +91,12 @@ class MLEstimateComponentBasedNormalisation //! Check if the data has been processed bool has_processed_data() const; - //! Get the efficiencies, nullptr if not calculated - std::shared_ptr get_efficiencies() const; - //! Get the geo data, nullptr if not calculated - std::shared_ptr get_geo_data() const; - //! Get the block data, nullptr if not calculated - std::shared_ptr get_block_data() const; + //! Get the efficiencies, requires process() to be called first + DetectorEfficiencies get_efficiencies() const; + //! Get the geo data, requires process() to be called first and do_geo to be true + GeoData3D get_geo_data() const; + //! Get the block data, requires process() to be called first and do_block to be true + BlockData3D get_block_data() const; BinNormalisationPETFromComponents construct_bin_normfactors_from_components() const; @@ -172,9 +172,9 @@ class MLEstimateComponentBasedNormalisation std::shared_ptr projdata_info; // Calculated variables - std::shared_ptr efficiencies_ptr; - std::shared_ptr norm_block_data_ptr; - std::shared_ptr norm_geo_data_ptr; + DetectorEfficiencies norm_efficiencies; + BlockData3D norm_block_data; + GeoData3D norm_geo_data; //! The threshold for the Kullback-Leibler divergence calculation float threshold_for_KL; FanProjData model_fan_data; diff --git a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx index baed4921d..8e8ad8573 100644 --- a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx @@ -60,18 +60,18 @@ ML_estimate_component_based_normalisation(const std::string& out_filename_prefix estimator.process(); } -MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation( - std::string out_filename_prefix_v, - const ProjData& measured_projdata_v, - const ProjData& model_projdata_v, - const int num_eff_iterations_v, - const int num_iterations_v, - const bool do_geo_v, - const bool do_block_v, - const bool do_symmetry_per_block_v, - const bool do_KL_v, - const bool do_display_v, - const bool do_save_to_file_v) +MLEstimateComponentBasedNormalisation:: +MLEstimateComponentBasedNormalisation(std::string out_filename_prefix_v, + const ProjData& measured_projdata_v, + const ProjData& model_projdata_v, + const int num_eff_iterations_v, + const int num_iterations_v, + const bool do_geo_v, + const bool do_block_v, + const bool do_symmetry_per_block_v, + const bool do_KL_v, + const bool do_display_v, + const bool do_save_to_file_v) : out_filename_prefix(std::move(out_filename_prefix_v)), num_eff_iterations(num_eff_iterations_v), num_iterations(num_iterations_v), @@ -119,20 +119,19 @@ MLEstimateComponentBasedNormalisation::MLEstimateComponentBasedNormalisation( // Setup the data structures given the PET scanner geometry data_fan_sums = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); - efficiencies_ptr = std::make_shared(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); + norm_efficiencies = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); measured_geo_data = GeoData3D(num_physical_axial_crystals_per_basic_unit, num_physical_transaxial_crystals_per_basic_unit / 2, num_physical_rings, num_physical_detectors_per_ring); - norm_geo_data_ptr = std::make_shared(num_physical_axial_crystals_per_basic_unit, - num_physical_transaxial_crystals_per_basic_unit / 2, - num_physical_rings, - num_physical_detectors_per_ring); + norm_geo_data = GeoData3D(num_physical_axial_crystals_per_basic_unit, + num_physical_transaxial_crystals_per_basic_unit / 2, + num_physical_rings, + num_physical_detectors_per_ring); measured_block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); - norm_block_data_ptr - = std::make_shared(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); + norm_block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); make_fan_data_remove_gaps(model_fan_data, model_projdata_v); make_fan_data_remove_gaps(measured_fan_data, measured_projdata_v); @@ -165,15 +164,15 @@ MLEstimateComponentBasedNormalisation::process() } // Initialize the efficiencies, geo data and block data to 1 - efficiencies_ptr->fill(sqrt(data_fan_sums.sum() / model_fan_data.sum())); - norm_geo_data_ptr->fill(1); - norm_block_data_ptr->fill(1); + norm_efficiencies.fill(sqrt(data_fan_sums.sum() / model_fan_data.sum())); + norm_geo_data.fill(1); + norm_block_data.fill(1); for (int iter_num = 1; iter_num <= std::max(num_iterations, 1); ++iter_num) { fan_data = model_fan_data; - apply_geo_norm(fan_data, *norm_geo_data_ptr); - apply_block_norm(fan_data, *norm_block_data_ptr); + apply_geo_norm(fan_data, norm_geo_data); + apply_block_norm(fan_data, norm_block_data); if (do_display) { display(fan_data, "model*geo*block"); @@ -210,34 +209,42 @@ MLEstimateComponentBasedNormalisation::has_processed_data() const return data_processed; } -std::shared_ptr +DetectorEfficiencies MLEstimateComponentBasedNormalisation::get_efficiencies() const { if (!data_processed) { - return nullptr; + error("MLEstimateComponentBasedNormalisation::get_efficiencies: data has not been processed yet"); } - return efficiencies_ptr; + return norm_efficiencies; } -std::shared_ptr +GeoData3D MLEstimateComponentBasedNormalisation::get_geo_data() const { if (!data_processed) { - return nullptr; + error("MLEstimateComponentBasedNormalisation::get_geo_data: data has not been processed yet"); } - return norm_geo_data_ptr; + if (!do_geo) + { + error("MLEstimateComponentBasedNormalisation::get_geo_data: geo data was not calculated"); + } + return norm_geo_data; } -std::shared_ptr +BlockData3D MLEstimateComponentBasedNormalisation::get_block_data() const { if (!data_processed) { - return nullptr; + error("MLEstimateComponentBasedNormalisation::get_block_data: data has not been processed yet"); + } + if (!do_block) + { + error("MLEstimateComponentBasedNormalisation::get_block_data: block data was not calculated"); } - return norm_block_data_ptr; + return norm_block_data; } BinNormalisationPETFromComponents @@ -249,14 +256,14 @@ MLEstimateComponentBasedNormalisation::construct_bin_normfactors_from_components } auto bin_norm = BinNormalisationPETFromComponents(); bin_norm.allocate(projdata_info, true, do_geo, do_block, do_symmetry_per_block); - bin_norm.set_crystal_efficiencies(*efficiencies_ptr); + bin_norm.set_crystal_efficiencies(norm_efficiencies); if (do_geo) { - bin_norm.set_geometric_factors(*norm_geo_data_ptr); + bin_norm.set_geometric_factors(norm_geo_data); } if (do_block) { - bin_norm.set_block_factors(*norm_block_data_ptr); + bin_norm.set_block_factors(norm_block_data); } return bin_norm; } @@ -264,14 +271,14 @@ MLEstimateComponentBasedNormalisation::construct_bin_normfactors_from_components void MLEstimateComponentBasedNormalisation::efficiency_iteration(const int iter_num, const int eff_iter_num) { - iterate_efficiencies(*efficiencies_ptr, data_fan_sums, fan_data); + iterate_efficiencies(norm_efficiencies, data_fan_sums, fan_data); if (do_save_to_file) { write_efficiencies_to_file(iter_num, eff_iter_num); } if (do_KL) { - apply_efficiencies(fan_data, *efficiencies_ptr); + apply_efficiencies(fan_data, norm_efficiencies); std::cerr << "measured*norm min " << measured_fan_data.find_min() << " ,max " << measured_fan_data.find_max() << std::endl; std::cerr << "model*norm min " << fan_data.find_min() << " ,max " << fan_data.find_max() << std::endl; if (do_display) @@ -279,18 +286,18 @@ MLEstimateComponentBasedNormalisation::efficiency_iteration(const int iter_num, info(boost::format("KL %1%") % KL(measured_fan_data, fan_data, threshold_for_KL)); // now restore for further iterations fan_data = model_fan_data; - apply_geo_norm(fan_data, *norm_geo_data_ptr); - apply_block_norm(fan_data, *norm_block_data_ptr); + apply_geo_norm(fan_data, norm_geo_data); + apply_block_norm(fan_data, norm_block_data); } if (do_display) { fan_data.fill(1); - apply_efficiencies(fan_data, *efficiencies_ptr); + apply_efficiencies(fan_data, norm_efficiencies); display(fan_data, "eff norm"); // now restore for further iterations fan_data = model_fan_data; - apply_geo_norm(fan_data, *norm_geo_data_ptr); - apply_block_norm(fan_data, *norm_block_data_ptr); + apply_geo_norm(fan_data, norm_geo_data); + apply_block_norm(fan_data, norm_block_data); } } @@ -298,10 +305,10 @@ void MLEstimateComponentBasedNormalisation::geo_normalization_iteration(int iter_num) { - fan_data = model_fan_data; // Reset fan_data to model_data - apply_efficiencies(fan_data, *efficiencies_ptr); // Apply efficiencies - apply_block_norm(fan_data, *norm_block_data_ptr); // Apply block norm - iterate_geo_norm(*norm_geo_data_ptr, measured_geo_data, fan_data); + fan_data = model_fan_data; // Reset fan_data to model_data + apply_efficiencies(fan_data, norm_efficiencies); // Apply efficiencies + apply_block_norm(fan_data, norm_block_data); // Apply block norm + iterate_geo_norm(norm_geo_data, measured_geo_data, fan_data); if (do_save_to_file) { @@ -309,13 +316,13 @@ MLEstimateComponentBasedNormalisation::geo_normalization_iteration(int iter_num) } if (do_KL) { - apply_geo_norm(fan_data, *norm_geo_data_ptr); + apply_geo_norm(fan_data, norm_geo_data); info(boost::format("KL %1%") % KL(measured_fan_data, fan_data, threshold_for_KL)); } if (do_display) { fan_data.fill(1); - apply_geo_norm(fan_data, *norm_geo_data_ptr); + apply_geo_norm(fan_data, norm_geo_data); display(fan_data, "geo norm"); } } @@ -324,10 +331,10 @@ void MLEstimateComponentBasedNormalisation::block_normalization_iteration(const int iter_num) { - fan_data = model_fan_data; // Reset fan_data to model_data - apply_efficiencies(fan_data, *efficiencies_ptr); // Apply efficiencies - apply_geo_norm(fan_data, *norm_geo_data_ptr); // Apply geo norm - iterate_block_norm(*norm_block_data_ptr, measured_block_data, fan_data); // Iterate block norm calculation + fan_data = model_fan_data; // Reset fan_data to model_data + apply_efficiencies(fan_data, norm_efficiencies); // Apply efficiencies + apply_geo_norm(fan_data, norm_geo_data); // Apply geo norm + iterate_block_norm(norm_block_data, measured_block_data, fan_data); // Iterate block norm calculation if (do_save_to_file) { @@ -335,14 +342,14 @@ MLEstimateComponentBasedNormalisation::block_normalization_iteration(const int i } if (do_KL) { - apply_block_norm(fan_data, *norm_block_data_ptr); + apply_block_norm(fan_data, norm_block_data); info(boost::format("KL %1%") % KL(measured_fan_data, fan_data, threshold_for_KL)); } if (do_display) { fan_data.fill(1); - apply_block_norm(fan_data, *norm_block_data_ptr); - display(*norm_block_data_ptr, "raw block norm"); + apply_block_norm(fan_data, norm_block_data); + display(norm_block_data, "raw block norm"); display(fan_data, "block norm"); } } @@ -353,7 +360,7 @@ MLEstimateComponentBasedNormalisation::write_efficiencies_to_file(const int iter char* out_filename = new char[out_filename_prefix.size() + 30]; sprintf(out_filename, "%s_%s_%d_%d.out", out_filename_prefix.c_str(), "eff", iter_num, eff_iter_num); std::ofstream out(out_filename); - out << *efficiencies_ptr; + out << norm_efficiencies; } void @@ -362,7 +369,7 @@ MLEstimateComponentBasedNormalisation::write_geo_data_to_file(const int iter_num char* out_filename = new char[out_filename_prefix.size() + 30]; sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "geo", iter_num); std::ofstream out(out_filename); - out << *norm_geo_data_ptr; + out << norm_geo_data; } void @@ -371,7 +378,7 @@ MLEstimateComponentBasedNormalisation::write_block_data_to_file(const int iter_n char* out_filename = new char[out_filename_prefix.size() + 30]; sprintf(out_filename, "%s_%s_%d.out", out_filename_prefix.c_str(), "block", iter_num); std::ofstream out(out_filename); - out << *norm_block_data_ptr; + out << norm_block_data; } float From b0e05c9f7a4dab9910451fe4658dd26f809c6ec9 Mon Sep 17 00:00:00 2001 From: rts Date: Tue, 3 Sep 2024 14:46:03 -0700 Subject: [PATCH 19/34] Documentation update. --- .../MLEstimateComponentBasedNormalisation.h | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h index 53748b82b..94c8d60d7 100644 --- a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h +++ b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h @@ -60,8 +60,8 @@ class MLEstimateComponentBasedNormalisation /*! \brief Constructor \param out_filename_prefix The prefix for the output files - \param measured_data The measured projection data - \param model_data The model projection data + \param measured_data_v The measured projection data + \param model_data_v The model projection data \param num_eff_iterations_v The number of (sub-)efficiency iterations to perform per iteration of the algorithm \param num_iterations_v The number of algorithm iterations to perform \param do_geo_v Whether to perform geo normalization calculations @@ -98,6 +98,7 @@ class MLEstimateComponentBasedNormalisation //! Get the block data, requires process() to be called first and do_block to be true BlockData3D get_block_data() const; + //! Construct a BinNormalisationPETFromComponents from the calculated efficiencies, geo data and block data BinNormalisationPETFromComponents construct_bin_normfactors_from_components() const; private: @@ -168,13 +169,17 @@ class MLEstimateComponentBasedNormalisation //! Whether to save the each iteration of the efficiencies, geo data and block data to file. bool do_save_to_file; - // The projdata info of the measured data + //! The projdata info of the measured data std::shared_ptr projdata_info; // Calculated variables + //! The efficiencies calculated by the algorithm DetectorEfficiencies norm_efficiencies; + //! The geo data calculated by the algorithm, if do_geo is true BlockData3D norm_block_data; + //! The block data calculated by the algorithm, if do_block is true GeoData3D norm_geo_data; + //! The threshold for the Kullback-Leibler divergence calculation float threshold_for_KL; FanProjData model_fan_data; From 1d6ca9bc5ad730a312589a29368c2825c4c91e21 Mon Sep 17 00:00:00 2001 From: rts Date: Wed, 4 Sep 2024 13:24:45 -0700 Subject: [PATCH 20/34] Fix variable names --- .../MLEstimateComponentBasedNormalisation.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h index 94c8d60d7..399c38d5e 100644 --- a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h +++ b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h @@ -59,9 +59,9 @@ class MLEstimateComponentBasedNormalisation public: /*! \brief Constructor - \param out_filename_prefix The prefix for the output files - \param measured_data_v The measured projection data - \param model_data_v The model projection data + \param out_filename_prefix_v The prefix for the output files + \param measured_projdata_v The measured projection data + \param model_projdata_v The model projection data \param num_eff_iterations_v The number of (sub-)efficiency iterations to perform per iteration of the algorithm \param num_iterations_v The number of algorithm iterations to perform \param do_geo_v Whether to perform geo normalization calculations @@ -71,9 +71,9 @@ class MLEstimateComponentBasedNormalisation \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, - const ProjData& measured_data_v, - const ProjData& model_data_v, + MLEstimateComponentBasedNormalisation(std::string out_filename_prefix_v, + const ProjData& measured_projdata_v, + const ProjData& model_projdata_v, int num_eff_iterations_v, int num_iterations_v, bool do_geo_v, From 0aae32e0422747dee66bf79dc63667c13273da40 Mon Sep 17 00:00:00 2001 From: rts Date: Wed, 4 Sep 2024 13:25:10 -0700 Subject: [PATCH 21/34] Add guard to ensure measured and model data use the same pdi --- .../MLEstimateComponentBasedNormalisation.cxx | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx index 8e8ad8573..f74d57eb2 100644 --- a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx @@ -82,6 +82,11 @@ MLEstimateComponentBasedNormalisation(std::string out_filename_prefix_v, do_display(do_display_v), do_save_to_file(do_save_to_file_v) { + if (measured_projdata_v.get_proj_data_info_sptr() != model_projdata_v.get_proj_data_info_sptr()) + { + error("MLEstimateComponentBasedNormalisation: measured and model data have different ProjDataInfo"); + } + projdata_info = measured_projdata_v.get_proj_data_info_sptr(); const int num_transaxial_blocks = measured_projdata_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks(); From ca1140c3e077bd1e4806eb8f1e5bb08340e6e381 Mon Sep 17 00:00:00 2001 From: rts Date: Wed, 4 Sep 2024 13:27:42 -0700 Subject: [PATCH 22/34] Add test for MLEstimateComponentBasedNormalisation --- src/recon_test/CMakeLists.txt | 1 + ..._MLEstimateComponentBasedNormalisation.cxx | 109 ++++++++++++++++++ 2 files changed, 110 insertions(+) create mode 100644 src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx diff --git a/src/recon_test/CMakeLists.txt b/src/recon_test/CMakeLists.txt index 93a0eb29e..9e25fcf41 100644 --- a/src/recon_test/CMakeLists.txt +++ b/src/recon_test/CMakeLists.txt @@ -26,6 +26,7 @@ set(${dir_SIMPLE_TEST_EXE_SOURCES} test_FBP3DRP.cxx test_blocks_on_cylindrical_projectors.cxx test_geometry_blocks_on_cylindrical.cxx + test_MLEstimateComponentBasedNormalisation.cxx ) diff --git a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx new file mode 100644 index 000000000..8bd3b0121 --- /dev/null +++ b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx @@ -0,0 +1,109 @@ +/* +Copyright (C) 2024, Robert Twyman skelly + This file is part of STIR. + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup recon_test + \brief Test program to ensure MLEstimateComponentBasedNormalisation works as expected and + produces the correct normalization factors from components + \author Robert Twyman Skelly +*/ + +#include "stir/RunTests.h" +#include "stir/num_threads.h" +#include "stir/recon_buildblock/ForwardProjectorByBinUsingProjMatrixByBin.h" +#include "stir/ProjDataInMemory.h" +#include "stir/SeparableCartesianMetzImageFilter.h" +#include "stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h" + +START_NAMESPACE_STIR + +/*! + \ingroup test + \brief Test class for MLEstimateComponentBasedNormalisation +*/ +class MLEstimateComponentBasedNormalisationTest : public RunTests +{ +public: + //! Constructor that can take some input data to run the test with + MLEstimateComponentBasedNormalisationTest() = default; + + ~MLEstimateComponentBasedNormalisationTest() override = default; + + void run_tests() override + { + + const auto scanner = std::make_shared(*Scanner::get_scanner_from_name("ECAT EXACT")); + const auto exam_info = std::make_shared(ImagingModality::PT); + exam_info->patient_position = PatientPosition(PatientPosition::HFS); + + const shared_ptr projdata_info(ProjDataInfo::ProjDataInfoCTI( + /*scanner_ptr*/ scanner, + /*span*/ 1, + /*max_delta*/ 23, + /*views*/ scanner->get_num_detectors_per_ring() / 2, + /*tang_pos*/ scanner->get_default_num_arccorrected_bins(), + /*arc_corrected*/ false, + /*tof_mash_factor*/ -1)); + + const auto measured_projdata = std::make_shared(ProjDataInMemory(exam_info, projdata_info, false)); + measured_projdata->fill(1.F); + + const auto model_projdata = std::make_shared(*measured_projdata); + model_projdata->fill(2.F); + + constexpr int num_eff_iterations = 6; + constexpr int num_iterations = 2; + constexpr bool do_geo = false; + constexpr bool do_block = false; + constexpr bool do_symmetry_per_block = false; + constexpr bool do_KL = false; + constexpr bool do_display = false; + constexpr bool do_save_to_file = false; + auto ml_estimator = MLEstimateComponentBasedNormalisation("", + *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(); + + auto bin_normalization = ml_estimator.construct_bin_normfactors_from_components(); + bin_normalization.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_normalization.apply(normalization_projdata); + check_if_less( + normalization_projdata.find_max(), 2.f * 1.1f, "The max value of the normalization projdata is greater than expected"); + check_if_less( + 2.f * 0.9f, normalization_projdata.find_max(), "The max value of the normalization projdata is less than expected"); + check_if_less( + normalization_projdata.find_min(), 2.f * 0.9f, "The min value of the normalization projdata is greater than expected"); + check_if_less( + 2.f * 1.1f, normalization_projdata.find_min(), "The min value of the normalization projdata is less than expected"); + + // ASSERT_EQ() + } +}; +END_NAMESPACE_STIR + +USING_NAMESPACE_STIR + +int +main() +{ + Verbosity::set(1); + MLEstimateComponentBasedNormalisationTest tests; + tests.run_tests(); + return tests.main_return_value(); +} From 2935b99b57bd595c8c198521df70ea85c798c0bc Mon Sep 17 00:00:00 2001 From: rts Date: Wed, 4 Sep 2024 13:59:41 -0700 Subject: [PATCH 23/34] Correct comparison to objects, not sptr --- src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx index f74d57eb2..0e83e5a07 100644 --- a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx @@ -82,7 +82,7 @@ MLEstimateComponentBasedNormalisation(std::string out_filename_prefix_v, do_display(do_display_v), do_save_to_file(do_save_to_file_v) { - if (measured_projdata_v.get_proj_data_info_sptr() != model_projdata_v.get_proj_data_info_sptr()) + if (*measured_projdata_v.get_proj_data_info_sptr() != *model_projdata_v.get_proj_data_info_sptr()) { error("MLEstimateComponentBasedNormalisation: measured and model data have different ProjDataInfo"); } From 3c115026df469f338f5c25f7cbaac26a171e330b Mon Sep 17 00:00:00 2001 From: rts Date: Wed, 4 Sep 2024 15:08:16 -0700 Subject: [PATCH 24/34] [ci skip] Fix min/max tolerance check --- .../test_MLEstimateComponentBasedNormalisation.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx index 8bd3b0121..eae7ea081 100644 --- a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx @@ -82,15 +82,16 @@ class MLEstimateComponentBasedNormalisationTest : public RunTests // Compute the projdata ProjDataInMemory normalization_projdata(*measured_projdata); normalization_projdata.fill(1.F); + bin_normalization.apply(normalization_projdata); check_if_less( normalization_projdata.find_max(), 2.f * 1.1f, "The max value of the normalization projdata is greater than expected"); check_if_less( 2.f * 0.9f, normalization_projdata.find_max(), "The max value of the normalization projdata is less than expected"); check_if_less( - normalization_projdata.find_min(), 2.f * 0.9f, "The min value of the normalization projdata is greater than expected"); + normalization_projdata.find_min(), 2.f * 1.1f, "The min value of the normalization projdata is greater than expected"); check_if_less( - 2.f * 1.1f, normalization_projdata.find_min(), "The min value of the normalization projdata is less than expected"); + 2.f * 0.9f, normalization_projdata.find_min(), "The min value of the normalization projdata is less than expected"); // ASSERT_EQ() } From 7c18bc742d5aa00932ef00d1f2514ed8d712ad03 Mon Sep 17 00:00:00 2001 From: robbietuk Date: Thu, 12 Sep 2024 13:10:36 -0700 Subject: [PATCH 25/34] Change scanner to Discovery This allows the test to pass and ignores the to ignore even num_tangential_poss issue --- src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx index eae7ea081..e0561d488 100644 --- a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx @@ -36,7 +36,7 @@ class MLEstimateComponentBasedNormalisationTest : public RunTests void run_tests() override { - const auto scanner = std::make_shared(*Scanner::get_scanner_from_name("ECAT EXACT")); + const auto scanner = std::make_shared(*Scanner::get_scanner_from_name("Discovery 690")); const auto exam_info = std::make_shared(ImagingModality::PT); exam_info->patient_position = PatientPosition(PatientPosition::HFS); From 42340f2d6a3df4e1b41485d7c4fb2c8d5c7db526 Mon Sep 17 00:00:00 2001 From: robbietuk Date: Thu, 12 Sep 2024 13:10:52 -0700 Subject: [PATCH 26/34] Add additional test on efficiency values --- .../test_MLEstimateComponentBasedNormalisation.cxx | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx index e0561d488..34a1ce7b5 100644 --- a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx @@ -76,6 +76,13 @@ class MLEstimateComponentBasedNormalisationTest : public RunTests do_save_to_file); ml_estimator.process(); + // Check the efficiencies, with measured data as uniform 1s and model data as uniform 2s, expect this to be ~0.707 + auto efficiencies = ml_estimator.get_efficiencies(); + check_if_less(efficiencies.find_max(), 0.75, "The max value of the efficiencies is greater than expected"); + check_if_less(0.65, efficiencies.find_max(), "The max value of the efficiencies is less than expected"); + check_if_less(efficiencies.find_min(), 0.75, "The min value of the efficiencies is greater than expected"); + check_if_less(0.65, efficiencies.find_min(), "The min value of the efficiencies is less than expected"); + auto bin_normalization = ml_estimator.construct_bin_normfactors_from_components(); bin_normalization.set_up(measured_projdata->get_exam_info_sptr(), measured_projdata->get_proj_data_info_sptr()); From f424ca80605538ff1d6c5361e6d53f2a3b6078aa Mon Sep 17 00:00:00 2001 From: robbietuk Date: Thu, 12 Sep 2024 13:11:38 -0700 Subject: [PATCH 27/34] Ensure normalization_projdata max and min values are about 2.0 --- ...test_MLEstimateComponentBasedNormalisation.cxx | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx index 34a1ce7b5..8b714ac16 100644 --- a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx @@ -89,18 +89,15 @@ class MLEstimateComponentBasedNormalisationTest : public RunTests // Compute the projdata ProjDataInMemory normalization_projdata(*measured_projdata); normalization_projdata.fill(1.F); - bin_normalization.apply(normalization_projdata); + + // Check the normalization factors, with measured data as uniform 1s and model data as uniform 2s, expect this to be 2.0f check_if_less( - normalization_projdata.find_max(), 2.f * 1.1f, "The max value of the normalization projdata is greater than expected"); - check_if_less( - 2.f * 0.9f, normalization_projdata.find_max(), "The max value of the normalization projdata is less than expected"); - check_if_less( - normalization_projdata.find_min(), 2.f * 1.1f, "The min value of the normalization projdata is greater than expected"); + normalization_projdata.find_max(), 2.2f, "The max value of the normalization projdata is greater than expected"); + check_if_less(1.8f, normalization_projdata.find_max(), "The max value of the normalization projdata is less than expected"); check_if_less( - 2.f * 0.9f, normalization_projdata.find_min(), "The min value of the normalization projdata is less than expected"); - - // ASSERT_EQ() + normalization_projdata.find_min(), 2.2f, "The min value of the normalization projdata is greater than expected"); + check_if_less(1.8f, normalization_projdata.find_min(), "The min value of the normalization projdata is less than expected"); } }; END_NAMESPACE_STIR From a5f4b9108596a6cfbc4a5272ed5e7f8edd6ca07b Mon Sep 17 00:00:00 2001 From: robbietuk Date: Fri, 13 Sep 2024 06:43:23 -0700 Subject: [PATCH 28/34] Make test use check_if_equal --- .../test_MLEstimateComponentBasedNormalisation.cxx | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx index 8b714ac16..5c059c459 100644 --- a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx @@ -92,12 +92,11 @@ class MLEstimateComponentBasedNormalisationTest : public RunTests bin_normalization.apply(normalization_projdata); // Check the normalization factors, with measured data as uniform 1s and model data as uniform 2s, expect this to be 2.0f - check_if_less( - normalization_projdata.find_max(), 2.2f, "The max value of the normalization projdata is greater than expected"); - check_if_less(1.8f, normalization_projdata.find_max(), "The max value of the normalization projdata is less than expected"); - check_if_less( - normalization_projdata.find_min(), 2.2f, "The min value of the normalization projdata is greater than expected"); - check_if_less(1.8f, normalization_projdata.find_min(), "The min value of the normalization projdata is less than expected"); + const auto prev_tolerance = get_tolerance(); + set_tolerance(1e-3); + check_if_equal(normalization_projdata.find_max(), 2.f, "The max value of the normalization projdata is not 2.0"); + check_if_equal(normalization_projdata.find_min(), 2.f, "The min value of the normalization projdata is not 0.0"); + set_tolerance(prev_tolerance); } }; END_NAMESPACE_STIR From 1eefeffde3a7673f6ccf91bba56d30aaacbb4454 Mon Sep 17 00:00:00 2001 From: robbietuk Date: Fri, 13 Sep 2024 06:44:02 -0700 Subject: [PATCH 29/34] Split test_normalization_calculation_with_efficiencies from the run_tests function Also documentation --- .../test_MLEstimateComponentBasedNormalisation.cxx | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx index 5c059c459..4c190a2a8 100644 --- a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx @@ -8,7 +8,8 @@ Copyright (C) 2024, Robert Twyman skelly \file \ingroup recon_test \brief Test program to ensure MLEstimateComponentBasedNormalisation works as expected and - produces the correct normalization factors from components + produces the correct normalization factors from components. + This is a very basic test that ensure the basic functionality. \author Robert Twyman Skelly */ @@ -33,9 +34,15 @@ class MLEstimateComponentBasedNormalisationTest : public RunTests ~MLEstimateComponentBasedNormalisationTest() override = default; - void run_tests() override - { + void run_tests() override { test_normalization_calculation_with_efficiencies(); } + /*! Runs a test to check the normalization calculation with efficiencies only. + * Uses two uniform projdata sets, one with 1s and the other with 2s to estimate the normalization factors. + * The efficiencies are then checked to ensure they are within the expected range. + * The normalization factors are then computed from PET components and applied to a uniform projdata set. + */ + void test_normalization_calculation_with_efficiencies() + { const auto scanner = std::make_shared(*Scanner::get_scanner_from_name("Discovery 690")); const auto exam_info = std::make_shared(ImagingModality::PT); exam_info->patient_position = PatientPosition(PatientPosition::HFS); From dcee2518e97644e418dc1db278b6c03d36a1f4a0 Mon Sep 17 00:00:00 2001 From: robbietuk Date: Fri, 13 Sep 2024 06:48:19 -0700 Subject: [PATCH 30/34] Improve usage of data_is_processed() --- .../MLEstimateComponentBasedNormalisation.h | 3 ++- .../MLEstimateComponentBasedNormalisation.cxx | 10 +++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h index 399c38d5e..2ec0a1fb0 100644 --- a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h +++ b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h @@ -89,7 +89,7 @@ class MLEstimateComponentBasedNormalisation void process(); //! Check if the data has been processed - bool has_processed_data() const; + bool data_is_processed() const; //! Get the efficiencies, requires process() to be called first DetectorEfficiencies get_efficiencies() const; @@ -188,6 +188,7 @@ class MLEstimateComponentBasedNormalisation BlockData3D measured_block_data; GeoData3D measured_geo_data; + //! Boolean to check if the data has been processed, see has_data_been_processed() bool data_processed = false; // do_KL specific varaibles diff --git a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx index 0e83e5a07..eb74265d5 100644 --- a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx @@ -209,7 +209,7 @@ MLEstimateComponentBasedNormalisation::process() data_processed = true; } bool -MLEstimateComponentBasedNormalisation::has_processed_data() const +MLEstimateComponentBasedNormalisation::data_is_processed() const { return data_processed; } @@ -217,7 +217,7 @@ MLEstimateComponentBasedNormalisation::has_processed_data() const DetectorEfficiencies MLEstimateComponentBasedNormalisation::get_efficiencies() const { - if (!data_processed) + if (!this->data_is_processed()) { error("MLEstimateComponentBasedNormalisation::get_efficiencies: data has not been processed yet"); } @@ -227,7 +227,7 @@ MLEstimateComponentBasedNormalisation::get_efficiencies() const GeoData3D MLEstimateComponentBasedNormalisation::get_geo_data() const { - if (!data_processed) + if (!this->data_is_processed()) { error("MLEstimateComponentBasedNormalisation::get_geo_data: data has not been processed yet"); } @@ -241,7 +241,7 @@ MLEstimateComponentBasedNormalisation::get_geo_data() const BlockData3D MLEstimateComponentBasedNormalisation::get_block_data() const { - if (!data_processed) + if (!this->data_is_processed()) { error("MLEstimateComponentBasedNormalisation::get_block_data: data has not been processed yet"); } @@ -255,7 +255,7 @@ MLEstimateComponentBasedNormalisation::get_block_data() const BinNormalisationPETFromComponents MLEstimateComponentBasedNormalisation::construct_bin_normfactors_from_components() const { - if (!data_processed) + if (!this->data_is_processed()) { error("MLEstimateComponentBasedNormalisation::construct_bin_normfactors_from_components: data has not been processed yet"); } From 63f85607e11156e593288f37150433d778c289d9 Mon Sep 17 00:00:00 2001 From: robbietuk Date: Fri, 13 Sep 2024 07:07:45 -0700 Subject: [PATCH 31/34] Add getters and setters --- .../MLEstimateComponentBasedNormalisation.h | 121 +++++++++++++++++- .../MLEstimateComponentBasedNormalisation.cxx | 74 ++++++++++- ..._MLEstimateComponentBasedNormalisation.cxx | 2 +- 3 files changed, 187 insertions(+), 10 deletions(-) diff --git a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h index 2ec0a1fb0..d7bb12a6d 100644 --- a/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h +++ b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h @@ -89,7 +89,7 @@ class MLEstimateComponentBasedNormalisation void process(); //! Check if the data has been processed - bool data_is_processed() const; + bool get_data_is_processed() const; //! Get the efficiencies, requires process() to be called first DetectorEfficiencies get_efficiencies() const; @@ -99,7 +99,121 @@ class MLEstimateComponentBasedNormalisation BlockData3D get_block_data() const; //! Construct a BinNormalisationPETFromComponents from the calculated efficiencies, geo data and block data - BinNormalisationPETFromComponents construct_bin_normfactors_from_components() const; + BinNormalisationPETFromComponents construct_bin_norm_from_pet_components() const; + + /*! + \brief Get the output filename prefix. + \return The output filename prefix. + */ + std::string get_output_filename_prefix() const { return out_filename_prefix; } + + /*! + \brief Set the output filename prefix. + \param out_filename_prefix The new output filename prefix. + */ + void set_output_filename_prefix(const std::string& out_filename_prefix); + + /*! + \brief Get the number of efficiency iterations. + There are the inner iterations of the algorithm, num_eff_iterations are performed per iteration. + \return The number of efficiency iterations. + */ + int get_num_eff_iterations() const { return num_eff_iterations; } + + /*! + \brief Set the number of efficiency iterations. + There are the inner iterations of the algorithm, num_eff_iterations are performed per iteration. + \param num_eff_iterations The new number of efficiency iterations. + */ + void set_num_eff_iterations(int num_eff_iterations); + + /*! + \brief Get the number of iterations. + These are the outer iterations of the algorithm. + \return The number of iterations. + */ + int get_num_iterations() const { return num_iterations; } + + /*! + \brief Set the number of iterations. + These are the outer iterations of the algorithm. + \param num_iterations The new number of iterations. + */ + void set_num_iterations(int num_iterations); + + /*! + \brief Check if geo normalization is enabled. + \return True if geo normalization is enabled, false otherwise. + */ + bool get_enable_geo_norm_calculation() const { return do_geo; } + + /*! + \brief Enable or disable geo normalization. + \param do_geo True to enable geo normalization, false to disable. + */ + void set_enable_geo_norm_calculation(bool do_geo); + + /*! + \brief Check if block normalization is enabled. + \return True if block normalization is enabled, false otherwise. + */ + bool get_enable_block_norm_calculation() const { return do_block; } + + /*! + \brief Enable or disable block normalization. + \param do_block True to enable block normalization, false to disable. + */ + void set_enable_block_norm_calculation(bool do_block); + + /*! + \brief Check if symmetry per block is enabled. + \return True if symmetry per block is enabled, false otherwise. + */ + bool get_enable_symmetry_per_block() const { return do_symmetry_per_block; } + + /*! + \brief Enable or disable symmetry per block. + \param do_symmetry_per_block True to enable symmetry per block, false to disable. + */ + void set_enable_symmetry_per_block(bool do_symmetry_per_block); + + /*! + \brief Check if KL divergence calculation is enabled. + This does not impact the calculation of the normalisation factors, only the display of the KL value. + \return True if KL divergence calculation is enabled, false otherwise. + */ + bool get_do_kl_calculation() const { return do_KL; } + + /*! + \brief Enable or disable KL divergence calculation. + This does not impact the calculation of the normalisation factors, it only prints the KL value to console. + \param do_kl True to enable KL divergence calculation, false to disable. + */ + void set_do_kl_calculation(bool do_kl); + + /*! + \brief Check if display is enabled. + \return True if display is enabled, false otherwise. + */ + bool get_write_display_data() const { return do_display; } + + /*! + \brief Enable or disable display. + \param do_display True to enable display, false to disable. + */ + void set_write_display_data(bool do_display); + + /*! + \brief Check if saving to file is enabled. + \return True if saving to file is enabled, false otherwise. + */ + bool get_write_intermediates_to_file() const { return do_save_to_file; } + + /*! + \brief Enable or disable saving to file. + \param do_save_to_file True to enable saving to file, false to disable. + */ + void set_write_intermediates_to_file(bool do_save_to_file); private: /*! @@ -166,7 +280,8 @@ class MLEstimateComponentBasedNormalisation bool do_KL; //! Whether to display the progress of the algorithm. bool do_display; - //! Whether to save the each iteration of the efficiencies, geo data and block data to file. + //! Whether to save each iteration of the efficiencies, geo data and block data to file. + //! This will only work if the out_filename_prefix is set. bool do_save_to_file; //! The projdata info of the measured data diff --git a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx index eb74265d5..e9e6fa844 100644 --- a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx @@ -209,7 +209,7 @@ MLEstimateComponentBasedNormalisation::process() data_processed = true; } bool -MLEstimateComponentBasedNormalisation::data_is_processed() const +MLEstimateComponentBasedNormalisation::get_data_is_processed() const { return data_processed; } @@ -217,7 +217,7 @@ MLEstimateComponentBasedNormalisation::data_is_processed() const DetectorEfficiencies MLEstimateComponentBasedNormalisation::get_efficiencies() const { - if (!this->data_is_processed()) + if (!this->get_data_is_processed()) { error("MLEstimateComponentBasedNormalisation::get_efficiencies: data has not been processed yet"); } @@ -227,7 +227,7 @@ MLEstimateComponentBasedNormalisation::get_efficiencies() const GeoData3D MLEstimateComponentBasedNormalisation::get_geo_data() const { - if (!this->data_is_processed()) + if (!this->get_data_is_processed()) { error("MLEstimateComponentBasedNormalisation::get_geo_data: data has not been processed yet"); } @@ -241,7 +241,7 @@ MLEstimateComponentBasedNormalisation::get_geo_data() const BlockData3D MLEstimateComponentBasedNormalisation::get_block_data() const { - if (!this->data_is_processed()) + if (!this->get_data_is_processed()) { error("MLEstimateComponentBasedNormalisation::get_block_data: data has not been processed yet"); } @@ -253,9 +253,9 @@ MLEstimateComponentBasedNormalisation::get_block_data() const } BinNormalisationPETFromComponents -MLEstimateComponentBasedNormalisation::construct_bin_normfactors_from_components() const +MLEstimateComponentBasedNormalisation::construct_bin_norm_from_pet_components() const { - if (!this->data_is_processed()) + if (!this->get_data_is_processed()) { error("MLEstimateComponentBasedNormalisation::construct_bin_normfactors_from_components: data has not been processed yet"); } @@ -273,6 +273,68 @@ MLEstimateComponentBasedNormalisation::construct_bin_normfactors_from_components return bin_norm; } +void +MLEstimateComponentBasedNormalisation::set_output_filename_prefix(const std::string& out_filename_prefix) +{ + this->out_filename_prefix = out_filename_prefix; +} + +void +MLEstimateComponentBasedNormalisation::set_num_eff_iterations(int num_eff_iterations) +{ + data_processed = false; + this->num_eff_iterations = num_eff_iterations; +} + +void +MLEstimateComponentBasedNormalisation::set_num_iterations(int num_iterations) +{ + data_processed = false; + this->num_iterations = num_iterations; +} + +void +MLEstimateComponentBasedNormalisation::set_enable_geo_norm_calculation(bool do_geo) +{ + data_processed = false; + this->do_geo = do_geo; +} + +void +MLEstimateComponentBasedNormalisation::set_enable_block_norm_calculation(bool do_block) +{ + data_processed = false; + this->do_block = do_block; +} + +void +MLEstimateComponentBasedNormalisation::set_enable_symmetry_per_block(bool do_symmetry_per_block) +{ + data_processed = false; + this->do_symmetry_per_block = do_symmetry_per_block; +} + +void +MLEstimateComponentBasedNormalisation::set_do_kl_calculation(bool do_kl) +{ + data_processed = false; + do_KL = do_kl; +} + +void +MLEstimateComponentBasedNormalisation::set_write_display_data(bool do_display) +{ + data_processed = false; + this->do_display = do_display; +} + +void +MLEstimateComponentBasedNormalisation::set_write_intermediates_to_file(bool do_save_to_file) +{ + data_processed = false; + this->do_save_to_file = do_save_to_file; +} + void MLEstimateComponentBasedNormalisation::efficiency_iteration(const int iter_num, const int eff_iter_num) { diff --git a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx index 4c190a2a8..2d6535695 100644 --- a/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_test/test_MLEstimateComponentBasedNormalisation.cxx @@ -90,7 +90,7 @@ class MLEstimateComponentBasedNormalisationTest : public RunTests check_if_less(efficiencies.find_min(), 0.75, "The min value of the efficiencies is greater than expected"); check_if_less(0.65, efficiencies.find_min(), "The min value of the efficiencies is less than expected"); - auto bin_normalization = ml_estimator.construct_bin_normfactors_from_components(); + auto bin_normalization = ml_estimator.construct_bin_norm_from_pet_components(); bin_normalization.set_up(measured_projdata->get_exam_info_sptr(), measured_projdata->get_proj_data_info_sptr()); // Compute the projdata From 3d7e9ff09a1e4e178823947ed730d46484689f66 Mon Sep 17 00:00:00 2001 From: robbietuk Date: Fri, 13 Sep 2024 09:26:00 -0700 Subject: [PATCH 32/34] [ci skip] Add const to args --- .../MLEstimateComponentBasedNormalisation.cxx | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx index e9e6fa844..5f4eb26b6 100644 --- a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx @@ -151,7 +151,7 @@ MLEstimateComponentBasedNormalisation(std::string out_filename_prefix_v, display(measured_block_data, "raw block data from measurements"); } - // Compute the do_KL specific varaibles from the measured data + // Compute the do_KL specific variables from the measured data fan_sums = DetectorEfficiencies(IndexRange2D(num_physical_rings, num_physical_detectors_per_ring)); geo_data = GeoData3D(num_physical_axial_crystals_per_basic_unit, num_physical_transaxial_crystals_per_basic_unit / 2, @@ -280,56 +280,56 @@ MLEstimateComponentBasedNormalisation::set_output_filename_prefix(const std::str } void -MLEstimateComponentBasedNormalisation::set_num_eff_iterations(int num_eff_iterations) +MLEstimateComponentBasedNormalisation::set_num_eff_iterations(const int num_eff_iterations) { data_processed = false; this->num_eff_iterations = num_eff_iterations; } void -MLEstimateComponentBasedNormalisation::set_num_iterations(int num_iterations) +MLEstimateComponentBasedNormalisation::set_num_iterations(const int num_iterations) { data_processed = false; this->num_iterations = num_iterations; } void -MLEstimateComponentBasedNormalisation::set_enable_geo_norm_calculation(bool do_geo) +MLEstimateComponentBasedNormalisation::set_enable_geo_norm_calculation(const bool do_geo) { data_processed = false; this->do_geo = do_geo; } void -MLEstimateComponentBasedNormalisation::set_enable_block_norm_calculation(bool do_block) +MLEstimateComponentBasedNormalisation::set_enable_block_norm_calculation(const bool do_block) { data_processed = false; this->do_block = do_block; } void -MLEstimateComponentBasedNormalisation::set_enable_symmetry_per_block(bool do_symmetry_per_block) +MLEstimateComponentBasedNormalisation::set_enable_symmetry_per_block(const bool do_symmetry_per_block) { data_processed = false; this->do_symmetry_per_block = do_symmetry_per_block; } void -MLEstimateComponentBasedNormalisation::set_do_kl_calculation(bool do_kl) +MLEstimateComponentBasedNormalisation::set_do_kl_calculation(const bool do_kl) { data_processed = false; do_KL = do_kl; } void -MLEstimateComponentBasedNormalisation::set_write_display_data(bool do_display) +MLEstimateComponentBasedNormalisation::set_write_display_data(const bool do_display) { data_processed = false; this->do_display = do_display; } void -MLEstimateComponentBasedNormalisation::set_write_intermediates_to_file(bool do_save_to_file) +MLEstimateComponentBasedNormalisation::set_write_intermediates_to_file(const bool do_save_to_file) { data_processed = false; this->do_save_to_file = do_save_to_file; @@ -369,7 +369,7 @@ MLEstimateComponentBasedNormalisation::efficiency_iteration(const int iter_num, } void -MLEstimateComponentBasedNormalisation::geo_normalization_iteration(int iter_num) +MLEstimateComponentBasedNormalisation::geo_normalization_iteration(const int iter_num) { fan_data = model_fan_data; // Reset fan_data to model_data From 55897ee202edeeed1ae3e4fc8749159ec08b21c0 Mon Sep 17 00:00:00 2001 From: robbietuk Date: Fri, 13 Sep 2024 09:32:20 -0700 Subject: [PATCH 33/34] [ci skip] Update release 6.3 notes --- documentation/release_6.3.htm | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/documentation/release_6.3.htm b/documentation/release_6.3.htm index 54a9c63bc..752ab841e 100644 --- a/documentation/release_6.3.htm +++ b/documentation/release_6.3.htm @@ -35,6 +35,11 @@

New functionality

Changed functionality

+
  • + Rework MLEstimateComponentBasedNormalisation into a class. + This allows for more flexibility in the use of the class without calling CLI executables.
    + PR # 1499 +

    • Bug fixes

    @@ -68,6 +73,12 @@

    Test changes

    C++ tests

    +
  • + Added MLEstimateComponentBasedNormalisationTest that tests the basic functionality of the + MLEstimateComponentBasedNormalisation class and creation of normalisation factors from + the calculated PET components.
    + PR # 1499 +

    • recon_test_pack

    From cbe81ec73b2d2f63aedb4bedd8bf4071d38940c5 Mon Sep 17 00:00:00 2001 From: robbietuk Date: Mon, 16 Sep 2024 14:52:30 -0700 Subject: [PATCH 34/34] [ci skip] Fix error message --- src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx index 5f4eb26b6..51e8e5ae9 100644 --- a/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx +++ b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx @@ -257,7 +257,7 @@ MLEstimateComponentBasedNormalisation::construct_bin_norm_from_pet_components() { if (!this->get_data_is_processed()) { - error("MLEstimateComponentBasedNormalisation::construct_bin_normfactors_from_components: data has not been processed yet"); + error("MLEstimateComponentBasedNormalisation::construct_bin_norm_from_pet_components: data has not been processed yet"); } auto bin_norm = BinNormalisationPETFromComponents(); bin_norm.allocate(projdata_info, true, do_geo, do_block, do_symmetry_per_block);