Bumped precision from float to double

[bbbxyz]

Calculations now performed in double precision

Performing the forwards and inverse FFT while storing values as floats is too inaccurate. (e.g. taking a volume, adding some noise on the order of 10^-6 and blurring the two volumes makes the difference jump up to 10^-4 in some cases).

[gdevenyi]

Do you have an example case to show how this fixes your described issue?

[bbbxyz]

You can try with those two volumes: blur_ex.zip

Comparing them with minccmp: ssq: 1.262429495e-07 rmse: 1.773632244e-06 maxdiff: 6.103515625e-05

After blurring them with a fwhm of 4 (in float): ssq: 0.0008500526486 rmse: 0.0001455401846 maxdiff: 0.0009765625

With calculations in double: ssq: 4.162689736e-06 rmse: 1.018467161e-05 maxdiff: 0.00048828125

Subtracting the volumes also helps to visualize how the noise gets amplified: Left column is the difference between the two original volumes, middle column is blurring in single precision and third column is in double precision.

[gdevenyi]

Interesting.

Often these types of error amplifications can be a result of the format of the underlying mathematical expression, where precision loss can cause such issues, moving to a double is often just a band-aid solution. I think perhaps the underlying implementation should be checked before just upgrading all the floating representations.

[bbbxyz]

I don't see anything clearly wrong with the FFT or vector multiplication implementations.
I guess another option would be to use FFTW (see https://github.com/BIC-MNI/mincfft).

Thoughts?

[vfonov]

https://github.com/BIC-MNI/EZminc/blob/ITK4/tools/fast_blur.cpp uses https://github.com/BIC-MNI/EZminc/blob/ITK4/image_proc/fftw_blur.cpp

[vfonov]

fast_blur --fwhm 4.0 volume1.mnc volume1_blur.mnc
fast_blur --fwhm 4.0 volume2.mnc volume2_blur.mnc
minccmp volume1_blur.mnc volume2_blur.mnc 

ssq: 4.824499886e-05 rmse: 3.467258004e-05 zscore: 7.913100942e-08

Calculations are done with float

[bbbxyz]

Thanks! That seems like a significant improvement compared to mincblur. I'lll see if that resolves the issues I've been encountering.

[gdevenyi]

@bbbxyz If this is something you're actually running into, there's also SmoothImage from the ANTs package which uses an ITK implementation (in float) and seems to work fine (very similar internals to fastblur.

[bbbxyz]

Closing this