STIR
STIR copied to clipboard
UCL
Scanner number of axial components is not logical for BlocksOnCylindrical
Issue
I have have been experiencing an issue related to the number of axial crystals in a scanner when the wrong num_rings value is input. I had the rough configuration
num_axial_crystals_per_block = 1
num_axial_blocks_per_bucket = 2
num_rings = 2 * 3 // Attempting 3 axial buckets
The issue I found is that everything is created cleanly but when back projecting a uniform 1s to a generated volume, the back projection only occurred in 1/3 of the volume. Note, this is only an issue with BlocksOnCylindrical.
CTest Example
Created a test to demonstrate the issue:
void
BlocksTests::run_my_test()
{
// create projadata info
auto scannerBlocks_sptr = std::make_shared<Scanner>(Scanner::SAFIRDualRingPrototype);
{
scannerBlocks_sptr->set_average_depth_of_interaction(5);
scannerBlocks_sptr->set_num_axial_crystals_per_block(1);
scannerBlocks_sptr->set_axial_block_spacing(scannerBlocks_sptr->get_axial_crystal_spacing()
* scannerBlocks_sptr->get_num_axial_crystals_per_block());
scannerBlocks_sptr->set_transaxial_block_spacing(scannerBlocks_sptr->get_transaxial_crystal_spacing()
* scannerBlocks_sptr->get_num_transaxial_crystals_per_block());
scannerBlocks_sptr->set_num_axial_blocks_per_bucket(2);
scannerBlocks_sptr->set_num_rings(2 * 3); // Attempting 3 axial buckets
scannerBlocks_sptr->set_scanner_geometry("BlocksOnCylindrical");
scannerBlocks_sptr->set_up();
}
auto proj_data_info_blocks_sptr = std::make_shared<ProjDataInfoBlocksOnCylindricalNoArcCorr>();
proj_data_info_blocks_sptr = set_blocks_projdata_info<ProjDataInfoBlocksOnCylindricalNoArcCorr>(scannerBlocks_sptr);
auto projdata_sptr
= std::make_shared<ProjDataInMemory>(std::make_shared<ExamInfo>(ImagingModality::PT), proj_data_info_blocks_sptr);
auto origin = CartesianCoordinate3D<float>(0, 0, 0);
auto volume_dimensions = CartesianCoordinate3D<int>(-1, -1, -1);
auto volume = std::make_shared<VoxelsOnCartesianGrid<float>>(*proj_data_info_blocks_sptr, 1, origin, volume_dimensions);
// Now run the test
volume->fill(0.0);
projdata_sptr->fill(1.0);
auto PM = std::make_shared<ProjMatrixByBinUsingRayTracing>();
PM->enable_cache(false);
auto back_projector = std::make_shared<BackProjectorByBinUsingProjMatrixByBin>(PM);
back_projector->set_up(proj_data_info_blocks_sptr, volume);
back_projector->back_project(*volume, *projdata_sptr, 0, 144);
auto center_axial_values = std::vector<float>(volume->get_z_size());
for (int z = volume->get_min_z(); z <= volume->get_max_z(); z++)
{
auto value = volume->get_plane(z).at(0).at(0);
center_axial_values[z] = value;
std::cout << "Central voxel value at z = " << z << "/" << volume->get_max_z() << " is " << value << std::endl;
}
}
Output
WARNING: FactoryRegistry:: overwriting previous value of key in registry.
key: None
INFO: Determined voxel size by dividing default_bin_size (1.1) by zoom
WARNING: Disabling all symmetries except for symmtery_z since they are not implemented in block geometry yet.
WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (x,y)=(1.1,1.1) that is smaller than the central bin size (3.00691) divided by num_tangential_LORs (1).
This matrix will completely miss some voxels for some (or all) views. It is therefore to best to increase 'number of rays in tangential direction to trace for each bin'.
INFO: Processing view 0 of segment -5, TOF bin 0
INFO: Processing view 0 of segment -4, TOF bin 0
INFO: Processing view 0 of segment -3, TOF bin 0
INFO: Processing view 0 of segment -2, TOF bin 0
INFO: Processing view 0 of segment -1, TOF bin 0
INFO: Processing view 0 of segment 0, TOF bin 0
INFO: Processing view 0 of segment 1, TOF bin 0
INFO: Processing view 0 of segment 2, TOF bin 0
INFO: Processing view 0 of segment 3, TOF bin 0
INFO: Processing view 0 of segment 4, TOF bin 0
INFO: Processing view 0 of segment 5, TOF bin 0
Central voxel value at z = 0/10 is 9.00096
Central voxel value at z = 1/10 is 13.5011
Central voxel value at z = 2/10 is 9.0001
Central voxel value at z = 3/10 is 2.25
Central voxel value at z = 4/10 is 0
Central voxel value at z = 5/10 is 0
Central voxel value at z = 6/10 is 0
Central voxel value at z = 7/10 is 0
Central voxel value at z = 8/10 is 0
Central voxel value at z = 9/10 is 0
Central voxel value at z = 10/10 is 0
Note there is no warning or errors that indicate an issue with the block geometry. I think part of the issue is there is no validation to match the num_rings to the number of axial crystals. Could this be because of the complications of virtual crystals in other scanners?
I have been looking at GeometryBlocksOnCylindrical::build_crystal_maps and I am wondering if there is something wrong with the loop ranges.
https://github.com/UCL/STIR/blob/ee307c721cae562ab163aee06dbab8407d861414/src/buildblock/GeometryBlocksOnCylindrical.cxx#L102-L104
Should the second loop terminate at num_axial_blocks_per_bucket = 2, not num_axial_blocks=6. num_axial_buckets = 3 and is derived by
https://github.com/UCL/STIR/blob/ee307c721cae562ab163aee06dbab8407d861414/src/buildblock/GeometryBlocksOnCylindrical.cxx#L68
using the num_rings
https://github.com/UCL/STIR/blob/ee307c721cae562ab163aee06dbab8407d861414/src/include/stir/Scanner.inl#L161-L167
I agree, it looks like the second loop should be be looping over the blocks per bucket. I wonder why I never ran into any issues with this - will look at it when I have a moment.
Nonetheless, I agree it would be very nice to have some errors if the scanner definition contains conflicting parameters.
I wonder why I never ran into any issues with this
I believe that most of the tests use a single axial bucket and that doesn't encounter any issue.
I will create a PR with the aforementioned test and an attempted fix today.