allow using full FOV for Parallelproj projector

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

@@ -95,6 +95,9 @@ class BackProjectorByBinParallelproj : public RegisteredParsingObject<BackProjec
 
   BackProjectorByBinParallelproj* clone() const override;
 
+  bool get_restrict_to_cylindrical_FOV() const;
+  void set_restrict_to_cylindrical_FOV(bool val);
+
 protected:
   void actual_back_project(const RelatedViewgrams<float>&,
                            const int min_axial_pos_num,
@@ -111,6 +114,7 @@ class BackProjectorByBinParallelproj : public RegisteredParsingObject<BackProjec
   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 @@ class ForwardProjectorByBinParallelproj : public RegisteredParsingObject<Forward
   /// 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:
   void actual_forward_project(RelatedViewgrams<float>& viewgrams,
                               const int min_axial_pos_num,
@@ -103,7 +102,7 @@ class ForwardProjectorByBinParallelproj : public RegisteredParsingObject<Forward
   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

@@ -53,13 +53,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

@@ -63,6 +63,7 @@ BackProjectorByBinParallelproj::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);
 }
 
@@ -71,6 +72,19 @@ BackProjectorByBinParallelproj::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
@@ -289,15 +303,12 @@ BackProjectorByBinParallelproj::get_output(DiscretisedDensity<3, float>& density
       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)
-  {
-    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;
-    truncate_rim(density, static_cast<int>(std::max((image_radius - radius) / _helper->voxsize[2], 0.F)));
-  }
+  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;
+      truncate_rim(density, static_cast<int>(std::max((image_radius - radius) / _helper->voxsize[2], 0.F)));
+    }
 }
 
 void

src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx

@@ -43,7 +43,7 @@ const char* const ForwardProjectorByBinParallelproj::registered_name = "Parallel
 
 ForwardProjectorByBinParallelproj::ForwardProjectorByBinParallelproj()
     : _cuda_verbosity(true),
-      _use_truncation(true),
+      _restrict_to_cylindrical_FOV(true),
       _num_gpu_chunks(1)
 {
   this->_already_set_up = false;
@@ -59,17 +59,30 @@ ForwardProjectorByBinParallelproj::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);
 }
 
 void
 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)
 {
@@ -144,15 +157,12 @@ ForwardProjectorByBinParallelproj::set_input(const DiscretisedDensity<3, float>&
 {
   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)
-  {
-    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;
-    truncate_rim(*_density_sptr, static_cast<int>(std::max((image_radius - radius) / _helper->voxsize[2], 0.F)));
-  }
+  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;
+      truncate_rim(*_density_sptr, static_cast<int>(std::max((image_radius - radius) / _helper->voxsize[2], 0.F)));
+    }
 
   std::vector<float> image_vec;
   float* image_ptr;

src/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.cxx

@@ -34,13 +34,16 @@ 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);
 }
 
 void
 ProjectorByBinPairUsingParallelproj::set_defaults()
 {
   base_type::set_defaults();
   this->set_verbosity(true);
+  this->set_restrict_to_cylindrical_FOV(true);
+  this->_already_set_up = false;
 }
 
 bool
@@ -60,6 +63,19 @@ ProjectorByBinPairUsingParallelproj::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;
+}
+
 BackProjectorByBinParallelproj*
 BackProjectorByBinParallelproj::clone() const
 {
@@ -70,15 +86,26 @@ 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;
 }