Added an additional test that compares sinogram bins for LORs that do not change during scanner downsampling.

Author: Markus Jehl

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()
 #
 

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"
 
@@ -57,6 +58,7 @@ class InterpolationTests : public RunTests
   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<float>& segment);
   void compare_segment(const SegmentBySinogram<float>& segment1, const SegmentBySinogram<float>& segment2, float maxDiff);
@@ -754,7 +756,7 @@ InterpolationTests::scatter_interpolation_test_blocks_downsampled()
   auto proj_data = ProjDataInMemory(std::make_shared<ExamInfo>(exam_info), proj_data_info);
   auto downsampled_proj_data = ProjDataInMemory(std::make_shared<ExamInfo>(exam_info), downsampled_proj_data_info);
 
-  // define a cylinder precisely in the middle of the FOV, such that symmetry can be used for validation
+  // define a cylinder and a box that are off-centre, such that the shapes in the sinogram can be compared
   auto emission_map = VoxelsOnCartesianGrid<float>(*proj_data_info, 1);
   auto cyl_map = VoxelsOnCartesianGrid<float>(*proj_data_info, 1);
   auto cylinder = EllipsoidalCylinder(40, 40, 20, CartesianCoordinate3D<float>(80, 100, 0));
@@ -763,7 +765,7 @@ InterpolationTests::scatter_interpolation_test_blocks_downsampled()
   box.construct_volume(emission_map, CartesianCoordinate3D<int>(1, 1, 1));
   emission_map += cyl_map;
 
-  // project the cylinder onto the full-scale scanner proj data
+  // 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);
@@ -801,6 +803,155 @@ InterpolationTests::scatter_interpolation_test_blocks_downsampled()
   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<ProjDataInfo>(
+      std::move(ProjDataInfo::construct_proj_data_info(std::make_shared<Scanner>(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<ExamInfo>(exam_info), proj_data_info);
+  auto downsampled_proj_data = ProjDataInMemory(std::make_shared<ExamInfo>(exam_info), downsampled_proj_data_info);
+
+  // define a cylinder precisely in the middle of the FOV
+  auto emission_map = VoxelsOnCartesianGrid<float>(*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<ProjMatrixByBinUsingRayTracing>(pm));
+  forw_proj.set_up(proj_data_info, std::make_shared<VoxelsOnCartesianGrid<float>>(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<float>(*downsampled_proj_data_info, 1);
+  make_symmetric_object(emission_map);
+  forw_proj.set_up(downsampled_proj_data_info, std::make_shared<VoxelsOnCartesianGrid<float>>(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<float> sinogram_downsampled = downsampled_proj_data.get_empty_segment_by_sinogram(0, false, 0);
+  SegmentBySinogram<float> sinogram_full_size = proj_data.get_empty_segment_by_sinogram(0, false, 0);
+  const auto pdi_downsampled
+      = dynamic_cast<const ProjDataInfoGenericNoArcCorr*>(&(*downsampled_proj_data.get_proj_data_info_sptr()));
+  const auto pdi_full_size = dynamic_cast<const ProjDataInfoGenericNoArcCorr*>(&(*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()
 {
@@ -809,6 +960,7 @@ InterpolationTests::run_tests()
   scatter_interpolation_test_blocks_asymmetric();
   scatter_interpolation_test_cyl_asymmetric();
   scatter_interpolation_test_blocks_downsampled();
+  transaxial_upsampling_interpolation_test_blocks();
 }
 
 END_NAMESPACE_STIR