Expansion to Buckets on Cylindrical geometries

[robbietuk]

Primarily an extension of issue #1374 and PR #1375. The motivation for this was detailed here with an initial plan laid out. This issue serves as an extension of that comment.

STIR does not support multiple axial buckets of blocks of crystals on a ring. #1375 disabled the ability to use multiple axial buckets due to a bug in an assumptive implementation.

I require the use of gaps_between blocks != gaps_between_buckets therefore I am planning the implementation of BucketsOnCylindrical using the BlocksOnCylindrical current frame work. I do not plan to change the name, just the implementation.

Below are some of the main points worth discussion.

---

What do are a set of parallel buckets called?

I will refer to them as modules in this issue but the name can change via a refactor.

Modules are rotated by angle alpha around the z axis of the scanner. It is important a suitable name is given to the collection of buckets

Some other ideas:

• In GATE these are called rsectors for CylindricalPET (see GATE example .mac file)
  

  STIR/documentation/STIR-UsersGuide.tex

  Lines 779 to 797 in 37f799f

  	An OpenGATE scanner geometry is modular. For CylindricalPET, \textit{STIR} uses the conversion
  	scheme detailed in Table \ref{tab:STIR-OpenGATE_terminology}.
  	\begin{table}[h]
  	\centering
  	\begin{tabular}{ |c | c| }
  	\hline
  	\textbf{GATE} & \textbf{STIR} \\
  	\hline
  	Crystal & Crystal \\
  	Submodule & Block \\
  	Module & Bucket \\
  	Rsector & Buckets around the ring \\
  	\hline
  	\end{tabular}
  	\caption{A guide on the terminology equivalents between STIR and OpenGATE for CylindricalPET
  	scanner geometries.}
  	\label{tab:STIR-OpenGATE_terminology}
  	\end{table}

  
  STIR could adopt rsectors, sectors or sections
• Alternatively, these can be refereed to as a modules. This would might make the translation between STIR and GATE even more annoying but it is a good fit.
• @KrisThielemans I don't know exactly what other manufacturers refer to them as.

Suggestions welcome

---

Transaxial Buckets are maybe not supported?

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.

#1375 (comment). I am not sure this is the case upon further review.

STIR/src/buildblock/GeometryBlocksOnCylindrical.cxx

Lines 128 to 129 in 37f799f

	float alpha = scanner.get_intrinsic_azimuthal_tilt() + trans_bucket_num * (2 * _PI) / num_transaxial_buckets
	+ csi_minus_csiGaps;

It seems trans_bucket_num is only used in the computation of alpha and in the tangential_coord (similar to axial_bucket_num in calculation of axial_coord). The alpha is used to determine the rotation_matrix, i.e. alpha is an angle around the z axis.

I understand this as no two buckets can be parallel.

In this case, additional variables is needed throughout: num_axial_buckets_per_module and num_transaxial_buckets_per_module.

• These variables inevitably need to be everywhere but will default to 1.
• GeometryBlocksOnCylindrical::build_crystal_maps loops will need to be extended.
• More methods added to get_num_modules(), get_num_axial_crystals_per_module() and get_num_transaxial_crystals_per_module()

---

Spacing between buckets (Axial & Transaxial)

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.

#1375 (comment) . Not too much to add. The value should probably default to num_axial_blocks_per_bucket*get_axial_block_spacing(), otherwise a user defined variable > 0.

---

More tests

#1375 adds a test requiring the axial crystals position to increase monotonically with the axial crystal coordinate.

Additional tests should be added for the transaxial. As an idea, the angle sequential coordinate vectors in the detector map should never be greater than 90 degrees, i.e. cos(A -> B, B->C) >= 0. Except maybe in the instance of only two modules.

---

Virtual Crystals

#776 by @NikEfth tries to add virtual crystals. I am not 100% sure on the status of this and neither do I want to delve into it but I assume any issues encountered with adding virtual crystals to BucketsOnCylindrical would also be present in the BlocksOnCylindrical, so no significant change there.

---

Data Correction

Will there be any other issues with data correction methods etc. that we have not taken into account? e.g. scatter correction and normalization? I expect the same issues, maybe further complicated by this proposal.

@KrisThielemans @markus-jehl @danieldeidda

[robbietuk]

Notes for myself.

It is currently assumed that a scanner is defined from two parts. Consider Scanner::set_params() method. A "module", as dubbed above, is defined from the following transaxial variables:

STIR/src/buildblock/Scanner.cxx

Lines 1428 to 1431 in 79a1e83

	int num_axial_blocks_per_bucket_v,
	int num_transaxial_blocks_per_bucket_v,
	int num_axial_crystals_per_block_v,
	int num_transaxial_crystals_per_block_v,

The number of these "modules", incrementally rotated by angle alpha about the z axis, is derived from this transaxial bucket information and the num_detectors_per_ring_v parameter. I am not sure there is a check on this...

The axial information is derived from the above axial variables in the code quote and the num_rings_v.

To expand to BucketsOnCylindrical, "modules" need to be derived from the bucket per module details too. This will require an expansion in the scanner definition to add num_axial_buckets_per_module_v and num_transaxial_buckets_per_module_v. Checks should be run on all of these.

---

I should also mention that in #1375 (comment) I suggested replacing num_rings (and by extension num_detectors_per_ring):

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 method is similar to how GATE handles PET geometries and it doesn't require the user to pre-calculate num_rings and num_detectors_per_ring in the scanner definition. It would also remove the potential issues with getting the value wrong... However, that depreciates functionality and is a rather large change. I propose using the current tooling and then make the switch later if desired, maybe STIR v7?

[KrisThielemans]

What are a set of parallel buckets called?

Do you mean "What are a set of parallel blocks called"? Maybe not...

I'm still confused if you want to have a difference between "buckets" and "modules".

These days, PET scanners can have a whole hierarchy of things. I know of at least

• crystals (potentially with layers)
• mini-blocks
• blocks
• a collection of blocks (usually (always?) parallel) handled by the same electronics (relevant for dead-time)
• a collection of parallel blocks with equal spacing

Other notes:

• For all scanners I know, the last 2 are the same, but obviously there is no strict need for this.
• STIR doesn't encode mini-blocks, but they are not related to geometry (but could be to dead-time)
• STIR doc:

  • STIR/documentation/STIR-glossary.tex

    Lines 59 to 66 in 79a1e83

    	\item[Detector]
    	Sometimes called \textbf{detector crystal}. Geometrically associated
    	to the inner face of a detector element. The \textbf{scanner} is then
    	considered as a tessellation of \textbf{detectors} constructing adjacent \textbf{rings}.
    	For many scanners, detectors are organized in a block. For instance,
    	on the HR+ scanner, a detector block consists of 8x8 detectors. Current scanners
    	have mini-blocks, related to the read-out.
    	Blocks are often organised in \textbf{modules}, currently called \textbf{buckets} in STIR.

  • STIR/src/include/stir/Scanner.h

    Lines 66 to 85 in 79a1e83

    	At present, these functions follow CTI terminology, but the concepts
    	are similar for other scanners.
    	\li \c crystal the smallest detection unit
    	\li \c block several crystals are grouped in a block, this can be in
    	3 dimensions (see layer). This information could be useful for finding the
    	geometry of the scanner, but we would plenty more
    	size info presumably.
    	\li \c layer Some scanners have multiple layers of detectors to give
    	Depth Of Interaction information
    	\li \c bucket several \c blocks send detected events to one \c bucket.
    	This has consequences for the dead-time modelling. For instance,
    	one bucket could have much more singles than another, and hence
    	presumably higher singles-dead-time.
    	\li \c singles_unit (non-standard terminology)
    	Most scanners report the singles detected during the acquisition.
    	Some scanners (such as GE scanners) report singles for every crystal,
    	while others (such as CTI scanners) give only singles for a
    	collection of blocks. A \c singles_unit is then a set of crystals
    	for which we can get singles rates.

• Siemens uses "buckets" in their interfile headers, but they have been told to switched internally to "modules". GE uses "modules"

Given all this, I think that at present it would be enough to stick to 3 levels as we currently do, and the highest level can be called buckets or modules, but they mean the same thing in STIR. If this is not enough for your scanner, we need to do more thinking (and code modification). If you can, a description with enough detail of what you need would be necessary.

Transaxial Buckets are maybe not supported?

I understand this as no two buckets can be parallel.

That's correct for transaxial buckets (obviously, in axial direction, they still can be)

So, do you need more?

[robbietuk]

The geometry I am trying to model is the following :

• n x n crystals in a block
• m x m blocks in a bucket
• i axial buckets.
  where n>1, m>1 and i>1. Currently i is limited to 1 in STIR. This with i>1, the scanner would have a set of parallel buckets, which is what I have refereed to as a module. The reason I cannot restructure to put all axial blocks into the same axial bucket this is dead time and the gaps/pitch between buckets is not the same as the gaps between blocks.

Note, the extension I am proposing would also allow >1 transaxial buckets.

[robbietuk]

I understand this as no two buckets can be parallel.
That's correct for transaxial buckets (obviously, in axial direction, they still can be)
So, do you need more?

Perhaps not... I might have got ahead of myself with the num_axial_buckets_per_module_v and num_transaxial_buckets_per_module_v additional variable requirements.

[KrisThielemans]

So, I see no need for a module level then. We just say there are i x m x n rings, n crystals per axial block, m blocks per axial bucket. We need a "bucket-spacing` capability. done.

Note, the extension I am proposing would also allow >1 transaxial buckets.

now you've confused me again. I think the current code does what you need already on the transaxial front.

[robbietuk]

So, I see no need for a module level then. We just say there are i x m x n rings, n crystals per axial block, m blocks per axial bucket. We need a "bucket-spacing` capability. done.

Yup! I over complicated it for myself.

Note, the extension I am proposing would also allow >1 transaxial buckets.

now you've confused me again. I think the current code does what you need already on the transaxial front.

I disagree.

STIR/src/buildblock/GeometryBlocksOnCylindrical.cxx

Lines 128 to 131 in 37f799f

	float alpha = scanner.get_intrinsic_azimuthal_tilt() + trans_bucket_num * (2 * _PI) / num_transaxial_buckets
	+ csi_minus_csiGaps;
	stir::Array<2, float> rotation_matrix = get_rotation_matrix(alpha);

alpha is proportional to the rotation_matrix and increases with trans_bucket_num. Sequential transaxial buckets rotate.

[KrisThielemans]

alpha is proportional to the rotation_matrix and increases with trans_bucket_num. Sequential transaxial buckets rotate.

alpha increases with trans_bucket_num, so all blocks in one trans_bucket are parallel. But isn't that what you want? (alpha is obviously not proportional to the matrix, but it must have been late).

However, see #1390

[robbietuk]

Is there an equation that explains

STIR/src/buildblock/GeometryBlocksOnCylindrical.cxx

Lines 89 to 91 in 37f799f

	float start_z = -(axial_block_spacing * (num_axial_blocks_per_bucket)*num_axial_buckets - axial_crystal_spacing
	- ax_blocks_gap * (num_axial_blocks_per_bucket * num_axial_buckets - 1))
	/ 2;

I believe it finds the distance from the first and last axial crystals and then divides by 2 but the exact methodology here is hard to understand. This should be refactored to make this clear.

[robbietuk]

I am pretty sure this can be replaced by

  float crystals_in_block_offset = (scanner.get_num_axial_crystals_per_block() - 1) * scanner.get_axial_crystal_spacing();
  float blocks_in_bucket_offset = (scanner.get_num_axial_blocks_per_bucket() - 1) * scanner.get_axial_block_spacing();
  float bucket_offset = (scanner.get_num_axial_buckets() - 1) * scanner.get_axial_bucket_spacing();
  return -0.5 * (crystals_in_block_offset + blocks_in_bucket_offset + bucket_offset);

The idea here is if the crystals in a block are centered, and the blocks in a bucket are centered, then the buckets are centered in the z axis, the scanner should be centered. Current tests look good.

[robbietuk]

I have an issue with the use of non-descriptive variable names and code readability in the BlocksOnCylindrical project. I don't understand half of the shorthand variables in https://github.com/UCL/STIR/blob/master/examples/python/plot_scanner_LORs.py and trying to comprehend how positions of detectors are calculated is challenging and likely prone to errors, see #1390 and the cause

STIR/src/buildblock/GeometryBlocksOnCylindrical.cxx

Lines 81 to 98 in 37f799f

	// estimate the angle covered by half bucket, csi
	float csi = _PI / num_transaxial_buckets;
	float trans_blocks_gap = transaxial_block_spacing - num_transaxial_crystals_per_block * transaxial_crystal_spacing;
	float ax_blocks_gap = axial_block_spacing - num_axial_crystals_per_block * axial_crystal_spacing;
	float csi_minus_csiGaps = csi - (csi / transaxial_block_spacing * 2) * (transaxial_crystal_spacing / 2 + trans_blocks_gap);
	// distance between the center of the scannner and the first crystal in the bucket, r=Reffective/cos(csi)
	float r = scanner.get_effective_ring_radius() / cos(csi_minus_csiGaps);
	float start_z = -(axial_block_spacing * (num_axial_blocks_per_bucket)*num_axial_buckets - axial_crystal_spacing
	- ax_blocks_gap * (num_axial_blocks_per_bucket * num_axial_buckets - 1))
	/ 2;
	float start_y = -1 * scanner.get_effective_ring_radius();
	float start_x
	= -1 * r
	* sin(csi_minus_csiGaps); //(
	// ((num_transaxial_blocks_per_bucket-1)/2.)*transaxial_block_spacing
	// + ((num_transaxial_crystals_per_block-1)/2.)*transaxial_crystal_spacing
	// ); //the first crystal in the bucket

[markus-jehl]

I agree with Robbie that in transaxial direction having the module concept would be required to support different spacings between blocks and buckets. We would actually require this with our scanner, but chose to spread out the overall sum of block and bucket spacing evenly on the module. Then we rely on the geometric normalisation to correct for the resulting small modelling errors. It would obviously be nice to be able to model any sort of geometry without resorting to hand-made crystal maps, but it could complicate the code a lot and lead to potential additional bugs that would be very hard to find.