Enabling PET reconstruction from list mode data

examples/Python/PET/osem_lm_reconstruction.py

@@ -0,0 +1,80 @@
+'''OSEM reconstruction from listmode demo.
+
+
+Usage:
+  osem_lm_reconstruction [--help | options] <data_path>
+
+Options:
+  -e <engn>, --engine=<engn>   reconstruction engine [default: STIR]
+'''
+
+
+__version__ = '0.1.0'
+from docopt import docopt
+args = docopt(__doc__, version=__version__)
+
+
+# import engine module
+exec('from sirf.' + args['--engine'] + ' import *')
+
+
+data_path = args['<data_path>'] + '/'
+
+
+def main():
+
+    # direct all engine's information and warnings printing to files
+    msg_red = MessageRedirector('info.txt', 'warn.txt')
+
+    cache_path = data_path;
+    sens_filename = cache_path + "sens_0.hv";
+    #print(sens_filename)
+    tmpl_projdata_filename = data_path + "tmpl_scanner.hs";
+    #print(tmpl_projdata_filename)
+
+    acq_tmpl = AcquisitionData(tmpl_projdata_filename)
+    #print(acq_tmpl.dimensions())
+
+    img_data = ImageData(sens_filename)
+    #print(img_data.dimensions())
+
+    #obj_fun = PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin()
+    obj_fun = make_Poisson_loglikelihood( \
+        likelihood_type='LinearModelForMeanAndListModeDataWithProjMatrixByBin')
+    obj_fun.set_cache_path(cache_path, True)
+    obj_fun.set_skip_lm_input_file(True);
+    obj_fun.set_skip_balanced_subsets(True);
+    obj_fun.set_acquisition_data(acq_tmpl);
+    obj_fun.set_skip_balanced_subsets(True);
+    obj_fun.set_max_ring_difference(60);
+    obj_fun.set_cache_max_size(1500000000);
+
+    #print(obj_fun.get_cache_path())
+    #print(obj_fun.get_cache_max_size())
+    #print(obj_fun.get_subsensitivity_filenames())
+
+    num_subiterations = 11
+    recon = OSMAPOSLReconstructor()
+    recon.set_objective_function(obj_fun)
+    recon.set_num_subiterations(num_subiterations)
+    recon.set_output_filename_prefix(data_path + "/my_output");
+    recon.set_save_interval(1);
+
+    print('setting up, please wait...')
+    recon.set_up(img_data)
+
+    print('reconstructing, please wait...')
+    recon.reconstruct(img_data)
+
+    img_out = recon.get_output()
+    print(img_out.dimensions())
+
+# if anything goes wrong, an exception will be thrown 
+# (cf. Error Handling section in the spec)
+try:
+    main()
+    print('\n=== done with %s' % __file__)
+
+except error as err:
+    # display error information
+    print('%s' % err.value)

src/Registration/cReg/CMakeLists.txt

@@ -74,17 +74,7 @@ if (NIFTYREG_BUILT_WITH_OPENCL)
 endif()
 
 # Add boost library dependencies
-if((CMAKE_VERSION VERSION_LESS 3.5.0) OR (NOT _Boost_IMPORTED_TARGETS))
-  # This is harder than it should be on older CMake versions to be able to cope with
-  # spaces in filenames.
-  foreach(C SYSTEM FILESYSTEM)
-    target_link_libraries(Reg PUBLIC optimized "${Boost_${C}_LIBRARY_RELEASE}")
-    target_link_libraries(Reg PUBLIC debug  "${Boost_${C}_LIBRARY_DEBUG}")
-  endforeach()
-else()
-  # Nice and simple for recent CMake (which knows about your Boost version)
-  target_link_libraries(Reg PUBLIC Boost::system Boost::filesystem)
-endif()
+target_link_libraries(Reg PUBLIC Boost::system Boost::filesystem)
 
 # Link Matlab for SPM
 IF(SPM_FOUND)

src/xGadgetron/cGadgetron/CMakeLists.txt

@@ -104,17 +104,7 @@ target_include_directories(cgadgetron PUBLIC "${ISMRMRD_INCLUDE_DIR}")
 
 target_link_libraries(cgadgetron PUBLIC iutilities csirf)
 # Add boost library dependencies
-if((CMAKE_VERSION VERSION_LESS 3.5.0) OR (NOT _Boost_IMPORTED_TARGETS))
-  # This is harder than it should be on older CMake versions to be able to cope with
-  # spaces in filenames.
-  foreach(C SYSTEM FILESYSTEM THREAD DATE_TIME CHRONO)
-    target_link_libraries(cgadgetron PUBLIC optimized "${Boost_${C}_LIBRARY_RELEASE}")
-    target_link_libraries(cgadgetron PUBLIC debug  "${Boost_${C}_LIBRARY_DEBUG}")
-  endforeach()
-else()
-  # Nice and simple for recent CMake (which knows about your Boost version)
-  target_link_libraries(cgadgetron PUBLIC Boost::system Boost::filesystem Boost::thread Boost::date_time Boost::chrono)
-endif()
+target_link_libraries(cgadgetron PUBLIC Boost::system Boost::filesystem Boost::thread Boost::date_time Boost::chrono)
 
 target_link_libraries(cgadgetron PUBLIC ISMRMRD::ISMRMRD)
 target_link_libraries(cgadgetron PUBLIC "${FFTW3_LIBRARIES}")

src/xSTIR/cSTIR/CMakeLists.txt

@@ -45,17 +45,7 @@ target_include_directories(cstir PUBLIC "${STIR_INCLUDE_DIRS}")
 target_link_libraries(cstir csirf iutilities)
 target_link_libraries(cstir "${STIR_LIBRARIES}")
 # Add boost library dependencies
-if((CMAKE_VERSION VERSION_LESS 3.5.0) OR (NOT _Boost_IMPORTED_TARGETS))
-  # This is harder than it should be on older CMake versions to be able to cope with
-  # spaces in filenames.
-  foreach(C SYSTEM FILESYSTEM THREAD DATE_TIME CHRONO)
-    target_link_libraries(cstir optimized "${Boost_${C}_LIBRARY_RELEASE}")
-    target_link_libraries(cstir debug  "${Boost_${C}_LIBRARY_DEBUG}")
-  endforeach()
-else()
-  # Nice and simple for recent CMake (which knows about your Boost version)
-  target_link_libraries(cstir Boost::system Boost::filesystem Boost::thread Boost::date_time Boost::chrono)
-endif()
+target_link_libraries(cstir Boost::system Boost::filesystem Boost::thread Boost::date_time Boost::chrono)
 
 install(DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/include/"
   COMPONENT Development

src/xSTIR/cSTIR/cstir.cpp

@@ -127,6 +127,10 @@ void* cSTIR_newObject(const char* name)
 			"PoissonLogLikelihoodWithLinearModelForMeanAndProjData"))
 			return NEW_OBJECT_HANDLE
 			(xSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndProjData3DF);
+		if (sirf::iequals(name,
+			"PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin"))
+			return NEW_OBJECT_HANDLE
+			(xSTIR_PoissonLLhLinModMeanListDataProjMatBin3DF);
 		if (sirf::iequals(name, "AcqModUsingMatrix"))
 			return NEW_OBJECT_HANDLE(AcqModUsingMatrix3DF);
 #ifdef STIR_WITH_NiftyPET_PROJECTOR
@@ -206,6 +210,11 @@ void* cSTIR_setParameter
 			return
 			cSTIR_setPoissonLogLikelihoodWithLinearModelForMeanAndProjDataParameter
 			(hs, name, hv);
+        else if (sirf::iequals(obj,
+            "PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin"))
+            return
+            cSTIR_setPoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinParameter
+            (hs, name, hv);
 		else if (sirf::iequals(obj, "Reconstruction"))
 			return cSTIR_setReconstructionParameter(hs, name, hv);
 		else if (sirf::iequals(obj, "IterativeReconstruction"))
@@ -256,11 +265,16 @@ void* cSTIR_parameter(const void* ptr, const char* obj, const char* name)
 			return cSTIR_PLSPriorParameter(handle, name);
 		else if (sirf::iequals(obj, "GeneralisedObjectiveFunction"))
 			return cSTIR_generalisedObjectiveFunctionParameter(handle, name);
-		else if (sirf::iequals(obj,
-			"PoissonLogLikelihoodWithLinearModelForMeanAndProjData"))
-			return
-			cSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndProjDataParameter
-			(handle, name);
+        else if (sirf::iequals(obj,
+            "PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin"))
+            return
+            cSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinParameter
+            (handle, name);
+        else if (sirf::iequals(obj,
+            "PoissonLogLikelihoodWithLinearModelForMeanAndProjData"))
+            return
+            cSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndProjDataParameter
+            (handle, name);
 		else if (sirf::iequals(obj, "IterativeReconstruction"))
 			return cSTIR_iterativeReconstructionParameter(handle, name);
 		else if (sirf::iequals(obj, "OSMAPOSL"))
@@ -307,8 +321,13 @@ void* cSTIR_objectFromFile(const char* name, const char* filename)
                 sptr.reset(new PETAcquisitionDataInMemory(filename));
 			return newObjectHandle(sptr);
 		}
+		if (sirf::iequals(name, "ListmodeData")) {
+			std::shared_ptr<ListmodeData>
+				sptr(new ListmodeData(filename));
+			return newObjectHandle(sptr);
+		}
 		if (sirf::iequals(name, "ListmodeToSinograms")) {
-                        std::shared_ptr<ListmodeToSinograms>
+			std::shared_ptr<ListmodeToSinograms>
 				sptr(new ListmodeToSinograms(filename));
 			return newObjectHandle(sptr);
 		}

src/xSTIR/cSTIR/cstir_p.cpp

@@ -39,6 +39,9 @@ using namespace sirf;
 extern "C"
 char* charDataFromHandle(const void* ptr);
 
+extern "C"
+void* charDataHandle(const char* ptr);
+
 static void*
 handle_error(const std::string& error_string, const char* file, int line) 
 {
@@ -660,6 +663,50 @@ sirf::cSTIR_setPoissonLogLikelihoodWithLinearModelForMeanAndProjDataParameter
 	return new DataHandle;
 }
 
+void*
+sirf::cSTIR_setPoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinParameter
+(DataHandle* hp, const char* name, const DataHandle* hv)
+{
+    xSTIR_PoissonLLhLinModMeanListDataProjMatBin3DF&
+        obj_fun = objectFromHandle
+        <xSTIR_PoissonLLhLinModMeanListDataProjMatBin3DF>
+        (hp);
+    std::string arg(name);
+    const char* prefix = arg.substr(0, 10).c_str();
+    std::string cp("cache_path");
+    if (sirf::iequals(prefix, cp)) {
+        std::string s(charDataFromDataHandle(hv));
+        auto found = s.find(",");
+        if (found == std::string::npos)
+            obj_fun.set_cache_path(s, false);
+        else
+            obj_fun.set_cache_path(s.substr(0, found), true);
+    }
+    else if (sirf::iequals(name, "acquisition_data")) {
+        SPTR_FROM_HANDLE(PETAcquisitionData, sptr_ad, hv);
+        obj_fun.set_acquisition_data(sptr_ad);
+    }
+//    else if (sirf::iequals(name, "acquisition_model")) {
+//        SPTR_FROM_HANDLE(PETAcquisitionModelUsingMatrix, sptr_pm, hv);
+//        obj_fun.set_acquisition_model(sptr_pm);
+//    }
+    else if (sirf::iequals(name, "skip_lm_input_file")) {
+        obj_fun.set_skip_lm_input_file(dataFromHandle<int>(hv));
+    }
+    else if (sirf::iequals(name, "skip_balanced_subsets")) {
+        obj_fun.set_skip_balanced_subsets(dataFromHandle<int>(hv));
+    }
+    else if (sirf::iequals(name, "max_ring_difference")) {
+        obj_fun.set_max_ring_difference(dataFromHandle<int>(hv));
+    }
+    else if (sirf::iequals(name, "cache_max_size")) {
+        obj_fun.set_cache_max_size(dataFromHandle<int>(hv));
+    }
+    else
+        return parameterNotFound(name, __FILE__, __LINE__);
+    return new DataHandle;
+}
+
 void*
 sirf::cSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndProjDataParameter
 (const DataHandle* handle, const char* name)
@@ -675,6 +722,24 @@ sirf::cSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndProjDataParameter
 	return parameterNotFound(name, __FILE__, __LINE__);
 }
 
+void*
+sirf::cSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinParameter
+(DataHandle* handle, const char* name)
+{
+    xSTIR_PoissonLLhLinModMeanListDataProjMatBin3DF&
+        obj_fun = objectFromHandle
+        <xSTIR_PoissonLLhLinModMeanListDataProjMatBin3DF>
+        (handle);
+    if (sirf::iequals(name, "cache_path")) {
+        return charDataHandle(obj_fun.get_cache_path().c_str());
+    }
+    if (sirf::iequals(name, "cache_max_size"))
+        return dataHandle<int>(obj_fun.get_cache_max_size());
+    if (sirf::iequals(name, "subsensitivity_filenames"))
+        return charDataHandle(obj_fun.get_subsensitivity_filenames().c_str());
+    return parameterNotFound(name, __FILE__, __LINE__);
+}
+
 void*
 sirf::cSTIR_setReconstructionParameter
 (DataHandle* hp, const char* name, const DataHandle* hv)
@@ -684,8 +749,10 @@ sirf::cSTIR_setReconstructionParameter
 	if (sirf::iequals(name, "output_filename_prefix"))
 		recon.set_output_filename_prefix(charDataFromDataHandle(hv));
 	else if (sirf::iequals(name, "input_data")) {
-		SPTR_FROM_HANDLE(PETAcquisitionData, sptr_ad, hv);
-		recon.set_input_data(sptr_ad->data());
+		//SPTR_FROM_HANDLE(PETAcquisitionData, sptr_ad, hv);
+		//recon.set_input_data(sptr_ad->data());
+		SPTR_FROM_HANDLE(PETScanData, sptr_ad, hv);
+		recon.set_input_data(sptr_ad->data_sptr());
 	}
 	else if (sirf::iequals(name, "disable_output")) {
 		recon.set_disable_output(true);

src/xSTIR/cSTIR/include/sirf/STIR/cstir_p.h

@@ -133,6 +133,14 @@ namespace sirf {
 		cSTIR_setPoissonLogLikelihoodWithLinearModelForMeanAndProjDataParameter
 		(DataHandle* hp, const char* name, const DataHandle* hv);
 
+    void*
+        cSTIR_setPoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinParameter
+        (DataHandle* hp, const char* name, const DataHandle* hv);
+
+    void*
+        cSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinParameter
+        (DataHandle* handle, const char* name);
+
 	void*
 		cSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndProjDataParameter
 		(const DataHandle* handle, const char* name);