Fast reconstruction for LAFOV scanners

* Building upon the cached lm reconstruction, we further parallelize the
calculation of the sensitity image to muliple passes when the scanner is
long.

* Now, a scanner is TOF ready even if it does not have set timing
resolution. If the timing resolution is 0 ps or if the TOF mashing leads
to a scanner with 1 TOF bin then we apply boundaries on the LOR's
length, equal to the coinc window. The probabilities do not change as we
set them to 1.

* cxx17 is nessesary for newer ROOT.

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/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/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/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,77 @@ 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;
+  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;
 
-            RelatedViewgrams<float> viewgrams = this->sens_proj_data_info_sptr->get_empty_related_viewgrams(
-                view_segment_num, symmetries_used, false); //, timing_pos_num);
+              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());
 
-            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();
+                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();
 
-              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>

src/recon_buildblock/ProjMatrixByBin.cxx

@@ -68,7 +68,7 @@ ProjMatrixByBin::enable_tof(const shared_ptr<const ProjDataInfo>& _proj_data_inf
 {
   if (v)
     {
-      if (proj_data_info_sptr->get_num_tof_poss() == 1) // This is a special case that we do a nonTOF backprojection with coincidence window.
+      if (proj_data_info_sptr->get_num_tof_poss() == 1 || proj_data_info_sptr->get_scanner_ptr()->get_timing_resolution() == 0 ) // This is a special case that we do a nonTOF backprojection with coincidence window.
         {
           tof_enabled = false;
         }