STIR
STIR copied to clipboard
issues with (P)OSMAPOSL sensitivity division
spotted by @AnderBiguri
We use the divide
function to divide by sensitivity
https://github.com/UCL/STIR/blob/40545147179fbe2b1131679672375cc91a39ab55/src/iterative/OSMAPOSL/OSMAPOSLReconstruction.cxx#L462-L465
with small_num
initialised here
https://github.com/UCL/STIR/blob/40545147179fbe2b1131679672375cc91a39ab55/src/iterative/OSMAPOSL/OSMAPOSLReconstruction.cxx#L436
actual division code is in https://github.com/UCL/STIR/blob/master/src/include/stir/numerics/divide.inl
This finds a threshold based on the max in the estimate. This is fine for normal images, but is dangerous for the parametric case, as the images can have very different scale. Therefore, the threshold could be inappropriate for the one of the images, resulting incorrect zeroes in the image.
We never see this in the logs due to #905
Of course, not so clear how to change the code (while still remaining generic), i.e. how do you find what is "small"?
the divide
strategy seems in any case dangerous as it assumes that numerator and denominator are the same scale. (it should arguably use different small
values for both of them)
I think that for the sensitivity division, we could/should set small_num
to zero. (the main question here is if the backprojection updates a voxel a not. If not, image is undefined there, and we set it to 0.)
Less clear what to do in the OSL code.
In my hand-made version, indeed what I do to avoid this is to remove NaNs and Infs after the division, and nothing else. This makes the bug dissapear. That would work, not sure if its the ideal solution though.
the
divide
strategy seems in any case dangerous as it assumes that numerator and denominator are the same scale. (it should arguably use differentsmall
values for both of them)
The divide
strategy is used to divide the backproj(quotient)
by backproj(1)
(i.e. subset-sensitivity), so these should (after a while) have similar scales, so this is probably ok.
However, it sets the 0/0 to 0, but this is inappropriate for OSEM. It should set it to 1 (as we should not update voxel values which are not affected by the subset).
I still have to think what we need to do with OSL (which is a whole load of other difficulties and heuristics)