Merge pull request #1448 from KrisThielemans/ProjDataMaths

add numerical operations to ProjData

Author: KrisThielemans

documentation/release_6.2.htm

@@ -44,9 +44,11 @@ <h2> Summary for end users (also to be read by developers)</h2>
 <h3>New functionality</h3>
 <ul>
   <li>
-    <code>ProjDataInMemory</code> now has same methods for numerical operations as <code>Array</code>,
-    i.e. =,-,*,/,+=,-=,*=,/=,<code>find_max()</code>,<code>find_min()</code>,<code>sum()</code>.<br>
-    <a href=https://github.com/UCL/STIR/pull/1439>PR #1439</a>
+    <code>ProjData</code> now has most of the methods for numerical operations as <code>Array</code>,
+    i.e. +=,-=,*=,/=,<code>find_max()</code>,<code>find_min()</code>,<code>sum()</code>.
+    <code>ProjDataInMemory</code> adds =,-,*,/ (as well as overloads that are faster than the implementations
+    in <code>ProjData</code>).<br>
+    <a href=https://github.com/UCL/STIR/pull/1439>PR #1439</a> and <a href=https://github.com/UCL/STIR/pull/1448>PR #1448</a>
   </li>
 </ul>
 

src/buildblock/ProjData.cxx

@@ -48,6 +48,7 @@
 #include "stir/IO/stir_ecat7.h"
 #include "stir/ViewgramIndices.h"
 #include "stir/is_null_ptr.h"
+#include "stir/numerics/norm.h"
 #include <cstring>
 #include <fstream>
 #include <algorithm>
@@ -478,6 +479,167 @@ ProjData::sapyb(const ProjData& a, const ProjData& y, const ProjData& b)
   this->xapyb(*this, a, y, b);
 }
 
+float
+ProjData::sum() const
+{
+  double t = 0.;
+  for (int timing_pos_num = this->get_min_tof_pos_num(); timing_pos_num <= this->get_max_tof_pos_num(); ++timing_pos_num)
+    for (int s = this->get_min_segment_num(); s <= this->get_max_segment_num(); ++s)
+      {
+        const SegmentIndices ind(s, timing_pos_num);
+        t += this->get_segment_by_sinogram(ind).sum();
+      }
+  return static_cast<float>(t);
+}
+
+float
+ProjData::find_max() const
+{
+  float t = 0;
+  bool init = true;
+  for (int timing_pos_num = this->get_min_tof_pos_num(); timing_pos_num <= this->get_max_tof_pos_num(); ++timing_pos_num)
+    for (int s = this->get_min_segment_num(); s <= this->get_max_segment_num(); ++s)
+      {
+        const SegmentIndices ind(s, timing_pos_num);
+        const auto t_seg = this->get_segment_by_sinogram(ind).find_max();
+        if (init)
+          {
+            init = false;
+            t = t_seg;
+          }
+        else
+          {
+            t = std::max(t, t_seg);
+          }
+      }
+  return t;
+}
+
+float
+ProjData::find_min() const
+{
+  float t = 0;
+  bool init = true;
+  for (int timing_pos_num = this->get_min_tof_pos_num(); timing_pos_num <= this->get_max_tof_pos_num(); ++timing_pos_num)
+    for (int s = this->get_min_segment_num(); s <= this->get_max_segment_num(); ++s)
+      {
+        const SegmentIndices ind(s, timing_pos_num);
+        const auto t_seg = this->get_segment_by_sinogram(ind).find_min();
+        if (init)
+          {
+            init = false;
+            t = t_seg;
+          }
+        else
+          {
+            t = std::min(t, t_seg);
+          }
+      }
+  return t;
+}
+
+double
+ProjData::norm_squared() const
+{
+  double t = 0.;
+  for (int timing_pos_num = this->get_min_tof_pos_num(); timing_pos_num <= this->get_max_tof_pos_num(); ++timing_pos_num)
+    for (int s = this->get_min_segment_num(); s <= this->get_max_segment_num(); ++s)
+      {
+        const SegmentIndices ind(s, timing_pos_num);
+        const auto seg = this->get_segment_by_sinogram(ind);
+        t += stir::norm_squared(seg.begin_all(), seg.end_all());
+      }
+  return t;
+}
+
+double
+ProjData::norm() const
+{
+  return std::sqrt(this->norm_squared());
+}
+
+// static helper functions to iterate over segment and apply a function
+
+// func(s1, s2) is supposed to modify s1
+template <typename Func>
+static ProjData&
+apply_func(ProjData& self, const ProjData& arg, Func func)
+{
+  for (int timing_pos_num = self.get_min_tof_pos_num(); timing_pos_num <= self.get_max_tof_pos_num(); ++timing_pos_num)
+    for (int s = self.get_min_segment_num(); s <= self.get_max_segment_num(); ++s)
+      {
+        const SegmentIndices ind(s, timing_pos_num);
+        auto seg = self.get_segment_by_sinogram(ind);
+        func(seg, arg.get_segment_by_sinogram(ind));
+        self.set_segment(seg);
+      }
+  return self;
+}
+
+// func(s1) is supposed to modify s1
+template <typename Func>
+static ProjData&
+apply_func(ProjData& self, Func func)
+{
+  for (int timing_pos_num = self.get_min_tof_pos_num(); timing_pos_num <= self.get_max_tof_pos_num(); ++timing_pos_num)
+    for (int s = self.get_min_segment_num(); s <= self.get_max_segment_num(); ++s)
+      {
+        const SegmentIndices ind(s, timing_pos_num);
+        auto seg = self.get_segment_by_sinogram(ind);
+        func(seg);
+        self.set_segment(seg);
+      }
+  return self;
+}
+
+ProjData&
+ProjData::operator+=(const ProjData& arg)
+{
+  return apply_func(*this, arg, [](SegmentBySinogram<float>& s, const SegmentBySinogram<float>& s_arg) { s += s_arg; });
+}
+
+ProjData&
+ProjData::operator-=(const ProjData& arg)
+{
+  return apply_func(*this, arg, [](SegmentBySinogram<float>& s, const SegmentBySinogram<float>& s_arg) { s -= s_arg; });
+}
+
+ProjData&
+ProjData::operator*=(const ProjData& arg)
+{
+  return apply_func(*this, arg, [](SegmentBySinogram<float>& s, const SegmentBySinogram<float>& s_arg) { s *= s_arg; });
+}
+
+ProjData&
+ProjData::operator/=(const ProjData& arg)
+{
+  return apply_func(*this, arg, [](SegmentBySinogram<float>& s, const SegmentBySinogram<float>& s_arg) { s /= s_arg; });
+}
+
+ProjData&
+ProjData::operator+=(float arg)
+{
+  return apply_func(*this, [arg](SegmentBySinogram<float>& s) { s += arg; });
+}
+
+ProjData&
+ProjData::operator-=(float arg)
+{
+  return apply_func(*this, [arg](SegmentBySinogram<float>& s) { s -= arg; });
+}
+
+ProjData&
+ProjData::operator*=(float arg)
+{
+  return apply_func(*this, [arg](SegmentBySinogram<float>& s) { s *= arg; });
+}
+
+ProjData&
+ProjData::operator/=(float arg)
+{
+  return apply_func(*this, [arg](SegmentBySinogram<float>& s) { s /= arg; });
+}
+
 std::vector<int>
 ProjData::standard_segment_sequence(const ProjDataInfo& pdi)
 {

src/buildblock/ProjDataInMemory.cxx

@@ -393,30 +393,44 @@ ProjDataInMemory::norm_squared() const
 }
 
 ProjDataInMemory&
-ProjDataInMemory::operator+=(const ProjDataInMemory& v)
+ProjDataInMemory::operator+=(const base_type& v)
 {
-  this->buffer += v.buffer;
+  if (auto vp = dynamic_cast<const ProjDataInMemory*>(&v))
+    this->buffer += vp->buffer;
+  else
+    base_type::operator+=(v);
+
   return *this;
 }
 
 ProjDataInMemory&
-ProjDataInMemory::operator-=(const ProjDataInMemory& v)
+ProjDataInMemory::operator-=(const base_type& v)
 {
-  this->buffer -= v.buffer;
+  if (auto vp = dynamic_cast<const ProjDataInMemory*>(&v))
+    this->buffer -= vp->buffer;
+  else
+    base_type::operator-=(v);
   return *this;
 }
 
 ProjDataInMemory&
-ProjDataInMemory::operator*=(const ProjDataInMemory& v)
+ProjDataInMemory::operator*=(const base_type& v)
 {
-  this->buffer *= v.buffer;
+  if (auto vp = dynamic_cast<const ProjDataInMemory*>(&v))
+    this->buffer *= vp->buffer;
+  else
+    base_type::operator*=(v);
   return *this;
 }
 
 ProjDataInMemory&
-ProjDataInMemory::operator/=(const ProjDataInMemory& v)
+ProjDataInMemory::operator/=(const base_type& v)
 {
-  this->buffer /= v.buffer;
+  if (auto vp = dynamic_cast<const ProjDataInMemory*>(&v))
+    this->buffer /= vp->buffer;
+  else
+    base_type::operator/=(v);
+
   return *this;
 }
 

src/include/stir/ProjData.h

@@ -102,6 +102,11 @@ class ProjDataInMemory;
 */
 class ProjData : public ExamData
 {
+#ifdef SWIG
+ public: // SWIG needs self_type exposed
+#endif
+  typedef ProjData self_type;
+
 public:
   //! A static member to get the projection data from a file
   static shared_ptr<ProjData> read_from_file(const std::string& filename, const std::ios::openmode open_mode = std::ios::in);
@@ -371,6 +376,47 @@ class ProjData : public ExamData
   //! writes data to a file in Interfile format
   Succeeded write_to_file(const std::string& filename) const;
 
+  //! @name arithmetic operations
+  ///@{
+  //! return sum of all elements
+  virtual float sum() const;
+
+  //! return maximum value of all elements
+  virtual float find_max() const;
+
+  //! return minimum value of all elements
+  virtual float find_min() const;
+
+  //! return L2-norm (sqrt of sum of squares)
+  virtual double norm() const;
+
+  //! return L2-norm squared (sum of squares)
+  virtual double norm_squared() const;
+
+  //! adding elements of \c v to the current data
+  virtual self_type& operator+=(const self_type& v);
+
+  //! subtracting elements of \c v from the current data
+  virtual self_type& operator-=(const self_type& v);
+
+  //! multiplying elements of the current data with elements of \c v
+  virtual self_type& operator*=(const self_type& v);
+
+  //! dividing all elements of the current data by elements of \c v
+  virtual self_type& operator/=(const self_type& v);
+
+  //! adding an \c float to the elements of the current data
+  virtual self_type& operator+=(const float v);
+
+  //! subtracting an \c float from the elements of the current data
+  virtual self_type& operator-=(const float v);
+
+  //! multiplying the elements of the current data with an \c float
+  virtual self_type& operator*=(const float v);
+
+  //! dividing the elements of the current data by an \c float
+  virtual self_type& operator/=(const float v);
+
   //! \deprecated a*x+b*y (use xapyb)
   STIR_DEPRECATED virtual void axpby(const float a, const ProjData& x, const float b, const ProjData& y);
 
@@ -385,6 +431,7 @@ class ProjData : public ExamData
 
   //! set values of the array to self*a+y*b where a, b and y are ProjData
   virtual void sapyb(const ProjData& a, const ProjData& y, const ProjData& b);
+  ///@}
 
 protected:
   shared_ptr<const ProjDataInfo> proj_data_info_sptr;

src/include/stir/ProjDataInMemory.h

@@ -42,6 +42,7 @@ class ProjDataInMemory : public ProjData
  public:
 #endif
   typedef ProjDataInMemory self_type;
+  typedef ProjData base_type;
 
 public:
   //! constructor with only info, but no data
@@ -106,19 +107,19 @@ class ProjDataInMemory : public ProjData
   //! @name arithmetic operations
   ///@{
   //! return sum of all elements
-  float sum() const;
+  float sum() const override;
 
   //! return maximum value of all elements
-  float find_max() const;
+  float find_max() const override;
 
   //! return minimum value of all elements
-  float find_min() const;
+  float find_min() const override;
 
   //! return L2-norm (sqrt of sum of squares)
-  double norm() const;
+  double norm() const override;
 
   //! return L2-norm squared (sum of squares)
-  double norm_squared() const;
+  double norm_squared() const override;
 
   //! elem by elem addition
   self_type operator+(const self_type& iv) const;
@@ -146,29 +147,29 @@ class ProjDataInMemory : public ProjData
 
   // corresponding assignment operators
 
-  //! adding elements of \c v to the current vector
-  self_type& operator+=(const self_type& v);
+  //! adding elements of \c v to the current data
+  self_type& operator+=(const base_type& v) override;
 
-  //! subtracting elements of \c v from the current vector
-  self_type& operator-=(const self_type& v);
+  //! subtracting elements of \c v from the current data
+  self_type& operator-=(const base_type& v) override;
 
-  //! multiplying elements of the current vector with elements of \c v
-  self_type& operator*=(const self_type& v);
+  //! multiplying elements of the current data with elements of \c v
+  self_type& operator*=(const base_type& v) override;
 
-  //! dividing all elements of the current vector by elements of \c v
-  self_type& operator/=(const self_type& v);
+  //! dividing all elements of the current data by elements of \c v
+  self_type& operator/=(const base_type& v) override;
 
-  //! adding an \c float to the elements of the current vector
-  self_type& operator+=(const float v);
+  //! adding an \c float to the elements of the current data
+  self_type& operator+=(const float v) override;
 
-  //! subtracting an \c float from the elements of the current vector
-  self_type& operator-=(const float v);
+  //! subtracting an \c float from the elements of the current data
+  self_type& operator-=(const float v) override;
 
-  //! multiplying the elements of the current vector with an \c float
-  self_type& operator*=(const float v);
+  //! multiplying the elements of the current data with an \c float
+  self_type& operator*=(const float v) override;
 
-  //! dividing the elements of the current vector by an \c float
-  self_type& operator/=(const float v);
+  //! dividing the elements of the current data by an \c float
+  self_type& operator/=(const float v) override;
 
   //! \deprecated a*x+b*y (use xapyb)
   STIR_DEPRECATED void axpby(const float a, const ProjData& x, const float b, const ProjData& y) override;