Continue working on moving Challenge data preparation utilities to SI…

…RF. (#1266)

Moved data preparation utilities from Challenge24 to SIRF.

* started moving Challenge nema-data utilities to SIRF

* reimplemented two nema data processing functions in SIRF (C++)

* moved the functions of the previous commit to more sensible places

* added C++ code for testing ScatterEstimator in test7.cpp

* rewrote all my Challenge functions in C++, hit incompatible proj data bug [ci skip]

* fixed incompatible proj data bug

* added simpler set up for PETAcquisitionSensitivityModel

* implemented Python interface for sinograms_and_randoms_from_listmode

* implemented Python interface to compute_ac_factors (not tested yet)

* fixed typos in STIR.py lines 1740 and 1741

* interfaced all C++ nema-data utilities in SIRF into Python

* made some corrections/amendments suggested by KT

* attended to Codacy issues

* removed hardwired prompts prefix from prompts and randoms computation code

* documented the type of objects returned by compute_attenuation_factors() [ci skip]

* removed unused import of existing_filepath from test_data_preparation.py

* compute_ac_factors now gets PETAcquisitionModel as an argument

* removed commented-out member function of ListmodeToSinograms

* attended to reviewers remarks on computing attenuation factors

* simplified compute_ac_factors methods (C++/C/Python)

* attended to reviewer's comments and suggestions

* updated CHANGES.md

* added methods for computing prompts only from list-mode data

* corrected the return of prompts_from_listmpde [ci skip]

* added modality() method to AcquisitionData objects

* commented out no longer used functions

* commented out no longer used functions in stir_x.h

* removed commented out lines confusing Codacy

* removed outdated PET/test_data_preparation.py

* updated CHANGES.md

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.

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

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()

src/iUtilities/include/sirf/iUtilities/DataHandle.h

@@ -279,6 +279,23 @@ getObjectSptrFromHandle(const void* h, std::shared_ptr<Object>& sptr) {
 	sptr = *ptr_sptr;
 }
 
+template<class Object>
+void
+setHandleObjectSptr(void* h, std::shared_ptr<Object>& sptr) {
+	auto handle = reinterpret_cast<ObjectHandle<Object>*>(h);
+	//ObjectHandle<Object>* handle = (ObjectHandle<Object>*)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<Object>, ptr_sptr, ptr);
+	//delete ptr_sptr;
+	*ptr_sptr = sptr;
+}
+
 #if defined(USE_BOOST)
 template<class Object>
 void

src/xSTIR/cSTIR/cstir.cpp

@@ -39,6 +39,8 @@ using namespace sirf;
 #define NEW_OBJECT_HANDLE(T) new ObjectHandle<T >(std::shared_ptr<T >(new T))
 #define SPTR_FROM_HANDLE(Object, X, H) \
   std::shared_ptr<Object> X; getObjectSptrFromHandle<Object>(H, X);
+#define HANDLE_FROM_SPTR(Object, X, H) \
+  setHandleObjectSptr<Object>(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<ListmodeToSinograms>(ptr_lm2s);
+		STIRListmodeData& lm_data = objectFromHandle<STIRListmodeData>(ptr_lmdata);
+		STIRAcquisitionData& templ = objectFromHandle<STIRAcquisitionData>(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<ListmodeToSinograms>(ptr_lm2s);
+		STIRListmodeData& lm_data = objectFromHandle<STIRListmodeData>(ptr_lmdata);
+		STIRAcquisitionData& templ = objectFromHandle<STIRAcquisitionData>(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<STIRAcquisitionData>(ptr_sino);
+		PETAttenuationModel& att = objectFromHandle<PETAttenuationModel>(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)

src/xSTIR/cSTIR/cstir_p.cpp

@@ -103,6 +103,8 @@ sirf::cSTIR_AcquisitionDataParameter(void* hp, const char* name)
 	STIRAcquisitionData& ad = objectFromHandle<STIRAcquisitionData>(hp);
 	if (sirf::iequals(name, "tof_mash_factor"))
 		return dataHandle<int>(ad.get_tof_mash_factor());
+	if (sirf::iequals(name, "modality"))
+		return charDataHandleFromCharData(ad.modality().c_str());
 	else
 		return parameterNotFound(name, __FILE__, __LINE__);
 }

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);

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(); }