Relative Difference Prior - Computing gradient can return image with Nan values

[robbietuk]

I have observed an issue that objective function gradients, from images with voxel values of 0 or very small, having Nan values. Setting the penalty_strength to 0 resolves this issue. Therefore the issue likely lies with Relative Difference Prior (RDP). The issue does not occur with the Quadratic Prior. I have only observed this issue on a Linux system but could not replicate it on MacOS using the same data.

The gradient computation involves divisions. To prevent the division by zeros there is a check if both compared voxels and epsilon have zero value. This issue might be caused by very small voxel values (passing the >0 check) being squared, which may result in in the denominator being set to zero and hence the Nan values.

Relevant code: https://github.com/UCL/STIR/blob/621154460e7f1bf6e0b82a08418663baacba945f/src/recon_buildblock/RelativeDifferencePrior.cxx#L387-L398

A resolution to this is to use an epsilon > 0