Merge pull request #1028 from SyneRBI/subset_object

Subset object in SIRF

examples/Python/PET/acquisition_data.py

@@ -86,6 +86,15 @@ def main():
         acq_data.show(range(dim[1]//4))
     acq_array = acq_data.as_array()
 
+    if storage[0] == 'm':  # for now, we can only subset acquisition data stored in memory
+        nv = dim[2]//2
+        views = numpy.arange(nv)
+        acq_subset = acq_data.get_subset(views)
+        dim_subset = acq_subset.dimensions()
+        print('subset dimensions: %d x %d x %d x %d' % dim_subset)
+        if show_plot:
+            acq_subset.show(range(dim[1]//4), title='Sinograms of a subset of views')
+
     # rebin the acquisition data
     new_acq_data = acq_data.rebin(3)
     rdim = new_acq_data.dimensions()

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

@@ -268,8 +268,10 @@ template<class Object>
 void
 getObjectSptrFromHandle(const void* h, std::shared_ptr<Object>& sptr) {
 	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");

src/xSTIR/cSTIR/cstir.cpp

@@ -19,6 +19,8 @@ limitations under the License.
 
 */
 
+#include <vector>
+
 #include "sirf/common/iequals.h"
 #include "sirf/STIR/stir_types.h"
 #include "sirf/iUtilities/DataHandle.h"
@@ -891,6 +893,19 @@ void* cSTIR_get_ProjDataInfo(void* ptr_acq)
 	CATCH;
 }
 
+extern "C"
+void* cSTIR_get_subset(void* ptr_acq, int nv, size_t ptr_views)
+{
+	try {
+		SPTR_FROM_HANDLE(PETAcquisitionData, sptr_ad, ptr_acq);
+		int* ptr_v = (int*)ptr_views;
+		std::vector<int> v(ptr_v, ptr_v + nv);
+		std::shared_ptr<PETAcquisitionData> sptr = std::move(sptr_ad->get_subset(v));
+		return newObjectHandle(sptr);
+	}
+	CATCH;
+}
+
 extern "C"
 void* cSTIR_setupFBP2DReconstruction(void* ptr_r, void* ptr_i)
 {

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

@@ -118,6 +118,7 @@ extern "C" {
 		(void* ptr_acq, const void * ptr_from);
 	void* cSTIR_writeAcquisitionData(void* ptr_acq, const char* filename);
 	void* cSTIR_get_ProjDataInfo(void* ptr_acq);
+	void* cSTIR_get_subset(void* ptr_acq, int nv, size_t ptr_views);
 
 	// Reconstruction methods
 	void* cSTIR_setupFBP2DReconstruction(void* ptr_r, void* ptr_i);

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

@@ -94,7 +94,7 @@ namespace sirf {
 			_filename(filename),
 			_owns_file(owns_file)
 		{}
-		ProjDataFile(stir::shared_ptr<stir::ExamInfo> sptr_exam_info,
+		ProjDataFile(stir::shared_ptr<const stir::ExamInfo> sptr_exam_info,
 			stir::shared_ptr<stir::ProjDataInfo> sptr_proj_data_info,
 			const std::string& filename, bool owns_file = true) :
 			stir::ProjDataInterfile(SPTR_WRAP(sptr_exam_info), SPTR_WRAP(sptr_proj_data_info),
@@ -160,6 +160,8 @@ namespace sirf {
 			return false;
 		}
 
+		virtual std::unique_ptr<PETAcquisitionData> get_subset(const std::vector<int>& views) const = 0;
+
 		//! rebin the data to lower resolution by adding
 		/*!
 		  \param num_segments_to_combine combines multiple oblique 'segments' together. If set to the
@@ -446,6 +448,24 @@ namespace sirf {
 			ptr->fill(0.0f);
 			_data.reset(ptr);
 		}
+		PETAcquisitionDataInFile(std::unique_ptr<stir::ProjData> uptr_pd) : _owns_file(true)
+		{
+//			auto *pd_ptr = dynamic_cast<stir::ProjDataInterfile*>(uptr_pd.get());
+			auto pd_ptr = dynamic_cast<stir::ProjDataInterfile*>(uptr_pd.get());
+			if (pd_ptr)
+				_data = std::move(uptr_pd);
+			else {
+				stir::ProjData& pd = *uptr_pd;
+				stir::shared_ptr<const stir::ExamInfo> sptr_exam_info =
+					pd.get_exam_info_sptr();
+				stir::shared_ptr<stir::ProjDataInfo> sptr_proj_data_info =
+					pd.get_proj_data_info_sptr()->create_shared_clone();
+				_data.reset(new ProjDataFile
+                                    (MAKE_SHARED<stir::ExamInfo>(*sptr_exam_info), sptr_proj_data_info,
+					_filename = SIRFUtilities::scratch_file_name()));
+				_data->fill(pd);
+			}
+		}
 		std::shared_ptr<PETAcquisitionData> new_acquisition_data(std::string filename)
 		{
 			std::shared_ptr<PETAcquisitionDataInFile> sptr_ad(new PETAcquisitionDataInFile);
@@ -493,6 +513,7 @@ namespace sirf {
 				(_template->same_acquisition_data(this->get_exam_info_sptr(),
 				this->get_proj_data_info_sptr()->create_shared_clone()));
 		}
+		virtual std::unique_ptr<PETAcquisitionData> get_subset(const std::vector<int>& views) const;
 
 	private:
 		bool _owns_file;
@@ -536,6 +557,22 @@ namespace sirf {
 			ptr->fill(0.0f);
 			_data.reset(ptr);
 		}
+		PETAcquisitionDataInMemory(std::unique_ptr<stir::ProjData> uptr_pd)
+		{
+//			auto *pd_ptr = dynamic_cast<stir::ProjDataInMemory*>(uptr_pd.get());
+			auto pd_ptr = dynamic_cast<stir::ProjDataInMemory*>(uptr_pd.get());
+			if (pd_ptr)
+				_data = std::move(uptr_pd);
+			else {
+				std::cout << "copying ProjData to ProjDataInMemory...\n";
+				const stir::ProjData& pd = *uptr_pd;
+				auto exam_info_sptr = SPTR_WRAP(pd.get_exam_info_sptr());
+				auto proj_data_info_sptr =
+					SPTR_WRAP(pd.get_proj_data_info_sptr()->create_shared_clone());
+				_data.reset(new stir::ProjDataInMemory(exam_info_sptr, proj_data_info_sptr));
+				_data->fill(pd);
+			}
+		}
         /// Constructor for PETAcquisitionDataInMemory from filename
         PETAcquisitionDataInMemory(const char* filename)
         {
@@ -556,10 +593,8 @@ namespace sirf {
 				(new stir::ProjDataInMemory(*pd_sptr));
         }
 
-		static void init() 
-		{ 
-			PETAcquisitionDataInFile::init(); 
-		}
+		static void init();
+
 		static void set_as_template()
 		{
 			init();
@@ -592,6 +627,8 @@ namespace sirf {
 				(this->get_exam_info_sptr(),
                                  this->get_proj_data_info_sptr()->create_shared_clone()));
 		}
+		virtual std::unique_ptr<PETAcquisitionData> get_subset(const std::vector<int>& views) const;
+
         /// fill with single value
         virtual void fill(const float v)
         {

src/xSTIR/cSTIR/stir_data_containers.cpp

@@ -245,6 +245,28 @@ PETAcquisitionData::binary_op_(
 	}
 }
 
+std::unique_ptr<PETAcquisitionData>
+PETAcquisitionDataInFile::get_subset(const std::vector<int>& views) const
+{
+	auto ptr_ad = new PETAcquisitionDataInFile(std::move(_data->get_subset(views)));
+//	auto ptr_ad = new PETAcquisitionDataInMemory(std::move(_data->get_subset(views)));
+	return std::unique_ptr<PETAcquisitionData>(ptr_ad);
+}
+
+std::unique_ptr<PETAcquisitionData>
+PETAcquisitionDataInMemory::get_subset(const std::vector<int>& views) const
+{
+	auto ptr_ad = new PETAcquisitionDataInMemory(std::move(_data->get_subset(views)));
+	return std::unique_ptr<PETAcquisitionData>(ptr_ad);
+}
+
+void
+PETAcquisitionDataInMemory::init()
+{
+	PETAcquisitionDataInFile::init();
+}
+
+
 STIRImageData::STIRImageData(const ImageData& id)
 {
     throw std::runtime_error("TODO - create STIRImageData from general SIRFImageData.");

src/xSTIR/pSTIR/STIR.py

@@ -1106,7 +1106,7 @@ def rebin(self, num_segments_to_combine,
         return ad
 
     def show(self, sino=None, tof=0, title=None):
-        '''Displays interactively selected sinograms.'''
+        '''Displays selected sinograms.'''
         if self.handle is None:
             raise AssertionError()
         if not HAVE_PYLAB:
@@ -1176,6 +1176,19 @@ def get_info(self):
         pyiutil.deleteDataHandle(handle)
         return info
 
+    def get_subset(self, views):
+        """Returns the subset of self data formed by specified views
+
+        views: array of views (will be converted to numpy ndarray)
+        """
+        # Ensure the array passed to C++ is a contiguous array of C++ int's
+        v = cpp_int_array(views)
+        n = len(views)
+        subset = AcquisitionData()
+        subset.handle = pystir.cSTIR_get_subset(self.handle, n, v.ctypes.data)
+        check_status(subset.handle)
+        return subset
+
     @property
     def shape(self):
         return self.dimensions()