diff --git a/CHANGES.md b/CHANGES.md index eb0419098..04b829a73 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -21,11 +21,6 @@ - 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. - -## 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/test_data_preparation.py b/examples/Python/PET/test_data_preparation.py deleted file mode 100644 index d5fd38815..000000000 --- a/examples/Python/PET/test_data_preparation.py +++ /dev/null @@ -1,74 +0,0 @@ -import os -from sirf.Utilities import examples_data_path - -import importlib -pet = importlib.import_module('sirf.STIR') - -data_path = examples_data_path('PET') + '/mMR' -print('Finding files in %s...' % data_path) - -#save_path = "~/tmp/data" -#print('Saving files in %s...' % save_path) - -f_listmode = os.path.join(data_path, "list.l.hdr"); -f_template = os.path.join(data_path, "mMR_template_span11_small.hs"); -f_attn = os.path.join(data_path, "mu_map.hv"); -f_norm = os.path.join(data_path, "norm.n.hdr"); - -# engine's messages go to files, except error messages, which go to stdout -_ = pet.MessageRedirector('info.txt', 'warn.txt') - -# select acquisition data storage scheme -pet.AcquisitionData.set_storage_scheme('memory') - -# read acquisition data template -acq_data_template = pet.AcquisitionData(f_template) - -listmode_data = pet.ListmodeData(f_listmode) - -# create listmode-to-sinograms converter object -lm2sino = pet.ListmodeToSinograms() -lm_data = pet.ListmodeData(f_listmode) -prompts, randoms = lm2sino.prompts_and_randoms_from_listmode(lm_data, 0, 10, acq_data_template) -print('data shape: %s' % repr(prompts.shape)) -print('prompts norm: %f' % prompts.norm()) -print('randoms norm: %f' % randoms.norm()) -#prompts.write(os.path.join(save_path, 'prompts.hs')) -#randoms.write(os.path.join(save_path, 'randoms.hs')) - -attn_image = pet.ImageData(f_attn) -attn, acf, iacf = pet.AcquisitionSensitivityModel.compute_attenuation_factors(prompts, attn_image) -print('norm of the attenuation correction factor: %f' % acf.norm()) -print('norm of the inverse of the attenuation correction factor: %f' % iacf.norm()) -#acf.write(os.path.join(save_path, 'acf.hs')) -#iacf.write(os.path.join(save_path, 'iacf.hs')) -#exit() - -asm = pet.AcquisitionSensitivityModel(f_norm) -se = pet.ScatterEstimator() -se.set_input(prompts) -se.set_attenuation_image(attn_image) -se.set_randoms(randoms) -se.set_asm(asm) -se.set_attenuation_correction_factors(iacf) -se.set_num_iterations(4) -se.set_OSEM_num_subsets(7) -se.set_output_prefix("scatter") -se.set_up() -se.process() -scatter = se.get_output() -print('norm of the scatter estimate: %f' % scatter.norm()) - -multfact = acf.clone() -asm.set_up(acf) -asm.unnormalise(multfact) -print(multfact.norm()) - -background = randoms + scatter -print('norm of the backgrount term: %f' % background.norm()) - -asm_mf = pet.AcquisitionSensitivityModel(multfact) -asm_mf.set_up(background) -asm_mf.normalise(background) -print('norm of the additive term: %f' % background.norm()) - diff --git a/src/iUtilities/include/sirf/iUtilities/DataHandle.h b/src/iUtilities/include/sirf/iUtilities/DataHandle.h index 5b2215fa3..0b6c65808 100644 --- a/src/iUtilities/include/sirf/iUtilities/DataHandle.h +++ b/src/iUtilities/include/sirf/iUtilities/DataHandle.h @@ -279,23 +279,6 @@ 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 081adf441..29680cb34 100644 --- a/src/xSTIR/cSTIR/cstir.cpp +++ b/src/xSTIR/cSTIR/cstir.cpp @@ -39,8 +39,6 @@ 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) @@ -493,42 +491,6 @@ 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) @@ -679,23 +641,6 @@ 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/include/sirf/STIR/cstir.h b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h index e78a428ff..7433da844 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h @@ -63,12 +63,6 @@ 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); @@ -82,8 +76,6 @@ 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); @@ -100,7 +92,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_x.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h index 0dd89785d..f3fd7f715 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h @@ -81,11 +81,175 @@ 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() {} @@ -103,11 +267,6 @@ 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 @@ -375,426 +534,28 @@ The actual algorithm is described in //shared_ptr sptr_normalisation_; }; - /*! - \ingroup PET - \brief Ray tracing matrix implementation of the PET acquisition model. + /*! + \ingroup PET - 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). - */ + \brief Class for simulating the scatter contribution to PET data. - 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); - } + This class uses the STIR Single Scatter simulation, taking as input an + activity and attenuation image, and a acquisition data template. - 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) - {} + 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) @@ -1050,6 +811,187 @@ 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 55fcd697d..8c185147b 100644 --- a/src/xSTIR/cSTIR/tests/CMakeLists.txt +++ b/src/xSTIR/cSTIR/tests/CMakeLists.txt @@ -29,15 +29,9 @@ 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 deleted file mode 100644 index 9fa4f4e2e..000000000 --- a/src/xSTIR/cSTIR/tests/test7.cpp +++ /dev/null @@ -1,142 +0,0 @@ -/* -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.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 bae99fb30..4d0ab9b7f 100644 --- a/src/xSTIR/pSTIR/STIR.py +++ b/src/xSTIR/pSTIR/STIR.py @@ -1599,29 +1599,6 @@ def get_time_at_which_num_prompts_exceeds_threshold(self, threshold): pyiutil.deleteDataHandle(h) return v - def prompts_from_listmode(self, lm_data, start, stop, templ, prefix="prompts"): - """Returns proampts computed from listmode raw data - - """ - assert_validity(lm_data, ListmodeData) - assert_validity(templ, AcquisitionData) - sino = AcquisitionData(templ) - try_calling(pystir.cSTIR_promptsFromListmode(self.handle, lm_data.handle, \ - start, stop, templ.handle, sino.handle, prefix)) - return sino - - def prompts_and_randoms_from_listmode(self, lm_data, start, stop, templ, prefix="prompts"): - """Returns proampts and randoms' estimates computed from listmode raw data - - """ - assert_validity(lm_data, ListmodeData) - assert_validity(templ, AcquisitionData) - sino = AcquisitionData(templ) - rand = AcquisitionData(templ) - try_calling(pystir.cSTIR_promptsAndRandomsFromListmode(self.handle, lm_data.handle, \ - start, stop, templ.handle, sino.handle, rand.handle, prefix)) - return sino, rand - class AcquisitionSensitivityModel(object): """ @@ -1746,20 +1723,6 @@ 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: