fix RFS for TOF but change get_all_det_pos_pairs_for_bin

- get_all_det_pos_pairs_for_bin now defaults to ignore TOF. This is
probably what is expected by most people.
- fix multiply_crystal_factors for TOF (in the previous version, it
was unexpectedly running over all unmashed TOF bins, giving the wrong
scale factor, and taking a long time

Author: KrisThielemans

src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx

@@ -336,31 +336,42 @@ initialise_det1det2_to_uncompressed_view_tangpos() const
 
 unsigned int
 ProjDataInfoCylindricalNoArcCorr::
-get_num_det_pos_pairs_for_bin(const Bin& bin) const
+get_num_det_pos_pairs_for_bin(const Bin& bin, bool ignore_non_spatial_dimensions) const
 {
   return
     get_num_ring_pairs_for_segment_axial_pos_num(bin.segment_num(),
 						 bin.axial_pos_num())*
     get_view_mashing_factor()*
-    std::max(1,get_tof_mash_factor());
+    (ignore_non_spatial_dimensions ? 1 : std::max(1,get_tof_mash_factor()));
 }
 
 void
 ProjDataInfoCylindricalNoArcCorr::
 get_all_det_pos_pairs_for_bin(vector<DetectionPositionPair<> >& dps,
-			      const Bin& bin) const
+			      const Bin& bin,
+                              bool ignore_non_spatial_dimensions) const
 {
   this->initialise_uncompressed_view_tangpos_to_det1det2_if_not_done_yet();
 
-  dps.resize(get_num_det_pos_pairs_for_bin(bin));
+  dps.resize(get_num_det_pos_pairs_for_bin(bin, ignore_non_spatial_dimensions));
 
   const ProjDataInfoCylindrical::RingNumPairs& ring_pairs =
     get_all_ring_pairs_for_segment_axial_pos_num(bin.segment_num(),
 						 bin.axial_pos_num());
   // not sure how to handle mashing with non-zero view offset...
   assert(get_min_view_num()==0);
-  // not sure how to handle even tof mashing
-  assert(!is_tof_data() || (get_tof_mash_factor() % 2 == 1));
+
+  int min_timing_pos_num = 0;
+  int max_timing_pos_num = 0;
+  if (!ignore_non_spatial_dimensions)
+    {
+      // not sure how to handle even tof mashing
+      assert(!is_tof_data() || (get_tof_mash_factor() % 2 == 1)); // TODOTOF
+      // we will need to add all (unmashed) timing_pos for the current bin
+      min_timing_pos_num = bin.timing_pos_num()*get_tof_mash_factor() - (get_tof_mash_factor() / 2);
+      max_timing_pos_num = bin.timing_pos_num()*get_tof_mash_factor() + (get_tof_mash_factor() / 2);
+    }
+
   unsigned int current_dp_num=0;
   for (int uncompressed_view_num=bin.view_num()*get_view_mashing_factor();
        uncompressed_view_num<(bin.view_num()+1)*get_view_mashing_factor();
@@ -374,25 +385,16 @@ get_all_det_pos_pairs_for_bin(vector<DetectionPositionPair<> >& dps,
         rings_iter != ring_pairs.end();
         ++rings_iter)
         {
-          for (int uncompressed_timing_pos_num = bin.timing_pos_num()*get_tof_mash_factor() - (get_tof_mash_factor() / 2);
-               uncompressed_timing_pos_num <= bin.timing_pos_num()*get_tof_mash_factor() + (get_tof_mash_factor() / 2);
+          for (int uncompressed_timing_pos_num = min_timing_pos_num;
+               uncompressed_timing_pos_num <= max_timing_pos_num;
                ++uncompressed_timing_pos_num)
             {
               assert(current_dp_num < get_num_det_pos_pairs_for_bin(bin));
               dps[current_dp_num].pos1().tangential_coord() = det1_num;
               dps[current_dp_num].pos1().axial_coord() = rings_iter->first;
               dps[current_dp_num].pos2().tangential_coord() = det2_num;
               dps[current_dp_num].pos2().axial_coord() = rings_iter->second;
-              // need to keep dp.timing_pos positive
-              if (uncompressed_timing_pos_num > 0)
-                {
-                  dps[current_dp_num].timing_pos() = static_cast<unsigned>(uncompressed_timing_pos_num);
-                }
-              else
-                {
-                  std::swap(dps[current_dp_num].pos1(), dps[current_dp_num].pos2());
-                  dps[current_dp_num].timing_pos() = static_cast<unsigned>(-uncompressed_timing_pos_num);
-                }
+              dps[current_dp_num].timing_pos() = uncompressed_timing_pos_num;
               ++current_dp_num;
             }
         }

src/buildblock/multiply_crystal_factors.cxx

@@ -24,6 +24,7 @@
 #include "stir/Bin.h"
 #include "stir/Sinogram.h"
 #include "stir/error.h"
+#include <memory>
 
 START_NAMESPACE_STIR
 
@@ -33,8 +34,9 @@ void multiply_crystal_factors_help(ProjData& proj_data,
                                    const TProjDataInfo& proj_data_info,
                                    const Array<2,float>& efficiencies, const float global_factor)
 {
-    const auto non_tof_proj_data_info_sptr = proj_data_info.create_non_tof_clone();
-    Bin bin;
+  const auto non_tof_proj_data_info_sptr =
+    std::dynamic_pointer_cast<TProjDataInfo>(proj_data_info.create_non_tof_clone());
+  Bin bin;
 
     for (bin.segment_num() = proj_data.get_min_segment_num();
      bin.segment_num() <= proj_data.get_max_segment_num();
@@ -59,8 +61,8 @@ void multiply_crystal_factors_help(ProjData& proj_data,
                  view_num <= proj_data.get_max_view_num();
                  ++ view_num)
               {
-                for (int tangential_pos_num = proj_data_info.get_min_tangential_pos_num();
-                     tangential_pos_num <= proj_data_info.get_max_tangential_pos_num();
+                for (int tangential_pos_num = proj_data.get_min_tangential_pos_num();
+                     tangential_pos_num <= proj_data.get_max_tangential_pos_num();
                      ++tangential_pos_num)
                   {
                     // Construct bin with appropriate values
@@ -70,7 +72,7 @@ void multiply_crystal_factors_help(ProjData& proj_data,
                     parallel_bin.tangential_pos_num() = tangential_pos_num;
 
                     std::vector<DetectionPositionPair<> > det_pos_pairs;
-                    proj_data_info.get_all_det_pos_pairs_for_bin(det_pos_pairs, parallel_bin);
+                    non_tof_proj_data_info_sptr->get_all_det_pos_pairs_for_bin(det_pos_pairs, parallel_bin); // using the default argument to ignore TOF here
                     float result = 0.F;
                     for (unsigned int i=0; i<det_pos_pairs.size(); ++i)
                       {

src/include/stir/ProjDataInfoCylindricalNoArcCorr.h

@@ -208,11 +208,17 @@ class ProjDataInfoCylindricalNoArcCorr : public ProjDataInfoCylindrical
 	                 const Bin&) const;
 
   //! This routine returns the number of detector pairs that correspond to a bin
+  /*!
+    \see get_all_det_pos_pairs_for_bin(). This function return the size of its output.
+  */
   unsigned int
-    get_num_det_pos_pairs_for_bin(const Bin&) const;
+  get_num_det_pos_pairs_for_bin(const Bin&, bool ignore_non_spatial_dimensions = true) const;
 
   //! This routine fills a vector with all the detector pairs that correspond to a bin.
-  /*!     
+  /*!
+    \param[in] ignore_non_spatial_dimensions if \c true, only 1 \c timing_pos_num
+       will be set, i.e. only the detection positions are returned, not TOF (nor energy)
+
     \see ProjDataInfoCylindrical::get_all_ring_pairs_for_segment_axial_pos_num
     for restrictions. 
     \todo It might be possible to return some weight factors in case there is 
@@ -222,7 +228,8 @@ class ProjDataInfoCylindricalNoArcCorr : public ProjDataInfoCylindrical
   */
   void
     get_all_det_pos_pairs_for_bin(std::vector<DetectionPositionPair<> >&,
-				  const Bin&) const;
+				  const Bin&,
+                                  bool ignore_non_spatial_dimensions = true) const;
 
  private:
   // old function, now private. Use get_bin_for_det_pos_pair instead.

src/test/test_proj_data_info.cxx

@@ -1408,7 +1408,7 @@ ProjDataInfoCylindricalNoArcCorrTests::test_proj_data_info(ProjDataInfoCylindric
                 bin.view_num() = view_num;
                 bin.tangential_pos_num() = tangential_pos_num;
                 std::vector<DetectionPositionPair<>> det_pos_pairs;
-                proj_data_info.get_all_det_pos_pairs_for_bin(det_pos_pairs, bin);
+                proj_data_info.get_all_det_pos_pairs_for_bin(det_pos_pairs, bin, false); // include TOF
                 Bin new_bin;
                 // set value for comparison with bin
                 new_bin.set_bin_value(0);