improved single precision support in skimage.transform

[grlee77]

Description

The warp, pyramid and radon transform functions all already supported float32, but this PR expands the tests a bit and makes sure float16 gets promoted to float32 rather than failing.

For integral_image, I changed it to promote any floating-point types to double precision for better accuracy with numpy.cumsum. (As far as I know, cumsum does not use the higher accuracy pairwise summation that some NumPy sums use). We had seen issues with limited precision of integral images in the past with non-local means and had changed the Cython code internal to that function to always use doubles for the integral image computation.

Checklist

• Docstrings for all functions
• Gallery example in ./doc/examples (new features only)
• Benchmark in ./benchmarks, if your changes aren't covered by an
  existing benchmark
• Unit tests
• Clean style in the spirit of PEP8

For reviewers

• Check that the PR title is short, concise, and will make sense 1 year
  later.
• Check that new functions are imported in corresponding __init__.py.
• Check that new features, API changes, and deprecations are mentioned in
  doc/release/release_dev.rst.

[pep8speaks]

Hello @grlee77! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2021-06-07 15:15:14 UTC

[rfezzani]

Thank you @grlee77.

[rfezzani]

@grlee77, the type promotion to double in integral is opposed to the strategy of preserving input image dtype 😕

[grlee77]

@grlee77, the type promotion to double in integral is opposed to the strategy of preserving input image dtype

It currently only preserves float dtypes. Apparently np.cumsum already uses int64 for, e.g. uint8 inputs. Otherwise you would be pretty much guaranteed to get integer overflows! I am just arguing that for accuracy sake, we should be promoting float inputs to float64. We could still restore the final output to the input float type, though, as a compromise.

[grlee77]

Another simple example. Suppose the image is somewhat large, but has float16 dtype. In that case it is also likely to encounter floating point overflow (if code is as in current main branch where cumsums would be in float16):

integral_image(np.ones((2048, 2048), dtype=np.float16))

gives

array([[1.000e+00, 2.000e+00, 3.000e+00, ..., 2.046e+03, 2.047e+03,
        2.048e+03],
       [2.000e+00, 4.000e+00, 6.000e+00, ..., 4.092e+03, 4.094e+03,
        4.096e+03],
       [3.000e+00, 6.000e+00, 9.000e+00, ..., 7.500e+03, 7.504e+03,
        7.508e+03],
       ...,
       [2.046e+03, 4.092e+03, 6.136e+03, ...,       inf,       inf,
              inf],
       [2.047e+03, 4.094e+03, 6.140e+03, ...,       inf,       inf,
              inf],
       [2.048e+03, 4.096e+03, 6.144e+03, ...,       inf,       inf,
              inf]], dtype=float16)

[rfezzani]

I understand the overflow problem, we can then consider this is a bug fix. May be should we document this behaviour...

[grlee77]

I understand the overflow problem, we can then consider this is a bug fix. May be should we document this behaviour...

I reverted the change to default dtype in integral_image. We can revisit in a separate PR.

[mkcor]

Thanks, @grlee77! I wasn't sure why float16 wasn't systematically included in the floating-point types when testing with parameters, so please let me know for my personal understanding, but this PR is good to go!

[rfezzani]

Merging 🎉 Thank you @grlee77