Repository: scikit-image/scikit-image

Issue:
rescale and other warp functions should preserve input range and dtype
## Description

See this StackOverflow question for the motivation:

https://stackoverflow.com/questions/66585291/skimage-rescale-segmentation-mask

For higher-order interpolation, we need to do floating point arithmetic so it makes sense to convert to float and change the intensity scale as per our data type rules. For order=0, though, I think it makes most sense to preserve the input dtype and range.

```python
In [1]: import numpy as np
In [2]: from skimage import transform
In [3]: seg_mask = np.random.randint(0, 2**16, size=(16, 16)).astype(np.uint16)
In [4]: transform.rescale(seg_mask, 0.5, order=0)
Out[4]:
array([[0.5222858 , 0.55330739, 0.76211185, 0.61666285, 0.45897612,
        0.19795529, 0.52283513, 0.56203555],
       [0.57460899, 0.81582361, 0.61560998, 0.68281071, 0.3169604 ,
        0.19830625, 0.40459297, 0.48615244],
       [0.47759213, 0.62240024, 0.76609445, 0.333196  , 0.88154421,
        0.43877317, 0.50019074, 0.65618372],
       [0.84681468, 0.58448157, 0.56137942, 0.2824445 , 0.76746777,
        0.75156786, 0.47310597, 0.66036469],
       [0.18850996, 0.60334173, 0.25194171, 0.83747616, 0.56156252,
        0.61664759, 0.80070192, 0.48683909],
       [0.29718471, 0.3195697 , 0.32796216, 0.58196384, 0.81673915,
        0.50057221, 0.48458076, 0.27048142],
       [0.52285039, 0.35194934, 0.84243534, 0.68207828, 0.66150912,
        0.46347753, 0.23147936, 0.84484627],
       [0.51114672, 0.35536736, 0.87023728, 0.44734874, 0.6835584 ,
        0.54543374, 0.43472953, 0.78928817]])

In [5]: transform.rescale(seg_mask, 0.5, order=0, preserve_range=True)
Out[5]:
array([[34228., 36261., 49945., 40413., 30079., 12973., 34264., 36833.],
       [37657., 53465., 40344., 44748., 20772., 12996., 26515., 31860.],
       [31299., 40789., 50206., 21836., 57772., 28755., 32780., 43003.],
       [55496., 38304., 36790., 18510., 50296., 49254., 31005., 43277.],
       [12354., 39540., 16511., 54884., 36802., 40412., 52474., 31905.],
       [19476., 20943., 21493., 38139., 53525., 32805., 31757., 17726.],
       [34265., 23065., 55209., 44700., 43352., 30374., 15170., 55367.],
       [33498., 23289., 57031., 29317., 44797., 35745., 28490., 51726.]])
```

## Way to reproduce
```python
import numpy as np
from skimage import transform
seg_mask = np.random.randint(0, 2**16, size=(16, 16)).astype(np.uint16)
print(transform.rescale(seg_mask, 0.5, order=0))
print(transform.rescale(seg_mask, 0.5, order=0, preserve_range=True))
```


Discussion:

Patch:
diff --git a/pyproject.toml b/pyproject.toml
index 4dee5c765..4b438079d 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -3,11 +3,28 @@ requires = [
     "wheel",
     "setuptools<=51.0.0",
     "Cython>=0.29.18",
+
     # We follow scipy for much of these pinnings
     # https://github.com/scipy/scipy/blob/master/pyproject.toml
-    "numpy==1.16.5; python_version=='3.7'",
-    "numpy==1.17.3; python_version=='3.8'",
-    "numpy==1.19.3; python_version=='3.9'",
-    # do not pin numpy on future versions of python to avoid incompatible numpy and python versions
+    # numpy 1.19 was the first minor release to provide aarch64 wheels, but
+    # wheels require fixes contained in numpy 1.19.2
+    "numpy==1.19.2; python_version=='3.6' and platform_machine=='aarch64'",
+    "numpy==1.19.2; python_version=='3.7' and platform_machine=='aarch64'",
+
+    # default numpy requirements
+    "numpy==1.16.5; python_version=='3.6' and platform_machine!='aarch64' and platform_python_implementation != 'PyPy'",
+    "numpy==1.16.5; python_version=='3.7' and platform_machine!='aarch64' and platform_python_implementation != 'PyPy'",
+    "numpy==1.17.3; python_version=='3.8' and platform_machine!='aarch64' and platform_python_implementation != 'PyPy'",
+    "numpy==1.19.3; python_version=='3.9' and platform_python_implementation != 'PyPy'",
+
+    # First PyPy versions for which there are numpy wheels
+    "numpy==1.19.0; python_version=='3.6' and platform_python_implementation=='PyPy'",
+    "numpy==1.20.0; python_version=='3.7' and platform_python_implementation=='PyPy'",
+
+    # For Python versions which aren't yet officially supported,
+    # we specify an unpinned NumPy which allows source distributions
+    # to be used and allows wheels to be used as soon as they
+    # become available.
     "numpy; python_version>='3.10'",
+    "numpy; python_version>='3.8' and platform_python_implementation=='PyPy'",
 ]
diff --git a/skimage/filters/_unsharp_mask.py b/skimage/filters/_unsharp_mask.py
index 129cd09f1..92ecaa39c 100644
--- a/skimage/filters/_unsharp_mask.py
+++ b/skimage/filters/_unsharp_mask.py
@@ -1,5 +1,5 @@
 import numpy as np
-from scipy.ndimage.filters import gaussian_filter
+from scipy.ndimage import gaussian_filter
 from skimage import img_as_float
 
 
diff --git a/skimage/measure/_regionprops.py b/skimage/measure/_regionprops.py
index f1ba93335..fa79c8a2a 100644
--- a/skimage/measure/_regionprops.py
+++ b/skimage/measure/_regionprops.py
@@ -255,7 +255,14 @@ class RegionProperties:
             # determine whether func requires intensity image
             if n_args == 2:
                 if self._intensity_image is not None:
-                    return func(self.image, self.intensity_image)
+                    if self._multichannel:
+                        multichannel_list = [func(self.image,
+                                                  self.intensity_image[..., i])
+                                             for i in range(
+                            self.intensity_image.shape[-1])]
+                        return np.stack(multichannel_list, axis=-1)
+                    else:
+                        return func(self.image, self.intensity_image)
                 else:
                     raise AttributeError(
                         f"intensity image required to calculate {attr}"
diff --git a/skimage/morphology/selem.py b/skimage/morphology/selem.py
index e9f172e97..b4648d282 100644
--- a/skimage/morphology/selem.py
+++ b/skimage/morphology/selem.py
@@ -363,4 +363,4 @@ def _default_selem(ndim):
         are 1 and 0 otherwise.
 
     """
-    return ndi.morphology.generate_binary_structure(ndim, 1)
+    return ndi.generate_binary_structure(ndim, 1)
diff --git a/skimage/registration/_optical_flow.py b/skimage/registration/_optical_flow.py
index b16c15e95..1d465b89d 100644
--- a/skimage/registration/_optical_flow.py
+++ b/skimage/registration/_optical_flow.py
@@ -275,7 +275,7 @@ def _ilk(reference_image, moving_image, flow0, radius, num_warp, gaussian,
 
     for _ in range(num_warp):
         if prefilter:
-            flow = ndi.filters.median_filter(flow, (1, ) + ndim * (3, ))
+            flow = ndi.median_filter(flow, (1, ) + ndim * (3, ))
 
         moving_image_warp = warp(moving_image, get_warp_points(grid, flow),
                                  mode='edge')
diff --git a/skimage/restoration/inpaint.py b/skimage/restoration/inpaint.py
index 201f9c396..ddcce13f6 100644
--- a/skimage/restoration/inpaint.py
+++ b/skimage/restoration/inpaint.py
@@ -3,7 +3,7 @@ import numpy as np
 from scipy import sparse
 from scipy.sparse.linalg import spsolve
 import scipy.ndimage as ndi
-from scipy.ndimage.filters import laplace
+from scipy.ndimage import laplace
 import skimage
 from ..measure import label
 
@@ -127,8 +127,8 @@ def inpaint_biharmonic(image, mask, multichannel=False):
     mask = mask.astype(bool)
 
     # Split inpainting mask into independent regions
-    kernel = ndi.morphology.generate_binary_structure(mask.ndim, 1)
-    mask_dilated = ndi.morphology.binary_dilation(mask, structure=kernel)
+    kernel = ndi.generate_binary_structure(mask.ndim, 1)
+    mask_dilated = ndi.binary_dilation(mask, structure=kernel)
     mask_labeled, num_labels = label(mask_dilated, return_num=True)
     mask_labeled *= mask
 
diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py
index a9bac7404..02a6b4ae6 100644
--- a/skimage/transform/_warps.py
+++ b/skimage/transform/_warps.py
@@ -139,7 +139,8 @@ def resize(image, output_shape, order=None, mode='reflect', cval=0, clip=True,
         # The grid_mode kwarg was introduced in SciPy 1.6.0
         order = _validate_interpolation_order(image.dtype, order)
         zoom_factors = [1 / f for f in factors]
-        image = convert_to_float(image, preserve_range)
+        if order > 0:
+            image = convert_to_float(image, preserve_range)
         out = ndi.zoom(image, zoom_factors, order=order, mode=ndi_mode,
                        cval=cval, grid_mode=True)
         _clip_warp_output(image, out, order, mode, cval, clip)
@@ -827,7 +828,8 @@ def warp(image, inverse_map, map_args={}, output_shape=None, order=None,
 
     order = _validate_interpolation_order(image.dtype, order)
 
-    image = convert_to_float(image, preserve_range)
+    if order > 0:
+        image = convert_to_float(image, preserve_range)
 
     input_shape = np.array(image.shape)
 
@@ -849,7 +851,7 @@ def warp(image, inverse_map, map_args={}, output_shape=None, order=None,
              "skimage's implementation is fixed, we recommend "
              "to use bi-linear or bi-cubic interpolation instead.")
 
-    if order in (0, 1, 3) and not map_args:
+    if order in (1, 3) and not map_args:
         # use fast Cython version for specific interpolation orders and input
 
         matrix = None


Test patch:
diff --git a/skimage/measure/tests/test_regionprops.py b/skimage/measure/tests/test_regionprops.py
index e83b614dc..83cba272a 100644
--- a/skimage/measure/tests/test_regionprops.py
+++ b/skimage/measure/tests/test_regionprops.py
@@ -722,9 +722,16 @@ def test_multichannel():
     labels = slic(astro.astype(float), start_label=1)
 
     segment_idx = np.max(labels) // 2
-    region = regionprops(labels, astro_green)[segment_idx]
-    region_multi = regionprops(labels, astro)[segment_idx]
-    for prop in PROPS:
+    region = regionprops(labels,
+                         astro_green,
+                         extra_properties=[median_intensity]
+                         )[segment_idx]
+    region_multi = regionprops(labels,
+                               astro,
+                               extra_properties=[median_intensity]
+                               )[segment_idx]
+
+    for prop in list(PROPS.keys()) + ["median_intensity"]:
         p = region[prop]
         p_multi = region_multi[prop]
         if np.shape(p) == np.shape(p_multi):
diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py
index c044910c5..2d3dc64a5 100644
--- a/skimage/transform/tests/test_warps.py
+++ b/skimage/transform/tests/test_warps.py
@@ -1,8 +1,9 @@
+import pytest
 import numpy as np
 from scipy.ndimage import map_coordinates
 
 from skimage.data import checkerboard, astronaut
-from skimage.util.dtype import img_as_float
+from skimage.util.dtype import img_as_float, _convert
 from skimage.color.colorconv import rgb2gray
 from skimage.draw.draw import circle_perimeter_aa
 from skimage.feature.peak import peak_local_max
@@ -160,7 +161,8 @@ def test_rotate_resize_center():
     ref_x45[6, 0] = 1
     ref_x45[7, 0] = 1
 
-    x45 = rotate(x, 45, resize=True, center=(3, 3), order=0)
+    x45 = rotate(x, 45, resize=True, center=(3, 3), order=0,
+                 mode='reflect')
     # new dimension should be d = sqrt(2 * (10/2)^2)
     assert x45.shape == (14, 14)
     assert_equal(x45, ref_x45)
@@ -349,8 +351,8 @@ def test_resize_dtype():
     assert resize(x, (10, 10), preserve_range=True).dtype == x.dtype
     assert resize(x_u8, (10, 10), preserve_range=False).dtype == np.double
     assert resize(x_u8, (10, 10), preserve_range=True).dtype == np.double
-    assert resize(x_b, (10, 10), preserve_range=False).dtype == np.double
-    assert resize(x_b, (10, 10), preserve_range=True).dtype == np.double
+    assert resize(x_b, (10, 10), preserve_range=False).dtype == bool
+    assert resize(x_b, (10, 10), preserve_range=True).dtype == bool
     assert resize(x_f32, (10, 10), preserve_range=False).dtype == x_f32.dtype
     assert resize(x_f32, (10, 10), preserve_range=True).dtype == x_f32.dtype
 
@@ -535,7 +537,7 @@ def test_keep_range():
                   mode='constant', multichannel=False, anti_aliasing=False,
                   clip=True, order=0)
     assert out.min() == 0
-    assert out.max() == 2 / 255.0
+    assert out.max() == 2
 
 
 def test_zero_image_size():
@@ -684,3 +686,15 @@ def test_boll_array_warnings():
 
     with expected_warnings(['Input image dtype is bool']):
         warp(img, np.eye(3), order=1)
+
+
+@pytest.mark.parametrize('dtype', [np.uint8, bool, np.float32, np.float64])
+def test_order_0_warp_dtype(dtype):
+
+    img = _convert(astronaut()[:10, :10, 0], dtype)
+
+    assert resize(img, (12, 12), order=0).dtype == dtype
+    assert rescale(img, 0.5, order=0).dtype == dtype
+    assert rotate(img, 45, order=0).dtype == dtype
+    assert warp_polar(img, order=0).dtype == dtype
+    assert swirl(img, order=0).dtype == dtype
diff --git a/skimage/viewer/tests/test_utils.py b/skimage/viewer/tests/test_utils.py
index a6f5d63e7..878c209d7 100644
--- a/skimage/viewer/tests/test_utils.py
+++ b/skimage/viewer/tests/test_utils.py
@@ -20,7 +20,6 @@ def test_format_filename():
     assert fname is None
 
 
-@testing.skipif(True, reason="Can't automatically close window. See #3081.")
 @testing.skipif(not has_qt, reason="Qt not installed")
 def test_open_file_dialog():
     QApp = utils.init_qtapp()
@@ -30,7 +29,6 @@ def test_open_file_dialog():
     assert filename is None
 
 
-@testing.skipif(True, reason="Can't automatically close window. See #3081.")
 @testing.skipif(not has_qt, reason="Qt not installed")
 def test_save_file_dialog():
     QApp = utils.init_qtapp()
diff --git a/skimage/viewer/tests/test_widgets.py b/skimage/viewer/tests/test_widgets.py
index 962c89cc1..d9bedc73f 100644
--- a/skimage/viewer/tests/test_widgets.py
+++ b/skimage/viewer/tests/test_widgets.py
@@ -87,7 +87,6 @@ def test_slider_float():
     assert_almost_equal(sld.val, 2.5, 2)
 
 
-@testing.skipif(True, reason="Can't automatically close window. See #3081.")
 @testing.skipif(not has_qt, reason="Qt not installed")
 def test_save_buttons():
     viewer = get_image_viewer()