issues with (P)OSMAPOSL sensitivity division

[KrisThielemans]

spotted by @AnderBiguri

We use the divide function to divide by sensitivity

STIR/src/iterative/OSMAPOSL/OSMAPOSLReconstruction.cxx

Lines 462 to 465 in 4054514

	divide(multiplicative_update_image_ptr->begin_all(),
	multiplicative_update_image_ptr->end_all(),
	sensitivity.begin_all(),
	small_num);

with small_num initialised here

STIR/src/iterative/OSMAPOSL/OSMAPOSLReconstruction.cxx

Line 436 in 4054514

static const float small_num = 0.000001F;

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.

[KrisThielemans]

We never see this in the logs due to #905

[KrisThielemans]

Of course, not so clear how to change the code (while still remaining generic), i.e. how do you find what is "small"?

[KrisThielemans]

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)

[KrisThielemans]

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.

[AnderBiguri]

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.

[KrisThielemans]

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)

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)