diff --git a/documentation/release_6.3.htm b/documentation/release_6.3.htm index 54a9c63bc2..752ab841ed 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

    diff --git a/src/buildblock/ML_norm.cxx b/src/buildblock/ML_norm.cxx index bc8c109109..986ed3b7ab 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/BinNormalisationPETFromComponents.h b/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h index 714ffca9b9..67f6e0a182 100644 --- a/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h +++ b/src/include/stir/recon_buildblock/BinNormalisationPETFromComponents.h @@ -34,7 +34,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 @@ -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 new file mode 100644 index 0000000000..d7bb12a6db --- /dev/null +++ b/src/include/stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h @@ -0,0 +1,316 @@ +/* + Copyright (C) 2022, University College London + 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_buildblock + \brief Declaration of stir::ML_estimate_component_based_normalisation + \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 + +/*! + \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. +*/ +void ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, + const ProjData& measured_projdata, + const ProjData& model_projdata, + 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 = true); + +/*! + \brief Find normalisation factors using a maximum likelihood approach + \ingroup recon_buildblock +*/ +class MLEstimateComponentBasedNormalisation +{ +public: + /*! + \brief Constructor + \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 + \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_v, + const ProjData& measured_projdata_v, + const ProjData& model_projdata_v, + 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(); + + //! Check if the data has been processed + bool get_data_is_processed() 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; + + //! Construct a BinNormalisationPETFromComponents from the calculated efficiencies, geo data and block data + 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: + /*! + \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); + + /*! + \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 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 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 + 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; + FanProjData fan_data; + DetectorEfficiencies data_fan_sums; + 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 + FanProjData measured_fan_data; + DetectorEfficiencies fan_sums; + GeoData3D geo_data; + BlockData3D block_data; +}; + +END_NAMESPACE_STIR 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 deleted file mode 100644 index 9fcbc0650e..0000000000 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ /dev/null @@ -1,38 +0,0 @@ -/* - Copyright (C) 2022, University College London - This file is part of STIR. - - SPDX-License-Identifier: Apache-2.0 - - See STIR/LICENSE.txt for details -*/ -/*! - \file - \ingroup recon_buildblock - \brief Declaration of stir::ML_estimate_component_based_normalisation - \author Kris Thielemans - */ -#include "stir/common.h" -#include - -START_NAMESPACE_STIR - -class ProjData; - -/*! - \brief Find normalisation factors using a maximum likelihood approach - - \ingroup recon_buildblock -*/ -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); - -END_NAMESPACE_STIR diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index 61286275f7..42948aec66 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/MLEstimateComponentBasedNormalisation.cxx b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx new file mode 100644 index 0000000000..51e8e5ae9c --- /dev/null +++ b/src/recon_buildblock/MLEstimateComponentBasedNormalisation.cxx @@ -0,0 +1,488 @@ +/* + Copyright (C) 2001- 2008, Hammersmith Imanet Ltd + Copyright (C) 2019-2020, 2022, University College London + Copyright (C) 2016-2017, PETsys Electronics + Copyright (C) 2021, Gefei Chen + 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_buildblock + \brief Implementation of ML_estimate_component_based_normalisation + + \author Kris Thielemans + \author Tahereh Niknejad + \author Gefei Chen + \author Robert Twyman Skelly + */ +#include "stir/recon_buildblock/MLEstimateComponentBasedNormalisation.h" +#include "stir/ML_norm.h" +#include "stir/Scanner.h" +#include "stir/display.h" +#include "stir/info.h" +#include "stir/ProjData.h" +#include "stir/stream.h" +#include +#include + +START_NAMESPACE_STIR + +void +ML_estimate_component_based_normalisation(const std::string& out_filename_prefix, + const ProjData& measured_projdata, + const ProjData& model_projdata, + 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_projdata, + model_projdata, + num_eff_iterations, + num_iterations, + do_geo, + do_block, + do_symmetry_per_block, + do_KL, + do_display, + do_save_to_file); + 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) + : out_filename_prefix(std::move(out_filename_prefix_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) +{ + 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(); + 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_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_virtual_axial_crystals_per_block(); + const int virtual_transaxial_crystals + = 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_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_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_v.get_proj_data_info_sptr()->get_scanner_sptr()->get_num_transaxial_blocks_per_bucket(); + const int 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_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_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) + { + 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 + data_fan_sums = DetectorEfficiencies(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 = 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 = 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); + + 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 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, + num_physical_rings, + num_physical_detectors_per_ring); + block_data = BlockData3D(num_axial_blocks, num_transaxial_blocks, num_axial_blocks - 1, num_transaxial_blocks - 1); +} + +void +MLEstimateComponentBasedNormalisation::process() +{ + if (do_display) + { + display(model_fan_data, "model"); + } + + // Initialize the efficiencies, geo data and block data to 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); + apply_block_norm(fan_data, norm_block_data); + if (do_display) + { + display(fan_data, "model*geo*block"); + } + + // 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); + } + + if (do_geo) + { + geo_normalization_iteration(iter_num); // Calculate geo norm iteration + } + if (do_block) + { + block_normalization_iteration(iter_num); // Calculate block norm iteration + } + 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); + 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::get_data_is_processed() const +{ + return data_processed; +} + +DetectorEfficiencies +MLEstimateComponentBasedNormalisation::get_efficiencies() const +{ + if (!this->get_data_is_processed()) + { + error("MLEstimateComponentBasedNormalisation::get_efficiencies: data has not been processed yet"); + } + return norm_efficiencies; +} + +GeoData3D +MLEstimateComponentBasedNormalisation::get_geo_data() const +{ + if (!this->get_data_is_processed()) + { + error("MLEstimateComponentBasedNormalisation::get_geo_data: data has not been processed yet"); + } + if (!do_geo) + { + error("MLEstimateComponentBasedNormalisation::get_geo_data: geo data was not calculated"); + } + return norm_geo_data; +} + +BlockData3D +MLEstimateComponentBasedNormalisation::get_block_data() const +{ + if (!this->get_data_is_processed()) + { + 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; +} + +BinNormalisationPETFromComponents +MLEstimateComponentBasedNormalisation::construct_bin_norm_from_pet_components() const +{ + if (!this->get_data_is_processed()) + { + 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); + bin_norm.set_crystal_efficiencies(norm_efficiencies); + if (do_geo) + { + bin_norm.set_geometric_factors(norm_geo_data); + } + if (do_block) + { + bin_norm.set_block_factors(norm_block_data); + } + 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(const int num_eff_iterations) +{ + data_processed = false; + this->num_eff_iterations = num_eff_iterations; +} + +void +MLEstimateComponentBasedNormalisation::set_num_iterations(const int num_iterations) +{ + data_processed = false; + this->num_iterations = num_iterations; +} + +void +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(const bool do_block) +{ + data_processed = false; + this->do_block = do_block; +} + +void +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(const bool do_kl) +{ + data_processed = false; + do_KL = do_kl; +} + +void +MLEstimateComponentBasedNormalisation::set_write_display_data(const bool do_display) +{ + data_processed = false; + this->do_display = do_display; +} + +void +MLEstimateComponentBasedNormalisation::set_write_intermediates_to_file(const 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) +{ + 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, 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) + 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, 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); + apply_block_norm(fan_data, norm_block_data); + } +} + +void +MLEstimateComponentBasedNormalisation::geo_normalization_iteration(const int iter_num) +{ + + 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) + { + write_geo_data_to_file(iter_num); + } + 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"); + } +} + +void +MLEstimateComponentBasedNormalisation::block_normalization_iteration(const int iter_num) +{ + + 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) + { + write_block_data_to_file(iter_num); + } + 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 +MLEstimateComponentBasedNormalisation::write_efficiencies_to_file(const int iter_num, const int eff_iter_num) const +{ + 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 << norm_efficiencies; +} + +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; +} + +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; +} + +float +MLEstimateComponentBasedNormalisation::compute_threshold_for_KL() +{ + // Set the max found value to -inf + float max_elem = -std::numeric_limits::infinity(); + + 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(); + + 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) + { + 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 diff --git a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx b/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx deleted file mode 100644 index a22ed6b0b9..0000000000 --- a/src/recon_buildblock/ML_estimate_component_based_normalisation.cxx +++ /dev/null @@ -1,308 +0,0 @@ -/* - Copyright (C) 2001- 2008, Hammersmith Imanet Ltd - Copyright (C) 2019-2020, 2022, University College London - Copyright (C) 2016-2017, PETsys Electronics - Copyright (C) 2021, Gefei Chen - This file is part of STIR. - - SPDX-License-Identifier: Apache-2.0 - - See STIR/LICENSE.txt for details -*/ -/*! - - \file - \ingroup recon_buildblock - \brief Implementation of ML_estimate_component_based_normalisation - - \author Kris Thielemans - \author Tahereh Niknejad - \author Gefei Chen - */ -#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 - -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) -{ - - 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; - 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; - } - } - } - - 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"); - - /* { - 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; - } - */ - } - - // 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 - { - 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); - } -#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); - - if (do_geo) - iterate_geo_norm(norm_geo_data, measured_geo_data, fan_data); - -#if 0 - - { // check - for (int a=0; a(*Scanner::get_scanner_from_name("Discovery 690")); + 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(); + + // 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_norm_from_pet_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 the normalization factors, with measured data as uniform 1s and model data as uniform 2s, expect this to be 2.0f + 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 + +USING_NAMESPACE_STIR + +int +main() +{ + Verbosity::set(1); + MLEstimateComponentBasedNormalisationTest tests; + tests.run_tests(); + return tests.main_return_value(); +} diff --git a/src/utilities/apply_normfactors3D.cxx b/src/utilities/apply_normfactors3D.cxx index a235aea2b1..0be120d051 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); diff --git a/src/utilities/find_ML_normfactors3D.cxx b/src/utilities/find_ML_normfactors3D.cxx index 2e78c0cbf8..f7bdcd998c 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" @@ -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());