From a2691acaa7564ba38c3d33258eae79e3d4be2692 Mon Sep 17 00:00:00 2001 From: Evgueni Ovtchinnikov Date: Tue, 16 Apr 2024 13:26:22 +0000 Subject: [PATCH] reimplemented two nema data processing functions in SIRF (C++) --- src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h | 802 +++++++++++---------- src/xSTIR/cSTIR/tests/test7.cpp | 14 +- 2 files changed, 423 insertions(+), 393 deletions(-) diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h index bbf0c257b..87b72ccdd 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h @@ -81,193 +81,11 @@ The actual algorithm is described in > [Online](http://dx.doi.org/10.1109/nssmic.2002.1239610). */ - class ListmodeToSinograms : public stir::LmToProjData { - public: - //! Constructor. - /*! Takes an optional text string argument with - the name of a STIR parameter file defining the conversion options. - If no argument is given, default settings apply except - for the names of input raw data file, template file and - output filename prefix, which must be set by the user by - calling respective methods. - - By default, `store_prompts` is `true` and `store_delayeds` is `false`. - */ - //ListmodeToSinograms(const char* const par) : stir::LmToProjData(par) {} - ListmodeToSinograms(const char* par) : stir::LmToProjData(par) {} - ListmodeToSinograms() : stir::LmToProjData() - { - set_defaults(); - fan_size = -1; - store_prompts = true; - store_delayeds = false; - delayed_increment = 0; - num_iterations = 10; - display_interval = 1; - KL_interval = 1; - save_interval = -1; - //num_events_to_store = -1; - } - void set_input(const STIRListmodeData& lm_data_v) - { - input_filename = "UNKNOWN"; - // call stir::LmToProjData::set_input_data - this->set_input_data(lm_data_v.data()); - exam_info_sptr_.reset(new ExamInfo(lm_data_ptr->get_exam_info())); - proj_data_info_sptr_.reset(lm_data_ptr->get_proj_data_info_sptr()->clone()); - } - void set_input(std::string lm_file) - { - this->set_input(STIRListmodeData(lm_file)); - this->input_filename = lm_file; - } - //! Specifies the prefix for the output file(s), - /*! This will be appended by `_g1f1d0b0.hs`. - */ - void set_output(std::string proj_data_file) - { - output_filename_prefix = proj_data_file; - } - void set_template(std::string proj_data_file) - { - STIRAcquisitionDataInFile acq_data_template(proj_data_file.c_str()); - set_template(acq_data_template); - } - void set_template(const STIRAcquisitionData& acq_data_template) - { - template_proj_data_info_ptr = - acq_data_template.get_proj_data_info_sptr()->create_shared_clone(); - } - void set_time_interval(double start, double stop) - { - std::pair interval(start, stop); - std::vector < std::pair > intervals; - intervals.push_back(interval); - frame_defs = stir::TimeFrameDefinitions(intervals); - do_time_frame = true; - } - int set_flag(const char* flag, bool value) - { - if (sirf::iequals(flag, "store_prompts")) - store_prompts = value; - else if (sirf::iequals(flag, "store_delayeds")) - store_delayeds = value; -#if 0 - else if (sirf::iequals(flag, "do_pre_normalisation")) - do_pre_normalisation = value; - else if (sirf::iequals(flag, "do_time_frame")) - do_time_frame = value; -#endif - else if (sirf::iequals(flag, "interactive")) - interactive = value; - else - return -1; - return 0; - } - bool get_store_prompts() const - { - return store_prompts; - } - bool get_store_delayeds() const - { - return store_delayeds; - } - virtual stir::Succeeded set_up() - { - if (LmToProjData::set_up() == Succeeded::no) - THROW("LmToProjData setup failed"); - fan_size = -1; -#if STIR_VERSION < 060000 - const auto max_fan_size = - lm_data_ptr->get_scanner_ptr()->get_max_num_non_arccorrected_bins(); -#else - const auto max_fan_size = - lm_data_ptr->get_scanner().get_max_num_non_arccorrected_bins(); -#endif - if (fan_size == -1) - fan_size = max_fan_size; - else - fan_size = - std::min(fan_size, max_fan_size); - half_fan_size = fan_size / 2; - fan_size = 2 * half_fan_size + 1; - - exam_info_sptr_->set_time_frame_definitions(frame_defs); - const float h = proj_data_info_sptr_->get_bed_position_horizontal(); - const float v = proj_data_info_sptr_->get_bed_position_vertical(); - stir::shared_ptr temp_proj_data_info_sptr(template_proj_data_info_ptr->clone()); - temp_proj_data_info_sptr->set_bed_position_horizontal(h); - temp_proj_data_info_sptr->set_bed_position_vertical(v); - randoms_sptr.reset(new STIRAcquisitionDataInMemory(exam_info_sptr_, temp_proj_data_info_sptr)); - - return stir::Succeeded::yes; - } - int estimate_randoms(); - void save_randoms() - { - std::string filename = "randoms_f1g1d0b0.hs"; - randoms_sptr->write(filename.c_str()); - } - std::shared_ptr get_output() - { - std::string filename = output_filename_prefix + "_f1g1d0b0.hs"; - return std::shared_ptr - (new STIRAcquisitionDataInFile(filename.c_str())); - } - std::shared_ptr get_randoms_sptr() - { - return randoms_sptr; - } - /// Get the time at which the number of prompts exceeds a certain threshold. - /// Returns -1 if not found. - float get_time_at_which_num_prompts_exceeds_threshold(const unsigned long threshold) const; - - void sinograms_and_randoms_from_listmode( - STIRListmodeData& lm_data, double start, double stop, - STIRAcquisitionData& acq_data_template, - std::shared_ptr& sinograms_sptr, - std::shared_ptr& randoms_sptr) - { - ListmodeToSinograms converter; - converter.set_input(lm_data); - converter.set_output("sinograms"); - converter.set_template(acq_data_template); - converter.set_time_interval(start, stop); - converter.set_up(); - converter.process_data(); - sinograms_sptr = converter.get_output(); - converter.estimate_randoms(); - randoms_sptr = converter.get_randoms_sptr(); - } - - protected: - // variables for ML estimation of singles/randoms - int fan_size; - int half_fan_size; - int max_ring_diff_for_fansums; - int num_iterations; - int display_interval; - int KL_interval; - int save_interval; - stir::shared_ptr exam_info_sptr_; - stir::shared_ptr proj_data_info_sptr_; - stir::shared_ptr > > fan_sums_sptr; - stir::shared_ptr det_eff_sptr; - std::shared_ptr randoms_sptr; - void compute_fan_sums_(bool prompt_fansum = false); - int compute_singles_(); -// void estimate_randoms_(); - static unsigned long compute_num_bins_(const int num_rings, - const int num_detectors_per_ring, - const int max_ring_diff, const int half_fan_size); - }; - /*! \ingroup PET \brief Class for PET scanner detector efficiencies model. */ - class PETAcquisitionSensitivityModel { public: PETAcquisitionSensitivityModel() {} @@ -552,39 +370,420 @@ The actual algorithm is described in //shared_ptr sptr_normalisation_; }; - /*! - \ingroup PET - - \brief Class for simulating the scatter contribution to PET data. + /*! + \ingroup PET + \brief Ray tracing matrix implementation of the PET acquisition model. - This class uses the STIR Single Scatter simulation, taking as input an - activity and attenuation image, and a acquisition data template. + In this implementation \e x and \e y are essentially vectors and \e G + a matrix. Each row of \e G corresponds to a line-of-response (LOR) + between two detectors (there may be more than one line for each pair). + The only non-zero elements of each row are those corresponding to + voxels through which LOR passes, so the matrix is very sparse. + Furthermore, owing to symmetries, many rows have the same values only + in different order, and thus only one set of values needs to be computed + and stored (see STIR documentation for details). + */ - WARNING: Currently this class does not use the low-resolution sampling - mechanism of STIR. This means that if you give it a full resolution acq_data, - you will likely run out of memory and/or time. - */ - class PETSingleScatterSimulator : public stir::SingleScatterSimulation - { - public: - //! Default constructor - PETSingleScatterSimulator() : stir::SingleScatterSimulation() - {} - //! Overloaded constructor which takes the parameter file - PETSingleScatterSimulator(std::string filename) : - stir::SingleScatterSimulation(filename) - {} + class PETAcquisitionModelUsingMatrix : public PETAcquisitionModel { + public: + PETAcquisitionModelUsingMatrix() + { + this->sptr_projectors_.reset(new ProjectorPairUsingMatrix); + } + void set_matrix(stir::shared_ptr sptr_matrix) + { + sptr_matrix_ = sptr_matrix; + ((ProjectorPairUsingMatrix*)this->sptr_projectors_.get())-> + set_proj_matrix_sptr(sptr_matrix); + } + stir::shared_ptr matrix_sptr() + { + return sptr_matrix_; + //return ((ProjectorPairUsingMatrix*)this->sptr_projectors_.get())-> + // get_proj_matrix_sptr(); + } + virtual void set_up( + std::shared_ptr sptr_acq, + std::shared_ptr sptr_image) + { + if (!sptr_matrix_.get()) + THROW("PETAcquisitionModelUsingMatrix setup failed - matrix not set"); + PETAcquisitionModel::set_up(sptr_acq, sptr_image); + } + + //! Enables or disables the caching mechanism. + void enable_cache(bool v = true) + { + sptr_matrix_->enable_cache(v); + } - void set_up(std::shared_ptr sptr_acq_template, - std::shared_ptr sptr_act_image_template) - { - this->sptr_acq_template_ = sptr_acq_template; + private: + stir::shared_ptr sptr_matrix_; + }; - stir::SingleScatterSimulation::set_template_proj_data_info( - *sptr_acq_template_->get_proj_data_info_sptr()); - stir::SingleScatterSimulation::set_exam_info( - *sptr_acq_template_->get_exam_info_sptr()); - // check if attenuation image is set + class PETAcquisitionModelUsingRayTracingMatrix : + public PETAcquisitionModelUsingMatrix { + public: + PETAcquisitionModelUsingRayTracingMatrix(int num_LORs = 2) : + PETAcquisitionModelUsingMatrix() + { + stir::shared_ptr matrix_sptr(new RayTracingMatrix); + matrix_sptr->set_num_tangential_LORs(num_LORs); + set_matrix(matrix_sptr); + } + void set_num_tangential_LORs(int num_LORs) + { + //RayTracingMatrix& matrix = (RayTracingMatrix&)*matrix_sptr(); + auto matrix = dynamic_cast(*matrix_sptr()); + //std::cout << matrix.get_num_tangential_LORs() << '\n'; + matrix.set_num_tangential_LORs(num_LORs); + //std::cout << get_num_tangential_LORs() << '\n'; + } + //!@ + int get_num_tangential_LORs() + { + auto matrix = dynamic_cast(*matrix_sptr()); + return matrix.get_num_tangential_LORs(); + } + //! Enables or disables using a circular axial FOV (vs rectangular) + void set_restrict_to_cylindrical_FOV(bool v = true) + { + auto matrix = dynamic_cast(*matrix_sptr()); + matrix.set_restrict_to_cylindrical_FOV(v); + } + //! \name Which symmetries will be used + //!@{ + //bool get_do_symmetry_90degrees_min_phi() const; + void set_do_symmetry_90degrees_min_phi(bool v = true) + { + auto matrix = dynamic_cast(*matrix_sptr()); + matrix.set_do_symmetry_90degrees_min_phi(v); + } + //bool get_do_symmetry_180degrees_min_phi() const; + void set_do_symmetry_180degrees_min_phi(bool v = true) + { + auto matrix = dynamic_cast(*matrix_sptr()); + matrix.set_do_symmetry_180degrees_min_phi(v); + } + //bool get_do_symmetry_swap_segment() const; + void set_do_symmetry_swap_segment(bool v = true) + { + auto matrix = dynamic_cast(*matrix_sptr()); + matrix.set_do_symmetry_swap_segment(v); + } + //bool get_do_symmetry_swap_s() const; + void set_do_symmetry_swap_s(bool v = true) + { + auto matrix = dynamic_cast(*matrix_sptr()); + matrix.set_do_symmetry_swap_s(v); + } + //bool get_do_symmetry_shift_z() const; + void set_do_symmetry_shift_z(bool v = true) + { + auto matrix = dynamic_cast(*matrix_sptr()); + matrix.set_do_symmetry_shift_z(v); + } + }; + + typedef PETAcquisitionModel AcqMod3DF; + typedef PETAcquisitionModelUsingMatrix AcqModUsingMatrix3DF; + typedef std::shared_ptr sptrAcqMod3DF; + +#ifdef STIR_WITH_NiftyPET_PROJECTOR + /*! + \ingroup PET + \brief NiftyPET implementation of the PET acquisition model. + */ + + class PETAcquisitionModelUsingNiftyPET : public PETAcquisitionModel { + public: + PETAcquisitionModelUsingNiftyPET() + { + _NiftyPET_projector_pair_sptr.reset(new ProjectorPairUsingNiftyPET); + this->sptr_projectors_ = _NiftyPET_projector_pair_sptr; + // Set verbosity to 0 by default + _NiftyPET_projector_pair_sptr->set_verbosity(0); + } + void set_cuda_verbosity(const bool verbosity) const + { + _NiftyPET_projector_pair_sptr->set_verbosity(verbosity); + } + void set_use_truncation(const bool use_truncation) const + { + _NiftyPET_projector_pair_sptr->set_use_truncation(use_truncation); + } + protected: + stir::shared_ptr _NiftyPET_projector_pair_sptr; + }; + typedef PETAcquisitionModelUsingNiftyPET AcqModUsingNiftyPET3DF; +#endif + +#ifdef STIR_WITH_Parallelproj_PROJECTOR + /*! + \ingroup PET + \brief Parallelproj implementation of the PET acquisition model + (see https://github.com/gschramm/parallelproj). + */ + class PETAcquisitionModelUsingParallelproj : public PETAcquisitionModel { + public: + PETAcquisitionModelUsingParallelproj() + { + this->sptr_projectors_.reset(new ProjectorByBinPairUsingParallelproj); + } + }; + typedef PETAcquisitionModelUsingParallelproj AcqModUsingParallelproj; +#endif + + /*! + \ingroup PET + \brief Attenuation model. + + */ + + class PETAttenuationModel : public PETAcquisitionSensitivityModel { + public: + PETAttenuationModel(STIRImageData& id, PETAcquisitionModel& am); + // multiply by bin efficiencies + virtual void unnormalise(STIRAcquisitionData& ad) const; + // divide by bin efficiencies + virtual void normalise(STIRAcquisitionData& ad) const; + protected: + stir::shared_ptr sptr_forw_projector_; + }; + + + class ListmodeToSinograms : public stir::LmToProjData { + public: + //! Constructor. + /*! Takes an optional text string argument with + the name of a STIR parameter file defining the conversion options. + If no argument is given, default settings apply except + for the names of input raw data file, template file and + output filename prefix, which must be set by the user by + calling respective methods. + + By default, `store_prompts` is `true` and `store_delayeds` is `false`. + */ + //ListmodeToSinograms(const char* const par) : stir::LmToProjData(par) {} + ListmodeToSinograms(const char* par) : stir::LmToProjData(par) {} + ListmodeToSinograms() : stir::LmToProjData() + { + set_defaults(); + fan_size = -1; + store_prompts = true; + store_delayeds = false; + delayed_increment = 0; + num_iterations = 10; + display_interval = 1; + KL_interval = 1; + save_interval = -1; + //num_events_to_store = -1; + } + void set_input(const STIRListmodeData& lm_data_v) + { + input_filename = "UNKNOWN"; + // call stir::LmToProjData::set_input_data + this->set_input_data(lm_data_v.data()); + exam_info_sptr_.reset(new ExamInfo(lm_data_ptr->get_exam_info())); + proj_data_info_sptr_.reset(lm_data_ptr->get_proj_data_info_sptr()->clone()); + } + void set_input(std::string lm_file) + { + this->set_input(STIRListmodeData(lm_file)); + this->input_filename = lm_file; + } + //! Specifies the prefix for the output file(s), + /*! This will be appended by `_g1f1d0b0.hs`. + */ + void set_output(std::string proj_data_file) + { + output_filename_prefix = proj_data_file; + } + void set_template(std::string proj_data_file) + { + STIRAcquisitionDataInFile acq_data_template(proj_data_file.c_str()); + set_template(acq_data_template); + } + void set_template(const STIRAcquisitionData& acq_data_template) + { + template_proj_data_info_ptr = + acq_data_template.get_proj_data_info_sptr()->create_shared_clone(); + } + void set_time_interval(double start, double stop) + { + std::pair interval(start, stop); + std::vector < std::pair > intervals; + intervals.push_back(interval); + frame_defs = stir::TimeFrameDefinitions(intervals); + do_time_frame = true; + } + int set_flag(const char* flag, bool value) + { + if (sirf::iequals(flag, "store_prompts")) + store_prompts = value; + else if (sirf::iequals(flag, "store_delayeds")) + store_delayeds = value; +#if 0 + else if (sirf::iequals(flag, "do_pre_normalisation")) + do_pre_normalisation = value; + else if (sirf::iequals(flag, "do_time_frame")) + do_time_frame = value; +#endif + else if (sirf::iequals(flag, "interactive")) + interactive = value; + else + return -1; + return 0; + } + bool get_store_prompts() const + { + return store_prompts; + } + bool get_store_delayeds() const + { + return store_delayeds; + } + virtual stir::Succeeded set_up() + { + if (LmToProjData::set_up() == Succeeded::no) + THROW("LmToProjData setup failed"); + fan_size = -1; +#if STIR_VERSION < 060000 + const auto max_fan_size = + lm_data_ptr->get_scanner_ptr()->get_max_num_non_arccorrected_bins(); +#else + const auto max_fan_size = + lm_data_ptr->get_scanner().get_max_num_non_arccorrected_bins(); +#endif + if (fan_size == -1) + fan_size = max_fan_size; + else + fan_size = + std::min(fan_size, max_fan_size); + half_fan_size = fan_size / 2; + fan_size = 2 * half_fan_size + 1; + + exam_info_sptr_->set_time_frame_definitions(frame_defs); + const float h = proj_data_info_sptr_->get_bed_position_horizontal(); + const float v = proj_data_info_sptr_->get_bed_position_vertical(); + stir::shared_ptr temp_proj_data_info_sptr(template_proj_data_info_ptr->clone()); + temp_proj_data_info_sptr->set_bed_position_horizontal(h); + temp_proj_data_info_sptr->set_bed_position_vertical(v); + randoms_sptr.reset(new STIRAcquisitionDataInMemory(exam_info_sptr_, temp_proj_data_info_sptr)); + + return stir::Succeeded::yes; + } + int estimate_randoms(); + void save_randoms() + { + std::string filename = "randoms_f1g1d0b0.hs"; + randoms_sptr->write(filename.c_str()); + } + std::shared_ptr get_output() + { + std::string filename = output_filename_prefix + "_f1g1d0b0.hs"; + return std::shared_ptr + (new STIRAcquisitionDataInFile(filename.c_str())); + } + std::shared_ptr get_randoms_sptr() + { + return randoms_sptr; + } + /// Get the time at which the number of prompts exceeds a certain threshold. + /// Returns -1 if not found. + float get_time_at_which_num_prompts_exceeds_threshold(const unsigned long threshold) const; + + void sinograms_and_randoms_from_listmode( + const STIRListmodeData& lm_data, double start, double stop, + const STIRAcquisitionData& acq_data_template, + std::shared_ptr& sinograms_sptr, + std::shared_ptr& randoms_sptr) + { + ListmodeToSinograms converter; + converter.set_input(lm_data); + converter.set_output("sinograms"); + converter.set_template(acq_data_template); + converter.set_time_interval(start, stop); + converter.set_up(); + converter.process_data(); + sinograms_sptr = converter.get_output(); + converter.estimate_randoms(); + randoms_sptr = converter.get_randoms_sptr(); + } + void acquisition_sensitivity_from_attenuation( + std::shared_ptr acq_templ_sptr, + std::shared_ptr attn_image_sptr, + std::shared_ptr& acq_sens_sptr) + { + PETAcquisitionModelUsingRayTracingMatrix acq_mod; + acq_mod.set_up(acq_templ_sptr, attn_image_sptr); + std::shared_ptr + asm_sptr(new PETAttenuationModel(*attn_image_sptr, acq_mod)); + PETAttenuationModel acq_sens_mod = *asm_sptr; + acq_sens_mod.set_up(acq_templ_sptr->get_exam_info_sptr(), + acq_templ_sptr->get_proj_data_info_sptr()->create_shared_clone()); + acq_mod.set_asm(asm_sptr); + std::cout << "applying attenuation (please wait, may take a while)...\n"; + acq_sens_sptr = acq_templ_sptr->clone(); + acq_sens_sptr->fill(1.0); + acq_sens_mod.unnormalise(*acq_sens_sptr); + } + + protected: + // variables for ML estimation of singles/randoms + int fan_size; + int half_fan_size; + int max_ring_diff_for_fansums; + int num_iterations; + int display_interval; + int KL_interval; + int save_interval; + stir::shared_ptr exam_info_sptr_; + stir::shared_ptr proj_data_info_sptr_; + stir::shared_ptr > > fan_sums_sptr; + stir::shared_ptr det_eff_sptr; + std::shared_ptr randoms_sptr; + void compute_fan_sums_(bool prompt_fansum = false); + int compute_singles_(); +// void estimate_randoms_(); + static unsigned long compute_num_bins_(const int num_rings, + const int num_detectors_per_ring, + const int max_ring_diff, const int half_fan_size); + }; + + /*! + \ingroup PET + + \brief Class for simulating the scatter contribution to PET data. + + This class uses the STIR Single Scatter simulation, taking as input an + activity and attenuation image, and a acquisition data template. + + WARNING: Currently this class does not use the low-resolution sampling + mechanism of STIR. This means that if you give it a full resolution acq_data, + you will likely run out of memory and/or time. + */ + class PETSingleScatterSimulator : public stir::SingleScatterSimulation + { + public: + //! Default constructor + PETSingleScatterSimulator() : stir::SingleScatterSimulation() + {} + //! Overloaded constructor which takes the parameter file + PETSingleScatterSimulator(std::string filename) : + stir::SingleScatterSimulation(filename) + {} + + void set_up(std::shared_ptr sptr_acq_template, + std::shared_ptr sptr_act_image_template) + { + this->sptr_acq_template_ = sptr_acq_template; + + stir::SingleScatterSimulation::set_template_proj_data_info( + *sptr_acq_template_->get_proj_data_info_sptr()); + stir::SingleScatterSimulation::set_exam_info( + *sptr_acq_template_->get_exam_info_sptr()); + // check if attenuation image is set try { auto& tmp = stir::SingleScatterSimulation::get_attenuation_image(); @@ -829,187 +1028,6 @@ The actual algorithm is described in } }; - /*! - \ingroup PET - \brief Ray tracing matrix implementation of the PET acquisition model. - - In this implementation \e x and \e y are essentially vectors and \e G - a matrix. Each row of \e G corresponds to a line-of-response (LOR) - between two detectors (there may be more than one line for each pair). - The only non-zero elements of each row are those corresponding to - voxels through which LOR passes, so the matrix is very sparse. - Furthermore, owing to symmetries, many rows have the same values only - in different order, and thus only one set of values needs to be computed - and stored (see STIR documentation for details). - */ - - class PETAcquisitionModelUsingMatrix : public PETAcquisitionModel { - public: - PETAcquisitionModelUsingMatrix() - { - this->sptr_projectors_.reset(new ProjectorPairUsingMatrix); - } - void set_matrix(stir::shared_ptr sptr_matrix) - { - sptr_matrix_ = sptr_matrix; - ((ProjectorPairUsingMatrix*)this->sptr_projectors_.get())-> - set_proj_matrix_sptr(sptr_matrix); - } - stir::shared_ptr matrix_sptr() - { - return sptr_matrix_; - //return ((ProjectorPairUsingMatrix*)this->sptr_projectors_.get())-> - // get_proj_matrix_sptr(); - } - virtual void set_up( - std::shared_ptr sptr_acq, - std::shared_ptr sptr_image) - { - if (!sptr_matrix_.get()) - THROW("PETAcquisitionModelUsingMatrix setup failed - matrix not set"); - PETAcquisitionModel::set_up(sptr_acq, sptr_image); - } - - //! Enables or disables the caching mechanism. - void enable_cache(bool v = true) - { - sptr_matrix_->enable_cache(v); - } - - private: - stir::shared_ptr sptr_matrix_; - }; - - class PETAcquisitionModelUsingRayTracingMatrix : - public PETAcquisitionModelUsingMatrix { - public: - PETAcquisitionModelUsingRayTracingMatrix(int num_LORs = 2) : - PETAcquisitionModelUsingMatrix() - { - stir::shared_ptr matrix_sptr(new RayTracingMatrix); - matrix_sptr->set_num_tangential_LORs(num_LORs); - set_matrix(matrix_sptr); - } - void set_num_tangential_LORs(int num_LORs) - { - //RayTracingMatrix& matrix = (RayTracingMatrix&)*matrix_sptr(); - auto matrix = dynamic_cast(*matrix_sptr()); - //std::cout << matrix.get_num_tangential_LORs() << '\n'; - matrix.set_num_tangential_LORs(num_LORs); - //std::cout << get_num_tangential_LORs() << '\n'; - } - //!@ - int get_num_tangential_LORs() - { - auto matrix = dynamic_cast(*matrix_sptr()); - return matrix.get_num_tangential_LORs(); - } - //! Enables or disables using a circular axial FOV (vs rectangular) - void set_restrict_to_cylindrical_FOV(bool v = true) - { - auto matrix = dynamic_cast(*matrix_sptr()); - matrix.set_restrict_to_cylindrical_FOV(v); - } - //! \name Which symmetries will be used - //!@{ - //bool get_do_symmetry_90degrees_min_phi() const; - void set_do_symmetry_90degrees_min_phi(bool v = true) - { - auto matrix = dynamic_cast(*matrix_sptr()); - matrix.set_do_symmetry_90degrees_min_phi(v); - } - //bool get_do_symmetry_180degrees_min_phi() const; - void set_do_symmetry_180degrees_min_phi(bool v = true) - { - auto matrix = dynamic_cast(*matrix_sptr()); - matrix.set_do_symmetry_180degrees_min_phi(v); - } - //bool get_do_symmetry_swap_segment() const; - void set_do_symmetry_swap_segment(bool v = true) - { - auto matrix = dynamic_cast(*matrix_sptr()); - matrix.set_do_symmetry_swap_segment(v); - } - //bool get_do_symmetry_swap_s() const; - void set_do_symmetry_swap_s(bool v = true) - { - auto matrix = dynamic_cast(*matrix_sptr()); - matrix.set_do_symmetry_swap_s(v); - } - //bool get_do_symmetry_shift_z() const; - void set_do_symmetry_shift_z(bool v = true) - { - auto matrix = dynamic_cast(*matrix_sptr()); - matrix.set_do_symmetry_shift_z(v); - } - }; - - typedef PETAcquisitionModel AcqMod3DF; - typedef PETAcquisitionModelUsingMatrix AcqModUsingMatrix3DF; - typedef std::shared_ptr sptrAcqMod3DF; - -#ifdef STIR_WITH_NiftyPET_PROJECTOR - /*! - \ingroup PET - \brief NiftyPET implementation of the PET acquisition model. - */ - - class PETAcquisitionModelUsingNiftyPET : public PETAcquisitionModel { - public: - PETAcquisitionModelUsingNiftyPET() - { - _NiftyPET_projector_pair_sptr.reset(new ProjectorPairUsingNiftyPET); - this->sptr_projectors_ = _NiftyPET_projector_pair_sptr; - // Set verbosity to 0 by default - _NiftyPET_projector_pair_sptr->set_verbosity(0); - } - void set_cuda_verbosity(const bool verbosity) const - { - _NiftyPET_projector_pair_sptr->set_verbosity(verbosity); - } - void set_use_truncation(const bool use_truncation) const - { - _NiftyPET_projector_pair_sptr->set_use_truncation(use_truncation); - } - protected: - stir::shared_ptr _NiftyPET_projector_pair_sptr; - }; - typedef PETAcquisitionModelUsingNiftyPET AcqModUsingNiftyPET3DF; -#endif - -#ifdef STIR_WITH_Parallelproj_PROJECTOR - /*! - \ingroup PET - \brief Parallelproj implementation of the PET acquisition model - (see https://github.com/gschramm/parallelproj). - */ - class PETAcquisitionModelUsingParallelproj : public PETAcquisitionModel { - public: - PETAcquisitionModelUsingParallelproj() - { - this->sptr_projectors_.reset(new ProjectorByBinPairUsingParallelproj); - } - }; - typedef PETAcquisitionModelUsingParallelproj AcqModUsingParallelproj; -#endif - - /*! - \ingroup PET - \brief Attenuation model. - - */ - - class PETAttenuationModel : public PETAcquisitionSensitivityModel { - public: - PETAttenuationModel(STIRImageData& id, PETAcquisitionModel& am); - // multiply by bin efficiencies - virtual void unnormalise(STIRAcquisitionData& ad) const; - // divide by bin efficiencies - virtual void normalise(STIRAcquisitionData& ad) const; - protected: - stir::shared_ptr sptr_forw_projector_; - }; - /*! \ingroup PET \brief Accessor classes. diff --git a/src/xSTIR/cSTIR/tests/test7.cpp b/src/xSTIR/cSTIR/tests/test7.cpp index 230155d81..a3480879c 100644 --- a/src/xSTIR/cSTIR/tests/test7.cpp +++ b/src/xSTIR/cSTIR/tests/test7.cpp @@ -39,6 +39,8 @@ limitations under the License. #include "sirf/common/iequals.h" #include "sirf/common/utilities.h" +#include "object.h" + using namespace stir; using namespace ecat; using namespace sirf; @@ -53,13 +55,20 @@ int main() return 1; } std::string path = append_path(SIRF_data_path, "mMR"); + std::cout << path << '\n'; std::string f_listmode = append_path(path, "list.l.hdr"); std::string f_template = append_path(path, "mMR_template_span11_small.hs"); + std::string f_mu_map = append_path(path, "mu_map.hv"); STIRAcquisitionDataInMemory::set_as_template(); - STIRAcquisitionDataInFile acq_data_template(f_template.c_str()); + //STIRAcquisitionDataInFile acq_data_template(f_template.c_str()); + CREATE_OBJECT(STIRAcquisitionData, STIRAcquisitionDataInFile, + acq_data_template, templ_sptr, f_template.c_str()); STIRListmodeData lm_data(f_listmode); + //STIRImageData mu_map(f_mu_map); + CREATE_OBJ(STIRImageData, mu_map, mu_map_sptr, f_mu_map); + std::shared_ptr acq_sens_sptr; ListmodeToSinograms converter; std::shared_ptr sinograms_sptr; @@ -69,6 +78,9 @@ int main() std::cout << sinograms_sptr->norm() << '\n'; std::cout << randoms_sptr->norm() << '\n'; + + converter.acquisition_sensitivity_from_attenuation(templ_sptr, mu_map_sptr, acq_sens_sptr); + std::cout << acq_sens_sptr->norm() << '\n'; std::cout << "done with test7.cpp...\n"; return 0;