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

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

Start with src/recon_buildblock/RelativeDifferencePrior.cxx around lines 387-398, focusing on the gradient divisions and the zero/epsilon check. Reproduce the issue with zero or very small voxel values on Linux, then verify that the Relative Difference Prior gradient contains no NaN values when epsilon is greater than zero.