Scanner number of axial components is not logical for BlocksOnCylindrical

[robbietuk]

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?

[robbietuk]

I have been looking at GeometryBlocksOnCylindrical::build_crystal_maps and I am wondering if there is something wrong with the loop ranges.

STIR/src/buildblock/GeometryBlocksOnCylindrical.cxx

Lines 102 to 104 in ee307c7

	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_crys_num = 0; ax_crys_num < num_axial_crystals_per_block; ++ax_crys_num)

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

STIR/src/buildblock/GeometryBlocksOnCylindrical.cxx

Line 68 in ee307c7

int num_axial_buckets = scanner.get_num_axial_blocks() / num_axial_blocks_per_bucket;

using the num_rings

STIR/src/include/stir/Scanner.inl

Lines 161 to 167 in ee307c7

	int
	Scanner::get_num_axial_blocks() const
	{
	// when using virtual crystals between blocks, there won't be one at the end, so we
	// need to take this into account.
	return (num_rings + get_num_virtual_axial_crystals_per_block()) / num_axial_crystals_per_block;
	}

[KrisThielemans]

@danieldeidda @markus-jehl

[markus-jehl]

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.

[robbietuk]

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.

[robbietuk]

I will create a PR with the aforementioned test and an attempted fix today.

[markus-jehl]

Yep, that's it - I never use more than one bucket axially. Thanks!

[robbietuk]

The conversion has been moved to #1374