From 589ce455be665fd83fefc2f24bb5518a383f8f14 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Mon, 21 Aug 2023 14:04:41 +0100 Subject: [PATCH 1/3] fixes STIR TOF AcquisitionData algebra inserted TOF loops. However, this needs the STIR_TOF preprocessor variable to be defined, which is still to be merged on the STIR TOF branch. --- .../include/sirf/STIR/stir_data_containers.h | 34 +++++-- src/xSTIR/cSTIR/stir_data_containers.cpp | 93 +++++++++++-------- 2 files changed, 82 insertions(+), 45 deletions(-) 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..7aa11e529 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; @@ -404,16 +408,32 @@ namespace sirf { { return data()->get_max_segment_num(); } +#ifdef STIR_TOF 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); - } + return data()->get_segment_by_sinogram(segment_num, timing_pos_num); + } +#else stir::SegmentBySinogram - get_empty_segment_by_sinogram(const int segment_num) const + get_segment_by_sinogram(const int segment_num) const { - return data()->get_empty_segment_by_sinogram(segment_num); - } + return data()->get_segment_by_sinogram(segment_num); + } +#endif +#ifdef STIR_TOF + stir::SegmentBySinogram + get_empty_segment_by_sinogram(const int segment_num, const int timing_pos_num) const + { + return data()->get_empty_segment_by_sinogram(segment_num, false, timing_pos_num); + } +#else + stir::SegmentBySinogram + get_empty_segment_by_sinogram(const int segment_num) const + { + return data()->get_empty_segment_by_sinogram(segment_num); + } +#endif void set_segment(const stir::SegmentBySinogram& s) { if (data()->set_segment(s) != stir::Succeeded::yes) 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*/) From f659bca0027b47cc06ad34565c5da253f812e683 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Mon, 21 Aug 2023 21:46:54 +0100 Subject: [PATCH 2/3] change test_data_container_algebra - increase tolerance for norm/dot tests, as large data sizes give larger error - change tests calls to use check_if_*tolerance to be more informative when the test fails - remove some unnecessary as_array() calls, speeding it up a little bit --- src/common/Utilities.py | 87 +++++++++++++++-------------------------- 1 file changed, 31 insertions(+), 56 deletions(-) diff --git a/src/common/Utilities.py b/src/common/Utilities.py index 46a2153c4..1526f5e53 100644 --- a/src/common/Utilities.py +++ b/src/common/Utilities.py @@ -597,93 +597,91 @@ def test_data_container_algebra(test, x, eps=1e-5): s = x.norm() t = numpy.linalg.norm(ax) - test.check_if_equal(1, abs(t - s) <= eps * abs(t)) + # needs increased tolerance for large data size + test.check_if_equal_within_tolerance(t, s, 0, eps * 10); s = x.max() t = numpy.max(ax) - test.check_if_equal(1, abs(t - s) <= eps * abs(t)) + test.check_if_equal_within_tolerance(t, s, 0, eps); s = x.sum() t = numpy.sum(ax) r = numpy.sum(abs(ax)) - test.check_if_equal(1, abs(t - s) <= eps * r) + test.check_if_equal_within_tolerance(t, s, 0, eps); s = x.dot(y) t = numpy.vdot(ay, ax) - test.check_if_equal(1, abs(t - s) <= eps * abs(t)) + # needs increased tolerance for large data size + test.check_if_equal_within_tolerance(t, s, 0, eps * 10); x2 = x.multiply(2) ax2 = x2.as_array() s = numpy.linalg.norm(ax2 - 2*ax) t = numpy.linalg.norm(ax2) - test.check_if_equal(1, abs(s) <= eps * abs(t)) + test.check_if_zero_within_tolerance(s, eps * t) x2 *= 0 x.multiply(2, out=x2) ax2 = x2.as_array() s = numpy.linalg.norm(ax2 - 2*ax) t = numpy.linalg.norm(ax2) - test.check_if_equal(1, abs(s) <= eps * abs(t)) + test.check_if_zero_within_tolerance(s, eps * t) t = x2.norm() x2 -= x*2 s = x2.norm() - test.check_if_equal(1, abs(s) <= eps * abs(t)) + test.check_if_zero_within_tolerance(s, eps * t) y = x.multiply(x) - ax = x.as_array() ay = y.as_array() s = numpy.linalg.norm(ay - ax * ax) t = numpy.linalg.norm(ay) - test.check_if_equal(1, abs(s) <= eps * abs(t)) + test.check_if_zero_within_tolerance(s, eps * t) y *= 0 x.multiply(x, out=y) - ax = x.as_array() ay = y.as_array() s = numpy.linalg.norm(ay - ax * ax) t = numpy.linalg.norm(ay) - test.check_if_equal(1, abs(s) <= eps * abs(t)) + test.check_if_zero_within_tolerance(s, eps * t) z = x*y az = z.as_array() s = numpy.linalg.norm(az - ax * ay) t = numpy.linalg.norm(az) - test.check_if_equal(1, abs(s) <= eps * abs(t)) + test.check_if_zero_within_tolerance(s, eps * t) y = x + 1 - ax = x.as_array() ay = y.as_array() s = numpy.linalg.norm(ay - (ax + 1)) t = numpy.linalg.norm(ay) - test.check_if_equal(1, abs(s) <= eps * abs(t)) + test.check_if_zero_within_tolerance(s, eps * t) y *= 0 x.add(1, out=y) - ax = x.as_array() ay = y.as_array() s = numpy.linalg.norm(ay - (ax + 1)) t = numpy.linalg.norm(ay) - test.check_if_equal(1, abs(s) <= eps * abs(t)) + test.check_if_zero_within_tolerance(s, eps * t) z = x/y az = z.as_array() s = numpy.linalg.norm(az - ax/ay) t = numpy.linalg.norm(az) - test.check_if_equal(1, abs(s) <= eps * abs(t)) + test.check_if_zero_within_tolerance(s, eps * t) z = x/2 az = z.as_array() s = numpy.linalg.norm(az - ax/2) t = numpy.linalg.norm(az) - test.check_if_equal(1, abs(s) <= eps * abs(t)) + test.check_if_zero_within_tolerance(s, eps * t) z *= 0 x.divide(y, out=z) az = z.as_array() s = numpy.linalg.norm(az - ax/ay) t = numpy.linalg.norm(az) - test.check_if_equal(1, abs(s) <= eps * abs(t)) + test.check_if_zero_within_tolerance(s, eps * t) y = x.sapyb(1, x, -1) s = y.norm() @@ -723,7 +721,6 @@ def test_data_container_algebra(test, x, eps=1e-5): test.check_if_equal(0, s) y = x.maximum(z) - ax = x.as_array() ay = y.as_array() az = z.as_array() ay -= numpy.maximum(ax, az) @@ -732,7 +729,6 @@ def test_data_container_algebra(test, x, eps=1e-5): y *= 0 x.maximum(z, out=y) - ax = x.as_array() ay = y.as_array() az = z.as_array() ay -= numpy.maximum(ax, az) @@ -740,7 +736,6 @@ def test_data_container_algebra(test, x, eps=1e-5): test.check_if_equal(0, s) y = x.maximum(0) - ax = x.as_array() ay = y.as_array() ay -= numpy.maximum(ax, 0) s = numpy.linalg.norm(ay) @@ -748,14 +743,12 @@ def test_data_container_algebra(test, x, eps=1e-5): y *= 0 x.maximum(0, out=y) - ax = x.as_array() ay = y.as_array() ay -= numpy.maximum(ax, 0) s = numpy.linalg.norm(ay) test.check_if_equal(0, s) y = x.minimum(z) - ax = x.as_array() ay = y.as_array() az = z.as_array() ay -= numpy.minimum(ax, az) @@ -764,7 +757,6 @@ def test_data_container_algebra(test, x, eps=1e-5): y *= 0 x.minimum(z, out=y) - ax = x.as_array() ay = y.as_array() az = z.as_array() ay -= numpy.minimum(ax, az) @@ -772,7 +764,6 @@ def test_data_container_algebra(test, x, eps=1e-5): test.check_if_equal(0, s) y = x.minimum(0) - ax = x.as_array() ay = y.as_array() ay -= numpy.minimum(ax, 0) s = numpy.linalg.norm(ay) @@ -780,31 +771,27 @@ def test_data_container_algebra(test, x, eps=1e-5): y *= 0 x.minimum(0, out=y) - ax = x.as_array() ay = y.as_array() ay -= numpy.minimum(ax, 0) s = numpy.linalg.norm(ay) test.check_if_equal(0, s) y = x.exp() - ax = x.as_array() ay = y.as_array() ay -= numpy.exp(ax) s = numpy.linalg.norm(ay) t = y.norm() - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) y *= 0 x.exp(out=y) - ax = x.as_array() ay = y.as_array() ay -= numpy.exp(ax) s = numpy.linalg.norm(ay) t = y.norm() - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) y = x.log() - ax = x.as_array() ay = y.as_array() az = numpy.log(ax) numpy.nan_to_num(ay, copy=False, posinf=0.0, neginf=0.0) @@ -812,11 +799,10 @@ def test_data_container_algebra(test, x, eps=1e-5): ay -= az s = numpy.linalg.norm(ay) t = numpy.linalg.norm(az) - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) y *= 0 x.log(out=y) - ax = x.as_array() ay = y.as_array() az = numpy.log(ax) numpy.nan_to_num(ay, copy=False, posinf=0.0, neginf=0.0) @@ -824,86 +810,76 @@ def test_data_container_algebra(test, x, eps=1e-5): ay -= az s = numpy.linalg.norm(ay) t = numpy.linalg.norm(az) - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) y = x.sqrt() - ax = x.as_array() ay = y.as_array() ay -= numpy.sqrt(ax) s = numpy.linalg.norm(ay) t = y.norm() - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) y *= 0 x.sqrt(out=y) - ax = x.as_array() ay = y.as_array() ay -= numpy.sqrt(ax) s = numpy.linalg.norm(ay) t = y.norm() - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) y = x.sign() - ax = x.as_array() ay = y.as_array() ay -= numpy.sign(ax) s = numpy.linalg.norm(ay) t = y.norm() - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) y *= 0 x.sign(out=y) - ax = x.as_array() ay = y.as_array() ay -= numpy.sign(ax) s = numpy.linalg.norm(ay) t = y.norm() - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) y = x.abs() - ax = x.as_array() ay = y.as_array() ay -= numpy.abs(ax) s = numpy.linalg.norm(ay) t = y.norm() - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) y *= 0 x.abs(out=y) - ax = x.as_array() ay = y.as_array() ay -= numpy.abs(ax) s = numpy.linalg.norm(ay) t = y.norm() - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) p = -0.5 z = x.power(p) - ax = x.as_array() az = z.as_array() numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0) t = numpy.linalg.norm(az) az -= numpy.power(ax, p) numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0) s = numpy.linalg.norm(az) - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) z *= 0 x.power(p, out=z) - ax = x.as_array() az = z.as_array() numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0) t = numpy.linalg.norm(az) az -= numpy.power(ax, p) numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0) s = numpy.linalg.norm(az) - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) - ax = x.as_array() ay = -numpy.ones_like(ax)/2 y.fill(ay) z = x.power(y) - ax = x.as_array() ay = y.as_array() az = z.as_array() numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0) @@ -911,11 +887,10 @@ def test_data_container_algebra(test, x, eps=1e-5): az -= numpy.power(ax, ay) numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0) s = numpy.linalg.norm(az) - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) z *= 0 x.power(y, out=z) - ax = x.as_array() ay = y.as_array() az = z.as_array() numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0) @@ -923,7 +898,7 @@ def test_data_container_algebra(test, x, eps=1e-5): az -= numpy.power(ax, ay) numpy.nan_to_num(az, copy=False, posinf=0.0, neginf=0.0) s = numpy.linalg.norm(az) - test.check_if_equal(1, s <= eps * t) + test.check_if_zero_within_tolerance(s, eps * t) class DataContainerAlgebraTests(object): From abaa91061fa03bc84a97deee31e18b72666de6b0 Mon Sep 17 00:00:00 2001 From: Kris Thielemans Date: Tue, 22 Aug 2023 08:16:08 +0100 Subject: [PATCH 3/3] simplify PETAcquisitionData::norm() simpler loop --- src/xSTIR/cSTIR/stir_data_containers.cpp | 23 +++++++---------------- 1 file changed, 7 insertions(+), 16 deletions(-) diff --git a/src/xSTIR/cSTIR/stir_data_containers.cpp b/src/xSTIR/cSTIR/stir_data_containers.cpp index c417d0665..2aa072f97 100644 --- a/src/xSTIR/cSTIR/stir_data_containers.cpp +++ b/src/xSTIR/cSTIR/stir_data_containers.cpp @@ -23,6 +23,8 @@ limitations under the License. #include "stir/KeyParser.h" #include "stir/is_null_ptr.h" #include "stir/zoom.h" +#include "stir/CartesianCoordinate3D.h" +#include "stir/numerics/norm.h" using namespace stir; using namespace sirf; @@ -45,23 +47,12 @@ STIRAcquisitionData::norm() const { double t = 0.0; TOF_LOOP - for (int s = 0; s <= get_max_segment_num(); ++s) + for (int s = data()->get_min_segment_num(); s <= data()->get_max_segment_num(); ++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 TOF_ARG); - for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();) { - double r = *seg_iter++; - t += r*r; - } - } + const auto seg = get_segment_by_sinogram(s TOF_ARG); + t += stir::norm_squared(seg.begin_all(), seg.end_all()); } - return std::sqrt((float)t); + return static_cast(std::sqrt(t)); } void @@ -887,7 +878,7 @@ STIRImageData::set_up_geom_info() = vox_image->get_LPS_coordinates_for_indices(next_vox_along_axis); Coord3DF axis_direction = next_vox_along_axis_coord - first_vox_coord; - axis_direction /= stir::norm(axis_direction); + axis_direction /= stir::norm(axis_direction.begin(), axis_direction.end()); for (int dim = 0; dim < 3; dim++) direction[dim][axis] = axis_direction[3 - dim]; }