diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h index 6d7011e53..de26f0afb 100644 --- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h +++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h @@ -249,7 +249,11 @@ namespace sirf { if (_is_empty != -1) return _is_empty ? 0 : 1; try { - get_segment_by_sinogram(0); +#ifdef STIR_TOF + get_segment_by_sinogram(0,0); +#else + get_segment_by_sinogram(0); +#endif } catch (...) { _is_empty = 1; @@ -405,10 +409,23 @@ namespace sirf { return data()->get_max_segment_num(); } stir::SegmentBySinogram - get_segment_by_sinogram(const int segment_num) const + get_segment_by_sinogram(const int segment_num, const int timing_pos_num) const { - return data()->get_segment_by_sinogram(segment_num); - } +#ifdef STIR_TOF + return data()->get_segment_by_sinogram(segment_num, timing_pos_num); +#else + return data()->get_segment_by_sinogram(segment_num); +#endif + } + stir::SegmentBySinogram + get_empty_segment_by_sinogram(const int segment_num, const int timing_pos_num ) const + { +#ifdef STIR_TOF + return data()->get_empty_segment_by_sinogram(segment_num, false, timing_pos_num); +#else + return data()->get_empty_segment_by_sinogram(segment_num); +#endif + } stir::SegmentBySinogram get_empty_segment_by_sinogram(const int segment_num) const { diff --git a/src/xSTIR/cSTIR/stir_data_containers.cpp b/src/xSTIR/cSTIR/stir_data_containers.cpp index 2105b9276..c417d0665 100644 --- a/src/xSTIR/cSTIR/stir_data_containers.cpp +++ b/src/xSTIR/cSTIR/stir_data_containers.cpp @@ -1,7 +1,7 @@ /* SyneRBI Synergistic Image Reconstruction Framework (SIRF) -Copyright 2017 - 2019 Rutherford Appleton Laboratory STFC -Copyright 2018 - 2020 University College London +Copyright 2017 - 2023 Rutherford Appleton Laboratory STFC +Copyright 2018 - 2023 University College London This is software developed for the Collaborative Computational Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR) @@ -32,20 +32,29 @@ using namespace sirf; std::string STIRAcquisitionData::_storage_scheme; std::shared_ptr STIRAcquisitionData::_template; +#ifdef STIR_TOF +#define TOF_LOOP for (int k=data()->get_min_tof_pos_num(); k<=data()->get_max_tof_pos_num(); ++k) +#define TOF_ARG , k +#else +#define TOF_LOOP +#define TOF_ARG +#endif + float STIRAcquisitionData::norm() const { double t = 0.0; + TOF_LOOP for (int s = 0; s <= get_max_segment_num(); ++s) { - SegmentBySinogram seg = get_segment_by_sinogram(s); + SegmentBySinogram seg = get_segment_by_sinogram(s TOF_ARG); SegmentBySinogram::full_iterator seg_iter; for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();) { double r = *seg_iter++; t += r*r; } if (s != 0) { - seg = get_segment_by_sinogram(-s); + seg = get_segment_by_sinogram(-s TOF_ARG); for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();) { double r = *seg_iter++; t += r*r; @@ -60,14 +69,15 @@ STIRAcquisitionData::sum(void* ptr) const { int n = get_max_segment_num(); double t = 0; + TOF_LOOP for (int s = 0; s <= n; ++s) { - SegmentBySinogram seg = get_segment_by_sinogram(s); + SegmentBySinogram seg = get_segment_by_sinogram(s TOF_ARG); SegmentBySinogram::full_iterator seg_iter; for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();) t += *seg_iter++; if (s != 0) { - seg = get_segment_by_sinogram(-s); + seg = get_segment_by_sinogram(-s TOF_ARG); for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();) t += *seg_iter++; } @@ -81,14 +91,15 @@ STIRAcquisitionData::max(void* ptr) const { int n = get_max_segment_num(); float t = 0; + TOF_LOOP for (int s = 0; s <= n; ++s) { - SegmentBySinogram seg = get_segment_by_sinogram(s); + SegmentBySinogram seg = get_segment_by_sinogram(s TOF_ARG); SegmentBySinogram::full_iterator seg_iter; for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();) t = std::max(t, *seg_iter++); if (s != 0) { - seg = get_segment_by_sinogram(-s); + seg = get_segment_by_sinogram(-s TOF_ARG); for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();) t = std::max(t, *seg_iter++); } @@ -104,10 +115,11 @@ STIRAcquisitionData::dot(const DataContainer& a_x, void* ptr) const int n = get_max_segment_num(); int nx = x.get_max_segment_num(); double t = 0; + TOF_LOOP for (int s = 0; s <= n && s <= nx; ++s) { - SegmentBySinogram seg = get_segment_by_sinogram(s); - SegmentBySinogram sx = x.get_segment_by_sinogram(s); + SegmentBySinogram seg = get_segment_by_sinogram(s TOF_ARG); + SegmentBySinogram sx = x.get_segment_by_sinogram(s TOF_ARG); SegmentBySinogram::full_iterator seg_iter; SegmentBySinogram::full_iterator sx_iter; for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all(); @@ -116,8 +128,8 @@ STIRAcquisitionData::dot(const DataContainer& a_x, void* ptr) const t += (*seg_iter++) * double(*sx_iter++); } if (s != 0) { - seg = get_segment_by_sinogram(-s); - sx = x.get_segment_by_sinogram(-s); + seg = get_segment_by_sinogram(-s TOF_ARG); + sx = x.get_segment_by_sinogram(-s TOF_ARG); for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all(); seg_iter != seg.end_all() && sx_iter != sx.end_all(); /*empty*/) @@ -188,11 +200,12 @@ STIRAcquisitionData::inv(float amin, const DataContainer& a_x) SIRF_DYNAMIC_CAST(const STIRAcquisitionData, x, a_x); int n = get_max_segment_num(); int nx = x.get_max_segment_num(); + TOF_LOOP for (int s = 0; s <= n && s <= nx; ++s) { //std::cout << "processing segment " << s << std::endl; - SegmentBySinogram seg = get_empty_segment_by_sinogram(s); - SegmentBySinogram sx = x.get_segment_by_sinogram(s); + SegmentBySinogram seg = get_empty_segment_by_sinogram(s TOF_ARG); + SegmentBySinogram sx = x.get_segment_by_sinogram(s TOF_ARG); SegmentBySinogram::full_iterator seg_iter; SegmentBySinogram::full_iterator sx_iter; for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all(); @@ -202,8 +215,8 @@ STIRAcquisitionData::inv(float amin, const DataContainer& a_x) set_segment(seg); if (s != 0) { //std::cout << "processing segment " << -s << std::endl; - seg = get_empty_segment_by_sinogram(-s); - sx = x.get_segment_by_sinogram(-s); + seg = get_empty_segment_by_sinogram(-s TOF_ARG); + sx = x.get_segment_by_sinogram(-s TOF_ARG); for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all(); seg_iter != seg.end_all() && sx_iter != sx.end_all(); /*empty*/) { @@ -223,9 +236,10 @@ STIRAcquisitionData::unary_op( SIRF_DYNAMIC_CAST(const STIRAcquisitionData, x, a_x); int n = get_max_segment_num(); int nx = x.get_max_segment_num(); + TOF_LOOP for (int s = 0; s <= n && s <= nx; ++s) { - SegmentBySinogram seg = get_empty_segment_by_sinogram(s); - SegmentBySinogram sx = x.get_segment_by_sinogram(s); + SegmentBySinogram seg = get_empty_segment_by_sinogram(s TOF_ARG); + SegmentBySinogram sx = x.get_segment_by_sinogram(s TOF_ARG); SegmentBySinogram::full_iterator seg_iter; SegmentBySinogram::full_iterator sx_iter; for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all(); @@ -233,8 +247,8 @@ STIRAcquisitionData::unary_op( *seg_iter++ = f(*sx_iter++); set_segment(seg); if (s > 0) { - SegmentBySinogram seg = get_empty_segment_by_sinogram(-s); - SegmentBySinogram sx = x.get_segment_by_sinogram(-s); + SegmentBySinogram seg = get_empty_segment_by_sinogram(-s TOF_ARG); + SegmentBySinogram sx = x.get_segment_by_sinogram(-s TOF_ARG); SegmentBySinogram::full_iterator seg_iter; SegmentBySinogram::full_iterator sx_iter; for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all(); @@ -256,9 +270,10 @@ STIRAcquisitionData::semibinary_op( SIRF_DYNAMIC_CAST(const STIRAcquisitionData, x, a_x); int n = get_max_segment_num(); int nx = x.get_max_segment_num(); + TOF_LOOP for (int s = 0; s <= n && s <= nx; ++s) { - SegmentBySinogram seg = get_empty_segment_by_sinogram(s); - SegmentBySinogram sx = x.get_segment_by_sinogram(s); + SegmentBySinogram seg = get_empty_segment_by_sinogram(s TOF_ARG); + SegmentBySinogram sx = x.get_segment_by_sinogram(s TOF_ARG); SegmentBySinogram::full_iterator seg_iter; SegmentBySinogram::full_iterator sx_iter; for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all(); @@ -266,8 +281,8 @@ STIRAcquisitionData::semibinary_op( *seg_iter++ = f(*sx_iter++, y); set_segment(seg); if (s > 0) { - SegmentBySinogram seg = get_empty_segment_by_sinogram(-s); - SegmentBySinogram sx = x.get_segment_by_sinogram(-s); + SegmentBySinogram seg = get_empty_segment_by_sinogram(-s TOF_ARG); + SegmentBySinogram sx = x.get_segment_by_sinogram(-s TOF_ARG); SegmentBySinogram::full_iterator seg_iter; SegmentBySinogram::full_iterator sx_iter; for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all(); @@ -296,11 +311,12 @@ STIRAcquisitionData::binary_op( SegmentBySinogram::full_iterator seg_iter; SegmentBySinogram::full_iterator sx_iter; SegmentBySinogram::full_iterator sy_iter; + TOF_LOOP for (int s = 0; s <= n; ++s) { - SegmentBySinogram seg = get_empty_segment_by_sinogram(s); - SegmentBySinogram sx = x.get_segment_by_sinogram(s); - SegmentBySinogram sy = y.get_segment_by_sinogram(s); + SegmentBySinogram seg = get_empty_segment_by_sinogram(s TOF_ARG); + SegmentBySinogram sx = x.get_segment_by_sinogram(s TOF_ARG); + SegmentBySinogram sy = y.get_segment_by_sinogram(s TOF_ARG); if (seg.size_all() != sx.size_all() || seg.size_all() != sy.size_all()) throw std::runtime_error("binary_op error: operands sizes differ"); for (seg_iter = seg.begin_all(), @@ -309,9 +325,9 @@ STIRAcquisitionData::binary_op( *seg_iter++ = f(*sx_iter++, *sy_iter++); set_segment(seg); if (s != 0) { - seg = get_empty_segment_by_sinogram(-s); - sx = x.get_segment_by_sinogram(-s); - sy = y.get_segment_by_sinogram(-s); + seg = get_empty_segment_by_sinogram(-s TOF_ARG); + sx = x.get_segment_by_sinogram(-s TOF_ARG); + sy = y.get_segment_by_sinogram(-s TOF_ARG); for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all(), sy_iter = sy.begin_all(); seg_iter != seg.end_all(); /*empty*/) @@ -340,12 +356,13 @@ STIRAcquisitionData::xapyb( SegmentBySinogram::full_iterator sx_iter; SegmentBySinogram::full_iterator sy_iter; SegmentBySinogram::full_iterator sb_iter; + TOF_LOOP for (int s = 0; s <= n; ++s) { - SegmentBySinogram seg = get_empty_segment_by_sinogram(s); - SegmentBySinogram sx = x.get_segment_by_sinogram(s); - SegmentBySinogram sy = y.get_segment_by_sinogram(s); - SegmentBySinogram sb = b.get_segment_by_sinogram(s); + SegmentBySinogram seg = get_empty_segment_by_sinogram(s TOF_ARG); + SegmentBySinogram sx = x.get_segment_by_sinogram(s TOF_ARG); + SegmentBySinogram sy = y.get_segment_by_sinogram(s TOF_ARG); + SegmentBySinogram sb = b.get_segment_by_sinogram(s TOF_ARG); size_t size_seg = seg.size_all(); size_t size_sx = sx.size_all(); size_t size_sy = sy.size_all(); @@ -358,10 +375,10 @@ STIRAcquisitionData::xapyb( *seg_iter++ = (*sx_iter++) * a + (*sy_iter++) * (*sb_iter++); set_segment(seg); if (s != 0) { - seg = get_empty_segment_by_sinogram(-s); - sx = x.get_segment_by_sinogram(-s); - sy = y.get_segment_by_sinogram(-s); - sb = b.get_segment_by_sinogram(-s); + seg = get_empty_segment_by_sinogram(-s TOF_ARG); + sx = x.get_segment_by_sinogram(-s TOF_ARG); + sy = y.get_segment_by_sinogram(-s TOF_ARG); + sb = b.get_segment_by_sinogram(-s TOF_ARG); for (seg_iter = seg.begin_all(), sx_iter = sx.begin_all(), sy_iter = sy.begin_all(), sb_iter = sb.begin_all(); seg_iter != seg.end_all(); /*empty*/)