Allow even number of ToF bins

documentation/STIR-glossary.tex

@@ -135,6 +135,9 @@ \section*{Basics}
 Set of \textbf{merged} \textbf{sinograms} with a common average \textbf{ring 
 difference} as shown in Annex 1.
 
+\item[TOF number] 
+The line of response is divided $N_t$ TOF positions the numbering goes from 0 to $N_t -1$
+
 \item[Viewgram] 
 Set of equal azimuth \textbf{merged} \textbf{LORs} of a \textbf{segment}.
 

src/IO/stir_ecat_common.cxx

@@ -208,7 +208,7 @@ find_timing_poss_sequence(const ProjDataInfo& pdi)
   for (int timing_pos_num = 1; timing_pos_num <= max_timing_pos_num; ++timing_pos_num)
     {
       timing_pos_sequence[2 * timing_pos_num - 1] = timing_pos_num;
-      timing_pos_sequence[2 * timing_pos_num] = -timing_pos_num;
+      timing_pos_sequence[2 * timing_pos_num] = max_timing_pos_num - timing_pos_num;
     }
   return timing_pos_sequence;
 }

src/buildblock/ProjDataFromStream.cxx

@@ -151,6 +151,10 @@ ProjDataFromStream::activate_TOF()
 void
 ProjDataFromStream::set_timing_poss_sequence_in_stream(const std::vector<int>& seq)
 {
+    if(!(seq[0]==this->proj_data_info_sptr->get_min_tof_pos_num() &&
+         seq[seq.size()-1]==this->proj_data_info_sptr->get_max_tof_pos_num()))
+        error("the timing position sequence is different from STIR convention 0-max!");
+
   this->timing_poss_sequence = seq;
 }
 

src/buildblock/ProjDataInfo.cxx

@@ -68,10 +68,7 @@ START_NAMESPACE_STIR
 float
 ProjDataInfo::get_k(const Bin& bin) const
 {
-  if (!(num_tof_bins % 2))
-    return bin.timing_pos_num() * tof_increament_in_mm + tof_increament_in_mm / 2.f;
-  else
-    return (bin.timing_pos_num() * tof_increament_in_mm);
+  return (bin.timing_pos_num() - (num_tof_bins - 1) / 2.f) * tof_increament_in_mm;
 }
 
 double
@@ -212,14 +209,15 @@ ProjDataInfo::set_tof_mash_factor(const int new_num)
       tof_increament_in_mm = tof_delta_time_to_mm(tof_mash_factor * scanner_ptr->get_size_of_timing_pos());
 
       // TODO cope with even numbers!
-      min_tof_pos_num = -(scanner_ptr->get_max_num_timing_poss() / tof_mash_factor) / 2;
-      max_tof_pos_num = min_tof_pos_num + (scanner_ptr->get_max_num_timing_poss() / tof_mash_factor) - 1;
+      min_tof_pos_num = 0; //- (scanner_ptr->get_max_num_timing_poss() / tof_mash_factor)/2;
+      max_tof_pos_num = (scanner_ptr->get_max_num_timing_poss() / tof_mash_factor)
+                        - 1; // min_tof_pos_num + (scanner_ptr->get_max_num_timing_poss() / tof_mash_factor) -1;
 
       num_tof_bins = max_tof_pos_num - min_tof_pos_num + 1;
 
       // Ensure that we have a central tof bin.
-      if (num_tof_bins % 2 == 0)
-        error("ProjDataInfo: Number of TOF bins should be an odd number. Abort.");
+      //        if (num_tof_bins%2 == 0)
+      //            error("ProjDataInfo: Number of TOF bins should be an odd number. Abort.");
 
       // Upper and lower boundaries of the timing poss;
       tof_bin_boundaries_mm.grow(min_tof_pos_num, max_tof_pos_num);

src/buildblock/ProjDataInfoCylindrical.cxx

@@ -534,12 +534,14 @@ ProjDataInfoCylindrical::get_LOR(LORInAxialAndNoArcCorrSinogramCoordinates<float
   const float z1 = (m_in_mm - max_a * tantheta);
   const float z2 = (m_in_mm - min_a * tantheta);
 
+  const bool swap = false; // get_tof_delta_time(bin)<0;
   lor = LORInAxialAndNoArcCorrSinogramCoordinates<float>(z1,
                                                          z2,
+
                                                          phi,
                                                          asin(s_in_mm / get_ring_radius()),
                                                          get_ring_radius(),
-                                                         false); // needs to set "swapped" to false given above code
+                                                         swap); // needs to set "swapped" to false given above code
 }
 
 #if 0

src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx

@@ -346,8 +346,12 @@ ProjDataInfoCylindricalNoArcCorr::get_all_det_pos_pairs_for_bin(vector<Detection
       // 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);
+      min_timing_pos_num = bin.timing_pos_num()
+                           * get_tof_mash_factor(); // bin.timing_pos_num()*get_tof_mash_factor() - (get_tof_mash_factor() / 2);
+      if (get_tof_mash_factor() == 0)
+        max_timing_pos_num = bin.timing_pos_num(); // bin.timing_pos_num()*get_tof_mash_factor() + (get_tof_mash_factor() / 2);
+      else
+        max_timing_pos_num = bin.timing_pos_num() * get_tof_mash_factor() + (get_tof_mash_factor()) - 1;
     }
 
   unsigned int current_dp_num = 0;
@@ -368,6 +372,7 @@ ProjDataInfoCylindricalNoArcCorr::get_all_det_pos_pairs_for_bin(vector<Detection
               dps[current_dp_num].pos2().tangential_coord() = det2_num;
               dps[current_dp_num].pos2().axial_coord() = rings_iter->second;
               dps[current_dp_num].timing_pos() = uncompressed_timing_pos_num;
+
               ++current_dp_num;
             }
         }
@@ -470,8 +475,7 @@ ProjDataInfoCylindricalNoArcCorr::find_cartesian_coordinates_given_scanner_coord
 {
   const int num_detectors_per_ring = get_scanner_ptr()->get_num_detectors_per_ring();
 
-  int d1, d2, r1, r2;
-  int tpos = timing_pos_num;
+  int d1, d2, r1, r2, tpos;
 
   this->initialise_det1det2_to_uncompressed_view_tangpos_if_not_done_yet();
 
@@ -481,14 +485,15 @@ ProjDataInfoCylindricalNoArcCorr::find_cartesian_coordinates_given_scanner_coord
       d2 = det1;
       r1 = Ring_B;
       r2 = Ring_A;
-      tpos *= -1;
+      tpos = get_max_tof_pos_num() - timing_pos_num;
     }
   else
     {
       d1 = det1;
       d2 = det2;
       r1 = Ring_A;
       r2 = Ring_B;
+      tpos = timing_pos_num;
     }
 
 #if 0
@@ -519,7 +524,8 @@ ProjDataInfoCylindricalNoArcCorr::find_cartesian_coordinates_given_scanner_coord
   coord_2 = lor.p2();
 
 #endif
-  if (tpos < 0)
+
+  if (tpos - get_max_tof_pos_num() / 2.F < 0)
     std::swap(coord_1, coord_2);
 }
 

src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl

@@ -124,7 +124,7 @@ ProjDataInfoCylindricalNoArcCorr::get_bin_for_det_pair(
     }
   else
     {
-      bin.timing_pos_num() = -timing_pos_num;
+      bin.timing_pos_num() = get_max_tof_pos_num() - timing_pos_num;
       return get_segment_axial_pos_num_for_ring_pair(bin.segment_num(), bin.axial_pos_num(), ring_num2, ring_num1);
     }
 }
@@ -137,9 +137,8 @@ ProjDataInfoCylindricalNoArcCorr::get_bin_for_det_pos_pair(Bin& bin, const Detec
                               dp.pos1().axial_coord(),
                               dp.pos2().tangential_coord(),
                               dp.pos2().axial_coord(),
-                              this->get_tof_mash_factor() == 0
-                                  ? 0 // use timing_pos==0 in the nonTOF case
-                                  : stir::round((float)dp.timing_pos() / this->get_tof_mash_factor()));
+                              this->get_tof_mash_factor() == 0 ? 0 // use timing_pos==0 in the nonTOF case
+                                                               : (int)(dp.timing_pos() / this->get_tof_mash_factor()));
 }
 void
 ProjDataInfoCylindricalNoArcCorr::get_det_pair_for_bin(
@@ -166,7 +165,7 @@ ProjDataInfoCylindricalNoArcCorr::get_det_pos_pair_for_bin(DetectionPositionPair
   // lousy work around because types don't match (short/int). TODO remove!
   int t1, a1, t2, a2;
   get_det_pair_for_bin(t1, a1, t2, a2, bin);
-  if (bin.timing_pos_num() >= 0)
+  if (bin.timing_pos_num() - get_max_tof_pos_num() / 2.F >= 0)
     {
       dp.pos1().tangential_coord() = t1;
       dp.pos1().axial_coord() = a1;

src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl

@@ -259,7 +259,7 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_general_bin(int s, int segm
                     view180, axial_pos_shift, z_shift, transform_z); // s < 0
             }
           else // neg_plus90
-            /////
+               /////
             if (segment_num < 0)
               {
                 if (!do_symmetry_swap_s || s > 0)
@@ -316,8 +316,8 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_general_bin(int s, int segm
               else
                 return new SymmetryOperation_PET_CartesianGrid_swap_ymy_zq(view180, axial_pos_shift, z_shift, transform_z);
             } // segment_num == 0
-          // /*else{   if ( !do_symmetry_swap_s || s > 0 ) return new SymmetryOperation_PET_CartesianGrid_swap_xmx();	else
-          // return new SymmetryOperation_PET_CartesianGrid_swap_ymy(view180, axial_pos_shift, z_shift);}*/
+              // /*else{   if ( !do_symmetry_swap_s || s > 0 ) return new SymmetryOperation_PET_CartesianGrid_swap_xmx();
+              // else return new SymmetryOperation_PET_CartesianGrid_swap_ymy(view180, axial_pos_shift, z_shift);}*/
         }
       else
         {
@@ -326,7 +326,8 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_general_bin(int s, int segm
           if (!do_symmetry_swap_segment || segment_num > 0)
             {
               if (do_symmetry_swap_s && s < 0)
-                return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(view180, axial_pos_shift, z_shift, transform_z);
+                return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(
+                    view180, axial_pos_shift, z_shift, transform_z, proj_data_info_ptr->get_max_tof_pos_num());
               else
                 {
                   if (z_shift == 0)
@@ -341,14 +342,18 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_general_bin(int s, int segm
                 return new SymmetryOperation_PET_CartesianGrid_swap_zq(view180, axial_pos_shift, z_shift, transform_z);
               else*/
               if (do_symmetry_swap_s && s < 0)
-                return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift);
+                return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(
+                    view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_max_tof_pos_num());
               else
                 return new SymmetryOperation_PET_CartesianGrid_swap_zq(view180, axial_pos_shift, z_shift, transform_z); // s > 0
             }
           else // segment_num = 0
             {
+
               if (do_symmetry_swap_s && s < 0)
-                return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift);
+                return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(
+                    view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_max_tof_pos_num());
+
               else
                 {
                   if (z_shift == 0)
@@ -462,7 +467,7 @@ DataSymmetriesForBins_PET_CartesianGrid::find_basic_bin(
         {
           // when swap_s, must invert timing pos for lor probs. Symmetry operation should correct bin
           tangential_pos_num *= -1;
-          timing_pos_num *= -1;
+          timing_pos_num = proj_data_info_ptr->get_max_tof_pos_num() - timing_pos_num;
           change = true;
         }
       if (do_symmetry_shift_z && axial_pos_num != 0)

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h

@@ -174,7 +174,7 @@ class PoissonLogLikelihoodWithLinearModelForMeanAndProjData
   const ProjData& get_proj_data() const;
   const shared_ptr<ProjData>& get_proj_data_sptr() const;
   const int get_max_segment_num_to_process() const;
-  const int get_max_timing_pos_num_to_process() const;
+  const int get_max_num_central_timing_poss_to_process() const;
   const bool get_zero_seg0_end_planes() const;
   const ProjData& get_additive_proj_data() const;
   const shared_ptr<ProjData>& get_additive_proj_data_sptr() const;
@@ -196,7 +196,7 @@ class PoissonLogLikelihoodWithLinearModelForMeanAndProjData
   int set_num_subsets(const int num_subsets) override;
   void set_proj_data_sptr(const shared_ptr<ProjData>&);
   void set_max_segment_num_to_process(const int);
-  void set_max_timing_pos_num_to_process(const int);
+  void set_max_num_central_timing_poss_to_process(const int);
   void set_zero_seg0_end_planes(const bool);
   // N.E. Changed to ExamData
   void set_additive_proj_data_sptr(const shared_ptr<ExamData>&) override;
@@ -298,7 +298,8 @@ class PoissonLogLikelihoodWithLinearModelForMeanAndProjData
 
   //! the maximum absolute time-of-flight bin number to use in the reconstruction
   /*! convention: if -1, use get_max_tof_pos_num()*/
-  int max_timing_pos_num_to_process;
+  int max_num_central_timing_poss_to_process;
+  int tof_bin_shift;
 
   /**********************/
   ParseAndCreateFrom<TargetT, ProjData> target_parameter_parser;

src/include/stir/recon_buildblock/SymmetryOperations_PET_CartesianGrid.h

@@ -299,14 +299,14 @@ class SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq : public SymmetryOpera
   typedef SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq self;
 
 public:
-  SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(const int num_views,
-                                                      const int axial_pos_shift,
-                                                      const int z_shift,
-                                                      const int q)
+  SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq(
+      const int num_views, const int axial_pos_shift, const int z_shift, const int q, const int max_tof_pos_num)
       : view180(num_views),
         axial_pos_shift(axial_pos_shift),
         z_shift(z_shift),
-        q(q)
+        q(q),
+        timing_shift(max_tof_pos_num)
+
   {}
 
   inline void transform_bin_coordinates(Bin&) const override;
@@ -322,6 +322,7 @@ class SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq : public SymmetryOpera
   int axial_pos_shift;
   int z_shift;
   int q;
+  int timing_shift;
 };
 
 class SymmetryOperation_PET_CartesianGrid_swap_xy_ymx_zq : public SymmetryOperation
@@ -441,10 +442,14 @@ class SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy : public SymmetryOperatio
   typedef SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy self;
 
 public:
-  SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(const int num_views, const int axial_pos_shift, const int z_shift)
+  SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(const int num_views,
+                                                   const int axial_pos_shift,
+                                                   const int z_shift,
+                                                   const int max_tof_pos_num)
       : view180(num_views),
         axial_pos_shift(axial_pos_shift),
-        z_shift(z_shift)
+        z_shift(z_shift),
+        timing_shift(max_tof_pos_num)
   {}
 
   inline void transform_bin_coordinates(Bin&) const override;
@@ -459,6 +464,7 @@ class SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy : public SymmetryOperatio
   int view180;
   int axial_pos_shift;
   int z_shift;
+  int timing_shift;
 };
 
 class SymmetryOperation_PET_CartesianGrid_swap_xmy_ymx_zq : public SymmetryOperation

src/include/stir/recon_buildblock/SymmetryOperations_PET_CartesianGrid.inl

@@ -328,7 +328,7 @@ SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq::transform_bin_coordinates(B
 {
   b.axial_pos_num() += axial_pos_shift;
   b.tangential_pos_num() *= -1;
-  b.timing_pos_num() *= -1;
+  b.timing_pos_num() = timing_shift - b.timing_pos_num();
 }
 
 void
@@ -514,7 +514,7 @@ SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy::transform_bin_coordinates(Bin&
   b.axial_pos_num() += axial_pos_shift;
   b.segment_num() *= -1;
   b.tangential_pos_num() *= -1;
-  b.timing_pos_num() *= -1;
+  b.timing_pos_num() = timing_shift - b.timing_pos_num();
 }
 void
 SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy::transform_view_segment_indices(ViewSegmentNumbers& vs) const