allow using full FOV for Parallelproj projector

introduce restrict_to_cylindrical_FOV variable/keyword (defaulting to true),
identical to what is used in the ray-tracing matrix

This variable is a rename of the _use_truncation variable (which was essentially
disabled), but now made to work.

Author: KrisThielemans

src/include/stir/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.h

@@ -85,6 +85,9 @@ class BackProjectorByBinParallelproj :
    void set_num_gpu_chunks(int num_gpu_chunks) {_num_gpu_chunks = num_gpu_chunks; }
    int  get_num_gpu_chunks() { return _num_gpu_chunks; }
 
+  bool get_restrict_to_cylindrical_FOV() const;
+  void set_restrict_to_cylindrical_FOV(bool val);
+
 protected:
 
  virtual void actual_back_project(const RelatedViewgrams<float>&,
@@ -100,6 +103,7 @@ class BackProjectorByBinParallelproj :
   void set_helper(shared_ptr<detail::ParallelprojHelper>);
   bool _cuda_verbosity;
   int _num_gpu_chunks;
+  bool _restrict_to_cylindrical_FOV;
 };
 
 END_NAMESPACE_STIR

src/include/stir/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.h

@@ -80,14 +80,13 @@ virtual void set_up(
     /// Set verbosity
     void set_verbosity(const bool verbosity) { _cuda_verbosity = verbosity; }
 
-    /// Set use truncation - truncate before forward
-    /// projection and after back projection
-    void set_use_truncation(const bool use_truncation) { _use_truncation = use_truncation; }
-
     // set/get number of gpu chunks to use
     void set_num_gpu_chunks(int num_gpu_chunks) {_num_gpu_chunks = num_gpu_chunks; }
     int  get_num_gpu_chunks() { return _num_gpu_chunks; }
 
+  bool get_restrict_to_cylindrical_FOV() const;
+  void set_restrict_to_cylindrical_FOV(bool val);
+
 protected:
 
   virtual void actual_forward_project(RelatedViewgrams<float>& viewgrams,
@@ -102,7 +101,7 @@ virtual void set_up(
     friend class ProjectorByBinPairUsingParallelproj;
     void set_helper(shared_ptr<detail::ParallelprojHelper>);
     bool _cuda_verbosity;
-    bool _use_truncation;
+    bool _restrict_to_cylindrical_FOV;
     int _num_gpu_chunks;
 };
 

src/include/stir/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.h

@@ -57,13 +57,18 @@ class ProjectorByBinPairUsingParallelproj :
   /// Set verbosity
   void set_verbosity(const bool verbosity);
 
+  bool get_restrict_to_cylindrical_FOV() const;
+  void set_restrict_to_cylindrical_FOV(bool val);
+
 private:
   shared_ptr<detail::ParallelprojHelper> _helper;
 
   void set_defaults() override;
   void initialise_keymap() override;
   bool post_processing() override;
   bool _verbosity;
+  bool _restrict_to_cylindrical_FOV;
+  bool _already_set_up;
 };
 
 END_NAMESPACE_STIR

src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx

@@ -67,6 +67,7 @@ initialise_keymap()
   parser.add_start_key("Back Projector Using Parallelproj Parameters");
   parser.add_stop_key("End Back Projector Using Parallelproj Parameters");
   parser.add_key("verbosity", &_cuda_verbosity);
+  parser.add_key("restrict to cylindrical FOV", &_restrict_to_cylindrical_FOV);
   parser.add_key("num_gpu_chunks", &_num_gpu_chunks);
 }
 
@@ -76,8 +77,22 @@ set_defaults()
 {
   _cuda_verbosity = true;
   _num_gpu_chunks = 1;
+  _restrict_to_cylindrical_FOV = true;
 }
 
+bool
+BackProjectorByBinParallelproj::
+get_restrict_to_cylindrical_FOV() const
+{
+  return this->_restrict_to_cylindrical_FOV;
+}
+
+void
+BackProjectorByBinParallelproj::
+set_restrict_to_cylindrical_FOV(bool val)
+{
+  this->_restrict_to_cylindrical_FOV = val;
+}
 
 void
 BackProjectorByBinParallelproj::set_helper(shared_ptr<detail::ParallelprojHelper> helper)
@@ -91,6 +106,10 @@ BackProjectorByBinParallelproj::
 set_up(const shared_ptr<const ProjDataInfo>& proj_data_info_sptr,
        const shared_ptr<const DiscretisedDensity<3,float> >& density_info_sptr)
 {
+#ifdef STIR_TOF
+  if (proj_data_info_sptr->get_num_tof_poss() > 1)
+    error("STIR-ParallelProj interface does not support TOF data yet. Sorry!");
+#endif
     BackProjectorByBin::set_up(proj_data_info_sptr,density_info_sptr);
     check(*proj_data_info_sptr, *_density_sptr);
     _symmetries_sptr.reset(new TrivialDataSymmetriesForBins(proj_data_info_sptr));
@@ -201,10 +220,7 @@ get_output(DiscretisedDensity<3,float> &density) const
     // --------------------------------------------------------------- //
     std::copy(image_vec.begin(), image_vec.end(), density.begin_all());
 
-    // After the back projection, we enforce a truncation outside of the FOV.
-    // This is because the parallelproj projector seems to have some trouble at the edges and this
-    // could cause some voxel values to spiral out of control.
-    //if (_use_truncation)
+    if (this->_restrict_to_cylindrical_FOV)
       {
         const float radius = p.get_proj_data_info_sptr()->get_scanner_sptr()->get_inner_ring_radius();
         const float image_radius = _helper->voxsize[2]*_helper->imgdim[2]/2;

src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx

@@ -42,7 +42,7 @@ ForwardProjectorByBinParallelproj::registered_name =
   "Parallelproj";
 
 ForwardProjectorByBinParallelproj::ForwardProjectorByBinParallelproj() :
-    _cuda_verbosity(true), _use_truncation(true), _num_gpu_chunks(1)
+    _cuda_verbosity(true), _restrict_to_cylindrical_FOV(true), _num_gpu_chunks(1)
 {
     this->_already_set_up = false;
     this->_do_not_setup_helper = false;
@@ -59,6 +59,7 @@ initialise_keymap()
   parser.add_start_key("Forward Projector Using Parallelproj Parameters");
   parser.add_stop_key("End Forward Projector Using Parallelproj Parameters");
   parser.add_key("verbosity", &_cuda_verbosity);
+  parser.add_key("restrict to cylindrical FOV", &_restrict_to_cylindrical_FOV);
   parser.add_key("num_gpu_chunks", &_num_gpu_chunks);
 }
 
@@ -67,10 +68,24 @@ ForwardProjectorByBinParallelproj::
 set_defaults()
 {
   _cuda_verbosity = true;
-  _use_truncation = true;
+  _restrict_to_cylindrical_FOV = true;
   _num_gpu_chunks = 1;
 }
 
+bool
+ForwardProjectorByBinParallelproj::
+get_restrict_to_cylindrical_FOV() const
+{
+  return this->_restrict_to_cylindrical_FOV;
+}
+
+void
+ForwardProjectorByBinParallelproj::
+set_restrict_to_cylindrical_FOV(bool val)
+{
+  this->_restrict_to_cylindrical_FOV = val;
+}
+
 
 void
 ForwardProjectorByBinParallelproj::set_helper(shared_ptr<detail::ParallelprojHelper> helper)
@@ -84,6 +99,10 @@ ForwardProjectorByBinParallelproj::
 set_up(const shared_ptr<const ProjDataInfo>& proj_data_info_sptr,
        const shared_ptr<const DiscretisedDensity<3,float> >& density_info_sptr)
 {
+#ifdef STIR_TOF
+  if (proj_data_info_sptr->get_num_tof_poss() > 1)
+    error("STIR-ParallelProj interface does not support TOF data yet. Sorry!");
+#endif
     ForwardProjectorByBin::set_up(proj_data_info_sptr,density_info_sptr);
     check(*proj_data_info_sptr, *_density_sptr);
     _symmetries_sptr.reset(new TrivialDataSymmetriesForBins(proj_data_info_sptr));
@@ -135,10 +154,7 @@ set_input(const DiscretisedDensity<3,float> & density)
 {
     ForwardProjectorByBin::set_input(density);
 
-    // Before forward projection, we enforce a truncation outside of the FOV.
-    // This is because the parallelproj projector seems to have some trouble at the edges and this
-    // could cause some voxel values to spiral out of control.
-    //if (_use_truncation)
+    if (this->_restrict_to_cylindrical_FOV)
       {
         const float radius = this->_proj_data_info_sptr->get_scanner_sptr()->get_inner_ring_radius();
         const float image_radius = _helper->voxsize[2]*_helper->imgdim[2]/2;

src/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.cxx

@@ -39,6 +39,7 @@ ProjectorByBinPairUsingParallelproj::initialise_keymap()
   parser.add_start_key("Projector Pair Using Parallelproj Parameters");
   parser.add_stop_key("End Projector Pair Using Parallelproj Parameters");
   parser.add_key("verbosity",&_verbosity);
+  parser.add_key("restrict to cylindrical FOV", &_restrict_to_cylindrical_FOV);
 }
 
 
@@ -47,6 +48,8 @@ ProjectorByBinPairUsingParallelproj::set_defaults()
 {
   base_type::set_defaults();
   this->set_verbosity(true);
+  this->set_restrict_to_cylindrical_FOV(true);
+  this->_already_set_up = false;
 }
 
 bool
@@ -67,22 +70,48 @@ ProjectorByBinPairUsingParallelproj()
   set_defaults();
 }
 
+bool
+ProjectorByBinPairUsingParallelproj::
+get_restrict_to_cylindrical_FOV() const
+{
+  return this->_restrict_to_cylindrical_FOV;
+}
+
+void
+ProjectorByBinPairUsingParallelproj::
+set_restrict_to_cylindrical_FOV(bool val)
+{
+  this->_already_set_up = this->_already_set_up && (this->_restrict_to_cylindrical_FOV == val);
+  this->_restrict_to_cylindrical_FOV = val;
+}
+
 Succeeded
 ProjectorByBinPairUsingParallelproj::
 set_up(const shared_ptr<const ProjDataInfo>& proj_data_info_sptr,
        const shared_ptr<const DiscretisedDensity<3,float> >& image_info_sptr)
 {
-  _helper = std::make_shared<detail::ParallelprojHelper>(*proj_data_info_sptr, *image_info_sptr);
-  dynamic_pointer_cast<ForwardProjectorByBinParallelproj>(this->forward_projector_sptr)
-    ->set_helper(_helper);
-  dynamic_pointer_cast<BackProjectorByBinParallelproj>(this->back_projector_sptr)
-    ->set_helper(_helper);
+  auto fwd_prj_downcast_sptr =
+    dynamic_pointer_cast<ForwardProjectorByBinParallelproj>(this->forward_projector_sptr);
+  if (!fwd_prj_downcast_sptr)
+    error("internal error: forward projector should be ParallelProj");
+
+  auto bck_prj_downcast_sptr =
+    dynamic_pointer_cast<BackProjectorByBinParallelproj>(this->back_projector_sptr);
+  if (!bck_prj_downcast_sptr)
+    error("internal error: back projector should be ParallelProj");
+
+  bck_prj_downcast_sptr->set_restrict_to_cylindrical_FOV(this->_restrict_to_cylindrical_FOV);
+  fwd_prj_downcast_sptr->set_restrict_to_cylindrical_FOV(this->_restrict_to_cylindrical_FOV);
+  this->_helper = std::make_shared<detail::ParallelprojHelper>(*proj_data_info_sptr, *image_info_sptr);
+  fwd_prj_downcast_sptr->set_helper(this->_helper);
+  bck_prj_downcast_sptr->set_helper(this->_helper);
 
   // the forward_projector->set_up etc will be called in the base class
 
   if (base_type::set_up(proj_data_info_sptr, image_info_sptr) != Succeeded::yes)
     return Succeeded::no;
 
+  this->_already_set_up = true;
   return Succeeded::yes;
 }