Skip to content

Commit

Permalink
Merge pull request #1375 from robbietuk/axial_bucket_projection_issue
Browse files Browse the repository at this point in the history
Throw errors if incorrect/unsupported blocks scanner. See #1374
  • Loading branch information
KrisThielemans authored Feb 15, 2024
2 parents bc56456 + 1bb2af0 commit 79a1e83
Show file tree
Hide file tree
Showing 7 changed files with 257 additions and 22 deletions.
6 changes: 4 additions & 2 deletions documentation/release_6.1.htm
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down
10 changes: 6 additions & 4 deletions src/buildblock/GeometryBlocksOnCylindrical.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

Expand All @@ -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();
Expand Down Expand Up @@ -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)
Expand Down
44 changes: 31 additions & 13 deletions src/buildblock/Scanner.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@
#include <iostream>
#include <algorithm>
#include <sstream>
#include <boost/format.hpp>
#include "stir/warning.h"
#include "stir/error.h"

Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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;
}
}
}
}

Expand Down
1 change: 1 addition & 0 deletions src/recon_test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
)


Expand Down
65 changes: 63 additions & 2 deletions src/recon_test/test_blocks_on_cylindrical_projectors.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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. */
Expand Down Expand Up @@ -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()
{
Expand Down Expand Up @@ -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
Expand Down
151 changes: 151 additions & 0 deletions src/recon_test/test_geometry_blocks_on_cylindrical.cxx
Original file line number Diff line number Diff line change
@@ -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();
}
2 changes: 1 addition & 1 deletion src/test/test_interpolate_projdata.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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>(
Expand Down

0 comments on commit 79a1e83

Please sign in to comment.