Merge pull request #1375 from robbietuk/axial_bucket_projection_issue

Throw errors if incorrect/unsupported blocks scanner. See #1374

Author: KrisThielemans

documentation/release_6.1.htm

@@ -85,8 +85,10 @@ <H2>What's new for developers (aside from what should be obvious
 
 <h3>Backward incompatibities</h3>
 <ul>
-<li>
-</li>
+    <li>
+        Additional checks on <code>GeometryBlocksOnCylindrical</code> scanner configuration, which may lead to an
+        error being raised, while previously the code silently proceeded.
+    </li>
 </ul>
 
 <h3>New functionality</h3>

src/buildblock/GeometryBlocksOnCylindrical.cxx

@@ -39,11 +39,14 @@ limitations under the License.
 #include <iostream>
 #include <fstream>
 #include <iomanip>
+#include <boost/format.hpp>
 
 START_NAMESPACE_STIR
 
 GeometryBlocksOnCylindrical::GeometryBlocksOnCylindrical(const Scanner& scanner)
 {
+  if (scanner.check_consistency() == Succeeded::no)
+    error("Error in GeometryBlocksOnCylindrical: scanner configuration not accepted. Please check warnings.");
   build_crystal_maps(scanner);
 }
 
@@ -61,11 +64,10 @@ GeometryBlocksOnCylindrical::build_crystal_maps(const Scanner& scanner)
   // local variables to describe scanner
   int num_axial_crystals_per_block = scanner.get_num_axial_crystals_per_block();
   int num_transaxial_crystals_per_block = scanner.get_num_transaxial_crystals_per_block();
-  int num_axial_blocks = scanner.get_num_axial_blocks();
   int num_transaxial_blocks_per_bucket = scanner.get_num_transaxial_blocks_per_bucket();
   int num_axial_blocks_per_bucket = scanner.get_num_axial_blocks_per_bucket();
-  int num_transaxial_buckets = scanner.get_num_transaxial_blocks() / num_transaxial_blocks_per_bucket;
-  int num_axial_buckets = scanner.get_num_axial_blocks() / num_axial_blocks_per_bucket;
+  int num_transaxial_buckets = scanner.get_num_transaxial_buckets();
+  int num_axial_buckets = scanner.get_num_axial_buckets();
   int num_detectors_per_ring = scanner.get_num_detectors_per_ring();
   float axial_block_spacing = scanner.get_axial_block_spacing();
   float transaxial_block_spacing = scanner.get_transaxial_block_spacing();
@@ -100,7 +102,7 @@ GeometryBlocksOnCylindrical::build_crystal_maps(const Scanner& scanner)
   stir::CartesianCoordinate3D<float> start_point(start_z, start_y, start_x);
 
   for (int ax_bucket_num = 0; ax_bucket_num < num_axial_buckets; ++ax_bucket_num)
-    for (int ax_block_num = 0; ax_block_num < num_axial_blocks; ++ax_block_num)
+    for (int ax_block_num = 0; ax_block_num < num_axial_blocks_per_bucket; ++ax_block_num)
       for (int ax_crys_num = 0; ax_crys_num < num_axial_crystals_per_block; ++ax_crys_num)
         for (int trans_bucket_num = 0; trans_bucket_num < num_transaxial_buckets; ++trans_bucket_num)
           for (int trans_block_num = 0; trans_block_num < num_transaxial_blocks_per_bucket; ++trans_block_num)

src/buildblock/Scanner.cxx

@@ -42,6 +42,7 @@
 #include <iostream>
 #include <algorithm>
 #include <sstream>
+#include <boost/format.hpp>
 #include "stir/warning.h"
 #include "stir/error.h"
 
@@ -1751,6 +1752,23 @@ Scanner::check_consistency() const
 
   if (get_scanner_geometry() == "BlocksOnCylindrical")
     { //! Check consistency of axial and transaxial spacing for block geometry
+      if (get_num_axial_buckets() != 1)
+        {
+          warning(boost::format("BlocksOnCylindrical num_axial_buckets (%d) is greater than 1. This is not supported yet. "
+                                "Consider multiplying the number of axial_blocks_per_bucket by %d.")
+                  % get_num_axial_buckets() % get_num_axial_buckets());
+          return Succeeded::no;
+        }
+      { // Assert that each block contains an equal number of axial crystals
+        if (get_num_rings() % get_num_axial_crystals_per_bucket() != 0)
+          {
+            warning(boost::format("Error in GeometryBlocksOnCylindrical: number of rings (%d) is not a multiple of the "
+                                  "get_num_axial_crystals_per_bucket "
+                                  "(%d) = num_axial_crystals_per_block (%d) * num_axial_blocks_per_bucket (%d)")
+                    % get_num_rings() % get_num_axial_crystals_per_bucket() % get_num_axial_crystals_per_block()
+                    % get_num_axial_blocks_per_bucket());
+          }
+      }
       if (get_axial_crystal_spacing() * get_num_axial_crystals_per_block() > get_axial_block_spacing())
         {
           warning(
@@ -1786,22 +1804,22 @@ Scanner::check_consistency() const
                   get_inner_ring_radius());
           return Succeeded::no;
         }
-      else if (get_scanner_geometry() == "Generic")
-        { //! Check if the crystal map is correct and given
-          if (get_crystal_map_file_name() == "")
+    }
+  else if (get_scanner_geometry() == "Generic")
+    { //! Check if the crystal map is correct and given
+      if (get_crystal_map_file_name().empty())
+        {
+          warning("No crystal map is provided. The scanner geometry Generic needs it! Please provide one.");
+          return Succeeded::no;
+        }
+      else
+        {
+          std::ifstream crystal_map(get_crystal_map_file_name());
+          if (!crystal_map)
             {
-              warning("No crystal map is provided. The scanner geometry Generic needs it! Please provide one.");
+              warning("No correct crystal map provided. Please check the file name.");
               return Succeeded::no;
             }
-          else
-            {
-              std::ifstream crystal_map(get_crystal_map_file_name());
-              if (!crystal_map)
-                {
-                  warning("No correct crystal map provided. Please check the file name.");
-                  return Succeeded::no;
-                }
-            }
         }
     }
 

src/recon_test/CMakeLists.txt

@@ -26,6 +26,7 @@ set(${dir_SIMPLE_TEST_EXE_SOURCES}
         test_FBP3DRP.cxx
         test_priors.cxx
         test_blocks_on_cylindrical_projectors.cxx
+        test_geometry_blocks_on_cylindrical.cxx
 )
 
 

src/recon_test/test_blocks_on_cylindrical_projectors.cxx

@@ -6,9 +6,11 @@
   \brief Test program for back projection and forward projection using stir::ProjDataInfoBlockOnCylindrical
 
   \author Daniel Deidda
+  \author Robert Twyman
 
 */
 /*  Copyright (C) 2021-2022, National Physical Laboratory
+    Copyright (C) 2024, Prescient Imaging
     This file is part of STIR.
 
     SPDX-License-Identifier: Apache-2.0
@@ -52,8 +54,6 @@
 #endif
 #include "stir/recon_buildblock/BackProjectorByBinUsingProjMatrixByBin.h"
 #include "stir/IO/write_to_file.h"
-#include "stir/VoxelsOnCartesianGrid.h"
-//#include "stir/Shape/Shape3D.h"
 #include <cmath>
 
 START_NAMESPACE_STIR
@@ -80,6 +80,7 @@ class BlocksTests : public RunTests
   void run_map_orientation_test(ForwardProjectorByBin& forw_projector1, ForwardProjectorByBin& forw_projector2);
   void run_projection_test(ForwardProjectorByBin& forw_projector1, ForwardProjectorByBin& forw_projector2);
   void run_intersection_with_cylinder_test();
+  void run_back_projection_test_with_axial_buckets(BackProjectorByBin& back_projector);
 };
 
 /*! The following is a function to allow a projdata_info BlocksOnCylindrical to be created from the scanner. */
@@ -871,6 +872,64 @@ BlocksTests::run_intersection_with_cylinder_test()
     }
 }
 
+void
+BlocksTests::run_back_projection_test_with_axial_buckets(BackProjectorByBin& back_projector)
+
+{
+  int num_buckets = 1;
+  auto scanner_sptr = std::make_shared<Scanner>(Scanner::SAFIRDualRingPrototype);
+  { // create geometry
+    scanner_sptr->set_average_depth_of_interaction(5);
+    scanner_sptr->set_num_axial_crystals_per_block(1);
+    scanner_sptr->set_axial_block_spacing(scanner_sptr->get_axial_crystal_spacing()
+                                          * scanner_sptr->get_num_axial_crystals_per_block());
+    scanner_sptr->set_transaxial_block_spacing(scanner_sptr->get_transaxial_crystal_spacing()
+                                               * scanner_sptr->get_num_transaxial_crystals_per_block());
+    scanner_sptr->set_num_axial_blocks_per_bucket(2);
+    scanner_sptr->set_num_rings(scanner_sptr->get_num_axial_crystals_per_bucket() * num_buckets);
+    scanner_sptr->set_scanner_geometry("BlocksOnCylindrical");
+    scanner_sptr->set_up();
+  }
+
+  auto projdata_info_sptr = set_direct_projdata_info<ProjDataInfoBlocksOnCylindricalNoArcCorr>(scanner_sptr, 1);
+  auto exam_info_sptr = std::make_shared<ExamInfo>(ImagingModality::PT);
+  auto projdata_sptr = std::make_shared<ProjDataInMemory>(exam_info_sptr, projdata_info_sptr);
+
+  auto origin = CartesianCoordinate3D<float>(0, 0, 0);
+  auto volume_dimensions = CartesianCoordinate3D<int>(-1, -1, -1);
+  auto volume_sptr
+      = std::make_shared<VoxelsOnCartesianGrid<float>>(exam_info_sptr, *projdata_info_sptr, 1, origin, volume_dimensions);
+
+  // Now run the test
+  volume_sptr->fill(0.0);
+  projdata_sptr->fill(1.0);
+
+  back_projector.set_up(projdata_info_sptr, volume_sptr);
+  back_projector.back_project(*volume_sptr, *projdata_sptr, 0, projdata_info_sptr->get_num_views());
+
+  bool test_ok = true;
+  std::ostringstream oss;
+  oss << "BlocksTests::run_back_projection_test_with_axial_buckets: Central voxel values are:" << std::endl;
+
+  auto centre_axial_values = std::vector<float>(volume_sptr->get_z_size());
+  for (int z = volume_sptr->get_min_z(); z <= volume_sptr->get_max_z(); z++)
+    {
+      centre_axial_values[z] = (*volume_sptr)[z][0][0];
+      oss << "\tz = " << z << "/" << volume_sptr->get_max_z() << " is " << centre_axial_values[z] << std::endl;
+      if (test_ok)
+        {
+          test_ok = check(centre_axial_values[z] > 0,
+                          "BlocksTests::run_back_projection_test_with_axial_buckets: Central voxel value <= 0");
+          everything_ok = everything_ok && test_ok;
+        }
+    }
+
+  if (test_ok)
+    return;
+  // Something went wrong, log the error
+  std::cerr << oss.str();
+}
+
 void
 BlocksTests::run_tests()
 {
@@ -905,6 +964,8 @@ BlocksTests::run_tests()
   print_time("map orientation test took: ");
   run_intersection_with_cylinder_test();
   print_time("intersection with cylinder test took: ");
+  run_back_projection_test_with_axial_buckets(back_projector);
+  print_time("back projection test with axial buckets took: ");
 
 #ifdef STIR_WITH_Parallelproj_PROJECTOR
   // run the same tests with parallelproj, if available

src/recon_test/test_geometry_blocks_on_cylindrical.cxx

@@ -0,0 +1,151 @@
+/*!
+
+\file
+\ingroup recontest
+
+\brief Test program to ensure the axial coordinates of blocks on cylindrical are monotonic with axial indices
+
+\author Robert Twyman
+
+    Copyright (C) 2024, Prescient Imaging
+    This file is part of STIR.
+
+    SPDX-License-Identifier: Apache-2.0
+
+    See STIR/LICENSE.txt for details
+*/
+
+#include "stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h"
+#include "stir/ExamInfo.h"
+#include "stir/Verbosity.h"
+#include "stir/LORCoordinates.h"
+#include "stir/ProjDataInfoGenericNoArcCorr.h"
+#include "stir/Succeeded.h"
+#include "stir/RunTests.h"
+#include "stir/Scanner.h"
+#include "stir/HighResWallClockTimer.h"
+#include "stir/GeometryBlocksOnCylindrical.h"
+#include "stir/IO/write_to_file.h"
+#include <cmath>
+
+START_NAMESPACE_STIR
+
+/*!
+\ingroup test
+\brief Test class for BlocksOnCylindrical geometry
+*/
+class GeometryBlocksOnCylindricalTests : public RunTests
+{
+public:
+  void run_tests() override;
+
+private:
+  /*! \brief Tests multiple axial blocks/bucket configurations to ensure the detector map's axial indices and coordinates
+   * are monotonic
+   */
+  void run_monotonic_axial_coordinates_in_detector_map_test();
+  //! Tests the axial indices and coordinates are monotonic in the detector map
+  static Succeeded monotonic_axial_coordinates_in_detector_map_test(const shared_ptr<Scanner>& scanner_sptr);
+};
+
+void
+GeometryBlocksOnCylindricalTests::run_monotonic_axial_coordinates_in_detector_map_test()
+{
+  auto scanner_sptr = std::make_shared<Scanner>(Scanner::SAFIRDualRingPrototype);
+  scanner_sptr->set_scanner_geometry("BlocksOnCylindrical");
+  scanner_sptr->set_transaxial_block_spacing(scanner_sptr->get_transaxial_crystal_spacing()
+                                             * scanner_sptr->get_num_transaxial_crystals_per_block());
+  int num_axial_buckets = 1; // TODO add for loop when support is added
+
+  for (int num_axial_crystals_per_blocks = 1; num_axial_crystals_per_blocks < 3; ++num_axial_crystals_per_blocks)
+    for (int num_axial_blocks_per_bucket = 1; num_axial_blocks_per_bucket < 3; ++num_axial_blocks_per_bucket)
+      {
+        scanner_sptr->set_num_axial_crystals_per_block(num_axial_crystals_per_blocks);
+        scanner_sptr->set_num_axial_blocks_per_bucket(num_axial_blocks_per_bucket);
+        scanner_sptr->set_num_rings(scanner_sptr->get_num_axial_crystals_per_bucket() * num_axial_buckets);
+        scanner_sptr->set_axial_block_spacing(scanner_sptr->get_axial_crystal_spacing()
+                                              * (scanner_sptr->get_num_axial_crystals_per_block() + 0.5));
+
+        if (monotonic_axial_coordinates_in_detector_map_test(scanner_sptr) == Succeeded::no)
+          {
+            warning(boost::format("Monothonic axial coordinates test failed for:\n"
+                                  "\taxial_crystal_per_block =\t%1%\n"
+                                  "\taxial_blocks_per_bucket =\t%2%\n"
+                                  "\tnum_axial_buckets =\t\t\t%3%")
+                    % num_axial_crystals_per_blocks % num_axial_blocks_per_bucket % num_axial_buckets);
+            everything_ok = false;
+            return;
+          }
+      }
+}
+
+Succeeded
+GeometryBlocksOnCylindricalTests::monotonic_axial_coordinates_in_detector_map_test(const shared_ptr<Scanner>& scanner_sptr)
+{
+  if (scanner_sptr->get_scanner_geometry() != "BlocksOnCylindrical")
+    {
+      warning("monotonic_axial_coordinates_in_detector_map_test is only for the BlocksOnCylindrical geometry");
+      return Succeeded::no;
+    }
+
+  shared_ptr<DetectorCoordinateMap> detector_map_sptr;
+  try
+    {
+      detector_map_sptr.reset(new GeometryBlocksOnCylindrical(*scanner_sptr));
+    }
+  catch (const std::runtime_error& e)
+    {
+      warning(boost::format("Caught runtime_error while creating GeometryBlocksOnCylindrical: %1%\n"
+                            "Failing the test.")
+              % e.what());
+      return Succeeded::no;
+    }
+
+  unsigned min_axial_pos = 0;
+  float prev_min_axial_coord = -std::numeric_limits<float>::max();
+
+  for (unsigned axial_idx = 0; axial_idx < detector_map_sptr->get_num_axial_coords(); ++axial_idx)
+    for (unsigned tangential_idx = 0; tangential_idx < detector_map_sptr->get_num_tangential_coords(); ++tangential_idx)
+      for (unsigned radial_idx = 0; radial_idx < detector_map_sptr->get_num_radial_coords(); ++radial_idx)
+        {
+          const DetectionPosition<> det_pos = DetectionPosition<>(tangential_idx, axial_idx, radial_idx);
+          CartesianCoordinate3D<float> coord = detector_map_sptr->get_coordinate_for_det_pos(det_pos);
+          if (coord.z() > prev_min_axial_coord)
+            {
+              min_axial_pos = axial_idx;
+              prev_min_axial_coord = coord.z();
+            }
+          else if (coord.z() < prev_min_axial_coord)
+            {
+              float delta = coord.z() - prev_min_axial_coord;
+              warning(boost::format("Axial Coordinates are not monotonic.\n"
+                                    "Next axial index =\t\t%1%, Next axial coord (mm) =\t\t%2%  (%3%)\n"
+                                    "Previous axial index =\t%4%, Previous axial coord (mm) =\t%5%")
+                      % axial_idx % coord.z() % delta % min_axial_pos % prev_min_axial_coord);
+              return Succeeded::no;
+            }
+        }
+
+  return Succeeded::yes;
+}
+
+void
+GeometryBlocksOnCylindricalTests::run_tests()
+{
+  HighResWallClockTimer timer;
+  timer.start();
+  run_monotonic_axial_coordinates_in_detector_map_test();
+  timer.stop();
+}
+END_NAMESPACE_STIR
+
+USING_NAMESPACE_STIR
+
+int
+main()
+{
+  Verbosity::set(1);
+  GeometryBlocksOnCylindricalTests tests;
+  tests.run_tests();
+  return tests.main_return_value();
+}

src/test/test_interpolate_projdata.cxx

@@ -496,7 +496,7 @@ InterpolationTests::scatter_interpolation_test_blocks_asymmetric()
                                      "BlocksOnCylindrical",
                                      10.0,
                                      16.0,
-                                     60.0,
+                                     120.0,
                                      96.0);
 
   auto proj_data_info = shared_ptr<ProjDataInfo>(