Moving SyneRBI-Challenge nema-data utilities to SIRF

CHANGES.md

@@ -2,6 +2,11 @@
 
 ## vx.x.x
 
+* 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.

examples/Python/PET/test_data_preparation.py

@@ -0,0 +1,74 @@
+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())
+

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)
@@ -478,6 +480,42 @@ 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)
@@ -628,6 +666,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/include/sirf/STIR/cstir.h

@@ -63,6 +63,12 @@ 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);
 
@@ -75,6 +81,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);
@@ -91,7 +99,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);