Possible problem with multiSiddon projector random sampling


I think there is an error in the GetRdmPositionsAndOrientations function of the file iScannerPET.cc. In short, the vectors used to create a random position in a detector are not orthogonal between each other and two of them are not guaranteed to have a norm of one.

The long version:
mp_crystalOrientation gives us a vector that goes along the depth of the current detector. However, CASToR does not keep a similar vector for the detector's transaxial and axial axes. For simplicity, let's define (mp_crystalOrientationX, mp_crystalOrientationY, mp_crystalOrientationZ) as (uX, uY, uZ). Currently, the transaxial axis of the detector is defined as (in the function GetRdmPositionsAndOrientations):
(uY, uX, 0)
while the axial axis of the detector is defined as:
(uX*uZ, uY*uZ, sqrt(1-uZ**2))
From there, we can observe several problems. Two examples:
1) (uX, uY, uZ) (uY, uX, 0)' = uX*uY + uY*uX =/= 0 (Most of the time.)
2) || (uY, uX, 0) || =/= 1 if uZ =/= 0 (Since ||(uX, uY, uZ)|| = 1.)
Thus, the set of positions generated by that function does not correspond exactly to the detector.

A possible correction: (it needs to be applied on the computation of ap_Position1 and ap_Position2)
Replace (uY, uX, 0) by (uY, -uX, 0) / || (uY, -uX, 0) ||
Replace (uX*uZ, uY*uZ, sqrt(1-uZ**2)) by (-uX*uZ, -uY*uZ, 1-uZ**2) / || (-uX*uZ, -uY*uZ, 1-uZ**2) ||

Note: Since mp_crystalOrientationZ is 0.0 for all detectors in most scanner, this error might have negligible repercussions in practice. I am not sure, since I didn't test its repercussions in reconstruction.

I have another question. In the comments prior to the function GetRdmPositionsAndOrientations, there is the following line: "fix the possibility to draw LOR outside the actual crystal volume (if mp_meanDepthOfInteraction != 0.5)". From what I understand, that function never uses that parameter. Was that comment meant for the function GetPositionsAndOrientations or did I miss something?

Have a nice day,
Maxime Toussaint

Hello Maxime,

Thank you for your feedback, the current implementation of the Siddon multi-lines projector class (and the associated GetRdmPositionsAndOrientations() function from the iScannerPET.cc file) has indeed some limitations. It is still work in progress and is only intended to be used with usual systems whose detectors have an axial orientation equal to 0. I should have made it clearer in the description of the projector.
It has not been assessed with systems using detectors with non-zero axial angle, and as you noticed it contains some errors in this case. Regarding the comment about mp_meanDepthOfInteraction, it refers to a previous implementation in the event a DOI was provided.

This projector is actually a draft and we plan to revamp the management of multi-line projector to make it more flexible. In short, the estimation of the random position in the detector should be computed in the mother class (or inside an additional intermediate class), from which any line projector class could be called to compute the system matrix elements. This would avoid to create a multi-lines projector class specific to each projection algorithm. Besides, we plan to rework on the orientations in order to make it more flexible for less generic systems. All of this will require some substantial work though.

Thank again for your remark and suggestions. I will have a look this summer to correct the current implementation.


Hi Thibaut,

It would be interesting to have the multi-lines approach for any projector. Good luck with the refactoring!

Just to make sure, the problem I talked about in my previous mail does not affect only the axial axis, but also the transaxial axis. Both of these vector are not guaranteed to be orthonormal with the detector orientation (depth axis). In its current form, it could generate position outside the detector. Thus, multiSiddon does not work as intended (except when uX*uY + uY*uX = 0).

While on the subject, are there plans to implement a deterministic multi-projector? For example, uniform sampling of the detector? If I understand correctly, the stochastic approach is advantageous when sample rate is low. However, the positions sampled are not kept, which result in a "different" system matrix at each iteration. It might not change much since PET detectors are small. Did your team, or any other groups you are aware of, studied the impact of deterministic vs stochastic multi-projector on the reconstruction?

Have a nice day,
Maxime Toussaint

Hi Maxime,

There is no current plan to implement new projector at the moment (at least not before the reorganization of the current projector classes). However such implementation should be feasible with the current version, as the geometric information of a detector (positions and orientations) could be recovered from the projectionLine structure, and the dimensions of the detector could be retrieved from the iScannerPET functions. From there, one could implement a multi-lines projector which reads an uniform distribution, or any kind of detection distribution model provided by a configuration file, and use the same detection location for each iteration.
There have been several works in our group regarding multi-lines projectors whose detection positions in the detectors are recovered from a probability distribution generated with Monte-Carlo simulations, but this is also a stochastic approach. I didn't read about comparisons between deterministic vs stochastic multi-projectors though.

Thanks for your advices !