Modifications for psmr2024

CMakeLists.txt

@@ -33,7 +33,7 @@ include(src/cmake/SetC++Version.cmake)
 
 # minimum C++ version required (you can still ask for a more recent one
 # by setting CMAKE_CXX_STANDARD)
-UseCXX(14)
+UseCXX(17)
 
 # set default build-type to Release
 if(NOT CMAKE_BUILD_TYPE)

src/IO/InterfileHeader.cxx

@@ -612,19 +612,19 @@ InterfilePDFSHeader::InterfilePDFSHeader()
   reference_energy = -1.f;
   add_key("Reference energy (in keV)", &reference_energy);
 
-  max_num_timing_poss = -1;
+  max_num_timing_poss = 1;
   add_key("Maximum number of (unmashed) TOF time bins", &max_num_timing_poss);
 #if STIR_VERSION < 070000
   add_alias_key("Maximum number of (unmashed) TOF time bins", "Number of TOF time bins");
 #endif
   timing_poss_sequence.clear();
   add_key("TOF bin order", &timing_poss_sequence);
-  size_of_timing_pos = -1.f;
+  size_of_timing_pos = 0.0f;
   add_key("Size of unmashed TOF time bins (ps)", &size_of_timing_pos);
 #if STIR_VERSION < 070000
   add_alias_key("Size of unmashed TOF time bins (ps)", "Size of timing bin (ps)");
 #endif
-  timing_resolution = -1.f;
+  timing_resolution = 0.0f;
   add_key("TOF timing resolution (ps)", &timing_resolution);
 #if STIR_VERSION < 070000
   add_alias_key("TOF timing resolution (ps)", "timing resolution (ps)");
@@ -1430,6 +1430,7 @@ InterfilePDFSHeader::post_processing()
       data_info_sptr.reset(new ProjDataInfoGenericNoArcCorr(
           scanner_sptr_from_file, sorted_num_rings_per_segment, sorted_min_ring_diff, sorted_max_ring_diff, num_views, num_bins));
     }
+  std::cout << data_info_sptr->get_num_tof_poss()  << " " << num_timing_poss << std::endl;
   if (data_info_sptr->get_num_tof_poss() != num_timing_poss)
     error(boost::format("Interfile header parsing with TOF: inconsistency between number of TOF bins in data (%d), "
                         "TOF mashing factor (%d) and max number of TOF bins in scanner info (%d)")

src/buildblock/ProjDataInfo.cxx

@@ -177,9 +177,41 @@ ProjDataInfo::set_tof_mash_factor(const int new_num)
     {
       tof_mash_factor = new_num;
       if (tof_mash_factor > scanner_ptr->get_max_num_timing_poss())
+        {
         error("ProjDataInfo::set_tof_mash_factor: TOF mashing factor (" + std::to_string(tof_mash_factor)
               + +") must be smaller than or equal to the scanner's number of max timing bins ("
               + std::to_string(scanner_ptr->get_max_num_timing_poss()) + ").");
+        }
+      else if (tof_mash_factor == num_tof_bins)
+        // This is a special case that we just want boundaries for the coincidence window.
+        {
+
+          Bin bin;
+          bin.timing_pos_num() = min_tof_pos_num;
+          //Get lowest low
+          float cur_low = get_k(bin) - get_sampling_in_k(bin) / 2.f;
+          //Get highest high
+          bin.timing_pos_num() = max_tof_pos_num;
+          float cur_high = get_k(bin) + get_sampling_in_k(bin) / 2.f;
+
+
+          min_tof_pos_num = 0;
+          max_tof_pos_num = 0;
+          num_tof_bins = 1;
+          // Upper and lower boundaries of the timing poss;
+          tof_bin_boundaries_mm.recycle();
+          tof_bin_boundaries_mm.grow(min_tof_pos_num, max_tof_pos_num);
+          tof_bin_boundaries_ps.recycle();
+          tof_bin_boundaries_ps.grow(min_tof_pos_num, max_tof_pos_num);
+
+          tof_bin_boundaries_mm[0].low_lim = cur_low;
+          tof_bin_boundaries_mm[0].high_lim = cur_high;
+
+          tof_bin_boundaries_ps[0].low_lim = static_cast<float>(mm_to_tof_delta_time(tof_bin_boundaries_mm[0].low_lim));
+          tof_bin_boundaries_ps[0].high_lim = static_cast<float>(mm_to_tof_delta_time(tof_bin_boundaries_mm[0].high_lim));
+
+          return;
+        }
 
 #if 0
         // KT: code disabled as buggy but currently not needed
@@ -288,7 +320,7 @@ ProjDataInfo::ProjDataInfo(const shared_ptr<Scanner>& scanner_ptr_v,
   max_tof_pos_num = 0;
   tof_increament_in_mm = 0.f;
   tof_mash_factor = 0;
-  num_tof_bins = 1;
+  num_tof_bins = 0;
 }
 
 // TOF version.

src/buildblock/Scanner.cxx

@@ -1231,6 +1231,64 @@ Scanner::Scanner(Type scanner_type)
       );
       break;
 
+
+    case PSMR_3rings_scanner:
+      set_params(PSMR_3rings_scanner,
+                 string_list("PSMR_3rings_scanner"),
+                 3 * 14 * 6,// + 8,
+                 344,
+                 344,
+                 24 * 7 * 5,
+                 405.0 - 18/2,
+                 7.0F,
+                 2.85F,
+                 2.08626F,
+                 static_cast<float>(-0.11),
+                 14,//            int num_axial_blocks_per_bucket_v,
+                 5,  //            int num_transaxial_blocks_per_bucket_v,
+                 6, //            int num_axial_crystals_per_block_v,
+                 7,//            int num_transaxial_crystals_per_block_v,
+                 252, //            int num_axial_crystals_per_singles_unit_v,
+                 35,  //            int num_transaxial_crystals_per_singles_unit_v,
+                 1,
+                 0.0F, 511.F,
+#ifdef STIR_TOF
+                 1,
+                 2011.5F*2.F, // Size of coincidence window
+                 0.F
+#endif
+                 );
+      break;
+
+  case PSMR_3rings_scanner_TOF:
+    set_params(PSMR_3rings_scanner_TOF,
+               string_list("PSMR_3rings_scanner_TOF"),
+               3 * 14 * 6,// + 8,
+               344, 344,
+               24 * 7 * 5,
+               405.0 - 18/2,
+               7.0F,
+               2.85F,
+               2.08626F,
+               static_cast<float>(-0.11),
+               14,//            int num_axial_blocks_per_bucket_v,
+               5,  //            int num_transaxial_blocks_per_bucket_v,
+               6, //            int num_axial_crystals_per_block_v,
+               7,//            int num_transaxial_crystals_per_block_v,
+               252, //            int num_axial_crystals_per_singles_unit_v,
+               35,  //            int num_transaxial_crystals_per_singles_unit_v,
+               1,
+               0.0F, 511.F,
+#ifdef STIR_TOF
+               //               (short int)(2*3.51*1000 -1),
+               //               (float)(1.f),
+               (short int)(13), // 13 TOF bins
+               (float)(2011.5F*2.F/13), //
+               (float)(500.F)
+#endif
+               );
+    break;
+
     case User_defined_scanner: // zlong, 08-04-2004, Userdefined support
 
       set_params(User_defined_scanner,

src/include/stir/ProjDataInfo.h

@@ -150,7 +150,9 @@ class ProjDataInfo
 
   //! Similar to create_shared_clone() but returns a non-tof version of ProjDataInfo setting tof mashing factor = 0
   inline shared_ptr<ProjDataInfo> create_non_tof_clone() const;
-
+  //! Similar to create_non_tof_clone() but returns a single TOF version of the ProjDataInfo setting the TOF mashing factor sush that,
+  //! a single TOF bin remains. This is usefull to apply the coincidence window boundaries without affecting the probabilities.
+    inline shared_ptr<ProjDataInfo> create_single_tof_clone() const;
   //! Destructor
   virtual ~ProjDataInfo() {}
 

src/include/stir/ProjDataInfo.inl

@@ -44,6 +44,19 @@ ProjDataInfo::create_non_tof_clone() const
   return sptr;
 }
 
+shared_ptr<ProjDataInfo>
+ProjDataInfo::create_single_tof_clone() const
+{
+  shared_ptr<ProjDataInfo> sptr(this->clone());
+  int a  = sptr->get_num_tof_poss();
+  if (a > 0)
+    sptr->set_tof_mash_factor(a);
+  // - TODO: remove this later
+  // int b  = sptr->get_num_tof_poss();
+  // std::cout <<"Old number of TOF poss " << a << " new number of TOF poss " << b << std::endl;
+  return sptr;
+}
+
 int
 ProjDataInfo::get_num_segments() const
 {

src/include/stir/Scanner.h

@@ -83,6 +83,12 @@ class Succeeded;
           while others (such as CTI scanners) give only singles for a
           collection of blocks. A \c singles_unit is then a set of crystals
           for which we can get singles rates.
+    \li \c If the user set number for TOF positions to 1 and a bin size equal to the
+          coincidence window then the reconstruction will essential be nonTOF but the
+            projector will restrict the size of the LOR to the size of the coincidence window.
+      \li \c The scanner will be classified as TOF enabled when the numer of TOF bins and
+TOF bin size are  >1 and >0, respectively. If the energy resolution is not set that will be fine
+as long as the final TOF possitions is 1. Then we just restict the size of the LOR.
 
       A further complication is that some scanners (including many Siemens scanners)
       insert virtual crystals in the sinogram data (corresponding to gaps between
@@ -168,6 +174,8 @@ class Scanner
     HZLR,
     RATPET,
     PANDA,
+    PSMR_3rings_scanner,
+    PSMR_3rings_scanner_TOF,
     HYPERimage,
     nanoPET,
     HRRT,

src/include/stir/Scanner.inl

@@ -262,7 +262,7 @@ Scanner::get_timing_resolution() const
 bool
 Scanner::is_tof_ready() const
 {
-  return (max_num_of_timing_poss > 0 && size_timing_pos > 0.0f && timing_resolution > 0.0f);
+  return (max_num_of_timing_poss > 0 && size_timing_pos > 0.0f ); //&& timing_resolution > 0.0f);
 }
 
 std::string

src/include/stir/recon_buildblock/ProjMatrixByBin.inl

@@ -67,7 +67,7 @@ ProjMatrixByBin::get_proj_matrix_elems_for_one_bin(ProjMatrixElemsForOneBin& pro
 #ifndef NDEBUG
           probabilities.check_state();
 #endif
-          if (proj_data_info_sptr->is_tof_data() && this->tof_enabled)
+          if (proj_data_info_sptr->is_tof_data())
             { // Apply TOF kernel to basic bin
               apply_tof_kernel(probabilities);
             }
@@ -97,7 +97,7 @@ ProjMatrixByBin::get_proj_matrix_elems_for_one_bin(ProjMatrixElemsForOneBin& pro
 #ifndef NDEBUG
               probabilities.check_state();
 #endif
-              if (proj_data_info_sptr->is_tof_data() && this->tof_enabled)
+              if (proj_data_info_sptr->is_tof_data())
                 { // Apply TOF kernel to basic bin
                   apply_tof_kernel(probabilities);
                 }
@@ -149,7 +149,16 @@ ProjMatrixByBin::get_tof_value(const float d1, const float d2) const
   if ((d1_n >= 4.f && d2_n >= 4.f) || (d1_n <= -4.f && d2_n <= -4.f))
     return 0.F;
   else
-    return 0.5f * (erf_interpolation(d2_n) - erf_interpolation(d1_n));
+    {
+       if (this->tof_enabled)
+        {
+             return 0.5f * (erf_interpolation(d2_n) - erf_interpolation(d1_n));
+        }
+       else
+         {
+           return 1.0f;
+         }
+    }
 }
 
 END_NAMESPACE_STIR

src/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.cxx

@@ -340,7 +340,7 @@ DataSymmetriesForBins_PET_CartesianGrid::DataSymmetriesForBins_PET_CartesianGrid
         }
 
       // RT Disabling some symmetries due to tof data
-      if (proj_data_info_ptr->is_tof_data())
+      if (proj_data_info_ptr->is_tof_data() && proj_data_info_ptr->get_num_tof_poss() > 1)
         {
           if (this->do_symmetry_90degrees_min_phi || this->do_symmetry_180degrees_min_phi)
             {

src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx

@@ -280,14 +280,14 @@ PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin<Tar
           > 1)) // TODO this check needs to cover the case if we reconstruct only TOF bin 0
     {
       // sets non-tof backprojector for sensitivity calculation (clone of the back_projector + set projdatainfo to non-tof)
-      this->sens_proj_data_info_sptr = this->proj_data_info_sptr->create_non_tof_clone();
+      this->sens_proj_data_info_sptr = this->proj_data_info_sptr->create_single_tof_clone();
       // TODO disable caching of the matrix
       this->sens_backprojector_sptr.reset(projector_pair_sptr->get_back_projector_sptr()->clone());
       this->sens_backprojector_sptr->set_up(this->sens_proj_data_info_sptr, target_sptr);
     }
   else
     {
-      // just use the normal backprojector
+      // just use a signle TOF backprojector
       this->sens_proj_data_info_sptr = this->proj_data_info_sptr;
       this->sens_backprojector_sptr = projector_pair_sptr->get_back_projector_sptr();
     }
@@ -622,60 +622,79 @@ PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin<Tar
 {
   // TODO replace with call to distributable function
 
-  const int min_segment_num = this->proj_data_info_sptr->get_min_segment_num();
-  const int max_segment_num = this->proj_data_info_sptr->get_max_segment_num();
-
   info(boost::format("Calculating sensitivity for subset %1%") % subset_num);
 
   int min_timing_pos_num = use_tofsens ? this->proj_data_info_sptr->get_min_tof_pos_num() : 0;
   int max_timing_pos_num = use_tofsens ? this->proj_data_info_sptr->get_max_tof_pos_num() : 0;
-  if (min_timing_pos_num < 0 || max_timing_pos_num > 0)
+  if (min_timing_pos_num < 0 || max_timing_pos_num > 1)
     error("TOF code for sensitivity needs work");
 
   this->sens_backprojector_sptr->start_accumulating_in_new_target();
 
-  // warning: has to be same as subset scheme used as in distributable_computation
+  int runs = 1;
+
+  //TODO: For long scanners we should run better split this. - Let's discuss more.
+  // if((this->proj_data_info_sptr->get_max_segment_num() -
+  //      this->proj_data_info_sptr->get_min_segment_num()) > 100)
+  //   {
+  //     runs = ceil(this->proj_data_info_sptr->get_max_segment_num() -
+  //      this->proj_data_info_sptr->get_min_segment_num())/ 100;
+  //     }
+
+  info(boost::format("The  number of runs needed for the sensitivity image is %1%:  ") % runs);
+
+  for (int run = 0; run < runs; ++run)
+    {
+      const int min_segment_num = this->proj_data_info_sptr->get_min_segment_num();
+      const int max_segment_num = this->proj_data_info_sptr->get_max_segment_num();
+      info(boost::format("Current run %1% of %2%:  ") % run % runs);
+      // warning: has to be same as subset scheme used as in distributable_computation
 #ifdef STIR_OPENMP
 #  if _OPENMP < 201107
 #    pragma omp parallel for schedule(dynamic)
 #  else
 #    pragma omp parallel for collapse(2) schedule(dynamic)
 #  endif
 #endif
-  for (int segment_num = min_segment_num; segment_num <= max_segment_num; ++segment_num)
-    {
-      for (int view = this->sens_proj_data_info_sptr->get_min_view_num() + subset_num;
-           view <= this->sens_proj_data_info_sptr->get_max_view_num();
-           view += this->num_subsets)
+      for (int segment_num = min_segment_num; segment_num <= max_segment_num; ++segment_num)
         {
-          const ViewSegmentNumbers view_segment_num(view, segment_num);
+          for (int view = this->sens_proj_data_info_sptr->get_min_view_num() + subset_num;
+               view <= this->sens_proj_data_info_sptr->get_max_view_num();
+               view += this->num_subsets)
+            {
+              const ViewSegmentNumbers view_segment_num(view, segment_num);
 
-          if (!this->sens_backprojector_sptr->get_symmetries_used()->is_basic(view_segment_num))
-            continue;
-          // for (int timing_pos_num = min_timing_pos_num; timing_pos_num <= max_timing_pos_num; ++timing_pos_num)
-          {
-            shared_ptr<DataSymmetriesForViewSegmentNumbers> symmetries_used(
-                this->sens_backprojector_sptr->get_symmetries_used()->clone());
+              if (!this->sens_backprojector_sptr->get_symmetries_used()->is_basic(view_segment_num))
+                continue;
+
+              info(boost::format("Current seg %1%,  view: %1% ") % segment_num, view);
+              // for (int timing_pos_num = min_timing_pos_num; timing_pos_num <= max_timing_pos_num; ++timing_pos_num)
+              {
+                shared_ptr<DataSymmetriesForViewSegmentNumbers> symmetries_used(
+                    this->sens_backprojector_sptr->get_symmetries_used()->clone());
 
-            RelatedViewgrams<float> viewgrams = this->sens_proj_data_info_sptr->get_empty_related_viewgrams(
-                view_segment_num, symmetries_used, false); //, timing_pos_num);
+                RelatedViewgrams<float> viewgrams = this->sens_proj_data_info_sptr->get_empty_related_viewgrams(
+                    view_segment_num, symmetries_used, false); //, timing_pos_num);
 
-            viewgrams.fill(1.F);
-            // find efficiencies
-            {
-              this->normalisation_sptr->undo(viewgrams);
-            }
-            // backproject
-            {
-              const int min_ax_pos_num = viewgrams.get_min_axial_pos_num();
-              const int max_ax_pos_num = viewgrams.get_max_axial_pos_num();
+                viewgrams.fill(1.F);
+                // find efficiencies
+                {
+                  this->normalisation_sptr->undo(viewgrams);
+                }
+                // backproject
+                {
+                  const int min_ax_pos_num = viewgrams.get_min_axial_pos_num();
+                  const int max_ax_pos_num = viewgrams.get_max_axial_pos_num();
 
-              this->sens_backprojector_sptr->back_project(viewgrams, min_ax_pos_num, max_ax_pos_num);
+                  this->sens_backprojector_sptr->back_project(viewgrams, min_ax_pos_num, max_ax_pos_num);
+                }
+              }
             }
-          }
         }
+      this->sens_backprojector_sptr->get_output(sensitivity);
+      // Next run
     }
-  this->sens_backprojector_sptr->get_output(sensitivity);
+
 }
 
 template <typename TargetT>