diff --git a/CHANGES.md b/CHANGES.md index 04b829a73..979073c2a 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -21,6 +21,12 @@ - provided prior and objective function objects with methods for computing the product of the Hessian and a vector. +* PET: + - incorporated into SIRF data processing utilities from SyneRBI-Challenge. + - Added method modality() to AcquisitionData classes. + +## v3.7.0 + * CMake/building: - add `DISABLE_Gadgetron_TOOLBOXES` option (defaults to `OFF`) to be able to cope with compilation problems with older Gadgetron versions. diff --git a/examples/Python/PET/acquisition_model.py b/examples/Python/PET/acquisition_model.py index 91d1ee0b8..2486185ab 100644 --- a/examples/Python/PET/acquisition_model.py +++ b/examples/Python/PET/acquisition_model.py @@ -169,6 +169,7 @@ def main(): # simulated_data = acq_model.forward(image, 0, num_subsets) if output_file is not None: simulated_data.write(output_file) + print('simulated data modality: %s' % simulated_data.modality) if show_plot: # show simulated acquisition data diff --git a/examples/Python/SPECT/acquisition_model.py b/examples/Python/SPECT/acquisition_model.py index 5651c3954..af8f35eb4 100644 --- a/examples/Python/SPECT/acquisition_model.py +++ b/examples/Python/SPECT/acquisition_model.py @@ -83,6 +83,7 @@ def main(): # raw data to be used as a template for the acquisition model acq_template = sirf.STIR.AcquisitionData(raw_data_file) + print('template modality: %s' % acq_template.modality) # create image with suitable sizes image = acq_template.create_uniform_image() @@ -109,6 +110,7 @@ def main(): # simulated_data = acq_model.forward(image, 0, 4) if output_file is not None: simulated_data.write(output_file) + print('simulated data modality: %s' % simulated_data.modality) # show simulated acquisition data simulated_data_as_array = simulated_data.as_array() diff --git a/src/iUtilities/include/sirf/iUtilities/DataHandle.h b/src/iUtilities/include/sirf/iUtilities/DataHandle.h index 0b6c65808..5b2215fa3 100644 --- a/src/iUtilities/include/sirf/iUtilities/DataHandle.h +++ b/src/iUtilities/include/sirf/iUtilities/DataHandle.h @@ -279,6 +279,23 @@ getObjectSptrFromHandle(const void* h, std::shared_ptr& sptr) { sptr = *ptr_sptr; } +template +void +setHandleObjectSptr(void* h, std::shared_ptr& sptr) { + auto handle = reinterpret_cast*>(h); + //ObjectHandle* handle = (ObjectHandle*)h; +#if defined(USE_BOOST) + if (handle->uses_boost_sptr()) + THROW("cannot cast boost::shared_ptr to std::shared_ptr"); +#endif + void* ptr = handle->data(); + if (ptr == 0) + THROW("zero data pointer cannot be dereferenced"); + CAST_PTR(std::shared_ptr, ptr_sptr, ptr); + //delete ptr_sptr; + *ptr_sptr = sptr; +} + #if defined(USE_BOOST) template void diff --git a/src/xSTIR/cSTIR/cstir.cpp b/src/xSTIR/cSTIR/cstir.cpp index 29680cb34..c74d6126e 100644 --- a/src/xSTIR/cSTIR/cstir.cpp +++ b/src/xSTIR/cSTIR/cstir.cpp @@ -39,6 +39,8 @@ using namespace sirf; #define NEW_OBJECT_HANDLE(T) new ObjectHandle(std::shared_ptr(new T)) #define SPTR_FROM_HANDLE(Object, X, H) \ std::shared_ptr X; getObjectSptrFromHandle(H, X); +#define HANDLE_FROM_SPTR(Object, X, H) \ + setHandleObjectSptr(H, X); static void* unknownObject(const char* obj, const char* name, const char* file, int line) @@ -490,7 +492,43 @@ void* cSTIR_convertListmodeToSinograms(void* ptr) } CATCH; } +/* +extern "C" +void* cSTIR_promptsFromListmode(void* ptr_lm2s, void* ptr_lmdata, + const float start, const float stop, + void* ptr_templ, void* ptr_sino, const char* prefix) +{ + try { + ListmodeToSinograms& lm2s = objectFromHandle(ptr_lm2s); + STIRListmodeData& lm_data = objectFromHandle(ptr_lmdata); + STIRAcquisitionData& templ = objectFromHandle(ptr_templ); + SPTR_FROM_HANDLE(STIRAcquisitionData, sptr_sino, ptr_sino); + lm2s.prompts_from_listmode(lm_data, start, stop, templ, sptr_sino, prefix); + HANDLE_FROM_SPTR(STIRAcquisitionData, sptr_sino, ptr_sino); + return new DataHandle; + } + CATCH; +} +extern "C" +void* cSTIR_promptsAndRandomsFromListmode(void* ptr_lm2s, void* ptr_lmdata, + const float start, const float stop, + void* ptr_templ, void* ptr_sino, void* ptr_rand, const char* prefix) +{ + try { + ListmodeToSinograms& lm2s = objectFromHandle(ptr_lm2s); + STIRListmodeData& lm_data = objectFromHandle(ptr_lmdata); + STIRAcquisitionData& templ = objectFromHandle(ptr_templ); + SPTR_FROM_HANDLE(STIRAcquisitionData, sptr_sino, ptr_sino); + SPTR_FROM_HANDLE(STIRAcquisitionData, sptr_rand, ptr_rand); + lm2s.prompts_and_randoms_from_listmode(lm_data, start, stop, templ, sptr_sino, sptr_rand, prefix); + HANDLE_FROM_SPTR(STIRAcquisitionData, sptr_sino, ptr_sino); + HANDLE_FROM_SPTR(STIRAcquisitionData, sptr_rand, ptr_rand); + return new DataHandle; + } + CATCH; +} +*/ extern "C" void* cSTIR_scatterSimulatorFwd (void* ptr_am, void* ptr_im) @@ -641,6 +679,23 @@ void* cSTIR_createPETAttenuationModel(const void* ptr_img, const void* ptr_am) CATCH; } +extern "C" +void* cSTIR_computeACF(const void* ptr_sino, + const void* ptr_att, void* ptr_af, void* ptr_acf) +{ + try { + STIRAcquisitionData& sino = objectFromHandle(ptr_sino); + PETAttenuationModel& att = objectFromHandle(ptr_att); + SPTR_FROM_HANDLE(STIRAcquisitionData, sptr_af, ptr_af); + SPTR_FROM_HANDLE(STIRAcquisitionData, sptr_acf, ptr_acf); + PETAttenuationModel::compute_ac_factors(sino, att, sptr_af, sptr_acf); + HANDLE_FROM_SPTR(STIRAcquisitionData, sptr_af, ptr_af); + HANDLE_FROM_SPTR(STIRAcquisitionData, sptr_acf, ptr_acf); + return (void*) new DataHandle; + } + CATCH; +} + extern "C" void* cSTIR_chainPETAcquisitionSensitivityModels (const void* ptr_first, const void* ptr_second) diff --git a/src/xSTIR/cSTIR/cstir_p.cpp b/src/xSTIR/cSTIR/cstir_p.cpp index af5859855..bb145db55 100644 --- a/src/xSTIR/cSTIR/cstir_p.cpp +++ b/src/xSTIR/cSTIR/cstir_p.cpp @@ -103,6 +103,8 @@ sirf::cSTIR_AcquisitionDataParameter(void* hp, const char* name) STIRAcquisitionData& ad = objectFromHandle(hp); if (sirf::iequals(name, "tof_mash_factor")) return dataHandle(ad.get_tof_mash_factor()); + if (sirf::iequals(name, "modality")) + return charDataHandleFromCharData(ad.modality().c_str()); else return parameterNotFound(name, __FILE__, __LINE__); } diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h index 7433da844..d8498ce2a 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h @@ -63,6 +63,14 @@ extern "C" { (void* ptr_lm2s, const char* flag, int v); void* cSTIR_setupListmodeToSinogramsConverter(void* ptr); void* cSTIR_convertListmodeToSinograms(void* ptr); + /* + void* cSTIR_promptsFromListmode(void* ptr_lm2s, void* ptr_lmdata, + const float start, const float stop, + void* ptr_templ, void* ptr_sino, const char* prefix); + void* cSTIR_promptsAndRandomsFromListmode(void* ptr_lm2s, void* ptr_lmdata, + const float start, const float stop, + void* ptr_templ, void* ptr_sino, void* ptr_rand, const char* prefix); + */ void* cSTIR_computeRandoms(void* ptr); void* cSTIR_lm_num_prompts_exceeds_threshold(void* ptr, const float threshold); void* cSTIR_objFunListModeSetInterval(void* ptr_f, size_t ptr_data); @@ -76,6 +84,8 @@ extern "C" { (const void* ptr_src, const char* src); void* cSTIR_createPETAttenuationModel (const void* ptr_img, const void* ptr_am); + void* cSTIR_computeACF + (const void* ptr_sino, const void* ptr_att, void* ptr_acf, void* ptr_iacf); void* cSTIR_chainPETAcquisitionSensitivityModels (const void* ptr_first, const void* ptr_second); void* cSTIR_setupAcquisitionSensitivityModel(void* ptr_sm, void* ptr_ad); @@ -92,7 +102,7 @@ extern "C" { int subset_num, int num_subsets); void* cSTIR_acquisitionModelBwdReplace(void* ptr_am, void* ptr_ad, int subset_num, int num_subsets, void* ptr_im); - void* cSTIR_get_MatrixInfo(void* ptr); + void* cSTIR_get_MatrixInfo(void* ptr); // Acquisition Model Matrix void* cSTIR_setupSPECTUBMatrix(const void* h_smx, const void* h_acq, const void* h_img); diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h index c83ec7c38..00efb33cf 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h @@ -480,6 +480,11 @@ namespace sirf { { return data()->get_proj_data_info_sptr(); } + std::string modality() const + { + const ExamInfo& ex_info = *get_exam_info_sptr(); + return ex_info.imaging_modality.get_name(); + } // ProjData casts operator stir::ProjData&() { return *data(); } diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h index f3fd7f715..80ddb3e6d 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h @@ -81,175 +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 = output_filename_prefix + "_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; - - 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() {} @@ -267,6 +103,11 @@ The actual algorithm is described in void set_up(const stir::shared_ptr& exam_info_sptr, const stir::shared_ptr&); + void set_up(const STIRAcquisitionData& ad) + { + set_up(ad.get_exam_info_sptr(), ad.get_proj_data_info_sptr()->create_shared_clone()); + } + // multiply by bin efficiencies virtual void unnormalise(STIRAcquisitionData& ad) const; // divide by bin efficiencies @@ -534,52 +375,450 @@ The actual algorithm is described in //shared_ptr sptr_normalisation_; }; - /*! - \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) - {} + /*! + \ingroup PET + \brief Ray tracing matrix implementation of the PET acquisition model. - void set_up(std::shared_ptr sptr_acq_template, - std::shared_ptr sptr_act_image_template) - { - this->sptr_acq_template_ = sptr_acq_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). + */ - 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(); - } - catch (...) - { - THROW("Fatal error in PETSingleScatterSimulator::set_up: attenuation_image has not been set"); - } - this->set_activity_image_sptr(sptr_act_image_template); + 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); + } - if (stir::SingleScatterSimulation::set_up() == Succeeded::no) - THROW("Fatal error in PETSingleScatterSimulator::set_up() failed."); - } + 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 (here attenuation factors), i.e. attenuate data in \a ad + virtual void unnormalise(STIRAcquisitionData& ad) const; + // divide by bin efficiencies (here attenuation factors), i.e. correct data in \a ad for attenuatio + virtual void normalise(STIRAcquisitionData& ad) const; + /*! Convenience function computing attenuation factor using unnormalise + and its inverse (attenuation correction factor) using STIRImageData::inv + */ + static void compute_ac_factors( + // input arguments + const STIRAcquisitionData& acq_templ, + const PETAttenuationModel& acq_sens_mod, + // output arguments + std::shared_ptr& af_sptr, + std::shared_ptr& acf_sptr) + { + af_sptr = acq_templ.clone(); + af_sptr->fill(1.0); + acf_sptr = af_sptr->clone(); + acq_sens_mod.unnormalise(*af_sptr); + acf_sptr->inv(0, *af_sptr); + } + + 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 prompts_from_listmode( + const STIRListmodeData& lm_data, + double start, double stop, + const STIRAcquisitionData& acq_data_template, + std::shared_ptr& prompts_sptr, + const std::string prompts_prefix = "prompts") + { + set_input(lm_data); + set_output(prompts_prefix); + set_template(acq_data_template); + set_time_interval(start, stop); + set_up(); + process_data(); + prompts_sptr = get_output(); + } + + void prompts_and_randoms_from_listmode( + const STIRListmodeData& lm_data, + double start, double stop, + const STIRAcquisitionData& acq_data_template, + std::shared_ptr& prompts_sptr, + std::shared_ptr& randoms_sptr, + const std::string prompts_prefix = "prompts") + { + set_input(lm_data); + set_output(prompts_prefix); + set_template(acq_data_template); + set_time_interval(start, stop); + set_up(); + process_data(); + prompts_sptr = get_output(); + estimate_randoms(); + randoms_sptr = 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 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(); + } + catch (...) + { + THROW("Fatal error in PETSingleScatterSimulator::set_up: attenuation_image has not been set"); + } + this->set_activity_image_sptr(sptr_act_image_template); + + if (stir::SingleScatterSimulation::set_up() == Succeeded::no) + THROW("Fatal error in PETSingleScatterSimulator::set_up() failed."); + } void set_activity_image_sptr(std::shared_ptr arg) { @@ -811,187 +1050,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/CMakeLists.txt b/src/xSTIR/cSTIR/tests/CMakeLists.txt index 8c185147b..55fcd697d 100644 --- a/src/xSTIR/cSTIR/tests/CMakeLists.txt +++ b/src/xSTIR/cSTIR/tests/CMakeLists.txt @@ -29,9 +29,15 @@ target_link_libraries(cstir_test6 csirf cstir ${STIR_LIBRARIES}) INSTALL(TARGETS cstir_test6 DESTINATION bin) + add_executable(cstir_test7 test7.cpp ${STIR_REGISTRIES}) + target_link_libraries(cstir_test7 csirf cstir ${STIR_LIBRARIES}) + INSTALL(TARGETS cstir_test7 DESTINATION bin) + ADD_TEST(NAME PET_TESTS_CPLUSPLUS_1 COMMAND cstir_tests WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) add_SIRF_DATA_PATH(PET_TESTS_CPLUSPLUS_1) ADD_TEST(NAME PET_TESTS_CPLUSPLUS_4 COMMAND cstir_test4 WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) add_SIRF_DATA_PATH(PET_TESTS_CPLUSPLUS_4) ADD_TEST(NAME PET_TESTS_CPLUSPLUS_6 COMMAND cstir_test6 WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) add_SIRF_DATA_PATH(PET_TESTS_CPLUSPLUS_6) +ADD_TEST(NAME PET_TESTS_CPLUSPLUS_7 COMMAND cstir_test7 WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) +add_SIRF_DATA_PATH(PET_TESTS_CPLUSPLUS_7) diff --git a/src/xSTIR/cSTIR/tests/test7.cpp b/src/xSTIR/cSTIR/tests/test7.cpp new file mode 100644 index 000000000..29bb7ea14 --- /dev/null +++ b/src/xSTIR/cSTIR/tests/test7.cpp @@ -0,0 +1,151 @@ +/* +SyneRBI Synergistic Image Reconstruction Framework (SIRF) +Copyright 2020 - 2024 Rutherford Appleton Laboratory STFC +Copyright 2020 - 2024 University College London + +This is software developed for the Collaborative Computational +Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR) +(http://www.ccpsynerbi.ac.uk/). + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at +http://www.apache.org/licenses/LICENSE-2.0 +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. + +*/ + +/*! +\file +\ingroup PET + +\author Kris Thielemans +\author Evgueni Ovtchinnikov +\author SyneRBI +*/ + +#include +#include + +#include "stir/common.h" +#include "stir/IO/stir_ecat_common.h" +#include "stir/Verbosity.h" + +#include "sirf/STIR/stir_x.h" +#include "sirf/common/getenv.h" +#include "sirf/common/iequals.h" +#include "sirf/common/csirf.h" +#include "sirf/common/utilities.h" + +#include "object.h" + +using namespace stir; +using namespace ecat; +using namespace sirf; + +int main() +{ + TextWriter w; // create writer with no output + TextWriterHandle h; + h.set_information_channel(&w); // suppress STIR info output + + std::cout << "running test7.cpp...\n"; + try { + std::string SIRF_data_path = examples_data_path("PET"); + if (SIRF_data_path.length() < 1) { + std::cout << "cannot find data" << std::endl; + 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"); + std::string f_norm = append_path(path, "norm.n.hdr"); + + STIRAcquisitionDataInMemory::set_as_template(); + + //PETAcquisitionModelUsingRayTracingMatrix acq_mod; + CREATE_OBJECT(PETAcquisitionModel, PETAcquisitionModelUsingRayTracingMatrix, am, am_sptr,); + //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; + std::shared_ptr randoms_sptr; + converter.set_input(lm_data); + converter.set_output("prompts"); + converter.set_template(acq_data_template); + converter.set_time_interval(0, 10); + converter.set_up(); + converter.process_data(); + sinograms_sptr = converter.get_output(); + converter.estimate_randoms(); + randoms_sptr = converter.get_randoms_sptr(); +// converter.prompts_and_randoms_from_listmode +// (lm_data, 0, 10, acq_data_template, sinograms_sptr, randoms_sptr); + + std::cout << "===== sinograms norm: " << sinograms_sptr->norm() << '\n'; + std::cout << "===== randoms norm: " << randoms_sptr->norm() << '\n'; + + am.set_up(sinograms_sptr, mu_map_sptr); + STIRAcquisitionData& sinograms = *sinograms_sptr; + PETAttenuationModel att(mu_map, am); + att.set_up(sinograms.get_exam_info_sptr(), + sinograms.get_proj_data_info_sptr()->create_shared_clone()); + std::shared_ptr af_sptr; // attenuation factor + std::shared_ptr acf_sptr; // the inverse of the above + PETAttenuationModel::compute_ac_factors(sinograms, att, af_sptr, acf_sptr); + std::cout << "===== norm of the attenuation factor: " << af_sptr->norm() << '\n'; + std::cout << "===== norm of the attenuation correction factor: " << acf_sptr->norm() << '\n'; + + CREATE_OBJ(PETAcquisitionSensitivityModel, acq_sm, acq_sm_sptr, f_norm); + + PETScatterEstimator se; + se.set_input_sptr(sinograms_sptr); + se.set_attenuation_image_sptr(mu_map_sptr); + se.set_background_sptr(randoms_sptr); + se.set_asm(acq_sm_sptr); + se.set_attenuation_correction_factors_sptr(acf_sptr); + se.set_num_iterations(4); + std::cout << "===== number of scatter iterations that will be used: " << se.get_num_iterations() << '\n'; + se.set_OSEM_num_subsets(7); + se.set_output_prefix("scatter"); + se.set_up(); + se.process(); + std::shared_ptr scatter_sptr = se.get_output(); + //scatter_estimate.write(scatter_file + '.hs'); + std::cout << "===== scatter estimate norm: " << scatter_sptr->norm() << '\n'; + + acq_sm.set_up(*acf_sptr); + std::shared_ptr mf_sptr = af_sptr->clone(); + acq_sm.unnormalise(*mf_sptr); + std::cout << "===== multfactors norm: " << mf_sptr->norm() << '\n'; + + shared_ptr background_sptr(randoms_sptr->new_acquisition_data()); + float alpha = 1.0; + background_sptr->axpby(&alpha, *randoms_sptr, &alpha, *scatter_sptr); + std::cout << "===== background norm: " << background_sptr->norm() << '\n'; + + CREATE_OBJ(PETAcquisitionSensitivityModel, asm_, asm_sptr, *mf_sptr); + asm_.set_up(*background_sptr); + asm_.normalise(*background_sptr); + std::cout << "===== additive term norm: " << background_sptr->norm() << '\n'; + + std::cout << "done with test7.cpp...\n"; + return 0; + } + catch (...) + { + return 1; + } +} diff --git a/src/xSTIR/pSTIR/STIR.py b/src/xSTIR/pSTIR/STIR.py index 4d0ab9b7f..fcc6328b1 100644 --- a/src/xSTIR/pSTIR/STIR.py +++ b/src/xSTIR/pSTIR/STIR.py @@ -457,6 +457,7 @@ def get_info(self): pyiutil.deleteDataHandle(handle) return info + @property def modality(self): """Returns imaging modality as Python string.""" return parms.char_par(self.handle, 'ImageData', 'modality') @@ -1296,6 +1297,11 @@ def get_tof_mash_factor(self): '''Returns TOF mashing factor.''' return parms.int_par(self.handle, 'AcquisitionData', 'tof_mash_factor') + @property + def modality(self): + """Returns imaging modality as Python string.""" + return parms.char_par(self.handle, 'AcquisitionData', 'modality') + def as_array(self): """Returns bin values as ndarray. @@ -1723,6 +1729,20 @@ def invert(self, ad): check_status(fd.handle) return fd + @staticmethod + def compute_attenuation_factors(sinograms, mu_map): + '''Creates attenuation model and returns the attenuation factor (af) + and the attenuation correction factor (acf) as AcquisitionData objects + ''' + am = AcquisitionModelUsingRayTracingMatrix() + attn = AcquisitionSensitivityModel(mu_map, am) + af = AcquisitionData(sinograms) + acf = AcquisitionData(sinograms) + am.set_up(sinograms, mu_map) + attn.set_up(sinograms) + try_calling(pystir.cSTIR_computeACF(sinograms.handle, attn.handle, af.handle, acf.handle)) + return af, acf + def __del__(self): """del.""" if self.handle is not None: