RelativeDifferencePrior gives wrong results due to use of abs()

[KrisThielemans]

We use the abs() function, e.g.

STIR/src/recon_buildblock/RelativeDifferencePrior.cxx

Line 416 in 8ced2d7

* abs(current_image_estimate[z][y][x] - current_image_estimate[z + dz][y + dy][x + dx])

It turns out that with g++-8 (and presumably any version), this falls back to the C version of abs, which will convert its argument to an integer type.

To fix this, we have to explicitly use std::abs (which has an overload for float and double in <cmath>).

Impact` of this is rather dramatic. When running test_priors.cxx, I get the following for the gradient of the RDP potential:

RDP old 0.590537 RDP new 0.287209
RDP old 0.893955 RDP new 0.326926
RDP old 0.319332 RDP new 0.20883
RDP old 0.00598219 RDP new 0.00592888
RDP old 0.373463 RDP new 0.229309
RDP old 0.980504 RDP new 0.332455
RDP old -0.00389029 RDP new -0.00386772
RDP old 0.348276 RDP new 0.220149
RDP old 0 RDP new 0
RDP old 0.153586 RDP new 0.123653
RDP old 0.893202 RDP new 0.326867
RDP old 0.493659 RDP new 0.265453
RDP old 0.773317 RDP new 0.315314
RDP old 0.18338 RDP new 0.141967
RDP old -0.159962 RDP new -0.129018
RDP old 0.731704 RDP new 0.31015
RDP old 0.973774 RDP new 0.332113

This is despite the tests all working, including a test on numerical vs analytic gradient which is surprising. Possibly this is because we use abs() consistently...

[KrisThielemans]

@Imraj-Singh, you faithfully reproduce this bug in #1408

[KrisThielemans]

Plots for derivative of the RDP

[KrisThielemans]

Looks like this was introduced in the PR for the RDP #335. @robbietuk, making you aware.

It is possible that this bug only transpired once we moved to C++-11, or certain compiler versions. I don't know.

[KrisThielemans]

Other places that use abs for floats/doubles that could be affected:

Library: (worrying!)

• ProjDataInfo::get_sampling_in*: are not affected, but I don't know why
• inverse_SSRB
• extend_projdata

Utilities:

• shift_image (maybe not due to #include <math.h>)
• find_fwhm_in_image.cxx (possibly not due to using namespace std;)

Relatively harmless as error checking:

• ExamInfo::operator==, - Radionuclide::operator==
• TimeFrameDefinitions::operator==
• error checking in ProjMatrixByBinSPECTUB
• error checking in ProjMatrixByBinUsingRayTracing

Test code:

• test_interpolate_projdata.cxx
• test_blocks_on_cylindrical_projectors.cxx

Commented out code, hence harmless (but confusing):

• warp_image