From a4a7f5c0527d11da90d5bb5eb858ebf609a4e657 Mon Sep 17 00:00:00 2001 From: rts Date: Wed, 7 Feb 2024 11:49:56 -0800 Subject: [PATCH 01/10] run_back_projection_test_with_axial_buckets --- .../test_blocks_on_cylindrical_projectors.cxx | 88 +++++++++++++++++++ 1 file changed, 88 insertions(+) diff --git a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx index 2d07afcbd..75c167c72 100644 --- a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx +++ b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx @@ -53,6 +53,7 @@ #include "stir/recon_buildblock/BackProjectorByBinUsingProjMatrixByBin.h" #include "stir/IO/write_to_file.h" #include "stir/VoxelsOnCartesianGrid.h" +#include "stir/PixelsOnCartesianGrid.h" //#include "stir/Shape/Shape3D.h" #include @@ -80,6 +81,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(); }; /*! The following is a function to allow a projdata_info BlocksOnCylindrical to be created from the scanner. */ @@ -871,6 +873,90 @@ BlocksTests::run_intersection_with_cylinder_test() } } +void +BlocksTests::run_back_projection_test_with_axial_buckets() + +{ + // create projadata info + auto scanner_sptr = std::make_shared(Scanner::SAFIRDualRingPrototype); + { + 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(2 * 5); // More than 1 bucket + } + auto segments = scanner_sptr->get_num_rings() - 1; + VectorWithOffset num_axial_pos_per_segment(-segments, segments); + VectorWithOffset min_ring_diff_v(-segments, segments); + VectorWithOffset max_ring_diff_v(-segments, segments); + for (int i = -segments; i <= segments; i++) + { + min_ring_diff_v[i] = i; + max_ring_diff_v[i] = i; + num_axial_pos_per_segment[i] = scanner_sptr->get_num_rings() - abs(i); + } + + scanner_sptr->set_scanner_geometry("BlocksOnCylindrical"); + scanner_sptr->set_up(); + auto projdata_info_sptr + = std::make_shared(scanner_sptr, + num_axial_pos_per_segment, + min_ring_diff_v, + max_ring_diff_v, + scanner_sptr->get_max_num_views(), + scanner_sptr->get_max_num_non_arccorrected_bins()); + // scanner_sptr->set_scanner_geometry("Cylindrical"); + // scanner_sptr->set_up(); + // auto projdata_info_sptr = std::make_shared(scanner_sptr, + // num_axial_pos_per_segment, + // min_ring_diff_v, + // max_ring_diff_v, + // scanner_sptr->get_max_num_views(), + // scanner_sptr->get_max_num_non_arccorrected_bins()); + + auto projdata_sptr = std::make_shared(std::make_shared(ImagingModality::PT), projdata_info_sptr); + + auto origin = CartesianCoordinate3D(0, 0, 0); + auto volume_dimensions = CartesianCoordinate3D(-1, -1, -1); + auto volume = std::make_shared>(*projdata_info_sptr, 1, origin, volume_dimensions); + + // Now run the test + volume->fill(0.0); + projdata_sptr->fill(1.0); + + auto PM = std::make_shared(); + PM->enable_cache(false); + auto back_projector = std::make_shared(PM); + back_projector->set_up(projdata_info_sptr, volume); + back_projector->back_project(*volume, *projdata_sptr, 0, 144); + + bool test_ok = true; + std::ostringstream oss; + oss << "BlocksTests::run_back_projection_test_with_axial_buckets: Central voxel values are:" << std::endl; + + auto center_axial_values = std::vector(volume->get_z_size()); + for (int z = volume->get_min_z(); z <= volume->get_max_z(); z++) + { + center_axial_values[z] = volume->get_plane(z).at(0).at(0); + oss << "\tz = " << z << "/" << volume->get_max_z() << " is " << center_axial_values[z] << std::endl; + if (test_ok) + { + test_ok = check(center_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 +991,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(); + print_time("back projection test with axial buckets took: "); #ifdef STIR_WITH_Parallelproj_PROJECTOR // run the same tests with parallelproj, if available From d5e958eab865ba05477f4531418631f01a0d99a9 Mon Sep 17 00:00:00 2001 From: rts Date: Thu, 8 Feb 2024 06:31:17 -0800 Subject: [PATCH 02/10] Add a guard for axial crystals per bucket being a factor of num_rings --- .../GeometryBlocksOnCylindrical.cxx | 24 +++++++++++++++---- .../stir/GeometryBlocksOnCylindrical.h | 3 +++ 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/src/buildblock/GeometryBlocksOnCylindrical.cxx b/src/buildblock/GeometryBlocksOnCylindrical.cxx index 81f3fd014..548cd0954 100644 --- a/src/buildblock/GeometryBlocksOnCylindrical.cxx +++ b/src/buildblock/GeometryBlocksOnCylindrical.cxx @@ -44,6 +44,7 @@ START_NAMESPACE_STIR GeometryBlocksOnCylindrical::GeometryBlocksOnCylindrical(const Scanner& scanner) { + check_scanner_configuration(scanner); build_crystal_maps(scanner); } @@ -55,17 +56,32 @@ GeometryBlocksOnCylindrical::get_rotation_matrix(float alpha) const stir::make_1d_array(0.F, -1 * std::sin(alpha), std::cos(alpha))); } +void +GeometryBlocksOnCylindrical::check_scanner_configuration(const Scanner& scanner) +{ + { // Ensure the number of axial crystals per bucket is a multiple of the number of rings + int num_axial_crystals_per_bucket = scanner.get_num_axial_crystals_per_block() * scanner.get_num_axial_blocks_per_bucket(); + if (scanner.get_num_rings() % (num_axial_crystals_per_bucket) != 0) + { + error("Error in GeometryBlocksOnCylindrical: number of rings (%d) is not a multiple of the num_axial_crystals_per_bucket " + "(%d) = num_axial_crystals_per_block (%d) * num_axial_blocks_per_bucket (%d)", + scanner.get_num_rings(), + num_axial_crystals_per_bucket, + scanner.get_num_axial_crystals_per_block(), + scanner.get_num_axial_blocks_per_bucket()); + } + } +} void 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 +116,7 @@ GeometryBlocksOnCylindrical::build_crystal_maps(const Scanner& scanner) stir::CartesianCoordinate3D 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) diff --git a/src/include/stir/GeometryBlocksOnCylindrical.h b/src/include/stir/GeometryBlocksOnCylindrical.h index 84677de4f..0a55b3709 100644 --- a/src/include/stir/GeometryBlocksOnCylindrical.h +++ b/src/include/stir/GeometryBlocksOnCylindrical.h @@ -56,6 +56,9 @@ class GeometryBlocksOnCylindrical : public DetectorCoordinateMap //! Get rotation matrix for a given angle around z axis stir::Array<2, float> get_rotation_matrix(float alpha) const; + //! Checks the scanner configuration is valid for a crystal map + static void check_scanner_configuration(const Scanner& scanner) ; + //! Build crystal map in cartesian coordinate void build_crystal_maps(const Scanner& scanner); }; From 8cd047dce58937c376c928f1a549a36227c153fa Mon Sep 17 00:00:00 2001 From: Robert Twyman Skelly Date: Thu, 8 Feb 2024 11:41:31 -0800 Subject: [PATCH 03/10] [ci skip] Refactor and renaming --- .../GeometryBlocksOnCylindrical.cxx | 12 +++---- .../test_blocks_on_cylindrical_projectors.cxx | 34 +++++++++---------- 2 files changed, 22 insertions(+), 24 deletions(-) diff --git a/src/buildblock/GeometryBlocksOnCylindrical.cxx b/src/buildblock/GeometryBlocksOnCylindrical.cxx index 548cd0954..9104397ec 100644 --- a/src/buildblock/GeometryBlocksOnCylindrical.cxx +++ b/src/buildblock/GeometryBlocksOnCylindrical.cxx @@ -39,6 +39,7 @@ limitations under the License. #include #include #include +#include START_NAMESPACE_STIR @@ -63,12 +64,11 @@ GeometryBlocksOnCylindrical::check_scanner_configuration(const Scanner& scanner) int num_axial_crystals_per_bucket = scanner.get_num_axial_crystals_per_block() * scanner.get_num_axial_blocks_per_bucket(); if (scanner.get_num_rings() % (num_axial_crystals_per_bucket) != 0) { - error("Error in GeometryBlocksOnCylindrical: number of rings (%d) is not a multiple of the num_axial_crystals_per_bucket " - "(%d) = num_axial_crystals_per_block (%d) * num_axial_blocks_per_bucket (%d)", - scanner.get_num_rings(), - num_axial_crystals_per_bucket, - scanner.get_num_axial_crystals_per_block(), - scanner.get_num_axial_blocks_per_bucket()); + error(boost::format("Error in GeometryBlocksOnCylindrical: number of rings (%d) is not a multiple of the " + "num_axial_crystals_per_bucket " + "(%d) = num_axial_crystals_per_block (%d) * num_axial_blocks_per_bucket (%d)") + % scanner.get_num_rings() % num_axial_crystals_per_bucket % scanner.get_num_axial_crystals_per_block() + % scanner.get_num_axial_blocks_per_bucket()); } } } diff --git a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx index 75c167c72..cc72c41a1 100644 --- a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx +++ b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx @@ -52,9 +52,7 @@ #endif #include "stir/recon_buildblock/BackProjectorByBinUsingProjMatrixByBin.h" #include "stir/IO/write_to_file.h" -#include "stir/VoxelsOnCartesianGrid.h" #include "stir/PixelsOnCartesianGrid.h" -//#include "stir/Shape/Shape3D.h" #include START_NAMESPACE_STIR @@ -81,7 +79,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(); + 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. */ @@ -874,7 +872,7 @@ BlocksTests::run_intersection_with_cylinder_test() } void -BlocksTests::run_back_projection_test_with_axial_buckets() +BlocksTests::run_back_projection_test_with_axial_buckets(BackProjectorByBin& back_projector) { // create projadata info @@ -887,7 +885,7 @@ BlocksTests::run_back_projection_test_with_axial_buckets() 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(2 * 5); // More than 1 bucket + scanner_sptr->set_num_rings(2 * 5); // More than 1 bucket } auto segments = scanner_sptr->get_num_rings() - 1; VectorWithOffset num_axial_pos_per_segment(-segments, segments); @@ -918,34 +916,33 @@ BlocksTests::run_back_projection_test_with_axial_buckets() // scanner_sptr->get_max_num_views(), // scanner_sptr->get_max_num_non_arccorrected_bins()); - auto projdata_sptr = std::make_shared(std::make_shared(ImagingModality::PT), projdata_info_sptr); + auto exam_info_sptr = std::make_shared(ImagingModality::PT); + auto projdata_sptr = std::make_shared(exam_info_sptr, projdata_info_sptr); auto origin = CartesianCoordinate3D(0, 0, 0); auto volume_dimensions = CartesianCoordinate3D(-1, -1, -1); - auto volume = std::make_shared>(*projdata_info_sptr, 1, origin, volume_dimensions); + auto volume_sptr + = std::make_shared>(exam_info_sptr, *projdata_info_sptr, 1, origin, volume_dimensions); // Now run the test - volume->fill(0.0); + volume_sptr->fill(0.0); projdata_sptr->fill(1.0); - auto PM = std::make_shared(); - PM->enable_cache(false); - auto back_projector = std::make_shared(PM); - back_projector->set_up(projdata_info_sptr, volume); - back_projector->back_project(*volume, *projdata_sptr, 0, 144); + 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 center_axial_values = std::vector(volume->get_z_size()); - for (int z = volume->get_min_z(); z <= volume->get_max_z(); z++) + auto centre_axial_values = std::vector(volume_sptr->get_z_size()); + for (int z = volume_sptr->get_min_z(); z <= volume_sptr->get_max_z(); z++) { - center_axial_values[z] = volume->get_plane(z).at(0).at(0); - oss << "\tz = " << z << "/" << volume->get_max_z() << " is " << center_axial_values[z] << std::endl; + centre_axial_values[z] = volume_sptr->get_plane(z).at(0).at(0); + oss << "\tz = " << z << "/" << volume_sptr->get_max_z() << " is " << centre_axial_values[z] << std::endl; if (test_ok) { - test_ok = check(center_axial_values[z] > 0, + 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; } @@ -992,6 +989,7 @@ BlocksTests::run_tests() run_intersection_with_cylinder_test(); print_time("intersection with cylinder test took: "); run_back_projection_test_with_axial_buckets(); + run_back_projection_test_with_axial_buckets(back_projector); print_time("back projection test with axial buckets took: "); #ifdef STIR_WITH_Parallelproj_PROJECTOR From bab0cddbc7395390ba24d72e7b0ee674c317b350 Mon Sep 17 00:00:00 2001 From: rts Date: Mon, 12 Feb 2024 09:42:43 -0800 Subject: [PATCH 04/10] [ci skip] Remove duplicate test and access volume array without PixelsOnCartesianGrid --- src/include/stir/GeometryBlocksOnCylindrical.h | 2 +- src/recon_test/test_blocks_on_cylindrical_projectors.cxx | 4 +--- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/include/stir/GeometryBlocksOnCylindrical.h b/src/include/stir/GeometryBlocksOnCylindrical.h index 0a55b3709..32cf54c57 100644 --- a/src/include/stir/GeometryBlocksOnCylindrical.h +++ b/src/include/stir/GeometryBlocksOnCylindrical.h @@ -57,7 +57,7 @@ class GeometryBlocksOnCylindrical : public DetectorCoordinateMap stir::Array<2, float> get_rotation_matrix(float alpha) const; //! Checks the scanner configuration is valid for a crystal map - static void check_scanner_configuration(const Scanner& scanner) ; + static void check_scanner_configuration(const Scanner& scanner); //! Build crystal map in cartesian coordinate void build_crystal_maps(const Scanner& scanner); diff --git a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx index cc72c41a1..7e0fce195 100644 --- a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx +++ b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx @@ -52,7 +52,6 @@ #endif #include "stir/recon_buildblock/BackProjectorByBinUsingProjMatrixByBin.h" #include "stir/IO/write_to_file.h" -#include "stir/PixelsOnCartesianGrid.h" #include START_NAMESPACE_STIR @@ -938,7 +937,7 @@ BlocksTests::run_back_projection_test_with_axial_buckets(BackProjectorByBin& bac auto centre_axial_values = std::vector(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->get_plane(z).at(0).at(0); + 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) { @@ -988,7 +987,6 @@ 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(); run_back_projection_test_with_axial_buckets(back_projector); print_time("back projection test with axial buckets took: "); From 4c4c857f606f620a4f5979cb211fe3b01db03e8c Mon Sep 17 00:00:00 2001 From: rts Date: Mon, 12 Feb 2024 11:32:58 -0800 Subject: [PATCH 05/10] Require GeometryBlocksOnCylindrical to check scanner consistency and error if succeed is no --- .../GeometryBlocksOnCylindrical.cxx | 18 ++---------------- src/buildblock/Scanner.cxx | 19 +++++++++++++++++++ .../stir/GeometryBlocksOnCylindrical.h | 5 ++--- .../test_blocks_on_cylindrical_projectors.cxx | 6 ++++-- 4 files changed, 27 insertions(+), 21 deletions(-) diff --git a/src/buildblock/GeometryBlocksOnCylindrical.cxx b/src/buildblock/GeometryBlocksOnCylindrical.cxx index 9104397ec..eb13fb234 100644 --- a/src/buildblock/GeometryBlocksOnCylindrical.cxx +++ b/src/buildblock/GeometryBlocksOnCylindrical.cxx @@ -45,7 +45,8 @@ START_NAMESPACE_STIR GeometryBlocksOnCylindrical::GeometryBlocksOnCylindrical(const Scanner& scanner) { - check_scanner_configuration(scanner); + if (scanner.check_consistency() == Succeeded::no) + error("Error in GeometryBlocksOnCylindrical: scanner configuration not accepted. Please check warnings."); build_crystal_maps(scanner); } @@ -57,21 +58,6 @@ GeometryBlocksOnCylindrical::get_rotation_matrix(float alpha) const stir::make_1d_array(0.F, -1 * std::sin(alpha), std::cos(alpha))); } -void -GeometryBlocksOnCylindrical::check_scanner_configuration(const Scanner& scanner) -{ - { // Ensure the number of axial crystals per bucket is a multiple of the number of rings - int num_axial_crystals_per_bucket = scanner.get_num_axial_crystals_per_block() * scanner.get_num_axial_blocks_per_bucket(); - if (scanner.get_num_rings() % (num_axial_crystals_per_bucket) != 0) - { - error(boost::format("Error in GeometryBlocksOnCylindrical: number of rings (%d) is not a multiple of the " - "num_axial_crystals_per_bucket " - "(%d) = num_axial_crystals_per_block (%d) * num_axial_blocks_per_bucket (%d)") - % scanner.get_num_rings() % num_axial_crystals_per_bucket % scanner.get_num_axial_crystals_per_block() - % scanner.get_num_axial_blocks_per_bucket()); - } - } -} void GeometryBlocksOnCylindrical::build_crystal_maps(const Scanner& scanner) { diff --git a/src/buildblock/Scanner.cxx b/src/buildblock/Scanner.cxx index a66b5c3d4..abe2ef02e 100644 --- a/src/buildblock/Scanner.cxx +++ b/src/buildblock/Scanner.cxx @@ -42,6 +42,7 @@ #include #include #include +#include #include "stir/warning.h" #include "stir/error.h" @@ -1751,6 +1752,24 @@ 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; + } + { // Ensure the number of axial crystals per bucket is a multiple of the number of rings + // This asserts the number of axial + if (get_num_rings() % get_num_axial_crystals_per_bucket() != 0) + { + error(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( diff --git a/src/include/stir/GeometryBlocksOnCylindrical.h b/src/include/stir/GeometryBlocksOnCylindrical.h index 32cf54c57..f00636b6a 100644 --- a/src/include/stir/GeometryBlocksOnCylindrical.h +++ b/src/include/stir/GeometryBlocksOnCylindrical.h @@ -1,6 +1,7 @@ /* Copyright 2017 ETH Zurich, Institute of Particle Physics and Astrophysics + Copyright (C) 2024, Prescient Imaging Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -22,6 +23,7 @@ limitations under the License. \brief Declaration of class stir::GeometryBlocksOnCylindrical \author Parisa Khateri + \author Robert Twyman */ #ifndef __stir_GeometryBlocksOnCylindrical_H__ @@ -56,9 +58,6 @@ class GeometryBlocksOnCylindrical : public DetectorCoordinateMap //! Get rotation matrix for a given angle around z axis stir::Array<2, float> get_rotation_matrix(float alpha) const; - //! Checks the scanner configuration is valid for a crystal map - static void check_scanner_configuration(const Scanner& scanner); - //! Build crystal map in cartesian coordinate void build_crystal_maps(const Scanner& scanner); }; diff --git a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx index 7e0fce195..52f70809e 100644 --- a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx +++ b/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 @@ -987,8 +989,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: "); +// 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 From 815c75e4be6f87ee993f59b779bb6b9b05ede319 Mon Sep 17 00:00:00 2001 From: rts Date: Mon, 12 Feb 2024 11:33:32 -0800 Subject: [PATCH 06/10] Fix a logic bug with get_scanner_geometry() == "Generic" --- src/buildblock/Scanner.cxx | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/buildblock/Scanner.cxx b/src/buildblock/Scanner.cxx index abe2ef02e..27bfb495d 100644 --- a/src/buildblock/Scanner.cxx +++ b/src/buildblock/Scanner.cxx @@ -1805,22 +1805,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; - } - } } } From 65ae1f80778ae588cbdb4b5fc90cabb3affcb1e5 Mon Sep 17 00:00:00 2001 From: rts Date: Mon, 12 Feb 2024 11:34:25 -0800 Subject: [PATCH 07/10] Fix an axial_block_spacing issue with a downsampled_scanner test --- src/test/test_interpolate_projdata.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/test/test_interpolate_projdata.cxx b/src/test/test_interpolate_projdata.cxx index 46925eba9..5a4e27029 100644 --- a/src/test/test_interpolate_projdata.cxx +++ b/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( From 57fb97b8d58a8858df6ab0aa631c5623b8ea3998 Mon Sep 17 00:00:00 2001 From: rts Date: Tue, 13 Feb 2024 11:34:42 -0800 Subject: [PATCH 08/10] [ci skip] Allow test to run and reduce code duplication --- .../test_blocks_on_cylindrical_projectors.cxx | 41 ++++--------------- 1 file changed, 8 insertions(+), 33 deletions(-) diff --git a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx index 52f70809e..c68144069 100644 --- a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx +++ b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx @@ -876,9 +876,9 @@ void BlocksTests::run_back_projection_test_with_axial_buckets(BackProjectorByBin& back_projector) { - // create projadata info + int num_buckets = 1; auto scanner_sptr = std::make_shared(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() @@ -886,37 +886,12 @@ BlocksTests::run_back_projection_test_with_axial_buckets(BackProjectorByBin& bac 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(2 * 5); // More than 1 bucket + scanner_sptr->set_num_rings(scanner_sptr->get_num_axial_blocks_per_bucket() * num_buckets); + scanner_sptr->set_scanner_geometry("BlocksOnCylindrical"); + scanner_sptr->set_up(); } - auto segments = scanner_sptr->get_num_rings() - 1; - VectorWithOffset num_axial_pos_per_segment(-segments, segments); - VectorWithOffset min_ring_diff_v(-segments, segments); - VectorWithOffset max_ring_diff_v(-segments, segments); - for (int i = -segments; i <= segments; i++) - { - min_ring_diff_v[i] = i; - max_ring_diff_v[i] = i; - num_axial_pos_per_segment[i] = scanner_sptr->get_num_rings() - abs(i); - } - - scanner_sptr->set_scanner_geometry("BlocksOnCylindrical"); - scanner_sptr->set_up(); - auto projdata_info_sptr - = std::make_shared(scanner_sptr, - num_axial_pos_per_segment, - min_ring_diff_v, - max_ring_diff_v, - scanner_sptr->get_max_num_views(), - scanner_sptr->get_max_num_non_arccorrected_bins()); - // scanner_sptr->set_scanner_geometry("Cylindrical"); - // scanner_sptr->set_up(); - // auto projdata_info_sptr = std::make_shared(scanner_sptr, - // num_axial_pos_per_segment, - // min_ring_diff_v, - // max_ring_diff_v, - // scanner_sptr->get_max_num_views(), - // scanner_sptr->get_max_num_non_arccorrected_bins()); + auto projdata_info_sptr = set_direct_projdata_info(scanner_sptr, 1); auto exam_info_sptr = std::make_shared(ImagingModality::PT); auto projdata_sptr = std::make_shared(exam_info_sptr, projdata_info_sptr); @@ -989,8 +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: "); + 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 From 070c06a0de5f1669ea2551b65d7ce989f5394026 Mon Sep 17 00:00:00 2001 From: rts Date: Tue, 13 Feb 2024 14:28:19 -0800 Subject: [PATCH 09/10] Add GeometryBlocksOnCylindricalTests --- src/recon_test/CMakeLists.txt | 1 + .../test_geometry_blocks_on_cylindrical.cxx | 151 ++++++++++++++++++ 2 files changed, 152 insertions(+) create mode 100644 src/recon_test/test_geometry_blocks_on_cylindrical.cxx diff --git a/src/recon_test/CMakeLists.txt b/src/recon_test/CMakeLists.txt index 5f82248fd..f7fe11631 100644 --- a/src/recon_test/CMakeLists.txt +++ b/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 ) diff --git a/src/recon_test/test_geometry_blocks_on_cylindrical.cxx b/src/recon_test/test_geometry_blocks_on_cylindrical.cxx new file mode 100644 index 000000000..e88c4e08d --- /dev/null +++ b/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 + +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_sptr); +}; + +void +GeometryBlocksOnCylindricalTests::run_monotonic_axial_coordinates_in_detector_map_test() +{ + auto scanner_sptr = std::make_shared(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_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 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::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 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(); +} \ No newline at end of file From 1bb2af0e065d1ca9f9cd4b1f0f9953663c19e010 Mon Sep 17 00:00:00 2001 From: rts Date: Tue, 13 Feb 2024 14:28:56 -0800 Subject: [PATCH 10/10] [ci skip] Small fixes --- documentation/release_6.1.htm | 6 ++++-- src/buildblock/Scanner.cxx | 15 +++++++-------- src/include/stir/GeometryBlocksOnCylindrical.h | 2 -- .../test_blocks_on_cylindrical_projectors.cxx | 2 +- 4 files changed, 12 insertions(+), 13 deletions(-) diff --git a/documentation/release_6.1.htm b/documentation/release_6.1.htm index 69cc7a310..0531a4677 100644 --- a/documentation/release_6.1.htm +++ b/documentation/release_6.1.htm @@ -84,8 +84,10 @@

What's new for developers (aside from what should be obvious

Backward incompatibities

    -
  • -
    • +
  • + Additional checks on GeometryBlocksOnCylindrical scanner configuration, which may lead to an + error being raised, while previously the code silently proceeded. +

New functionality

diff --git a/src/buildblock/Scanner.cxx b/src/buildblock/Scanner.cxx index 27bfb495d..b66dee692 100644 --- a/src/buildblock/Scanner.cxx +++ b/src/buildblock/Scanner.cxx @@ -1754,20 +1754,19 @@ Scanner::check_consistency() const { //! 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." + 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; } - { // Ensure the number of axial crystals per bucket is a multiple of the number of rings - // This asserts the number of axial + { // Assert that each block contains an equal number of axial crystals if (get_num_rings() % get_num_axial_crystals_per_bucket() != 0) { - error(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()); + 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()) diff --git a/src/include/stir/GeometryBlocksOnCylindrical.h b/src/include/stir/GeometryBlocksOnCylindrical.h index f00636b6a..84677de4f 100644 --- a/src/include/stir/GeometryBlocksOnCylindrical.h +++ b/src/include/stir/GeometryBlocksOnCylindrical.h @@ -1,7 +1,6 @@ /* Copyright 2017 ETH Zurich, Institute of Particle Physics and Astrophysics - Copyright (C) 2024, Prescient Imaging Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -23,7 +22,6 @@ limitations under the License. \brief Declaration of class stir::GeometryBlocksOnCylindrical \author Parisa Khateri - \author Robert Twyman */ #ifndef __stir_GeometryBlocksOnCylindrical_H__ diff --git a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx index c68144069..773a37d97 100644 --- a/src/recon_test/test_blocks_on_cylindrical_projectors.cxx +++ b/src/recon_test/test_blocks_on_cylindrical_projectors.cxx @@ -886,7 +886,7 @@ BlocksTests::run_back_projection_test_with_axial_buckets(BackProjectorByBin& bac 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_blocks_per_bucket() * num_buckets); + 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(); }