Axial bucket projection issue

[robbietuk]

Changes in this pull request

In response to #1374, add an integration CTest that ensures BlocksOnCylindrical with more than 1 axial bucket back project into all z slices of a volume derived from the projection data.

Add a guard to ensure num_rings % num_crystals_per_block == 0.

Testing performed

Added a new CTest

Related issues

#1374

Checklist before requesting a review

• [] I have performed a self-review of my code
• [] I have added docstrings/doxygen in line with the guidance in the developer guide
• [] I have implemented unit tests that cover any new or modified functionality (if applicable)
• [] The code builds and runs on my machine
• [] documentation/release_XXX.md has been updated with any functionality change (if applicable)

[robbietuk]

@markus-jehl, the modification to use num_axial_blocks_per_bucket as the upper bound of the loop didn't resolve the issue. There must be something else. Investigating...

[KrisThielemans]

house-keeping comments, not about the actual methodology (which does look fine to me)

[robbietuk]

I believe the issue described in #1374 is caused by the axial bucket spacing is not taken into account in the z component of the transformation_matrix (should be translation_matrix):

STIR/src/buildblock/GeometryBlocksOnCylindrical.cxx

Lines 123 to 127 in ee307c7

	// calculate cartesian coordinate for a given detector
	stir::CartesianCoordinate3D<float> transformation_matrix(
	ax_block_num * axial_block_spacing + ax_crys_num * axial_crystal_spacing,
	0.,
	trans_block_num * transaxial_block_spacing + trans_crys_num * transaxial_crystal_spacing);

In my tests, the cart_coord.z() was always less than 0. Detectors at higher axial_coord are assigned CartesianCoordinates z components that do not account for the "ax_bucket_spacing". However, this variable does not exist.

---

Is this an oversight or a feature? Crystals > Blocks > Bucket(s). BlocksOnCylindrical does not support multiple buckets, transaxially or axially as there is no *_bucket_spacing variables in Scanner. This is not a previous problem for a cylindrical PET scanner because the ring spacing and num_rings determined the number of axial buckets.

---

Possible solutions

• Add *_bucket_spacing variables in Scanner and implement usage.
  • This is the simplest solution and uses the current methodology
• Disable the use of multiple buckets per "module" and ensure each parallel plane is a single bucket (axially and transaxially).
  • This is limiting and may break a few things?
• Transition away from the use of num_rings to define the geometry and implement entirely from num_crystals_per_block, num_blocks_per_bucket and num_buckets_per_(module?) (in each dimension)
  • I personally perfer this method as it is the simplest for the user as they no longer need to precalculate the num_rings, STIR does it for them.

@KrisThielemans @markus-jehl @danieldeidda thoughts?

[KrisThielemans]

First of all, "buckets" consist of a number of "blocks" which have a number of "crystals". BlocksOnCylindrical was developed without thinking about buckets at all (corresponding to "1 block per bucket"), as that's what was needed at the time, hence the bug. Unfortunately, most scanners actually would correspond to "BucketsOnCylindrical". Furthermore, for most scanners it would make more sense to have a larger axial gap between buckets than between blocks. We currently don't have the variables for that.

#776 by @NikEfth tried to address some of this by adding a facility for extra virtual crystals between buckets. Sadly, this get stuck in changing meaning of various variables with a long-ish discussion. I cannot recall how far this got, and it will need revision now since various merges. In any case, that isn't completely general anyway as the extra spacing can only be a crystal.

The whole issue of arbitrary geometries and virtual crystals is a very difficult one. Getting the detector-map sorted is relatively easy, but things like symmetries (ok, we just switch them off), component-based normalisation, and down/upsampling for scatter and even nomenclature are hard.

While this should be rather high on our priority list, I suggest we don't try to resolve that in this PR, or discuss it here either. I think for now we should just fix the code such that it does a sensible thing (and document it). Luckily, the fact that the code is currently broken for num_axial_blocks_per_bucket > 1 means we are just bug fixing, and don't introduce extra incompatibility! 😄

Looking at the current definition,the alpha definition is in terms of buckets. So it will put all blocks in a bucket next/parallel to eachother. That looks good to me. So, I think that currently BlocksOnCylindrical does support multiple blocks per bucket in transaxial direction. However, clearly in axial direction the bucket location is completely ignored.

On to the proposed solutions. Keep in mind there are 2 different things: how do we get information (get_num_rings() is all over the place), and how do we define the scanner. The former is independent of the latter. The latter is where things get quite hard due to existing Interfile headers, and loads of scanner definitions. Changing the Interfile headers would need versions and all that. Not so easy therefore.

• Add *_bucket_spacing variables in Scanner and implement usage.

  • This is the simplest solution and uses the current methodology

Yes, this is rather easy. In the first instance, I would suggest we do this initially with "gaps between blocks" = "gaps between buckets". So, I suggest to have only a get_bucket_spacing() that computes it as num_axial_blocks_per_bucket*get_axial_block_spacing(). i.e this is only an additional the get. In the next step, making the gap larger needs extra variables/parsing etc, so it's a bit more work. (not much in implementation, probably loads of work in testing).

• Disable the use of multiple buckets per "module" and ensure each parallel plane is a single bucket (axially and transaxially).

  • This is limiting and may break a few things?

"module" is not STIR terminology ATM, so I have no idea what you mean @robbietuk

• Transition away from the use of num_rings to define the geometry and implement entirely from num_crystals_per_block, num_blocks_per_bucket and num_buckets_per_(module?) (in each dimension)

  • I personally perfer this method as it is the simplest for the user as they no longer need to precalculate the num_rings, STIR does it for them.

This is a set change. This is a lot harder as discussed above.

So, I prefer the first option. It's a fairly easy fix, and the 3rd option can then be done later. (Thinking how we add virtual crystals to all this gives me a headache. They are a pain, but they do enable easy use of symmetries...)

[KrisThielemans]

extra comments, also, don't forget to add you as \author, and when done release_6.1

[robbietuk]

I agree with all of the above (#1375 (comment)). The third option requires reconsidering of "how do we define the scanner" which would require a "change to existing Interfile headers, and loads of scanner definitions". This is uncomfortable and I don't have time.

I agree the first option is the most useful.

• Add *_bucket_spacing variables in Scanner and implement usage.

  • This is the simplest solution and uses the current methodology

Yes, this is rather easy. In the first instance, I would suggest we do this initially with "gaps between blocks" = "gaps between buckets". So, I suggest to have only a get_bucket_spacing() that computes it as num_axial_blocks_per_bucket*get_axial_block_spacing(). i.e this is only an additional the get. In the next step, making the gap larger needs extra variables/parsing etc, so it's a bit more work. (not much in implementation, probably loads of work in testing).

For my own scanner development I need "gaps between blocks" != "gaps between buckets" but I will make this a separate PR. I suggest we hold off implementing "incorrectly" functioning get_*_bucket_spacing() methods until then.

I therefore I have two options.

1. Update the build_crystal_maps method to calculate the pitch of multiple axial buckets with the assumption that all the axial blocks have the same pitch/spacing, regardless of bucket index. I could introduce the get_*_bucket_spacing() methods here but their definitions would be "incorrect".
2. Disable the use of multiple axial buckets with this geometry via an error. num_axial_buckets > 1 does not work due to the current state of build_crystal_maps and computing axial_bucket_spacing = num_axial_blocks_per_bucket*get_axial_block_spacing() is assumptive about the geometry.

Really this is a case of trying to distinguish a line between two PRs, one a bug fix and the other an added feature. I would prefer to proceed with 2. and bug fix the current issue and then create a new PR with bucket_spacing interfile variables and more tests. Note, I would also have the disable the test added in this PR because it would fail.

[robbietuk]

Some comments on the latest changes. Running CI to ensure everything still works with GeometryBlocksOnCylindrical constructor changes

[KrisThielemans]

Really this is a case of trying to distinguish a line between two PRs, one a bug fix and the other an added feature. I would prefer to proceed with 2. and bug fix the current issue and then create a new PR with bucket_spacing interfile variables and more tests.

This is of course up to you. I guess in the second PR, you can still do it in stages (i.e. first commit(s), assume bucket_gaps==block_gap and see if tests work).

Note, I would also have the disable the test added in this PR because it would fail.

Which one? It seems a bit to introduce test code now that will fail soon. Or maybe you mean "disable the check on num_axial_blocks_per_bucket"? That'd of course be fine.

[KrisThielemans]

I'm a bit confused by the clang-format situation. Did you run it through? I guess so, as the pre-commit action is fine, but I guess I'm mildly surprised at the output.

[robbietuk]

Which one? It seems a bit to introduce test code now that will fail soon. Or maybe you mean "disable the check on num_axial_blocks_per_bucket"? That'd of course be fine.

Ignore this, I modified the tests to use num_axial_buckets=1 and allowed the tests (all pass).

I'm a bit confused by the clang-format situation. Did you run it through? I guess so, as the pre-commit action is fine, but I guess I'm mildly surprised at the output.

I am running clang-format on save. I don't see where the problem is...

[KrisThielemans]

I am running clang-format on save.

great

I don't see where the problem is...

no problem. I'm just surprised by some of the formatting choices it makes, but anyway.

[robbietuk]

@KrisThielemans b758eaf CI passed. Everything since is documentation or simple changes. I believe this is ready.

[KrisThielemans]

@markus-jehl can you have a look please?

[markus-jehl]

Looks good to me, apart from the release note comment that seems to refer to something that isn't there (anymore?).

[robbietuk]

Looks good to me, apart from the release note comment that seems to refer to something that isn't there (anymore?).

The release notes comment refers to the inclusion of Scanner::check_consistency() in the constructor of GeometryBlocksOnCylindrical that may result in an error if an invalid geometry.
https://github.com/UCL/STIR/pull/1375/files#diff-cdce2d1c9cb44f5a4ab3ef4228bd84366641b8a6e8fb077b459d1e6e68bf27aaR48-R49

[KrisThielemans]

thanks! Squash-merge ok?

[robbietuk]

Yes. Remember there is another bug fix in there with else if (get_scanner_geometry() == "Generic")

[KrisThielemans]

ok. well, you could update the release notes, and/or clean-up the commits, or just let me know if you prefer squash-merge or normal-merge...

[robbietuk]

I combine the past five commits into one. They were mostly documentation changes. I suggest a regular merge.

[KrisThielemans]

Thanks! Will merge after #1386 as #1383 broke the gcc12 job as it runs out of disk space (first noted'/fixed on #1236, but that PR will have to wait till 6.2)