Output of rescale different from 0.18.3 to 0.19.0

[dstansby]

Description

The output of skimage.transform.rescale has changed from version 0.18.3 to 0.19. The changes all appear to be around the edge of the output - see below. As far as I can tell this change isn't documented anywhere in the release notes (https://scikit-image.org/docs/stable/release_notes.html), so I was wondering if this was either a bug, an undocumented deliberate change, or if I've missed something in the release notes?

Way to reproduce

import skimage.data as images
from skimage import transform as tf
import skimage
import numpy as np

# Rescale and save a copy
original = images.camera().astype('float')
newim = tf.rescale(original / original.max(), 0.25, order=4,
                   mode='constant', anti_aliasing=False) * original.max()
np.savetxt(f'new_{skimage.__version__}.txt', newim)

# Load both copies and plot the difference
import matplotlib.pyplot as plt
new_18 = np.loadtxt('new_0.18.3.txt')
new_19 = np.loadtxt('new_0.19.0.txt')

plt.imshow(new_19 - new_18, cmap='RdBu', vmin=-15, vmax=15)
plt.colorbar(label='0.19 - 0.18')
plt.show()

Version information

# Paste the output of the following python commands
from __future__ import print_function
import sys; print(sys.version)
import platform; print(platform.platform())
import skimage; print(f'scikit-image version: {skimage.__version__}')
import numpy; print(f'numpy version: {numpy.__version__}')

3.9.7 (default, Sep 16 2021, 08:50:36) 
[Clang 10.0.0 ]
macOS-10.16-x86_64-i386-64bit
scikit-image version: 0.19.0
numpy version: 1.21.2

[grlee77]

I assume you also have SciPy >= 1.6.0? If so, then the following is likely the source of this change:

scikit-image/skimage/_shared/utils.py

Lines 699 to 705 in 75f532b

	def _fix_ndimage_mode(mode):
	# SciPy 1.6.0 introduced grid variants of constant and wrap which
	# have less surprising behavior for images. Use these when available
	grid_modes = {'constant': 'grid-constant', 'wrap': 'grid-wrap'}
	if NumpyVersion(scipy.__version__) >= '1.6.0':
	mode = grid_modes.get(mode, mode)
	return mode

Note, that the user-specified "constant" gets converted to SciPy's "grid-constant" and "wrap" gets converted to "grid-wrap". In the case of wrap, the SciPy code is not what one would expect, so I would consider "grid-wrap" a bug fix rather than unintentional change.

For "constant" vs. "grid-constant" it is more a matter of convention than a true bug fix so we can debate if this change should have been made. The difference between the two is shown in the top two rows of this figure from the SciPy docs:
https://docs.scipy.org/doc/scipy/reference/_images/plot_boundary_modes.png

[rfezzani]

For "constant" vs. "grid-constant" it is more a matter of convention than a true bug fix so we can debate if this change should have been made. The difference between the two is shown in the top two rows of this figure from the SciPy docs:
https://docs.scipy.org/doc/scipy/reference/_images/plot_boundary_modes.png

This decision introduced a silent breaking change with no easy fix for our users (please see sunpy/sunpy#5750). Maybe should we rollback this decision and simply propose grid-constant as a new mode???

[grlee77]

Yeah, that seems fine.

[grlee77]

Actually, on closer look, I am not sure we should change this on the v0.19.x branch. If I take the SciPy 1.6-based code path vs. commenting that out so that the alternative, legacy 2D code path is taken, the images are the same (to within ~1e-11 floating point precision). So the new scipy.ndimage 'grid-constant' mode is the one that is consistent with our own 2D-only Cython-based code as called in the other code path. If I change grid-constant back to constant then I get inconsistent values at the edges between these two code paths.

That said, I could reproduce the difference from 0.18.3 to 0.19.0 as reported here, although I am not yet sure of the cause. The "legacy" code path linked above looks the same as the v0.18.3 code.

[grlee77]

Okay, I understand now. Our own cython-based interpolation code is only called within warp for order 1 or 3. For order=4 as in the example here ndimage.map_coordinates is used. Thus the SciPy 1.6 code path with grid-constant is inconsistent with legacy 'constant'.

However, we should not change 'grid-constant' back to 'constant' in 0.19.x or we would make two other issues:
1.) order=1 results would become inconsistent across 0.18.3 and 0.19.0
2.) within 0.19.0 results would become inconsistent between SciPy >= 1.6 and SciPy < 1.6 code paths

given that, I will close this issue as I consider the current change to be a bug fix