warp/rotate: fixed a bug with clipping when cval is not in the input range

[ayshih]

Description

Fixes #6305

For warp/rotate, when specifying clip=True (which is the default), the behavior when cval is outside of the range of the input image gives slightly wrong results for orders 1 and 3 (which use native code) and very wrong results for orders 2, 4, and 5 (which fall back to SciPy code). This is due to the use of a strict equality check.

Example code:

import matplotlib.pyplot as plt
import numpy as np
from skimage.transform import rotate

image = np.ones((20, 20))
image[5:-5, 5:-5] = 2

fig, axs = plt.subplots(2, 6, figsize=(13, 4))

for i in range(6):
    axs[0, i].set_title(f'order={i}')
    axs[0, i].imshow(rotate(image, 30, order=i, cval=0, resize=True, clip=False), vmin=0, vmax=3)
    axs[1, i].imshow(rotate(image, 30, order=i, cval=0, resize=True, clip=True), vmin=0, vmax=3)
    
axs[0, 0].set_ylabel('clip=False')
axs[1, 0].set_ylabel('clip=True')

plt.show()

Before this PR:

After this PR:

This PR also probably closes #5599 since the second block in that issue would no longer appear inconsistent:

>>> import numpy as np
>>> from skimage.transform import AffineTransform, warp
>>> img = np.ones((5, 5))
>>> tmat = np.array([[1, 0, -2.3], [0, 1, 0], [0, 0, 1]])
>>> tform = AffineTransform(matrix=tmat)
>>> warp(img, tform, cval=np.nan)
array([[nan, nan, nan,  1.,  1.],
       [nan, nan, nan,  1.,  1.],
       [nan, nan, nan,  1.,  1.],
       [nan, nan, nan,  1.,  1.],
       [nan, nan, nan,  1.,  1.]])
>>> warp(img, tform, cval=0)
array([[0., 0., 0.7, 1., 1.],
       [0., 0., 0.7, 1., 1.],
       [0., 0., 0.7, 1., 1.],
       [0., 0., 0.7, 1., 1.],
       [0., 0., 0.7, 1., 1.]])

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
• Descriptive commit messages (see below)

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.
• There is a bot to help automate backporting a PR to an older branch. For
  example, to backport to v0.19.x after merging, add the following in a PR
  comment: @meeseeksdev backport to v0.19.x
• To run benchmarks on a PR, add the run-benchmark label. To rerun, the label
  can be removed and then added again. The benchmark output can be checked in
  the "Actions" tab.

[pep8speaks]

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

• In the file skimage/transform/_warps.py:

Line 735:80: E501 line too long (86 > 79 characters)

• In the file skimage/transform/tests/test_warps.py:

Line 174:41: E226 missing whitespace around arithmetic operator
Line 174:45: E226 missing whitespace around arithmetic operator

Comment last updated at 2022-04-16 02:21:40 UTC

[grlee77]

Thanks, @ayshih!

Can you add a test case that verifies this fix to the clipping bug to skimage/transform/tests/test_warps.py?

Please also see the suggestion on avoiding NaN-related overhead for the common case where no NaN values are present.

[ayshih]

Unit tests added

[lagru]

Thank you @ayshih for tackling this. I don't understand this code fully yet and might look deeper when I have more time but see below...

[lagru]

I'm a little bit confused by the (new) logic in this function. So I tested it locally and noticed that commenting these 3 lines has no effect on testing test_warps.py at all! So either the test are not covering this function fully or the logic does not work the way you expect it to.

Why is it relevant to preserve cval if it's within min_func(output_image) <= cval <= max_func(output_image)?

If not min_val <= cval <= max_val, then the calls to min() / max() within np.clip(...) ensure that it is preserved. So I think that is why this line has no effect.

For not np.isnan(cval) to matter, the test suite should use cval=np.nan.

[ayshih]

Sorry, yes, there is incomplete test coverage, so I can add more tests.

Why is it relevant to preserve cval if it's within min_func(output_image) <= cval <= max_func(output_image)?

My thinking here is that even if cval is not in the input range, we don't want to expand the effective input range to include cval unless we know that cval was actually used in the output image. That is, if cval is not in the output range, cval wasn't actually used, so we proceed with clipping to the input range of min_val to max_val.

Note that it's not possible to test this logic using rotate() because that invariably uses cval for the output image for non-trivial rotations, while trivial rotations would not result in any pixel values that would need to be clipped. The test may be able to get away with using resize(..., mode='constant'), but may require using warp().

[ayshih]

Tests added

It turns out that the explicit not np.isnan(cval) is redundant because the min_func(output_image) <= cval <= max_func(output_image) can only be true if cval is not NaN.

Surprisingly, it turns out that even if you remove both checks, the rest of the new code as written happens to still pass the test for when cval is NaN due to a later quirk. cval will be NaN for the clip() call, but that doesn't break the clipping because min(min_val, np.nan) evaluates to min_val. I say "surprisingly" because not only does np.min([min_val, np.nan]) evaluate to np.nan, min(np.nan, min_val) (i.e., swapping the order of the arguments) evaluates to np.nan as well. I wasn't aware of that consequence of comparison operators for np.nan always evaluating to False: min(min_val, np.nan) != min(np.nan, min_val).

[lagru]

Just to make sure that I understand the way clip and cval interact together: if cval is outside the range of the input image and clip=True, what takes precedence? I think the old behavior tried to give cval precedence, although I think that output_image == cval never worked as intended. When I inspect the output image, the "outside area" seems to be close to but not precisely cval... So I wonder if cval_mask ever evaluated to true?

[lagru]

Ah, reading the original issue might help with that question. 😅 So cval should take precedence.

[ayshih]

At last, all tests green

[grlee77]

Thanks for adding the tests. I made just one new suggestion here

[alexdesiqueira]

Thank you for taking care of this, @ayshih!
@grlee77 thanks for your reviews, too 🙂 How are we on this one?
Thanks again!

[grlee77]

Thanks @ayshih, it looks good now. sorry for taking so long to get back to this

[grlee77]

@meeseeksdev backport to v0.19.x