From e30cf184a764a13d968f5c2553e424fd51b64bd4 Mon Sep 17 00:00:00 2001
From: Kris Thielemans <k.thielemans@ucl.ac.uk>
Date: Mon, 18 Sep 2023 00:23:33 +0100
Subject: [PATCH] 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.
---
.../BackProjectorByBinParallelproj.h | 4 ++
.../ForwardProjectorByBinParallelproj.h | 9 ++---
.../ProjectorByBinPairUsingParallelproj.h | 5 +++
.../BackProjectorByBinParallelproj.cxx | 24 ++++++++++--
.../ForwardProjectorByBinParallelproj.cxx | 28 ++++++++++---
.../ProjectorByBinPairUsingParallelproj.cxx | 39 ++++++++++++++++---
6 files changed, 89 insertions(+), 20 deletions(-)
diff --git a/src/include/stir/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.h b/src/include/stir/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.h
index 9038cab98e..73b9148fb6 100644
--- a/src/include/stir/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.h
+++ b/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
diff --git a/src/include/stir/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.h b/src/include/stir/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.h
index fceeb1c534..94756b4685 100644
--- a/src/include/stir/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.h
+++ b/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;
};
diff --git a/src/include/stir/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.h b/src/include/stir/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.h
index b08a2e0589..a4d12e4801 100644
--- a/src/include/stir/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.h
+++ b/src/include/stir/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.h
@@ -57,6 +57,9 @@ 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;
@@ -64,6 +67,8 @@ class ProjectorByBinPairUsingParallelproj :
void initialise_keymap() override;
bool post_processing() override;
bool _verbosity;
+ bool _restrict_to_cylindrical_FOV;
+ bool _already_set_up;
};
END_NAMESPACE_STIR
diff --git a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx
index 77291a8d68..c9a2a79481 100644
--- a/src/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.cxx
+++ b/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;
diff --git a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx
index 020a84e255..73a0ae36f2 100644
--- a/src/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.cxx
+++ b/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;
diff --git a/src/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.cxx b/src/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.cxx
index 4627ccc847..22e9cf983c 100644
--- a/src/recon_buildblock/Parallelproj_projector/ProjectorByBinPairUsingParallelproj.cxx
+++ b/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;
}