diff --git a/documentation/release_6.3.htm b/documentation/release_6.3.htm
new file mode 100644
index 000000000..fc818f8db
--- /dev/null
+++ b/documentation/release_6.3.htm
@@ -0,0 +1,66 @@
+
+
+
+ Summary of changes in STIR release 6.3
+
+
+
+ Summary of changes in STIR release 6.3
+
+
+Overall summary
+
+
+Patch release info
+
+
+ Summary for end users (also to be read by developers)
+
+
+New functionality
+
+ -
+
ScatterSimulation can now downsample the scanner transaxially (crystals per ring) for BlocksOnCylindrical,
+ scanners, which speeds up ScatterEstimation considerably. By default, downsampling the detectors per reading
+ is disabled for backwards compatibility.
+ PR #1291
+
+
+
+
+Changed functionality
+
+
+Bug fixes
+
+
+Build system
+
+
+Known problems
+ See our issue tracker.
+
+
+What is new for developers (aside from what should be obvious from the above):
+
+
+Changed functionality
+
+
+Bug fixes
+
+
+Other code changes
+
+
+Test changes
+
+
+C++ tests
+
+
+recon_test_pack
+
+
+
+
diff --git a/src/buildblock/extend_projdata.cxx b/src/buildblock/extend_projdata.cxx
index 1ecd11c15..59fa14288 100644
--- a/src/buildblock/extend_projdata.cxx
+++ b/src/buildblock/extend_projdata.cxx
@@ -93,7 +93,7 @@ extend_segment(const SegmentBySinogram& segment,
if (extend_without_wrapping)
{
out[axial_pos][min_dim[2] + view_edge] = out[axial_pos][min_dim[2] + view_extension];
- out[axial_pos][max_dim[2] - view_extension] = out[axial_pos][max_dim[2] - view_extension];
+ out[axial_pos][max_dim[2] - view_edge] = out[axial_pos][max_dim[2] - view_extension];
}
else if (flip_views)
{
diff --git a/src/buildblock/interpolate_projdata.cxx b/src/buildblock/interpolate_projdata.cxx
index 060e8aa35..8e0eee6a8 100644
--- a/src/buildblock/interpolate_projdata.cxx
+++ b/src/buildblock/interpolate_projdata.cxx
@@ -23,6 +23,7 @@
//#include "stir/display.h"
#include "stir/ProjDataInfo.h"
#include "stir/ProjDataInfoCylindricalNoArcCorr.h"
+#include "stir/ProjDataInfoGenericNoArcCorr.h"
#include "stir/IndexRange.h"
#include "stir/BasicCoordinate.h"
#include "stir/Sinogram.h"
@@ -191,6 +192,10 @@ interpolate_projdata(ProjData& proj_data_out,
error("interpolate_projdata needs both projection to be of a scanner with the same ring radius");
}
+ // handle BlocksOnCylindrical interpolation manually
+ if (proj_data_in_info.get_scanner_sptr()->get_scanner_geometry() != "Cylindrical")
+ return interpolate_blocks_on_cylindrical_projdata(proj_data_out, proj_data_in, remove_interleaving);
+
// initialise interpolator
BSpline::BSplinesRegularGrid<3, float, float> proj_data_interpolator(these_types);
for (int k = proj_data_out_info.get_min_tof_pos_num(); k <= proj_data_out_info.get_max_tof_pos_num(); ++k)
@@ -200,96 +205,319 @@ interpolate_projdata(ProjData& proj_data_out,
proj_data_in.get_segment_by_sinogram(0, k))
: proj_data_in.get_segment_by_sinogram(0, k);
+ // for Cylindrical, spacing is regular in all directions, which makes mapping trivial
std::function(const BasicCoordinate<3, int>&)> index_converter;
- if (proj_data_in_info.get_scanner_sptr()->get_scanner_geometry() == "Cylindrical")
- { // for Cylindrical, spacing is regular in all directions, which makes mapping trivial
- // especially in view direction, extending by 5 leads to much smaller artifacts
- proj_data_interpolator.set_coef(extend_segment(segment, 5, 5, 5));
-
- BasicCoordinate<3, double> offset, step;
- // out_index * step + offset = in_index
- const float in_sampling_m = proj_data_in_info.get_sampling_in_m(Bin(0, 0, 0, 0));
- const float out_sampling_m = proj_data_out_info.get_sampling_in_m(Bin(0, 0, 0, 0));
- // offset in 'in' index units
- offset[1] = (proj_data_out_info.get_m(Bin(0, 0, 0, 0)) - proj_data_in_info.get_m(Bin(0, 0, 0, 0))) / in_sampling_m;
- step[1] = out_sampling_m / in_sampling_m;
-
- const float in_sampling_phi = (proj_data_in_info.get_phi(Bin(0, 1, 0, 0)) - proj_data_in_info.get_phi(Bin(0, 0, 0, 0)))
- / (remove_interleaving ? 2 : 1);
- const float out_sampling_phi
- = proj_data_out_info.get_phi(Bin(0, 1, 0, 0)) - proj_data_out_info.get_phi(Bin(0, 0, 0, 0));
- offset[2]
- = (proj_data_out_info.get_phi(Bin(0, 0, 0, 0)) - proj_data_in_info.get_phi(Bin(0, 0, 0, 0))) / in_sampling_phi;
- step[2] = out_sampling_phi / in_sampling_phi;
-
- const float in_sampling_s = proj_data_in_info.get_sampling_in_s(Bin(0, 0, 0, 0));
- const float out_sampling_s = proj_data_out_info.get_sampling_in_s(Bin(0, 0, 0, 0));
- offset[3] = (proj_data_out_info.get_s(Bin(0, 0, 0, 0)) - proj_data_in_info.get_s(Bin(0, 0, 0, 0))) / in_sampling_s;
- step[3] = out_sampling_s / in_sampling_s;
-
- // define a function to translate indices in the output proj data to indices in input proj data
- index_converter
- = [&proj_data_out_info, offset, step](const BasicCoordinate<3, int>& index_out) -> BasicCoordinate<3, double> {
- // translate to indices in input proj data
- BasicCoordinate<3, double> index_in;
- for (auto dim = 1; dim <= 3; dim++)
- index_in[dim] = index_out[dim] * step[dim] + offset[dim];
-
- return index_in;
- };
- }
- else
- { // for BlocksOnCylindrical, views and tangential positions are not subsampled and can be mapped 1:1
- if (proj_data_in_info.get_num_tangential_poss() != proj_data_out_info.get_num_tangential_poss())
- {
- error("Interpolation of BlocksOnCylindrical scanners assumes that number of tangential positions "
- "is the same in the downsampled scanner.");
- }
- if (proj_data_in_info.get_num_views() != proj_data_out_info.get_num_views())
- {
- error("Interpolation of BlocksOnCylindrical scanners assumes that number of views "
- "is the same in the downsampled scanner.");
- }
+ // especially in view direction, extending by 5 leads to much smaller artifacts
+ proj_data_interpolator.set_coef(extend_segment(segment, 5, 5, 5));
+
+ BasicCoordinate<3, double> offset, step;
+ // out_index * step + offset = in_index
+ const float in_sampling_m = proj_data_in_info.get_sampling_in_m(Bin(0, 0, 0, 0));
+ const float out_sampling_m = proj_data_out_info.get_sampling_in_m(Bin(0, 0, 0, 0));
+ // offset in 'in' index units
+ offset[1] = (proj_data_out_info.get_m(Bin(0, 0, 0, 0)) - proj_data_in_info.get_m(Bin(0, 0, 0, 0))) / in_sampling_m;
+ step[1] = out_sampling_m / in_sampling_m;
+
+ const float in_sampling_phi = (proj_data_in_info.get_phi(Bin(0, 1, 0, 0)) - proj_data_in_info.get_phi(Bin(0, 0, 0, 0)))
+ / (remove_interleaving ? 2 : 1);
+ const float out_sampling_phi = proj_data_out_info.get_phi(Bin(0, 1, 0, 0)) - proj_data_out_info.get_phi(Bin(0, 0, 0, 0));
+ offset[2] = (proj_data_out_info.get_phi(Bin(0, 0, 0, 0)) - proj_data_in_info.get_phi(Bin(0, 0, 0, 0))) / in_sampling_phi;
+ step[2] = out_sampling_phi / in_sampling_phi;
+
+ const float in_sampling_s = proj_data_in_info.get_sampling_in_s(Bin(0, 0, 0, 0));
+ const float out_sampling_s = proj_data_out_info.get_sampling_in_s(Bin(0, 0, 0, 0));
+ offset[3] = (proj_data_out_info.get_s(Bin(0, 0, 0, 0)) - proj_data_in_info.get_s(Bin(0, 0, 0, 0))) / in_sampling_s;
+ step[3] = out_sampling_s / in_sampling_s;
+
+ // define a function to translate indices in the output proj data to indices in input proj data
+ index_converter
+ = [&proj_data_out_info, offset, step](const BasicCoordinate<3, int>& index_out) -> BasicCoordinate<3, double> {
+ // translate to indices in input proj data
+ BasicCoordinate<3, double> index_in;
+ for (auto dim = 1; dim <= 3; dim++)
+ index_in[dim] = index_out[dim] * step[dim] + offset[dim];
+
+ return index_in;
+ };
- // only extending in axial direction - an extension of 2 was found to be sufficient
- proj_data_interpolator.set_coef(extend_segment(segment, 0, 2, 0));
+ SegmentBySinogram sino_3D_out = proj_data_out.get_empty_segment_by_sinogram(0, false, k);
+ sample_function_using_index_converter(sino_3D_out, proj_data_interpolator, index_converter);
- auto m_offset = proj_data_in_info.get_m(Bin(0, 0, 0, 0));
- auto m_sampling = proj_data_in_info.get_sampling_in_m(Bin(0, 0, 0, 0));
+ if (proj_data_out.set_segment(sino_3D_out) == Succeeded::no)
+ return Succeeded::no;
+ }
+ return Succeeded::yes;
+}
- // confirm that proj_data_in has equidistant sampling in m
- for (auto axial_pos = proj_data_in_info.get_min_axial_pos_num(0);
- axial_pos <= proj_data_in_info.get_max_axial_pos_num(0);
- axial_pos++)
- {
- if (abs(m_sampling - proj_data_in_info.get_sampling_in_m(Bin(0, 0, axial_pos, 0))) > 1E-4)
- error("input projdata to interpolate_projdata are not equidistantly sampled in m.");
- }
+//! This function interpolates BlocksOnCylindrical proj data taking bucket intersections and gaps into account.
+/*! The above proj data interpolation function works well for cylindrical scanners where there are no large discontinuities.
+For BlocksOnCylindrical scanners, the sudden change of angle from one bucket to the next and the gaps between blocks and buckets
+lead to interpolation artifacts if not taken into account. Therefore, this function does the interpolation in physical space
+rather than directly on the proj data bin values. For each bin in proj_data_out (the full size proj data), we find the four
+closest LORs in the downsampled proj_data_in. These are then weighted using bilinear interpolation based on the crystal positions
+of the two endpoints of the LOR.
+*/
+Succeeded
+interpolate_blocks_on_cylindrical_projdata(ProjData& proj_data_out, const ProjData& proj_data_in, bool remove_interleaving)
+{
+ const ProjDataInfo& proj_data_in_info = *proj_data_in.get_proj_data_info_sptr();
+ const ProjDataInfo& proj_data_out_info = *proj_data_out.get_proj_data_info_sptr();
- // define a function to translate indices in the output proj data to indices in input proj data
- index_converter = [&proj_data_out_info, m_offset, m_sampling](
- const BasicCoordinate<3, int>& index_out) -> BasicCoordinate<3, double> {
- // translate index on output to coordinate
- auto bin
- = Bin(0 /* segment */, index_out[2] /* view */, index_out[1] /* axial pos */, index_out[3] /* tangential pos */);
- auto out_m = proj_data_out_info.get_m(bin);
-
- // translate to indices in input proj data
- BasicCoordinate<3, double> index_in;
- index_in[1] = (out_m - m_offset) / m_sampling;
- index_in[2] = index_out[2];
- index_in[3] = index_out[3];
-
- return index_in;
- };
- }
+ auto m_offset = proj_data_in_info.get_m(Bin(0, 0, 0, 0));
+ auto m_sampling = proj_data_in_info.get_sampling_in_m(Bin(0, 0, 0, 0));
+ for (int k = proj_data_out_info.get_min_tof_pos_num(); k <= proj_data_out_info.get_max_tof_pos_num(); ++k)
+ {
+ SegmentBySinogram segment
+ = remove_interleaving ? make_non_interleaved_segment(*(make_non_interleaved_proj_data_info(proj_data_in_info)),
+ proj_data_in.get_segment_by_sinogram(0, k))
+ : proj_data_in.get_segment_by_sinogram(0, k);
SegmentBySinogram sino_3D_out = proj_data_out.get_empty_segment_by_sinogram(0, false, k);
- sample_function_using_index_converter(sino_3D_out, proj_data_interpolator, index_converter);
+ BasicCoordinate<3, int> min_out, max_out;
+ IndexRange<3> out_range = sino_3D_out.get_index_range();
+ if (!out_range.get_regular_range(min_out, max_out))
+ warning("Output must be regular range!");
+
+ // confirm that proj_data_in has equidistant sampling in m
+ for (auto axial_pos = proj_data_in_info.get_min_axial_pos_num(0); axial_pos <= proj_data_in_info.get_max_axial_pos_num(0);
+ axial_pos++)
+ {
+ if (abs(m_sampling - proj_data_in_info.get_sampling_in_m(Bin(0, 0, axial_pos, 0))) > 1E-4)
+ error("input projdata to interpolate_projdata are not equidistantly sampled in m.");
+ }
+
+ BasicCoordinate<3, int> index_out;
+ for (index_out[1] = min_out[1]; index_out[1] <= max_out[1]; ++index_out[1])
+ {
+ for (index_out[2] = min_out[2]; index_out[2] <= max_out[2]; ++index_out[2])
+ {
+ for (index_out[3] = min_out[3]; index_out[3] <= max_out[3]; ++index_out[3])
+ {
+ // translate index on output to coordinate
+ auto bin
+ = Bin(0 /* segment */, index_out[2] /* view */,
+ index_out[1] /* axial pos */, index_out[3] /* tangential pos
+ */);
+ auto out_m = proj_data_out_info.get_m(bin);
+ double axial_idx = (out_m - m_offset) / m_sampling;
+ int axial_floor = std::floor(axial_idx);
+ int axial_ceil = std::ceil(axial_idx);
+ if (axial_floor == axial_ceil)
+ {
+ if (axial_floor == 0)
+ axial_ceil++;
+ else
+ axial_floor--;
+ }
+
+ // find the two crystals for this bin and on which buckets (modules) they are
+ int det1_num_out, det2_num_out;
+ const auto proj_data_out_info_ptr = dynamic_cast(&proj_data_out_info);
+ proj_data_out_info_ptr->get_det_num_pair_for_view_tangential_pos_num(det1_num_out,
+ det2_num_out,
+ index_out[2], /*view*/
+ index_out[3] /* tangential pos */);
+ const int dets_per_module_out = proj_data_out_info.get_scanner_sptr()->get_num_transaxial_crystals_per_bucket();
+ const int det1_module = det1_num_out / dets_per_module_out;
+ const int det2_module = det2_num_out / dets_per_module_out;
+
+ // translate the crystal position on each module from the full size scanner to the downsampled scanner
+ const auto proj_data_in_info_ptr = dynamic_cast(&proj_data_in_info);
+ const int dets_per_module_in
+ = proj_data_in_info_ptr->get_scanner_sptr()->get_num_transaxial_crystals_per_bucket();
+ const int crystal1_out_module_idx = det1_num_out % dets_per_module_out;
+ const double crystal1_out_module_pos
+ = std::floor(static_cast(crystal1_out_module_idx)
+ / proj_data_out_info_ptr->get_scanner_sptr()->get_num_transaxial_crystals_per_block())
+ * proj_data_out_info_ptr->get_scanner_sptr()->get_transaxial_block_spacing()
+ + static_cast(
+ crystal1_out_module_idx
+ % proj_data_out_info_ptr->get_scanner_sptr()->get_num_transaxial_crystals_per_block())
+ * proj_data_out_info_ptr->get_scanner_sptr()->get_transaxial_crystal_spacing();
+ const double crystal1_num_in
+ = det1_module * dets_per_module_in
+ + crystal1_out_module_pos / proj_data_in_info.get_scanner_sptr()->get_transaxial_crystal_spacing();
+ const int crystal2_out_module_idx = det2_num_out % dets_per_module_out;
+ const double crystal2_out_module_pos
+ = std::floor(static_cast(crystal2_out_module_idx)
+ / proj_data_out_info_ptr->get_scanner_sptr()->get_num_transaxial_crystals_per_block())
+ * proj_data_out_info_ptr->get_scanner_sptr()->get_transaxial_block_spacing()
+ + static_cast(
+ crystal2_out_module_idx
+ % proj_data_out_info_ptr->get_scanner_sptr()->get_num_transaxial_crystals_per_block())
+ * proj_data_out_info_ptr->get_scanner_sptr()->get_transaxial_crystal_spacing();
+ const double crystal2_num_in
+ = det2_module * dets_per_module_in
+ + crystal2_out_module_pos / proj_data_in_info.get_scanner_sptr()->get_transaxial_crystal_spacing();
+ const auto crystal1_num_in_floor
+ = std::max(static_cast(std::floor(crystal1_num_in)), det1_module * dets_per_module_in);
+ const auto crystal1_num_in_ceil
+ = std::min(static_cast(std::ceil(crystal1_num_in)), (det1_module + 1) * dets_per_module_in - 1);
+ const auto crystal2_num_in_floor
+ = std::max(static_cast(std::floor(crystal2_num_in)), det2_module * dets_per_module_in);
+ const auto crystal2_num_in_ceil
+ = std::min(static_cast(std::ceil(crystal2_num_in)), (det2_module + 1) * dets_per_module_in - 1);
+
+ // translate to indices in input proj data
+ BasicCoordinate<3, int> index_in;
+
+ // in this case we can skip parts of the interpolation
+ if (crystal1_num_in_floor == crystal1_num_in_ceil)
+ {
+ if (crystal2_num_in_floor == crystal2_num_in_ceil)
+ {
+ int one_view, one_tang;
+ proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair(
+ one_view, one_tang, crystal1_num_in_floor, crystal2_num_in_floor);
+ one_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), one_tang),
+ proj_data_in_info_ptr->get_max_tangential_pos_num());
+
+ index_in[1] = std::max(axial_floor, proj_data_in_info_ptr->get_min_axial_pos_num(0));
+ index_in[2] = one_view;
+ index_in[3] = one_tang;
+ sino_3D_out[index_out] += segment[index_in] * (axial_ceil - axial_idx);
+
+ index_in[1] = std::min(axial_ceil, proj_data_in_info_ptr->get_max_axial_pos_num(0));
+ sino_3D_out[index_out] += segment[index_in] * (axial_idx - axial_floor);
+ }
+ else
+ {
+ int ff_view, fc_view;
+ int ff_tang, fc_tang;
+ proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair(
+ ff_view, ff_tang, crystal1_num_in_floor, crystal2_num_in_floor);
+ proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair(
+ fc_view, fc_tang, crystal1_num_in_floor, crystal2_num_in_ceil);
+
+ ff_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), ff_tang),
+ proj_data_in_info_ptr->get_max_tangential_pos_num());
+ fc_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), fc_tang),
+ proj_data_in_info_ptr->get_max_tangential_pos_num());
+
+ index_in[1] = std::max(axial_floor, proj_data_in_info_ptr->get_min_axial_pos_num(0));
+ index_in[2] = ff_view;
+ index_in[3] = ff_tang;
+ sino_3D_out[index_out]
+ += segment[index_in] * (axial_ceil - axial_idx) * (crystal2_num_in_ceil - crystal2_num_in);
+ index_in[2] = fc_view;
+ index_in[3] = fc_tang;
+ sino_3D_out[index_out]
+ += segment[index_in] * (axial_ceil - axial_idx) * (crystal2_num_in - crystal2_num_in_floor);
+ index_in[1] = std::min(axial_ceil, proj_data_in_info_ptr->get_max_axial_pos_num(0));
+ sino_3D_out[index_out]
+ += segment[index_in] * (axial_idx - axial_floor) * (crystal2_num_in - crystal2_num_in_floor);
+ index_in[2] = ff_view;
+ index_in[3] = ff_tang;
+ sino_3D_out[index_out]
+ += segment[index_in] * (axial_idx - axial_floor) * (crystal2_num_in_ceil - crystal2_num_in);
+ }
+ }
+ else if (crystal2_num_in_floor == crystal2_num_in_ceil)
+ {
+ int ff_view, cf_view;
+ int ff_tang, cf_tang;
+ proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair(
+ ff_view, ff_tang, crystal1_num_in_floor, crystal2_num_in_floor);
+ proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair(
+ cf_view, cf_tang, crystal1_num_in_ceil, crystal2_num_in_floor);
+
+ ff_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), ff_tang),
+ proj_data_in_info_ptr->get_max_tangential_pos_num());
+ cf_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), cf_tang),
+ proj_data_in_info_ptr->get_max_tangential_pos_num());
+
+ index_in[1] = std::max(axial_floor, proj_data_in_info_ptr->get_min_axial_pos_num(0));
+ index_in[2] = ff_view;
+ index_in[3] = ff_tang;
+ sino_3D_out[index_out]
+ += segment[index_in] * (axial_ceil - axial_idx) * (crystal1_num_in_ceil - crystal1_num_in);
+ index_in[2] = cf_view;
+ index_in[3] = cf_tang;
+ sino_3D_out[index_out]
+ += segment[index_in] * (axial_ceil - axial_idx) * (crystal1_num_in - crystal1_num_in_floor);
+ index_in[1] = std::min(axial_ceil, proj_data_in_info_ptr->get_max_axial_pos_num(0));
+ sino_3D_out[index_out]
+ += segment[index_in] * (axial_idx - axial_floor) * (crystal1_num_in - crystal1_num_in_floor);
+ index_in[2] = ff_view;
+ index_in[3] = ff_tang;
+ sino_3D_out[index_out]
+ += segment[index_in] * (axial_idx - axial_floor) * (crystal1_num_in_ceil - crystal1_num_in);
+ }
+ else // in this case we need to do a trilinear interpolation
+ {
+ int ff_view, fc_view, cf_view, cc_view;
+ int ff_tang, fc_tang, cf_tang, cc_tang;
+ proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair(
+ ff_view, ff_tang, crystal1_num_in_floor, crystal2_num_in_floor);
+ proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair(
+ fc_view, fc_tang, crystal1_num_in_floor, crystal2_num_in_ceil);
+ proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair(
+ cf_view, cf_tang, crystal1_num_in_ceil, crystal2_num_in_floor);
+ proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair(
+ cc_view, cc_tang, crystal1_num_in_ceil, crystal2_num_in_ceil);
+
+ // TODO: why can we get positions out that are not even in the proj data?!
+ ff_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), ff_tang),
+ proj_data_in_info_ptr->get_max_tangential_pos_num());
+ fc_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), fc_tang),
+ proj_data_in_info_ptr->get_max_tangential_pos_num());
+ cf_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), cf_tang),
+ proj_data_in_info_ptr->get_max_tangential_pos_num());
+ cc_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), cc_tang),
+ proj_data_in_info_ptr->get_max_tangential_pos_num());
+
+ index_in[1] = std::max(axial_floor, proj_data_in_info_ptr->get_min_axial_pos_num(0));
+ index_in[2] = ff_view;
+ index_in[3] = ff_tang;
+ sino_3D_out[index_out] += segment[index_in] * (axial_ceil - axial_idx)
+ * (crystal1_num_in_ceil - crystal1_num_in)
+ * (crystal2_num_in_ceil - crystal2_num_in);
+ index_in[2] = fc_view;
+ index_in[3] = fc_tang;
+ sino_3D_out[index_out] += segment[index_in] * (axial_ceil - axial_idx)
+ * (crystal1_num_in_ceil - crystal1_num_in)
+ * (crystal2_num_in - crystal2_num_in_floor);
+ index_in[2] = cc_view;
+ index_in[3] = cc_tang;
+ sino_3D_out[index_out] += segment[index_in] * (axial_ceil - axial_idx)
+ * (crystal1_num_in - crystal1_num_in_floor)
+ * (crystal2_num_in - crystal2_num_in_floor);
+ index_in[2] = cf_view;
+ index_in[3] = cf_tang;
+ sino_3D_out[index_out] += segment[index_in] * (axial_ceil - axial_idx)
+ * (crystal1_num_in - crystal1_num_in_floor)
+ * (crystal2_num_in_ceil - crystal2_num_in);
+
+ index_in[1] = std::min(axial_ceil, proj_data_in_info_ptr->get_max_axial_pos_num(0));
+ index_in[2] = ff_view;
+ index_in[3] = ff_tang;
+ sino_3D_out[index_out] += segment[index_in] * (axial_idx - axial_floor)
+ * (crystal1_num_in_ceil - crystal1_num_in)
+ * (crystal2_num_in_ceil - crystal2_num_in);
+ index_in[2] = fc_view;
+ index_in[3] = fc_tang;
+ sino_3D_out[index_out] += segment[index_in] * (axial_idx - axial_floor)
+ * (crystal1_num_in_ceil - crystal1_num_in)
+ * (crystal2_num_in - crystal2_num_in_floor);
+ index_in[2] = cc_view;
+ index_in[3] = cc_tang;
+ sino_3D_out[index_out] += segment[index_in] * (axial_idx - axial_floor)
+ * (crystal1_num_in - crystal1_num_in_floor)
+ * (crystal2_num_in - crystal2_num_in_floor);
+ index_in[2] = cf_view;
+ index_in[3] = cf_tang;
+ sino_3D_out[index_out] += segment[index_in] * (axial_idx - axial_floor)
+ * (crystal1_num_in - crystal1_num_in_floor)
+ * (crystal2_num_in_ceil - crystal2_num_in);
+ }
+ }
+ }
+ }
if (proj_data_out.set_segment(sino_3D_out) == Succeeded::no)
return Succeeded::no;
}
+
return Succeeded::yes;
}
diff --git a/src/include/stir/interpolate_projdata.h b/src/include/stir/interpolate_projdata.h
index 43889ca8a..0fe97060c 100644
--- a/src/include/stir/interpolate_projdata.h
+++ b/src/include/stir/interpolate_projdata.h
@@ -59,6 +59,8 @@ Succeeded interpolate_projdata(ProjData& proj_data_out,
const ProjData& proj_data_in,
const BasicCoordinate<3, BSpline::BSplineType>& these_types,
const bool remove_interleaving);
+Succeeded
+interpolate_blocks_on_cylindrical_projdata(ProjData& proj_data_out, const ProjData& proj_data_in, bool remove_interleaving);
//@}
END_NAMESPACE_STIR
diff --git a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl
index 16f92a2e7..661744d50 100644
--- a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl
+++ b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl
@@ -216,6 +216,7 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_bin0(int segment_num, int v
// No symmetry is implemented for generic scanner
return new TrivialSymmetryOperation();
}
+ return new TrivialSymmetryOperation();
}
// from symmetries
@@ -390,6 +391,8 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_general_bin(int s, int segm
// No symmetry is implemented for generic scanner
return new TrivialSymmetryOperation();
}
+
+ return new TrivialSymmetryOperation();
}
bool
diff --git a/src/include/stir/scatter/ScatterEstimation.h b/src/include/stir/scatter/ScatterEstimation.h
index 38df94df4..785cd3b25 100644
--- a/src/include/stir/scatter/ScatterEstimation.h
+++ b/src/include/stir/scatter/ScatterEstimation.h
@@ -174,6 +174,12 @@ class ScatterEstimation : public ParsingObject
inline void set_num_iterations(int);
+ inline unsigned int get_half_filter_width() const;
+ inline void set_half_filter_width(unsigned int);
+
+ inline void
+ set_downsample_scanner(bool downsample_scanner, int downsampled_number_of_rings = -1, int downsampled_detectors_per_ring = -1);
+
void set_output_scatter_estimate_prefix(const std::string&);
void set_export_scatter_estimates_of_each_iteration(bool);
@@ -382,6 +388,9 @@ class ScatterEstimation : public ParsingObject
float min_scale_value;
bool downsample_scanner_bool;
+ int downsampled_number_of_rings;
+ int downsampled_detectors_per_ring;
+
//!
unsigned int half_filter_width;
diff --git a/src/include/stir/scatter/ScatterEstimation.inl b/src/include/stir/scatter/ScatterEstimation.inl
index 73336a24d..dc8bd9f2b 100644
--- a/src/include/stir/scatter/ScatterEstimation.inl
+++ b/src/include/stir/scatter/ScatterEstimation.inl
@@ -66,4 +66,26 @@ ScatterEstimation::set_num_iterations(int arg)
this->num_scatter_iterations = arg;
}
+unsigned int
+ScatterEstimation::get_half_filter_width() const
+{
+ return this->half_filter_width;
+}
+
+void
+ScatterEstimation::set_half_filter_width(unsigned int arg)
+{
+ this->half_filter_width = arg;
+}
+
+void
+ScatterEstimation::set_downsample_scanner(bool downsample_scanner,
+ int downsampled_number_of_rings,
+ int downsampled_detectors_per_ring)
+{
+ this->downsample_scanner_bool = downsample_scanner;
+ this->downsampled_number_of_rings = downsampled_number_of_rings;
+ this->downsampled_detectors_per_ring = downsampled_detectors_per_ring;
+}
+
END_NAMESPACE_STIR
diff --git a/src/scatter_buildblock/ScatterEstimation.cxx b/src/scatter_buildblock/ScatterEstimation.cxx
index 3ceb4cf1d..07772ae99 100644
--- a/src/scatter_buildblock/ScatterEstimation.cxx
+++ b/src/scatter_buildblock/ScatterEstimation.cxx
@@ -23,6 +23,7 @@
#include "stir/recon_buildblock/ChainedBinNormalisation.h"
#include "stir/ProjDataInterfile.h"
#include "stir/ProjDataInMemory.h"
+#include "stir/inverse_SSRB.h"
#include "stir/ExamInfo.h"
#include "stir/ProjDataInfo.h"
#include "stir/ProjDataInfoCylindricalNoArcCorr.h"
@@ -77,6 +78,8 @@ ScatterEstimation::set_defaults()
this->override_scanner_template = true;
this->override_density_image = true;
this->downsample_scanner_bool = true;
+ this->downsampled_number_of_rings = -1;
+ this->downsampled_detectors_per_ring = -1;
this->remove_interleaving = true;
this->atten_image_filename = "";
this->atten_coeff_filename = "";
@@ -124,6 +127,8 @@ ScatterEstimation::initialise_keymap()
this->parser.add_parsing_key("Scatter Simulation type", &this->scatter_simulation_sptr);
this->parser.add_key("scatter simulation parameter filename", &this->scatter_sim_par_filename);
this->parser.add_key("use scanner downsampling in scatter simulation", &this->downsample_scanner_bool);
+ this->parser.add_key("override number of downsampled rings", &this->downsampled_number_of_rings);
+ this->parser.add_key("override number of downsampled detectors per ring", &this->downsampled_detectors_per_ring);
this->parser.add_key("override attenuation image", &this->override_density_image);
this->parser.add_key("override scanner template", &this->override_scanner_template);
@@ -597,7 +602,7 @@ ScatterEstimation::set_up()
}
if (this->downsample_scanner_bool)
- this->scatter_simulation_sptr->downsample_scanner();
+ this->scatter_simulation_sptr->downsample_scanner(this->downsampled_number_of_rings, this->downsampled_detectors_per_ring);
// Check if Load a mask proj_data
@@ -834,7 +839,7 @@ ScatterEstimation::process_data()
float local_min_scale_value = 0.5f;
float local_max_scale_value = 0.5f;
- stir::BSpline::BSplineType spline_type = stir::BSpline::quadratic;
+ stir::BSpline::BSplineType spline_type = stir::BSpline::linear;
// This has been set to 2D or 3D in the set_up()
shared_ptr unscaled_est_projdata_sptr(
@@ -1026,16 +1031,13 @@ ScatterEstimation::process_data()
shared_ptr normalisation_factors_3d_sptr
= this->get_normalisation_object_sptr(this->multiplicative_binnorm_sptr);
- upsample_and_fit_scatter_estimate(*scatter_estimate_sptr,
- *this->input_projdata_sptr,
- *temp_projdata,
- *normalisation_factors_3d_sptr,
- *this->input_projdata_sptr,
- 1.0f,
- 1.0f,
- 1,
- spline_type,
- false);
+ ProjDataInMemory interpolated_scatter(this->input_projdata_sptr->get_exam_info_sptr(),
+ this->input_projdata_sptr->get_proj_data_info_sptr()->create_shared_clone());
+ inverse_SSRB(interpolated_scatter, *temp_projdata);
+ normalisation_factors_3d_sptr->set_up(this->input_projdata_sptr->get_exam_info_sptr(),
+ this->input_projdata_sptr->get_proj_data_info_sptr()->create_shared_clone());
+ normalisation_factors_3d_sptr->undo(interpolated_scatter);
+ scatter_estimate_sptr->fill(interpolated_scatter);
}
else
{
diff --git a/src/scatter_buildblock/ScatterSimulation.cxx b/src/scatter_buildblock/ScatterSimulation.cxx
index 8411b3f53..522178b83 100644
--- a/src/scatter_buildblock/ScatterSimulation.cxx
+++ b/src/scatter_buildblock/ScatterSimulation.cxx
@@ -841,19 +841,51 @@ ScatterSimulation::downsample_scanner(int new_num_rings, int new_num_dets)
float approx_num_non_arccorrected_bins;
if (new_scanner_sptr->get_scanner_geometry() != "Cylindrical")
{
- new_num_dets = this->proj_data_info_sptr->get_scanner_ptr()->get_num_detectors_per_ring();
- approx_num_non_arccorrected_bins = this->proj_data_info_sptr->get_num_tangential_poss();
+ if (new_num_dets <= 0)
+ { // by default, do not downsample the detectors per ring for BlocksOnCylindrical
+ new_num_dets = this->proj_data_info_sptr->get_scanner_ptr()->get_num_detectors_per_ring();
+ }
+ // STIR does not like block spacings that are smaller than number of crystals times crystal spacing,
+ // therefore add a 1% extension on top for the downsampled scanner, to avoid running into floating point issues
+ const float block_spacing_factor = 1.01;
+
+ // extend the bins by a small amount to avoid edge-effects, at the expense of longer computation time
+ approx_num_non_arccorrected_bins = ceil(this->proj_data_info_sptr->get_num_tangential_poss() * float(new_num_dets)
+ / old_scanner_ptr->get_num_detectors_per_ring())
+ + 1;
// preserve the length of the scanner, but place crystals equidistantly
float scanner_length = old_scanner_ptr->get_axial_length();
float new_ring_spacing = scanner_length / (new_num_rings - 1);
-
new_scanner_sptr->set_num_axial_blocks_per_bucket(1);
+ new_scanner_sptr->set_num_transaxial_blocks_per_bucket(1);
+
new_scanner_sptr->set_num_rings(new_num_rings);
+ new_scanner_sptr->set_num_detectors_per_ring(new_num_dets);
+
new_scanner_sptr->set_axial_crystal_spacing(new_ring_spacing);
new_scanner_sptr->set_ring_spacing(new_ring_spacing);
new_scanner_sptr->set_num_axial_crystals_per_block(new_num_rings);
- new_scanner_sptr->set_axial_block_spacing(new_ring_spacing * new_scanner_sptr->get_num_axial_crystals_per_block());
+ new_scanner_sptr->set_axial_block_spacing(new_ring_spacing * new_num_rings * block_spacing_factor);
+
+ float transaxial_bucket_width
+ = old_scanner_ptr->get_transaxial_block_spacing() * (old_scanner_ptr->get_num_transaxial_blocks_per_bucket() - 1)
+ + old_scanner_ptr->get_transaxial_crystal_spacing() * (old_scanner_ptr->get_num_transaxial_crystals_per_block() - 1);
+
+ int num_trans_buckets = old_scanner_ptr->get_num_transaxial_buckets();
+ // get a new number of detectors that is a multiple of the number of buckets to preserve scanner shape
+ new_scanner_sptr->set_num_detectors_per_ring(new_num_dets);
+ int new_transaxial_dets_per_bucket = new_num_dets / num_trans_buckets;
+ float new_det_spacing = transaxial_bucket_width / (new_transaxial_dets_per_bucket - 1);
+
+ new_scanner_sptr->set_num_transaxial_blocks_per_bucket(1);
+ new_scanner_sptr->set_num_transaxial_crystals_per_block(new_transaxial_dets_per_bucket);
+ new_scanner_sptr->set_transaxial_crystal_spacing(new_det_spacing);
+ new_scanner_sptr->set_transaxial_block_spacing(new_transaxial_dets_per_bucket * new_det_spacing * block_spacing_factor);
+ // avoid problems with Scanner checks by setting singles_units to 1 crystal
+ // (only used for dead-time processing in ECAY norm)
+ new_scanner_sptr->set_num_axial_crystals_per_singles_unit(1);
+ new_scanner_sptr->set_num_transaxial_crystals_per_singles_unit(1);
}
else
{
diff --git a/src/test/CMakeLists.txt b/src/test/CMakeLists.txt
index 6b7ae72fd..e111be65b 100644
--- a/src/test/CMakeLists.txt
+++ b/src/test/CMakeLists.txt
@@ -74,7 +74,7 @@ set(buildblock_simple_tests
include(stir_test_exe_targets)
foreach(source ${buildblock_simple_tests})
- create_stir_test(${source} "buildblock;IO;buildblock;numerics_buildblock;display;IO;recon_buildblock;Shape_buildblock" "")
+ create_stir_test(${source} "buildblock;IO;buildblock;numerics_buildblock;display;IO;recon_buildblock;Shape_buildblock;scatter_buildblock" "")
endforeach()
#
diff --git a/src/test/test_interpolate_projdata.cxx b/src/test/test_interpolate_projdata.cxx
index 29c554ce7..3153678b3 100644
--- a/src/test/test_interpolate_projdata.cxx
+++ b/src/test/test_interpolate_projdata.cxx
@@ -41,6 +41,7 @@
#include "stir/Shape/Box3D.h"
#include "stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h"
#include "stir/recon_buildblock/ForwardProjectorByBinUsingProjMatrixByBin.h"
+#include "stir/scatter/SingleScatterSimulation.h"
#include "stir/RunTests.h"
@@ -56,6 +57,8 @@ class InterpolationTests : public RunTests
void scatter_interpolation_test_cyl();
void scatter_interpolation_test_blocks_asymmetric();
void scatter_interpolation_test_cyl_asymmetric();
+ void scatter_interpolation_test_blocks_downsampled();
+ void transaxial_upsampling_interpolation_test_blocks();
void check_symmetry(const SegmentBySinogram& segment);
void compare_segment(const SegmentBySinogram& segment1, const SegmentBySinogram& segment2, float maxDiff);
@@ -675,6 +678,280 @@ InterpolationTests::scatter_interpolation_test_cyl_asymmetric()
compare_segment_shape(full_size_model_sino.get_segment_by_sinogram(0), interpolated_proj_data.get_segment_by_sinogram(0), 2);
}
+void
+InterpolationTests::scatter_interpolation_test_blocks_downsampled()
+{
+ info("Performing downampled interpolation test for BlocksOnCylindrical scanner");
+ auto time_frame_def = TimeFrameDefinitions();
+ time_frame_def.set_num_time_frames(1);
+ time_frame_def.set_time_frame(1, 0, 1e9);
+ auto exam_info = ExamInfo();
+ exam_info.set_high_energy_thres(650);
+ exam_info.set_low_energy_thres(425);
+ exam_info.set_time_frame_definitions(time_frame_def);
+
+ // define the original scanner and a downsampled one, as it would be used for scatter simulation
+ auto scanner = Scanner(Scanner::User_defined_scanner,
+ "Some_BlocksOnCylindrical_Scanner",
+ 96,
+ 30,
+ int(150 * 96 / 192),
+ int(150 * 96 / 192),
+ 127,
+ 6.5,
+ 3.313,
+ 4.156,
+ -3.1091819,
+ 5,
+ 3,
+ 6,
+ 4,
+ 1,
+ 1,
+ 1,
+ 0.14,
+ 511,
+ 1,
+ 0,
+ 500,
+ "BlocksOnCylindrical",
+ 3.313,
+ 7.0,
+ 20.0,
+ 29.0);
+ auto downsampled_scanner = Scanner(Scanner::User_defined_scanner,
+ "Some_Downsampled_BlocksOnCylindrical_Scanner",
+ 64,
+ 8,
+ int(150 * 64 / 192 + 1),
+ int(150 * 64 / 192 + 1),
+ 127,
+ 6.5,
+ 16.652,
+ 6.234,
+ -3.1091819,
+ 1,
+ 1,
+ 8,
+ 8,
+ 1,
+ 1,
+ 1,
+ 0.17,
+ 511,
+ 1,
+ 0,
+ 500,
+ "BlocksOnCylindrical",
+ 13.795,
+ 15.4286,
+ 112.0,
+ 125.0);
+
+ auto proj_data_info = shared_ptr(std::move(
+ ProjDataInfo::construct_proj_data_info(std::make_shared(scanner), 1, 29, 48, int(150 * 96 / 192), false)));
+ auto downsampled_proj_data_info = shared_ptr(std::move(ProjDataInfo::construct_proj_data_info(
+ std::make_shared(downsampled_scanner), 1, 0, 32, int(150 * 64 / 192 + 1), false)));
+
+ auto proj_data = ProjDataInMemory(std::make_shared(exam_info), proj_data_info);
+ auto downsampled_proj_data = ProjDataInMemory(std::make_shared(exam_info), downsampled_proj_data_info);
+
+ // define a cylinder and a box that are off-centre, such that the shapes in the sinogram can be compared
+ auto emission_map = VoxelsOnCartesianGrid(*proj_data_info, 1);
+ auto cyl_map = VoxelsOnCartesianGrid(*proj_data_info, 1);
+ auto cylinder = EllipsoidalCylinder(40, 40, 20, CartesianCoordinate3D(80, 100, 0));
+ cylinder.construct_volume(cyl_map, CartesianCoordinate3D(1, 1, 1));
+ auto box = Box3D(20, 20, 20, CartesianCoordinate3D(30, -20, 70));
+ box.construct_volume(emission_map, CartesianCoordinate3D(1, 1, 1));
+ emission_map += cyl_map;
+
+ // project the emission map onto the full-scale scanner proj data
+ auto pm = ProjMatrixByBinUsingRayTracing();
+ pm.set_use_actual_detector_boundaries(true);
+ pm.enable_cache(false);
+ auto forw_proj = ForwardProjectorByBinUsingProjMatrixByBin(std::make_shared(pm));
+ forw_proj.set_up(proj_data_info, std::make_shared>(emission_map));
+ auto full_size_model_sino = ProjDataInMemory(proj_data);
+ full_size_model_sino.fill(0);
+ forw_proj.forward_project(full_size_model_sino, emission_map);
+
+ // also project onto the downsampled scanner
+ emission_map = VoxelsOnCartesianGrid(*downsampled_proj_data_info, 1);
+ cyl_map = VoxelsOnCartesianGrid(*downsampled_proj_data_info, 1);
+ cylinder.construct_volume(cyl_map, CartesianCoordinate3D(1, 1, 1));
+ box.construct_volume(emission_map, CartesianCoordinate3D(1, 1, 1));
+ emission_map += cyl_map;
+ forw_proj.set_up(downsampled_proj_data_info, std::make_shared>(emission_map));
+ auto downsampled_model_sino = ProjDataInMemory(downsampled_proj_data);
+ downsampled_model_sino.fill(0);
+ forw_proj.forward_project(downsampled_model_sino, emission_map);
+
+ // write the proj data to file
+ downsampled_model_sino.write_to_file("transaxially_downsampled_sino.hs");
+ full_size_model_sino.write_to_file("transaxially_full_size_sino.hs");
+
+ // interpolate the downsampled proj data to the original scanner size and fill in oblique sinograms
+ auto interpolated_direct_proj_data = ProjDataInMemory(proj_data);
+ interpolate_projdata(interpolated_direct_proj_data, downsampled_model_sino, BSpline::linear, false);
+ auto interpolated_proj_data = ProjDataInMemory(proj_data);
+ inverse_SSRB(interpolated_proj_data, interpolated_direct_proj_data);
+
+ // write the proj data to file
+ interpolated_proj_data.write_to_file("transaxially_interpolated_sino.hs");
+
+ // compare to ground truth
+ compare_segment_shape(full_size_model_sino.get_segment_by_sinogram(0), interpolated_proj_data.get_segment_by_sinogram(0), 3);
+}
+
+void
+InterpolationTests::transaxial_upsampling_interpolation_test_blocks()
+{
+ info("Performing transaxial downampled interpolation test for BlocksOnCylindrical scanner");
+ auto time_frame_def = TimeFrameDefinitions();
+ time_frame_def.set_num_time_frames(1);
+ time_frame_def.set_time_frame(1, 0, 1e9);
+ auto exam_info = ExamInfo();
+ exam_info.set_high_energy_thres(650);
+ exam_info.set_low_energy_thres(425);
+ exam_info.set_time_frame_definitions(time_frame_def);
+
+ // define the original scanner and a downsampled one, as it would be used for scatter simulation
+ auto scanner = Scanner(Scanner::User_defined_scanner,
+ "Some_BlocksOnCylindrical_Scanner",
+ 96,
+ 3,
+ 60,
+ 60,
+ 127,
+ 6.5,
+ 3.313,
+ 1.65,
+ -3.1091819,
+ 1,
+ 3,
+ 3,
+ 4,
+ 1,
+ 1,
+ 1,
+ 0.14,
+ 511,
+ 1,
+ 0,
+ 500,
+ "BlocksOnCylindrical",
+ 3.313,
+ 7.0,
+ 20.0,
+ 29.0);
+ auto proj_data_info = shared_ptr(
+ std::move(ProjDataInfo::construct_proj_data_info(std::make_shared(scanner), 1, 0, 48, 60, false)));
+
+ // use the code in scatter simulation to downsample the scanner
+ auto scatter_simulation = SingleScatterSimulation();
+ scatter_simulation.set_template_proj_data_info(*proj_data_info);
+ scatter_simulation.set_exam_info(exam_info);
+ scatter_simulation.downsample_scanner(-1, 96 / 4); // number of detectors per ring reduced by factor of four
+ auto downsampled_proj_data_info = scatter_simulation.get_template_proj_data_info_sptr();
+
+ auto proj_data = ProjDataInMemory(std::make_shared(exam_info), proj_data_info);
+ auto downsampled_proj_data = ProjDataInMemory(std::make_shared(exam_info), downsampled_proj_data_info);
+
+ // define a cylinder precisely in the middle of the FOV
+ auto emission_map = VoxelsOnCartesianGrid(*downsampled_proj_data_info, 1);
+ make_symmetric_object(emission_map);
+
+ // project the cylinder onto the full-scale scanner proj data
+ auto pm = ProjMatrixByBinUsingRayTracing();
+ pm.set_use_actual_detector_boundaries(true);
+ pm.enable_cache(false);
+ auto forw_proj = ForwardProjectorByBinUsingProjMatrixByBin(std::make_shared(pm));
+ forw_proj.set_up(proj_data_info, std::make_shared>(emission_map));
+ auto full_size_model_sino = ProjDataInMemory(proj_data);
+ full_size_model_sino.fill(0);
+ forw_proj.forward_project(full_size_model_sino, emission_map);
+
+ // also project onto the downsampled scanner
+ emission_map = VoxelsOnCartesianGrid(*downsampled_proj_data_info, 1);
+ make_symmetric_object(emission_map);
+ forw_proj.set_up(downsampled_proj_data_info, std::make_shared>(emission_map));
+ auto downsampled_model_sino = ProjDataInMemory(downsampled_proj_data);
+ downsampled_model_sino.fill(0);
+ forw_proj.forward_project(downsampled_model_sino, emission_map);
+
+ // write the proj data to file
+ downsampled_model_sino.write_to_file("transaxially_downsampled_sino_for_LOR.hs");
+
+ // interpolate the downsampled proj data to the original scanner size and fill in oblique sinograms
+ auto interpolated_direct_proj_data = ProjDataInMemory(proj_data);
+ interpolate_projdata(interpolated_direct_proj_data, downsampled_model_sino, BSpline::linear, false);
+ auto interpolated_proj_data = ProjDataInMemory(proj_data);
+ inverse_SSRB(interpolated_proj_data, interpolated_direct_proj_data);
+
+ // write the proj data to file
+ interpolated_proj_data.write_to_file("transaxially_interpolated_sino_for_LOR.hs");
+
+ // Identify the bins which should be identical between the downsampled and the interpolated sinogram:
+ // Each module has 96 / 8 = 12 crystal in the full size scanner, organised in 3 blocks of 4 crystals, while
+ // the downsampled scanner has 3 crystals per module. The idea is that the centre of the outer two
+ // is in exactly the same position than the centre of the first and last crystal in the full size scanner.
+ SegmentBySinogram sinogram_downsampled = downsampled_proj_data.get_empty_segment_by_sinogram(0, false, 0);
+ SegmentBySinogram sinogram_full_size = proj_data.get_empty_segment_by_sinogram(0, false, 0);
+ const auto pdi_downsampled
+ = dynamic_cast(&(*downsampled_proj_data.get_proj_data_info_sptr()));
+ const auto pdi_full_size = dynamic_cast(&(*sinogram_full_size.get_proj_data_info_sptr()));
+
+ int tested_LORs = 0;
+ for (int det1_downsampled = 0; det1_downsampled < 3 * 8; det1_downsampled++)
+ {
+ if (det1_downsampled % 3 == 1)
+ continue; // skip the central crystal of each module
+ for (int det2_downsampled = 0; det2_downsampled < 3 * 8; det2_downsampled++)
+ {
+ if (det2_downsampled % 3 == 1 || det1_downsampled == det2_downsampled)
+ continue; // skip the central crystal of each module
+ if (det1_downsampled / 3 == det2_downsampled / 3)
+ continue; // skip the LORs that lie on the same module
+
+ int view_ds, tang_pos_ds;
+ pdi_downsampled->get_view_tangential_pos_num_for_det_num_pair(view_ds, tang_pos_ds, det1_downsampled, det2_downsampled);
+ BasicCoordinate<3, int> index_downsampled;
+ index_downsampled[1] = 1; // looking at central slice
+ index_downsampled[2] = view_ds;
+ index_downsampled[3] = tang_pos_ds;
+
+ if (tang_pos_ds < pdi_downsampled->get_min_tangential_pos_num()
+ || tang_pos_ds > pdi_downsampled->get_max_tangential_pos_num())
+ continue;
+
+ int view_fs, tang_pos_fs;
+ pdi_full_size->get_view_tangential_pos_num_for_det_num_pair(
+ view_fs,
+ tang_pos_fs,
+ (det1_downsampled / 3) * 12 + ((det1_downsampled % 3) / 2) * 11,
+ (det2_downsampled / 3) * 12 + ((det2_downsampled % 3) / 2) * 11);
+
+ BasicCoordinate<3, int> index_full_size;
+ index_full_size[1] = 1; // looking at central slice
+ index_full_size[2] = view_fs;
+ index_full_size[3] = tang_pos_fs;
+
+ if (tang_pos_fs < pdi_full_size->get_min_tangential_pos_num()
+ || tang_pos_fs > pdi_full_size->get_max_tangential_pos_num())
+ continue;
+
+ // confirm that the difference is smaller than an empirically found value
+ check_if_less(std::abs(sinogram_downsampled[index_downsampled] - sinogram_full_size[index_full_size]),
+ 0.01,
+ "difference between sinogram bin is larger than expected");
+
+ tested_LORs++;
+ }
+ }
+
+ info(boost::format("A total of %1% LORs were compared between the downsampled and the interpolated sinogram.") % tested_LORs);
+}
+
void
InterpolationTests::run_tests()
{
@@ -682,6 +959,8 @@ InterpolationTests::run_tests()
scatter_interpolation_test_cyl();
scatter_interpolation_test_blocks_asymmetric();
scatter_interpolation_test_cyl_asymmetric();
+ scatter_interpolation_test_blocks_downsampled();
+ transaxial_upsampling_interpolation_test_blocks();
}
END_NAMESPACE_STIR