Allow even number of ToF bins

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
@@ -223,14 +220,14 @@ 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

@@ -582,13 +582,13 @@ get_LOR(LORInAxialAndNoArcCorrSinogramCoordinates<float>& lor,
   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

@@ -334,7 +334,7 @@ get_num_det_pos_pairs_for_bin(const Bin& bin, bool ignore_non_spatial_dimensions
 {
   return
     get_num_ring_pairs_for_segment_axial_pos_num(bin.segment_num(),
-						 bin.axial_pos_num())*
+                         bin.axial_pos_num())*
     get_view_mashing_factor()*
     (ignore_non_spatial_dimensions ? 1 : std::max(1,get_tof_mash_factor()));
 }
@@ -362,8 +362,12 @@ get_all_det_pos_pairs_for_bin(vector<DetectionPositionPair<> >& dps,
       // 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;
@@ -387,6 +391,7 @@ get_all_det_pos_pairs_for_bin(vector<DetectionPositionPair<> >& dps,
               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;
             }
         }
@@ -494,8 +499,8 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3D<floa
   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();
 
@@ -505,14 +510,15 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3D<floa
       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
@@ -543,8 +549,9 @@ find_cartesian_coordinates_given_scanner_coordinates (CartesianCoordinate3D<floa
   coord_2 = lor.p2();
 
 #endif
-  if (tpos < 0)
-    std::swap(coord_1, coord_2);
+
+  if (tpos -get_max_tof_pos_num()/2.F < 0)
+      std::swap(coord_1, coord_2);
 }
 
 void 

src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl

@@ -138,7 +138,7 @@ get_bin_for_det_pair(Bin& bin,
   }
   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);
   }
 }
@@ -156,7 +156,7 @@ get_bin_for_det_pos_pair(Bin& bin,
                          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()));
+                           : (int)(dp.timing_pos()/this->get_tof_mash_factor()));
 }
 void
 ProjDataInfoCylindricalNoArcCorr::
@@ -189,7 +189,7 @@ get_det_pos_pair_for_bin(
   //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

@@ -353,7 +353,7 @@ find_sym_op_general_bin(
     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()-1);
       else
       {
         if (z_shift==0)
@@ -369,13 +369,13 @@ find_sym_op_general_bin(
           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()-1);
           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);
+        if ( do_symmetry_swap_s && s < 0) return new SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy(view180, axial_pos_shift, z_shift, proj_data_info_ptr->get_max_tof_pos_num()-1);
         else
         {
           if (z_shift==0)
@@ -484,7 +484,7 @@ find_basic_bin(int &segment_num, int &view_num, int &axial_pos_num, int &tangent
     {
       //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    )  { axial_pos_num  =  0;     change = true; }   

src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h

@@ -179,7 +179,7 @@ public  RegisteredParsingObject<PoissonLogLikelihoodWithLinearModelForMeanAndPro
   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;
@@ -201,7 +201,7 @@ public  RegisteredParsingObject<PoissonLogLikelihoodWithLinearModelForMeanAndPro
   int set_num_subsets(const int num_subsets);
   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
   virtual void set_additive_proj_data_sptr(const shared_ptr<ExamData>&);
@@ -311,7 +311,8 @@ public  RegisteredParsingObject<PoissonLogLikelihoodWithLinearModelForMeanAndPro
 
   //! 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

@@ -339,8 +339,8 @@ class SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq : public SymmetryOpera
 private:
   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)
-    : view180(num_views), axial_pos_shift(axial_pos_shift), z_shift(z_shift), q(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), timing_shift(max_tof_pos_num)
   {}
 
   inline void 
@@ -364,6 +364,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
@@ -497,8 +498,8 @@ class SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy : public SymmetryOperatio
 private:
   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)
-    : view180(num_views), axial_pos_shift(axial_pos_shift), z_shift(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), timing_shift(max_tof_pos_num)
   {}
 
   inline void 
@@ -521,6 +522,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

@@ -351,7 +351,7 @@ SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy_zq::
 {
   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 
@@ -553,7 +553,7 @@ SymmetryOperation_PET_CartesianGrid_swap_xmx_ymy::
   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::

src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx

@@ -95,7 +95,7 @@ set_defaults()
 
   this->input_filename="";
   this->max_segment_num_to_process=-1;
-  this->max_timing_pos_num_to_process=0;
+  this->max_num_central_timing_poss_to_process=-1;
   // KT 20/06/2001 disabled
   //num_views_to_add=1;  
   this->proj_data_sptr.reset(); //MJ added
@@ -322,8 +322,8 @@ get_max_segment_num_to_process() const
 template <typename TargetT>
 const int
 PoissonLogLikelihoodWithLinearModelForMeanAndProjData<TargetT>::
-get_max_timing_pos_num_to_process() const
-{ return this->max_timing_pos_num_to_process; }
+get_max_num_central_timing_poss_to_process() const
+{ return this->max_num_central_timing_poss_to_process; }
 
 template <typename TargetT>
 const bool 
@@ -416,10 +416,10 @@ set_max_segment_num_to_process(const int arg)
 template<typename TargetT>
 void
 PoissonLogLikelihoodWithLinearModelForMeanAndProjData<TargetT>::
-set_max_timing_pos_num_to_process(const int arg)
+set_max_num_central_timing_poss_to_process(const int arg)
 {
-  this->already_set_up = this->already_set_up && (this->max_timing_pos_num_to_process == arg);
-  this->max_timing_pos_num_to_process = arg;
+  this->already_set_up = this->already_set_up && (this->max_num_central_timing_poss_to_process == arg);
+  this->max_num_central_timing_poss_to_process = arg;
 }
 
 template<typename TargetT>
@@ -614,8 +614,11 @@ set_up_before_sensitivity(shared_ptr<const TargetT > const& target_sptr)
       return Succeeded::no;
     }
 
-  this->max_timing_pos_num_to_process =
-    this->proj_data_sptr->get_max_tof_pos_num();
+    this->max_num_central_timing_poss_to_process =
+      this->proj_data_sptr->get_max_tof_pos_num();
+//  The following "shift"is used to estimate the maximum and minimum tof position number to process
+//  i.e. min_pos= num_tof_bins+shift   and   max_pos= num_tof_bins-shift   if we use all than shift =0
+  this->tof_bin_shift = (proj_data_sptr->get_num_tof_poss()-this->max_num_central_timing_poss_to_process)/2;
 
   shared_ptr<ProjDataInfo> proj_data_info_sptr(this->proj_data_sptr->get_proj_data_info_sptr()->clone());
 
@@ -745,8 +748,8 @@ actual_compute_subset_gradient_without_penalty(TargetT& gradient,
                                  this->additive_proj_data_sptr,
                                  this->normalisation_sptr,
                                  caching_info_ptr,
-                                   -this->max_timing_pos_num_to_process,
-                                   this->max_timing_pos_num_to_process,
+                                 proj_data_sptr->get_min_tof_pos_num() + tof_bin_shift,
+                                 proj_data_sptr->get_max_tof_pos_num() - tof_bin_shift,
                                  add_sensitivity);
 }
 
@@ -789,8 +792,8 @@ actual_compute_objective_function_without_penalty(const TargetT& current_estimat
                                          this->get_time_frame_definitions().get_start_time(this->get_time_frame_num()),
                                          this->get_time_frame_definitions().get_end_time(this->get_time_frame_num()),
                                          this->caching_info_ptr,
-										 -this->max_timing_pos_num_to_process,
-										 this->max_timing_pos_num_to_process
+                                         proj_data_sptr->get_min_tof_pos_num() + tof_bin_shift,
+                                         proj_data_sptr->get_max_tof_pos_num() - tof_bin_shift
                                          );
 
 
@@ -805,10 +808,11 @@ sum_projection_data() const
 {
 
   float counts=0.0F;
-
+  int min_timing_poss_to_process = proj_data_sptr->get_min_tof_pos_num() + tof_bin_shift;
+  int max_timing_poss_to_process = proj_data_sptr->get_max_tof_pos_num() - tof_bin_shift;
   for (int segment_num = -max_segment_num_to_process; segment_num <= max_segment_num_to_process; ++segment_num)
   {
-	  for (int timing_pos_num = -max_timing_pos_num_to_process; timing_pos_num <= max_timing_pos_num_to_process; ++timing_pos_num)
+      for (int timing_pos_num = min_timing_poss_to_process; timing_pos_num <= max_timing_poss_to_process; ++timing_pos_num)
 	  {
 		for (int view_num = proj_data_sptr->get_min_view_num();
 			 view_num <= proj_data_sptr->get_max_view_num();
@@ -901,8 +905,8 @@ add_subset_sensitivity(TargetT& sensitivity, const int subset_num) const
                                  this->get_time_frame_definitions().get_start_time(this->get_time_frame_num()),
                                  this->get_time_frame_definitions().get_end_time(this->get_time_frame_num()),
                                  this->caching_info_ptr,
-                                 use_tofsens ? -this->max_timing_pos_num_to_process : 0,
-                                 use_tofsens ? this->max_timing_pos_num_to_process : 0);
+                                 use_tofsens ? proj_data_sptr->get_min_tof_pos_num() + tof_bin_shift : 0,
+                                 use_tofsens ? proj_data_sptr->get_max_tof_pos_num() - tof_bin_shift : 0);
 
   std::transform(sensitivity.begin_all(), sensitivity.end_all(), 
                  sensitivity_this_subset_sptr->begin_all(), sensitivity.begin_all(), 

src/recon_test/test_DataSymmetriesForBins_PET_CartesianGrid.cxx

@@ -128,16 +128,8 @@ run_tests_2_proj_matrices_1_bin(const ProjMatrixByBin& proj_matrix_no_symm,
 	{
 	  // SYM const Bin bin=*bin_iter;
 
-      cerr << "\nCurrent bin: \tsegment = " << bin.segment_num()
-           << ", \taxial pos " << bin.axial_pos_num()
-           << ", \tview = " << bin.view_num()
-           << ", \ttangential_pos_num = " << bin.tangential_pos_num()
-           << ", timing position index = " << bin.timing_pos_num()
-           << "\nSymm bin: \t\tsegment = " << elems_with_sym.get_bin().segment_num()
-           << ", \taxial pos " << elems_with_sym.get_bin().axial_pos_num()
-           << ", \tview = " << elems_with_sym.get_bin().view_num()
-           << ", \ttangential_pos_num = " << elems_with_sym.get_bin().tangential_pos_num()
-           << ", timing position index = " << bin.timing_pos_num() << "\n";
+      cerr << "\nCurrent bin: " << elems_no_sym.get_bin()
+           << "\nSymm bin:    " << elems_with_sym.get_bin()<< "\n";
 
         if (elems_no_sym != elems_with_sym)
         {