From 47cb0bb970bff4a2e0592ec7421079f4e10f7162 Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Sun, 30 May 2021 23:50:10 +0200
Subject: [PATCH 01/24] Add dask_image.ndinterp.map_coordinates + some tests.
 Implementation: pre-map coordinates onto the chunks of 'image' and execute
 `ndimage.map_coordinates` for every chunk.

---
 dask_image/ndinterp/__init__.py               |  66 +++++++++++
 .../test_ndinterp/test_map_coordinates.py     | 106 ++++++++++++++++++
 2 files changed, 172 insertions(+)
 create mode 100644 tests/test_dask_image/test_ndinterp/test_map_coordinates.py

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 8d6d6837..f53b233e 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -11,6 +11,7 @@
 from dask.highlevelgraph import HighLevelGraph
 import scipy
 from scipy.ndimage import affine_transform as ndimage_affine_transform
+from scipy.ndimage import map_coordinates as ndimage_map_coordinates
 
 from ..dispatch._dispatch_ndinterp import (
     dispatch_affine_transform,
@@ -381,3 +382,68 @@ def spline_filter1d(
     )
 
     return result
+
+
+def map_coordinates(image, coordinates, order=3,
+                    mode='constant', cval=0.0, prefilter=False):
+    """
+    Pre-map coordinates onto the chunks of 'image' and execute
+    `ndimage.map_coordinates` for every chunk.
+    """
+
+    coordinates = np.asarray(coordinates)
+
+    coord_chunk_locations = [np.searchsorted(np.cumsum(image.chunks[dim]),
+                                             coordinates[dim])
+                             for dim in range(image.ndim)]
+
+    coord_chunk_locations = np.array(coord_chunk_locations).T
+
+    required_chunks, coord_rc_inds = np.unique(coord_chunk_locations,
+                                               axis=0, return_inverse=True)
+
+    rc_coord_inds = [[] for _ in range(len(required_chunks))]
+
+    rois = np.ones((2, len(required_chunks), image.ndim)) * np.nan
+    for ic, c in enumerate(coordinates.T):
+
+        rois[0][coord_rc_inds[ic]] = np.nanmin([rois[0][coord_rc_inds[ic]],
+                                                np.floor(c - order // 2)], 0)
+        rois[1][coord_rc_inds[ic]] = np.nanmax([rois[1][coord_rc_inds[ic]],
+                                                np.ceil(c + order // 2)], 0)
+
+        rc_coord_inds[coord_rc_inds[ic]] += [ic]
+
+    name = "map_coordinates-%s" % tokenize(image,
+                                           coordinates,
+                                           order,
+                                           mode,
+                                           cval,
+                                           prefilter)
+
+    keys = [(name, i) for i in range(len(required_chunks))]
+
+    values = []
+    for irc, rc in enumerate(required_chunks):
+        offset = np.min([np.max([[0] * image.ndim,
+                                 rois[0][irc].astype(np.int64)], 0),
+                         np.array(image.shape) - 1], 0)
+        sl = [slice(offset[dim],
+                    min(image.shape[dim], int(rois[1][irc][dim]) + 1))
+              for dim in range(image.ndim)]
+
+        values.append((ndimage_map_coordinates,
+                       image[tuple(sl)],
+                       coordinates[:, rc_coord_inds[irc]] - offset[:, None],
+                       None,
+                       order,
+                       mode,
+                       cval,
+                       prefilter))
+
+    dsk = dict(zip(keys, values))
+    ar = da.Array(dsk, name, tuple([[len(rc_ci) for rc_ci in rc_coord_inds]]), image.dtype)
+
+    orig_order = np.argsort([ic for rc_ci in rc_coord_inds for ic in rc_ci])
+
+    return ar[orig_order]
diff --git a/tests/test_dask_image/test_ndinterp/test_map_coordinates.py b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
new file mode 100644
index 00000000..29270728
--- /dev/null
+++ b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
@@ -0,0 +1,106 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+from packaging import version
+
+import dask.array as da
+import numpy as np
+import pytest
+import scipy
+import scipy.ndimage
+
+import dask_image.ndinterp
+
+# mode lists for the case with prefilter = False
+_supported_modes = ['constant', 'nearest']
+_unsupported_modes = ['wrap', 'reflect', 'mirror']
+
+# mode lists for the case with prefilter = True
+_supported_prefilter_modes = ['constant']
+_unsupported_prefilter_modes = _unsupported_modes + ['nearest']
+
+have_scipy16 = version.parse(scipy.__version__) >= version.parse('1.6.0')
+
+# additional modes are present in SciPy >= 1.6.0
+if have_scipy16:
+    _supported_modes += ['grid-constant']
+    _unsupported_modes += ['grid-mirror', 'grid-wrap']
+    _unsupported_prefilter_modes += ['grid-constant', 'grid-mirror',
+                                     'grid-wrap']
+
+
+def validate_map_coordinates_general(n=2,
+                                     interp_order=1,
+                                     interp_mode='constant',
+                                     coord_len=12,
+                                     coord_chunksize=6,
+                                     coord_offset=0.,
+                                     im_shape_per_dim=12,
+                                     im_chunksize_per_dim=6,
+                                     random_seed=0,
+                                     prefilter=False,
+                                     ):
+
+    if interp_order > 1 and interp_mode == 'nearest' and not have_scipy16:
+        # not clear on the underlying cause, but this fails on older SciPy
+        pytest.skip("requires SciPy >= 1.6.0")
+
+    # define test image
+    np.random.seed(random_seed)
+    image = np.random.random([im_shape_per_dim] * n)
+    image_da = da.from_array(image, chunks=im_chunksize_per_dim)
+
+    # define test coordinates
+    coords = np.random.random((n, coord_len)) * im_shape_per_dim + coord_offset
+
+    # ndimage result
+    ints_scipy = scipy.ndimage.map_coordinates(
+        image,
+        coords,
+        output=None,
+        order=interp_order,
+        mode=interp_mode,
+        cval=0.0,
+        prefilter=prefilter)
+
+    # dask-image result
+    ints_dask = dask_image.ndinterp.map_coordinates(
+        image_da,
+        coords,
+        order=interp_order,
+        mode=interp_mode,
+        cval=0.0,
+        prefilter=prefilter)
+
+    ints_dask_computed = ints_dask.compute()
+
+    assert np.allclose(ints_scipy, ints_dask_computed)
+
+
+@pytest.mark.parametrize("n",
+                         [1, 2, 3, 4])
+@pytest.mark.parametrize("random_seed",
+                         range(2))
+def test_map_coordinates_basic(n,
+                               random_seed,
+                               ):
+
+    kwargs = dict()
+    kwargs['n'] = n
+    kwargs['random_seed'] = random_seed
+
+    validate_map_coordinates_general(**kwargs)
+
+
+@pytest.mark.timeout(3)
+def test_map_coordinates_large_input():
+
+    # define large test image
+    image_da = da.random.random([1000] * 3, chunks=200)
+
+    # define sparse test coordinates
+    coords = np.random.random((3, 2)) * 1000
+
+    # dask-image result
+    dask_image.ndinterp.map_coordinates(
+        image_da,
+        coords).compute()

From 42be4a86c57e5a3acac35114bbadeaf60b93ff6d Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Mon, 31 May 2021 00:00:29 +0200
Subject: [PATCH 02/24] Tick function in coverage doc

---
 docs/coverage.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/coverage.rst b/docs/coverage.rst
index d849c485..1979d80d 100644
--- a/docs/coverage.rst
+++ b/docs/coverage.rst
@@ -145,7 +145,7 @@ This table shows which SciPy ndimage functions are supported by dask-image.
      - ✓
    * - ``map_coordinates``
      - ✓
-     -
+     - ✓
    * - ``maximum``
      - ✓
      - ✓

From 4c009766e414af5c7233d5cb521a63b99ad02e88 Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Sun, 9 Oct 2022 18:35:51 +0200
Subject: [PATCH 03/24] Added comments explaining the steps taken. Increased
 verbosity. Simplified code. Clarified tests.

---
 dask_image/ndinterp/__init__.py               | 77 ++++++++++++-------
 .../test_ndinterp/test_map_coordinates.py     |  5 ++
 2 files changed, 56 insertions(+), 26 deletions(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index f53b233e..686ab6b1 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -387,33 +387,56 @@ def spline_filter1d(
 def map_coordinates(image, coordinates, order=3,
                     mode='constant', cval=0.0, prefilter=False):
     """
-    Pre-map coordinates onto the chunks of 'image' and execute
-    `ndimage.map_coordinates` for every chunk.
+    Wraps ndimage.map_coordinates.
+
+    Covers the use case in which the coordinates are non-dask arrays
+    and the input image array consists of a (potentially large) dask array.
+
+    Implementation details and steps:
+
+    1) associate each coordinate in coordinates with the chunk
+       it maps to in the input image
+    2) for each input image chunk that has been associated to at least one
+       coordinate, calculate the minimal slice required to map
+       all coordinates that are associated to it (note that resulting slice
+       coordinates can lie outside of the coordinate's chunk)
+    3) for each previously obtained slice and its associated coordinates,
+       define a dask task and apply ndimage.map_coordinates
+    4) outputs of ndimage.map_coordinates are rearranged to match input order
     """
 
-    coordinates = np.asarray(coordinates)
-
-    coord_chunk_locations = [np.searchsorted(np.cumsum(image.chunks[dim]),
-                                             coordinates[dim])
-                             for dim in range(image.ndim)]
-
-    coord_chunk_locations = np.array(coord_chunk_locations).T
-
-    required_chunks, coord_rc_inds = np.unique(coord_chunk_locations,
-                                               axis=0, return_inverse=True)
-
-    rc_coord_inds = [[] for _ in range(len(required_chunks))]
-
-    rois = np.ones((2, len(required_chunks), image.ndim)) * np.nan
+    # STEP 1: Associate each coordinate in coordinates with
+    # the chunk it maps to in the input image
+
+    # determine which chunk each coordinate maps to
+    coord_chunk_locations = np.array([
+        np.searchsorted(np.cumsum(image.chunks[dim]), coordinates[dim])
+        for dim in range(image.ndim)]).T
+    # associate coordinates with the chunks they map into
+    required_input_chunks, coord_rc_inds, coord_rc_count = np.unique(
+        coord_chunk_locations, axis=0, return_inverse=True, return_counts=True)
+
+    # STEP 2: for each input image chunk that has been associated to at least
+    # one coordinate, calculate the minimal slice required to map all
+    # coordinates that are associated to it (note that resulting slice
+    # coordinates can lie outside of the coordinate's chunk)
+
+    input_slices = np.ones((2, len(required_input_chunks), image.ndim)) * np.nan
+    rc_coord_inds = [np.where(coord_rc_inds == irc)[0]
+                     for irc in range(len(required_input_chunks))]
     for ic, c in enumerate(coordinates.T):
 
-        rois[0][coord_rc_inds[ic]] = np.nanmin([rois[0][coord_rc_inds[ic]],
-                                                np.floor(c - order // 2)], 0)
-        rois[1][coord_rc_inds[ic]] = np.nanmax([rois[1][coord_rc_inds[ic]],
-                                                np.ceil(c + order // 2)], 0)
+        input_slices[0][coord_rc_inds[ic]]\
+            = np.nanmin([input_slices[0][coord_rc_inds[ic]],
+                         np.floor(c - order // 2)], 0)
+        input_slices[1][coord_rc_inds[ic]]\
+            = np.nanmax([input_slices[1][coord_rc_inds[ic]],
+                         np.ceil(c + order // 2)], 0)
 
-        rc_coord_inds[coord_rc_inds[ic]] += [ic]
+    # STEP 3: For each previously obtained slice and its associated
+    # coordinates, define a dask task and apply ndimage.map_coordinates
 
+    # prepare building of dask graph
     name = "map_coordinates-%s" % tokenize(image,
                                            coordinates,
                                            order,
@@ -421,15 +444,15 @@ def map_coordinates(image, coordinates, order=3,
                                            cval,
                                            prefilter)
 
-    keys = [(name, i) for i in range(len(required_chunks))]
+    keys = [(name, i) for i in range(len(required_input_chunks))]
 
     values = []
-    for irc, rc in enumerate(required_chunks):
+    for irc, _rc in enumerate(required_input_chunks):
         offset = np.min([np.max([[0] * image.ndim,
-                                 rois[0][irc].astype(np.int64)], 0),
+                                 input_slices[0][irc].astype(np.int64)], 0),
                          np.array(image.shape) - 1], 0)
         sl = [slice(offset[dim],
-                    min(image.shape[dim], int(rois[1][irc][dim]) + 1))
+                    min(image.shape[dim], int(input_slices[1][irc][dim]) + 1))
               for dim in range(image.ndim)]
 
         values.append((ndimage_map_coordinates,
@@ -442,8 +465,10 @@ def map_coordinates(image, coordinates, order=3,
                        prefilter))
 
     dsk = dict(zip(keys, values))
-    ar = da.Array(dsk, name, tuple([[len(rc_ci) for rc_ci in rc_coord_inds]]), image.dtype)
+    ar = da.Array(dsk, name, tuple([list(coord_rc_count)]), image.dtype)
 
+    # STEP 4: rearrange outputs of ndimage.map_coordinates
+    # to match input order
     orig_order = np.argsort([ic for rc_ci in rc_coord_inds for ic in rc_ci])
 
     return ar[orig_order]
diff --git a/tests/test_dask_image/test_ndinterp/test_map_coordinates.py b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
index 29270728..4c212102 100644
--- a/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
+++ b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
@@ -94,6 +94,11 @@ def test_map_coordinates_basic(n,
 @pytest.mark.timeout(3)
 def test_map_coordinates_large_input():
 
+    """
+    This test assesses whether relatively large
+    inputs are processed before timeout.
+    """
+
     # define large test image
     image_da = da.random.random([1000] * 3, chunks=200)
 

From 9b8b0574bef9416f74d07d80b22ed1a7c2f8e4a7 Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Sun, 8 Jan 2023 22:14:49 +0100
Subject: [PATCH 04/24] Extended `map_coordinates` to work also on chunked
 coordinate arrays, which requires two compute passes. Adapted tests.

---
 dask_image/ndinterp/__init__.py               | 83 +++++++++++++++++--
 .../test_ndinterp/test_map_coordinates.py     |  3 +-
 2 files changed, 77 insertions(+), 9 deletions(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 686ab6b1..b4ea3fc6 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -5,6 +5,7 @@
 from itertools import product
 import warnings
 
+from dask import compute, delayed
 import dask.array as da
 import numpy as np
 from dask.base import tokenize
@@ -384,16 +385,14 @@ def spline_filter1d(
     return result
 
 
-def map_coordinates(image, coordinates, order=3,
+def _map_single_coordinates_array_chunk(image, coordinates, order=3,
                     mode='constant', cval=0.0, prefilter=False):
     """
-    Wraps ndimage.map_coordinates.
+    Central helper function for implementing map_coordinates.
 
-    Covers the use case in which the coordinates are non-dask arrays
-    and the input image array consists of a (potentially large) dask array.
+    Receives 'image' as a dask array and computed coordinates.
 
     Implementation details and steps:
-
     1) associate each coordinate in coordinates with the chunk
        it maps to in the input image
     2) for each input image chunk that has been associated to at least one
@@ -412,7 +411,7 @@ def map_coordinates(image, coordinates, order=3,
     coord_chunk_locations = np.array([
         np.searchsorted(np.cumsum(image.chunks[dim]), coordinates[dim])
         for dim in range(image.ndim)]).T
-    # associate coordinates with the chunks they map into
+    # associate coordinates_chunk with the chunks they map into
     required_input_chunks, coord_rc_inds, coord_rc_count = np.unique(
         coord_chunk_locations, axis=0, return_inverse=True, return_counts=True)
 
@@ -437,7 +436,7 @@ def map_coordinates(image, coordinates, order=3,
     # coordinates, define a dask task and apply ndimage.map_coordinates
 
     # prepare building of dask graph
-    name = "map_coordinates-%s" % tokenize(image,
+    name = "map_coordinates_chunk-%s" % tokenize(image,
                                            coordinates,
                                            order,
                                            mode,
@@ -471,4 +470,72 @@ def map_coordinates(image, coordinates, order=3,
     # to match input order
     orig_order = np.argsort([ic for rc_ci in rc_coord_inds for ic in rc_ci])
 
-    return ar[orig_order]
+    return ar[orig_order].compute()
+
+
+def map_coordinates(image, coordinates, order=3,
+                    mode='constant', cval=0.0, prefilter=False):
+    """
+    Wraps ndimage.map_coordinates.
+
+    Both the image and coordinate arrays can be dask arrays.
+
+    For each chunk in the coordinates array, the coordinates are computed
+    and mapped onto the required slices of the input image. One task is
+    is defined for each input image chunk that has been associated to at
+    least one coordinate. The outputs of the tasks are then rearranged to
+    match the input order. For more details see the docstring of
+    '_map_single_coordinates_array_chunk'.
+
+    image : array_like
+        The input array.
+    coordinates : array_like
+        The coordinates at which to sample the input.
+    order : int, optional
+        The order of the spline interpolation, default is 3. The order has to
+        be in the range 0-5.
+    mode : boundary behavior mode, optional
+    cval : float, optional
+        Value to fill past edges of input if mode is 'constant'. Default is 0.0
+    prefilter : bool, optional
+
+    Comments:
+     - in case of a small coordinate array, it might make sense to rechunk
+       it into a single chunk
+    """
+
+    # if coordinate array is not a dask array, convert it into one
+    if type(coordinates) is not da.Array:
+        coordinates = da.from_array(coordinates, chunks=coordinates.shape)
+    else:
+        # make sure indices are not split across chunks, i.e. that there's
+        # no chunking along the first dimension 
+        if len(coordinates.chunks[0]) > 1:
+            coordinates = da.rechunk(coordinates,
+                                    (-1,) + coordinates.chunks[1:])
+
+    # if image array is not a dask array, convert it into one
+    if type(image) is not da.Array:
+        image = da.from_array(image, chunks=image.shape)
+
+    # Map each chunk of the coordinates array onto the entire input image.
+    # 'image' is passed to `_map_single_coordinates_array_chunk` using a
+    # dirty trick: it is split into its components and passed as a delayed
+    # object, which reconstructs the original array when the task is
+    # executed. Therefore two `compute` calls are required to obtain the
+    # final result, one of which is peformed by
+    # `_map_single_coordinates_array_chunk`
+    output = da.map_blocks(
+        _map_single_coordinates_array_chunk,
+        delayed(da.Array)(image.dask, image.name, image.chunks, image.dtype),
+        coordinates,
+        order=order,
+        mode=mode,
+        cval=cval,
+        prefilter=prefilter,
+        dtype=image.dtype,
+        chunks=coordinates.chunks[1:],
+        drop_axis=0,
+    )
+
+    return output
diff --git a/tests/test_dask_image/test_ndinterp/test_map_coordinates.py b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
index 4c212102..318e847f 100644
--- a/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
+++ b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
@@ -51,6 +51,7 @@ def validate_map_coordinates_general(n=2,
 
     # define test coordinates
     coords = np.random.random((n, coord_len)) * im_shape_per_dim + coord_offset
+    coords_da = da.from_array(coords, chunks=(n, coord_chunksize))
 
     # ndimage result
     ints_scipy = scipy.ndimage.map_coordinates(
@@ -65,7 +66,7 @@ def validate_map_coordinates_general(n=2,
     # dask-image result
     ints_dask = dask_image.ndinterp.map_coordinates(
         image_da,
-        coords,
+        coords_da,
         order=interp_order,
         mode=interp_mode,
         cval=0.0,

From 22cedf9056e0ba8be3601bbdfce130a1f52972b0 Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Sun, 8 Jan 2023 23:27:49 +0100
Subject: [PATCH 05/24] Rename image to input

---
 dask_image/ndinterp/__init__.py | 54 ++++++++++++++++-----------------
 1 file changed, 27 insertions(+), 27 deletions(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index b4ea3fc6..483a366a 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -385,17 +385,17 @@ def spline_filter1d(
     return result
 
 
-def _map_single_coordinates_array_chunk(image, coordinates, order=3,
+def _map_single_coordinates_array_chunk(input, coordinates, order=3,
                     mode='constant', cval=0.0, prefilter=False):
     """
     Central helper function for implementing map_coordinates.
 
-    Receives 'image' as a dask array and computed coordinates.
+    Receives 'input' as a dask array and computed coordinates.
 
     Implementation details and steps:
     1) associate each coordinate in coordinates with the chunk
-       it maps to in the input image
-    2) for each input image chunk that has been associated to at least one
+       it maps to in the input
+    2) for each input chunk that has been associated to at least one
        coordinate, calculate the minimal slice required to map
        all coordinates that are associated to it (note that resulting slice
        coordinates can lie outside of the coordinate's chunk)
@@ -405,22 +405,22 @@ def _map_single_coordinates_array_chunk(image, coordinates, order=3,
     """
 
     # STEP 1: Associate each coordinate in coordinates with
-    # the chunk it maps to in the input image
+    # the chunk it maps to in the input array
 
     # determine which chunk each coordinate maps to
     coord_chunk_locations = np.array([
-        np.searchsorted(np.cumsum(image.chunks[dim]), coordinates[dim])
-        for dim in range(image.ndim)]).T
+        np.searchsorted(np.cumsum(input.chunks[dim]), coordinates[dim])
+        for dim in range(input.ndim)]).T
     # associate coordinates_chunk with the chunks they map into
     required_input_chunks, coord_rc_inds, coord_rc_count = np.unique(
         coord_chunk_locations, axis=0, return_inverse=True, return_counts=True)
 
-    # STEP 2: for each input image chunk that has been associated to at least
+    # STEP 2: for each input chunk that has been associated to at least
     # one coordinate, calculate the minimal slice required to map all
     # coordinates that are associated to it (note that resulting slice
     # coordinates can lie outside of the coordinate's chunk)
 
-    input_slices = np.ones((2, len(required_input_chunks), image.ndim)) * np.nan
+    input_slices = np.ones((2, len(required_input_chunks), input.ndim)) * np.nan
     rc_coord_inds = [np.where(coord_rc_inds == irc)[0]
                      for irc in range(len(required_input_chunks))]
     for ic, c in enumerate(coordinates.T):
@@ -436,7 +436,7 @@ def _map_single_coordinates_array_chunk(image, coordinates, order=3,
     # coordinates, define a dask task and apply ndimage.map_coordinates
 
     # prepare building of dask graph
-    name = "map_coordinates_chunk-%s" % tokenize(image,
+    name = "map_coordinates_chunk-%s" % tokenize(input,
                                            coordinates,
                                            order,
                                            mode,
@@ -447,15 +447,15 @@ def _map_single_coordinates_array_chunk(image, coordinates, order=3,
 
     values = []
     for irc, _rc in enumerate(required_input_chunks):
-        offset = np.min([np.max([[0] * image.ndim,
+        offset = np.min([np.max([[0] * input.ndim,
                                  input_slices[0][irc].astype(np.int64)], 0),
-                         np.array(image.shape) - 1], 0)
+                         np.array(input.shape) - 1], 0)
         sl = [slice(offset[dim],
-                    min(image.shape[dim], int(input_slices[1][irc][dim]) + 1))
-              for dim in range(image.ndim)]
+                    min(input.shape[dim], int(input_slices[1][irc][dim]) + 1))
+              for dim in range(input.ndim)]
 
         values.append((ndimage_map_coordinates,
-                       image[tuple(sl)],
+                       input[tuple(sl)],
                        coordinates[:, rc_coord_inds[irc]] - offset[:, None],
                        None,
                        order,
@@ -464,7 +464,7 @@ def _map_single_coordinates_array_chunk(image, coordinates, order=3,
                        prefilter))
 
     dsk = dict(zip(keys, values))
-    ar = da.Array(dsk, name, tuple([list(coord_rc_count)]), image.dtype)
+    ar = da.Array(dsk, name, tuple([list(coord_rc_count)]), input.dtype)
 
     # STEP 4: rearrange outputs of ndimage.map_coordinates
     # to match input order
@@ -473,21 +473,21 @@ def _map_single_coordinates_array_chunk(image, coordinates, order=3,
     return ar[orig_order].compute()
 
 
-def map_coordinates(image, coordinates, order=3,
+def map_coordinates(input, coordinates, order=3,
                     mode='constant', cval=0.0, prefilter=False):
     """
     Wraps ndimage.map_coordinates.
 
-    Both the image and coordinate arrays can be dask arrays.
+    Both the input and coordinate arrays can be dask arrays.
 
     For each chunk in the coordinates array, the coordinates are computed
-    and mapped onto the required slices of the input image. One task is
+    and mapped onto the required slices of the input array. One task is
     is defined for each input image chunk that has been associated to at
     least one coordinate. The outputs of the tasks are then rearranged to
     match the input order. For more details see the docstring of
     '_map_single_coordinates_array_chunk'.
 
-    image : array_like
+    input : array_like
         The input array.
     coordinates : array_like
         The coordinates at which to sample the input.
@@ -515,25 +515,25 @@ def map_coordinates(image, coordinates, order=3,
                                     (-1,) + coordinates.chunks[1:])
 
     # if image array is not a dask array, convert it into one
-    if type(image) is not da.Array:
-        image = da.from_array(image, chunks=image.shape)
+    if type(input) is not da.Array:
+        input = da.from_array(input, chunks=input.shape)
 
-    # Map each chunk of the coordinates array onto the entire input image.
-    # 'image' is passed to `_map_single_coordinates_array_chunk` using a
-    # dirty trick: it is split into its components and passed as a delayed
+    # Map each chunk of the coordinates array onto the entire input array.
+    # 'input' is passed to `_map_single_coordinates_array_chunk` using a bit of
+    # a dirty trick: it is split into its components and passed as a delayed
     # object, which reconstructs the original array when the task is
     # executed. Therefore two `compute` calls are required to obtain the
     # final result, one of which is peformed by
     # `_map_single_coordinates_array_chunk`
     output = da.map_blocks(
         _map_single_coordinates_array_chunk,
-        delayed(da.Array)(image.dask, image.name, image.chunks, image.dtype),
+        delayed(da.Array)(input.dask, input.name, input.chunks, input.dtype),
         coordinates,
         order=order,
         mode=mode,
         cval=cval,
         prefilter=prefilter,
-        dtype=image.dtype,
+        dtype=input.dtype,
         chunks=coordinates.chunks[1:],
         drop_axis=0,
     )

From e42f6affbde495a303b89328ec66c2c3bd96f8ed Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Sun, 8 Jan 2023 23:28:52 +0100
Subject: [PATCH 06/24] Rename image to input consistently.

---
 dask_image/ndinterp/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 483a366a..04fa83a5 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -482,7 +482,7 @@ def map_coordinates(input, coordinates, order=3,
 
     For each chunk in the coordinates array, the coordinates are computed
     and mapped onto the required slices of the input array. One task is
-    is defined for each input image chunk that has been associated to at
+    is defined for each input array chunk that has been associated to at
     least one coordinate. The outputs of the tasks are then rearranged to
     match the input order. For more details see the docstring of
     '_map_single_coordinates_array_chunk'.

From 431de07f34d2391c957bd01a94df9863f1bb7699 Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Sun, 8 Jan 2023 23:46:01 +0100
Subject: [PATCH 07/24] Added tests for different combinations of numpy/dask
 inputs

---
 .../test_ndinterp/test_map_coordinates.py     | 33 ++++++++++---------
 1 file changed, 17 insertions(+), 16 deletions(-)

diff --git a/tests/test_dask_image/test_ndinterp/test_map_coordinates.py b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
index 318e847f..3bade121 100644
--- a/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
+++ b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
@@ -44,37 +44,38 @@ def validate_map_coordinates_general(n=2,
         # not clear on the underlying cause, but this fails on older SciPy
         pytest.skip("requires SciPy >= 1.6.0")
 
-    # define test image
+    # define test input
     np.random.seed(random_seed)
-    image = np.random.random([im_shape_per_dim] * n)
-    image_da = da.from_array(image, chunks=im_chunksize_per_dim)
+    input = np.random.random([im_shape_per_dim] * n)
+    input_da = da.from_array(input, chunks=im_chunksize_per_dim)
 
     # define test coordinates
     coords = np.random.random((n, coord_len)) * im_shape_per_dim + coord_offset
     coords_da = da.from_array(coords, chunks=(n, coord_chunksize))
 
     # ndimage result
-    ints_scipy = scipy.ndimage.map_coordinates(
-        image,
+    mapped_scipy = scipy.ndimage.map_coordinates(
+        input,
         coords,
-        output=None,
         order=interp_order,
         mode=interp_mode,
         cval=0.0,
         prefilter=prefilter)
 
-    # dask-image result
-    ints_dask = dask_image.ndinterp.map_coordinates(
-        image_da,
-        coords_da,
-        order=interp_order,
-        mode=interp_mode,
-        cval=0.0,
-        prefilter=prefilter)
+    # dask-image results
+    for input_array in [input, input_da]:
+        for coords_array in [coords, coords_da]:
+            mapped_dask = dask_image.ndinterp.map_coordinates(
+                input_array,
+                coords_array,
+                order=interp_order,
+                mode=interp_mode,
+                cval=0.0,
+                prefilter=prefilter)
 
-    ints_dask_computed = ints_dask.compute()
+            mapped_dask_computed = mapped_dask.compute()
 
-    assert np.allclose(ints_scipy, ints_dask_computed)
+            assert np.allclose(mapped_scipy, mapped_dask_computed)
 
 
 @pytest.mark.parametrize("n",

From e46dc310184065ddbe66b8ea45ae927df89e92c3 Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Fri, 13 Jan 2023 22:19:04 +0100
Subject: [PATCH 08/24] Vectorized coordinate-to-chunk association

---
 dask_image/ndinterp/__init__.py | 70 +++++++++++++++++++++------------
 1 file changed, 44 insertions(+), 26 deletions(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 04fa83a5..212760a8 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -13,6 +13,7 @@
 import scipy
 from scipy.ndimage import affine_transform as ndimage_affine_transform
 from scipy.ndimage import map_coordinates as ndimage_map_coordinates
+from scipy.ndimage import labeled_comprehension as ndimage_labeled_comprehension
 
 from ..dispatch._dispatch_ndinterp import (
     dispatch_affine_transform,
@@ -407,30 +408,47 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
     # STEP 1: Associate each coordinate in coordinates with
     # the chunk it maps to in the input array
 
-    # determine which chunk each coordinate maps to
-    coord_chunk_locations = np.array([
-        np.searchsorted(np.cumsum(input.chunks[dim]), coordinates[dim])
-        for dim in range(input.ndim)]).T
+    coord_chunk_locations = coordinates.T // np.array(input.chunksize)
+    
     # associate coordinates_chunk with the chunks they map into
-    required_input_chunks, coord_rc_inds, coord_rc_count = np.unique(
-        coord_chunk_locations, axis=0, return_inverse=True, return_counts=True)
+    input_chunk_locations = np.array([i for i in np.ndindex(input.numblocks)])
+    coordinate_chunk_indices = np.sum(coord_chunk_locations *\
+        np.array([np.product(input.numblocks[:dim])
+            for dim in range(input.ndim)])[::-1], axis=1, dtype=np.int64)
+    chunk_indices_count = np.bincount(coordinate_chunk_indices,
+                                      minlength=np.product(input.numblocks))
+    required_input_chunk_indices = np.where(chunk_indices_count > 0)[0]
+    required_input_chunks = input_chunk_locations[required_input_chunk_indices]
+    coord_rc_count = chunk_indices_count[required_input_chunk_indices]
 
     # STEP 2: for each input chunk that has been associated to at least
     # one coordinate, calculate the minimal slice required to map all
     # coordinates that are associated to it (note that resulting slice
     # coordinates can lie outside of the coordinate's chunk)
 
-    input_slices = np.ones((2, len(required_input_chunks), input.ndim)) * np.nan
-    rc_coord_inds = [np.where(coord_rc_inds == irc)[0]
-                     for irc in range(len(required_input_chunks))]
-    for ic, c in enumerate(coordinates.T):
-
-        input_slices[0][coord_rc_inds[ic]]\
-            = np.nanmin([input_slices[0][coord_rc_inds[ic]],
-                         np.floor(c - order // 2)], 0)
-        input_slices[1][coord_rc_inds[ic]]\
-            = np.nanmax([input_slices[1][coord_rc_inds[ic]],
-                         np.ceil(c + order // 2)], 0)
+    input_slices_lower = np.array([np.clip(
+            ndimage_labeled_comprehension(
+                np.floor(coordinates[dim] - order // 2),
+                coordinate_chunk_indices,
+                required_input_chunk_indices,
+                np.min, np.int64, 0),
+            0, input.shape[dim] - 1)
+        for dim in range(input.ndim)])
+
+    input_slices_upper = np.array([np.clip(
+            ndimage_labeled_comprehension(
+                np.ceil(coordinates[dim] + order // 2) + 1,
+                coordinate_chunk_indices,
+                required_input_chunk_indices,
+                np.max, np.int64, 0),
+            0, input.shape[dim])
+        for dim in range(input.ndim)])
+
+    input_slices = np.array([input_slices_lower, input_slices_upper])
+    input_slices = np.swapaxes(input_slices, 1, 2)
+
+    rc_coord_inds = [np.where(coordinate_chunk_indices == irc)[0]
+                     for irc in required_input_chunk_indices]
 
     # STEP 3: For each previously obtained slice and its associated
     # coordinates, define a dask task and apply ndimage.map_coordinates
@@ -447,16 +465,15 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
 
     values = []
     for irc, _rc in enumerate(required_input_chunks):
-        offset = np.min([np.max([[0] * input.ndim,
-                                 input_slices[0][irc].astype(np.int64)], 0),
-                         np.array(input.shape) - 1], 0)
-        sl = [slice(offset[dim],
-                    min(input.shape[dim], int(input_slices[1][irc][dim]) + 1))
-              for dim in range(input.ndim)]
+
+        ric_slice = [slice(input_slices[0][irc][dim],
+                           input_slices[1][irc][dim])
+                                for dim in range(input.ndim)]
+        ric_offset = input_slices[0][irc]
 
         values.append((ndimage_map_coordinates,
-                       input[tuple(sl)],
-                       coordinates[:, rc_coord_inds[irc]] - offset[:, None],
+                       input[tuple(ric_slice)],
+                       coordinates[:, rc_coord_inds[irc]] - ric_offset[:, None],
                        None,
                        order,
                        mode,
@@ -470,7 +487,8 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
     # to match input order
     orig_order = np.argsort([ic for rc_ci in rc_coord_inds for ic in rc_ci])
 
-    return ar[orig_order].compute()
+    # return ar[orig_order].compute()
+    return ar.compute()[orig_order]
 
 
 def map_coordinates(input, coordinates, order=3,

From 85331a25622e6552414043399cab75510305dc2d Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Fri, 20 Jan 2023 19:33:37 +0100
Subject: [PATCH 09/24] Improve comments and refactor

---
 dask_image/ndinterp/__init__.py | 53 ++++++++++++++++++++++-----------
 1 file changed, 36 insertions(+), 17 deletions(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 212760a8..8f89e8ac 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -408,28 +408,44 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
     # STEP 1: Associate each coordinate in coordinates with
     # the chunk it maps to in the input array
 
-    coord_chunk_locations = coordinates.T // np.array(input.chunksize)
+    # get the input chunks each coordinates maps onto
+    coords_input_chunk_locations = coordinates.T // np.array(input.chunksize)
     
-    # associate coordinates_chunk with the chunks they map into
+    # all input chunk locations
     input_chunk_locations = np.array([i for i in np.ndindex(input.numblocks)])
-    coordinate_chunk_indices = np.sum(coord_chunk_locations *\
+
+    # linearize input chunk locations
+    coords_input_chunk_locations_linear = np.sum(coords_input_chunk_locations *\
         np.array([np.product(input.numblocks[:dim])
             for dim in range(input.ndim)])[::-1], axis=1, dtype=np.int64)
-    chunk_indices_count = np.bincount(coordinate_chunk_indices,
+
+    # determine the input chunks that have coords associated and
+    # count how many coords map onto each input chunk
+    chunk_indices_count = np.bincount(coords_input_chunk_locations_linear,
                                       minlength=np.product(input.numblocks))
     required_input_chunk_indices = np.where(chunk_indices_count > 0)[0]
     required_input_chunks = input_chunk_locations[required_input_chunk_indices]
     coord_rc_count = chunk_indices_count[required_input_chunk_indices]
 
+    # inverse mapping: input chunks to coordinates
+    required_input_chunk_coords_indices = \
+        [np.where(coords_input_chunk_locations_linear == irc)[0]
+                     for irc in required_input_chunk_indices]
+
     # STEP 2: for each input chunk that has been associated to at least
     # one coordinate, calculate the minimal slice required to map all
     # coordinates that are associated to it (note that resulting slice
     # coordinates can lie outside of the coordinate's chunk)
 
+    # determine the slices of the input array that are required for
+    # mapping all coordinates associated to a given input chunk.
+    # Note that this slice can be larger than the given chunk when coords
+    # lie at chunk borders.
+    # (probably there's a more efficient way to do this)
     input_slices_lower = np.array([np.clip(
             ndimage_labeled_comprehension(
                 np.floor(coordinates[dim] - order // 2),
-                coordinate_chunk_indices,
+                coords_input_chunk_locations_linear,
                 required_input_chunk_indices,
                 np.min, np.int64, 0),
             0, input.shape[dim] - 1)
@@ -438,22 +454,20 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
     input_slices_upper = np.array([np.clip(
             ndimage_labeled_comprehension(
                 np.ceil(coordinates[dim] + order // 2) + 1,
-                coordinate_chunk_indices,
+                coords_input_chunk_locations_linear,
                 required_input_chunk_indices,
                 np.max, np.int64, 0),
             0, input.shape[dim])
         for dim in range(input.ndim)])
 
-    input_slices = np.array([input_slices_lower, input_slices_upper])
-    input_slices = np.swapaxes(input_slices, 1, 2)
-
-    rc_coord_inds = [np.where(coordinate_chunk_indices == irc)[0]
-                     for irc in required_input_chunk_indices]
+    input_slices = np.array([input_slices_lower, input_slices_upper])\
+        .swapaxes(1, 2)
 
     # STEP 3: For each previously obtained slice and its associated
     # coordinates, define a dask task and apply ndimage.map_coordinates
 
-    # prepare building of dask graph
+    # prepare building dask graph
+    # define one task per associated input chunk
     name = "map_coordinates_chunk-%s" % tokenize(input,
                                            coordinates,
                                            order,
@@ -463,8 +477,9 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
 
     keys = [(name, i) for i in range(len(required_input_chunks))]
 
+    # pair map_coordinates calls with input slices and mapped coordinates
     values = []
-    for irc, _rc in enumerate(required_input_chunks):
+    for irc in range(len(required_input_chunks)):
 
         ric_slice = [slice(input_slices[0][irc][dim],
                            input_slices[1][irc][dim])
@@ -473,21 +488,25 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
 
         values.append((ndimage_map_coordinates,
                        input[tuple(ric_slice)],
-                       coordinates[:, rc_coord_inds[irc]] - ric_offset[:, None],
+                       coordinates[:, required_input_chunk_coords_indices[irc]]
+                                      - ric_offset[:, None],
                        None,
                        order,
                        mode,
                        cval,
                        prefilter))
 
+    # build dask graph
     dsk = dict(zip(keys, values))
     ar = da.Array(dsk, name, tuple([list(coord_rc_count)]), input.dtype)
 
     # STEP 4: rearrange outputs of ndimage.map_coordinates
     # to match input order
-    orig_order = np.argsort([ic for rc_ci in rc_coord_inds for ic in rc_ci])
+    orig_order = np.argsort(
+        [ic for ric_ci in required_input_chunk_coords_indices for ic in ric_ci])
 
-    # return ar[orig_order].compute()
+    # compute result and reorder
+    # (ordering first would probably unnecessarily inflate the task graph)
     return ar.compute()[orig_order]
 
 
@@ -532,7 +551,7 @@ def map_coordinates(input, coordinates, order=3,
             coordinates = da.rechunk(coordinates,
                                     (-1,) + coordinates.chunks[1:])
 
-    # if image array is not a dask array, convert it into one
+    # if the input array is not a dask array, convert it into one
     if type(input) is not da.Array:
         input = da.from_array(input, chunks=input.shape)
 

From e4970e7cd2f5998760c3cc0001e22fbdc027671e Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Fri, 4 Jul 2025 14:47:33 +0200
Subject: [PATCH 10/24] Use np.prod instead of deprecated np.product

---
 dask_image/ndinterp/__init__.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 8f89e8ac..f88a3bc2 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -416,13 +416,13 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
 
     # linearize input chunk locations
     coords_input_chunk_locations_linear = np.sum(coords_input_chunk_locations *\
-        np.array([np.product(input.numblocks[:dim])
+        np.array([np.prod(input.numblocks[:dim])
             for dim in range(input.ndim)])[::-1], axis=1, dtype=np.int64)
 
     # determine the input chunks that have coords associated and
     # count how many coords map onto each input chunk
     chunk_indices_count = np.bincount(coords_input_chunk_locations_linear,
-                                      minlength=np.product(input.numblocks))
+                                      minlength=np.prod(input.numblocks))
     required_input_chunk_indices = np.where(chunk_indices_count > 0)[0]
     required_input_chunks = input_chunk_locations[required_input_chunk_indices]
     coord_rc_count = chunk_indices_count[required_input_chunk_indices]

From c35dde66d837883901dce1d8955394ea39e8acba Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Wed, 9 Jul 2025 00:24:26 +0200
Subject: [PATCH 11/24] Fix flake8 errors

---
 dask_image/ndinterp/__init__.py | 73 +++++++++++++++++++--------------
 1 file changed, 42 insertions(+), 31 deletions(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 51926e48..55918ab7 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -4,7 +4,7 @@
 import math
 from itertools import product
 
-from dask import compute, delayed
+from dask import delayed
 import dask.array as da
 import numpy as np
 from dask.base import tokenize
@@ -12,7 +12,8 @@
 import scipy
 from scipy.ndimage import affine_transform as ndimage_affine_transform
 from scipy.ndimage import map_coordinates as ndimage_map_coordinates
-from scipy.ndimage import labeled_comprehension as ndimage_labeled_comprehension
+from scipy.ndimage import labeled_comprehension as\
+    ndimage_labeled_comprehension
 from scipy.special import sindg, cosdg
 
 from ..dispatch._dispatch_ndinterp import (
@@ -541,8 +542,9 @@ def spline_filter1d(
     return result
 
 
-def _map_single_coordinates_array_chunk(input, coordinates, order=3,
-                    mode='constant', cval=0.0, prefilter=False):
+def _map_single_coordinates_array_chunk(
+        input, coordinates, order=3, mode='constant',
+        cval=0.0, prefilter=False):
     """
     Central helper function for implementing map_coordinates.
 
@@ -565,14 +567,16 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
 
     # get the input chunks each coordinates maps onto
     coords_input_chunk_locations = coordinates.T // np.array(input.chunksize)
-    
+
     # all input chunk locations
     input_chunk_locations = np.array([i for i in np.ndindex(input.numblocks)])
 
     # linearize input chunk locations
-    coords_input_chunk_locations_linear = np.sum(coords_input_chunk_locations *\
-        np.array([np.prod(input.numblocks[:dim])
-            for dim in range(input.ndim)])[::-1], axis=1, dtype=np.int64)
+    coords_input_chunk_locations_linear = np.sum(
+        coords_input_chunk_locations * np.array(
+            [np.prod(input.numblocks[:dim])
+                for dim in range(input.ndim)])[::-1],
+        axis=1, dtype=np.int64)
 
     # determine the input chunks that have coords associated and
     # count how many coords map onto each input chunk
@@ -585,7 +589,7 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
     # inverse mapping: input chunks to coordinates
     required_input_chunk_coords_indices = \
         [np.where(coords_input_chunk_locations_linear == irc)[0]
-                     for irc in required_input_chunk_indices]
+            for irc in required_input_chunk_indices]
 
     # STEP 2: for each input chunk that has been associated to at least
     # one coordinate, calculate the minimal slice required to map all
@@ -623,12 +627,14 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
 
     # prepare building dask graph
     # define one task per associated input chunk
-    name = "map_coordinates_chunk-%s" % tokenize(input,
-                                           coordinates,
-                                           order,
-                                           mode,
-                                           cval,
-                                           prefilter)
+    name = "map_coordinates_chunk-%s" % tokenize(
+        input,
+        coordinates,
+        order,
+        mode,
+        cval,
+        prefilter
+        )
 
     keys = [(name, i) for i in range(len(required_input_chunks))]
 
@@ -636,20 +642,23 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
     values = []
     for irc in range(len(required_input_chunks)):
 
-        ric_slice = [slice(input_slices[0][irc][dim],
-                           input_slices[1][irc][dim])
-                                for dim in range(input.ndim)]
+        ric_slice = [slice(
+            input_slices[0][irc][dim],
+            input_slices[1][irc][dim])
+            for dim in range(input.ndim)]
         ric_offset = input_slices[0][irc]
 
-        values.append((ndimage_map_coordinates,
-                       input[tuple(ric_slice)],
-                       coordinates[:, required_input_chunk_coords_indices[irc]]
-                                      - ric_offset[:, None],
-                       None,
-                       order,
-                       mode,
-                       cval,
-                       prefilter))
+        values.append((
+            ndimage_map_coordinates,
+            input[tuple(ric_slice)],
+            coordinates[:, required_input_chunk_coords_indices[irc]],
+            - ric_offset[:, None],
+            None,
+            order,
+            mode,
+            cval,
+            prefilter
+        ))
 
     # build dask graph
     dsk = dict(zip(keys, values))
@@ -658,7 +667,8 @@ def _map_single_coordinates_array_chunk(input, coordinates, order=3,
     # STEP 4: rearrange outputs of ndimage.map_coordinates
     # to match input order
     orig_order = np.argsort(
-        [ic for ric_ci in required_input_chunk_coords_indices for ic in ric_ci])
+        [ic for ric_ci in required_input_chunk_coords_indices
+            for ic in ric_ci])
 
     # compute result and reorder
     # (ordering first would probably unnecessarily inflate the task graph)
@@ -701,10 +711,11 @@ def map_coordinates(input, coordinates, order=3,
         coordinates = da.from_array(coordinates, chunks=coordinates.shape)
     else:
         # make sure indices are not split across chunks, i.e. that there's
-        # no chunking along the first dimension 
+        # no chunking along the first dimension
         if len(coordinates.chunks[0]) > 1:
-            coordinates = da.rechunk(coordinates,
-                                    (-1,) + coordinates.chunks[1:])
+            coordinates = da.rechunk(
+                coordinates,
+                (-1,) + coordinates.chunks[1:])
 
     # if the input array is not a dask array, convert it into one
     if type(input) is not da.Array:

From d259b432fe395c8e095fa0ce1914ac2ce927b0f3 Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Wed, 9 Jul 2025 00:39:10 +0200
Subject: [PATCH 12/24] Fix error introduced when fixing flake8

---
 dask_image/ndinterp/__init__.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 55918ab7..953779d1 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -651,7 +651,7 @@ def _map_single_coordinates_array_chunk(
         values.append((
             ndimage_map_coordinates,
             input[tuple(ric_slice)],
-            coordinates[:, required_input_chunk_coords_indices[irc]],
+            coordinates[:, required_input_chunk_coords_indices[irc]]
             - ric_offset[:, None],
             None,
             order,

From 3187e9fae2d35430ffa6fcaf5e0b176c58325860 Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Fri, 18 Jul 2025 13:01:10 +0200
Subject: [PATCH 13/24] Add test for out-of-bounds

---
 .../test_ndinterp/test_map_coordinates.py     | 21 +++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/tests/test_dask_image/test_ndinterp/test_map_coordinates.py b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
index 3bade121..633b1e2c 100644
--- a/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
+++ b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
@@ -111,3 +111,24 @@ def test_map_coordinates_large_input():
     dask_image.ndinterp.map_coordinates(
         image_da,
         coords).compute()
+
+
+def test_map_coordinates_out_of_bounds():
+    """
+    This test checks that an error is raised when out-of-bounds
+    coordinates are used.
+    """
+
+    image_da = da.random.random((128, 128))
+
+    coords = np.array([
+        [128] * 2,  # out of bounds
+        [-1] * 2,  # negative coordinates
+        [-1, 128],  # mixed
+    ]).T
+
+    dask_image.ndinterp.map_coordinates(
+        image_da,
+        coords
+        ).compute(scheduler='single-threaded')
+

From 95ac264cb66f326991eb6af52ae3595b720fe36b Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Fri, 18 Jul 2025 15:42:24 +0200
Subject: [PATCH 14/24] Improve test

---
 .../test_ndinterp/test_map_coordinates.py     | 22 ++++++++-----------
 1 file changed, 9 insertions(+), 13 deletions(-)

diff --git a/tests/test_dask_image/test_ndinterp/test_map_coordinates.py b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
index 633b1e2c..a6f55879 100644
--- a/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
+++ b/tests/test_dask_image/test_ndinterp/test_map_coordinates.py
@@ -113,22 +113,18 @@ def test_map_coordinates_large_input():
         coords).compute()
 
 
-def test_map_coordinates_out_of_bounds():
+@pytest.mark.parametrize("interp_mode",
+                         _supported_modes)
+def test_map_coordinates_out_of_bounds(interp_mode):
     """
     This test checks that an error is raised when out-of-bounds
     coordinates are used.
     """
 
-    image_da = da.random.random((128, 128))
-
-    coords = np.array([
-        [128] * 2,  # out of bounds
-        [-1] * 2,  # negative coordinates
-        [-1, 128],  # mixed
-    ]).T
-
-    dask_image.ndinterp.map_coordinates(
-        image_da,
-        coords
-        ).compute(scheduler='single-threaded')
+    kwargs = dict()
+    kwargs['random_seed'] = 0
+    kwargs['interp_mode'] = interp_mode
+    kwargs['im_shape_per_dim'] = 10
+    kwargs['coord_offset'] = 10  # coordinates will be out of bounds
 
+    validate_map_coordinates_general(**kwargs)

From d1aa154d9d3cec4940915434f7cc83455350f4d1 Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Fri, 18 Jul 2025 15:43:30 +0200
Subject: [PATCH 15/24] Map out-of-bounds chunk locations to closest valid
 chunks

---
 dask_image/ndinterp/__init__.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 953779d1..17a17c7c 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -565,9 +565,14 @@ def _map_single_coordinates_array_chunk(
     # STEP 1: Associate each coordinate in coordinates with
     # the chunk it maps to in the input array
 
-    # get the input chunks each coordinates maps onto
+    # get the input chunks each coordinate maps onto
     coords_input_chunk_locations = coordinates.T // np.array(input.chunksize)
 
+    # map out-of-bounds chunk locations to valid input chunks
+    coords_input_chunk_locations = np.clip(
+        coords_input_chunk_locations, 0, np.array(input.numblocks) - 1
+    )
+
     # all input chunk locations
     input_chunk_locations = np.array([i for i in np.ndindex(input.numblocks)])
 

From 2e728b43da8101d4f71d82edae6f46b6cccd22a4 Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Mon, 28 Jul 2025 15:33:36 +0200
Subject: [PATCH 16/24] Clarify prefilter default value in docstring

---
 dask_image/ndinterp/__init__.py | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 17a17c7c..7793577e 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -705,10 +705,17 @@ def map_coordinates(input, coordinates, order=3,
     cval : float, optional
         Value to fill past edges of input if mode is 'constant'. Default is 0.0
     prefilter : bool, optional
+        If True, prefilter the input before interpolation. Default is False.
+        Warning: prefilter is True by default in `scipy.ndimage.map_coordinates`.
+        Prefiltering here is performed on a chunk-by-chunk basis, which may lead
+        to different results than `scipy.ndimage.map_coordinates` in case of
+        chunked input arrays and order > 1.
 
     Comments:
-     - in case of a small coordinate array, it might make sense to rechunk
-       it into a single chunk
+      - in case of a small coordinate array, it might make sense to rechunk
+        it into a single chunk
+      - note the different default for `prefilter` compared to
+        `scipy.ndimage.map_coordinates`, which is True by default.
     """
 
     # if coordinate array is not a dask array, convert it into one

From b3471e865bb3b738c1acd4169eb9ee6977d5125f Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Mon, 28 Jul 2025 15:53:30 +0200
Subject: [PATCH 17/24] Fix line length

---
 dask_image/ndinterp/__init__.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 7793577e..6256cd6e 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -706,10 +706,11 @@ def map_coordinates(input, coordinates, order=3,
         Value to fill past edges of input if mode is 'constant'. Default is 0.0
     prefilter : bool, optional
         If True, prefilter the input before interpolation. Default is False.
-        Warning: prefilter is True by default in `scipy.ndimage.map_coordinates`.
-        Prefiltering here is performed on a chunk-by-chunk basis, which may lead
-        to different results than `scipy.ndimage.map_coordinates` in case of
-        chunked input arrays and order > 1.
+        Warning: prefilter is True by default in
+        `scipy.ndimage.map_coordinates`. Prefiltering here is performed on a
+        chunk-by-chunk basis, which may lead to different results than
+        `scipy.ndimage.map_coordinates` in case of chunked input arrays
+        and order > 1.
 
     Comments:
       - in case of a small coordinate array, it might make sense to rechunk

From a0102635b4f0151e6cfb6d6936a8abb0309b6d96 Mon Sep 17 00:00:00 2001
From: Marvin Albert <marvin.albert@gmail.com>
Date: Tue, 30 Sep 2025 19:18:25 +0200
Subject: [PATCH 18/24] Improve top level docstring

---
 dask_image/ndinterp/__init__.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 6256cd6e..bb4e0960 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -694,6 +694,13 @@ def map_coordinates(input, coordinates, order=3,
     match the input order. For more details see the docstring of
     '_map_single_coordinates_array_chunk'.
 
+    Using this function together with schedulers that support
+    parallelism (threads, processes, distributed) makes sense in the
+    case of either a large input array or a large coordinates array.
+    When both arrays are large, it is recommended to use the
+    single-threaded scheduler. A scheduler can be specified using e.g.
+    `with dask.config.set(scheduler='threads'): ...`.
+
     input : array_like
         The input array.
     coordinates : array_like

From 42bf4111d72171091bee55d9d7217abffb4a6225 Mon Sep 17 00:00:00 2001
From: GenevieveBuckley <30920819+GenevieveBuckley@users.noreply.github.com>
Date: Tue, 11 Nov 2025 16:02:27 +1100
Subject: [PATCH 19/24] Add map_coordinates to __all__, makes visible in public
 API function list

---
 dask_image/ndinterp/__init__.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index bb4e0960..b99dcde8 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -26,6 +26,7 @@
 
 __all__ = [
     "affine_transform",
+    "map_coordinates",
     "rotate",
     "spline_filter",
     "spline_filter1d",

From 8d7804736fa8e90934be332c6642b44b9051cc97 Mon Sep 17 00:00:00 2001
From: GenevieveBuckley <30920819+GenevieveBuckley@users.noreply.github.com>
Date: Tue, 11 Nov 2025 16:08:30 +1100
Subject: [PATCH 20/24] Comment add link to dask discourse with delayed compute
 discussion

---
 dask_image/ndinterp/__init__.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index b99dcde8..48156340 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -749,6 +749,7 @@ def map_coordinates(input, coordinates, order=3,
     # executed. Therefore two `compute` calls are required to obtain the
     # final result, one of which is peformed by
     # `_map_single_coordinates_array_chunk`
+    # Discussion https://dask.discourse.group/t/passing-dask-objects-to-delayed-computations-without-triggering-compute/1441
     output = da.map_blocks(
         _map_single_coordinates_array_chunk,
         delayed(da.Array)(input.dask, input.name, input.chunks, input.dtype),

From 148bdcf1435ed10f33d3ff5cabc8eea5855ce2d7 Mon Sep 17 00:00:00 2001
From: GenevieveBuckley <30920819+GenevieveBuckley@users.noreply.github.com>
Date: Tue, 11 Nov 2025 16:08:49 +1100
Subject: [PATCH 21/24] map_overlap more details about prefilter kwarg added to
 docstring

---
 dask_image/ndinterp/__init__.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 48156340..72ba4895 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -719,6 +719,11 @@ def map_coordinates(input, coordinates, order=3,
         chunk-by-chunk basis, which may lead to different results than
         `scipy.ndimage.map_coordinates` in case of chunked input arrays
         and order > 1.
+        Note: prefilter is not necessary when:
+          - You are using nearest neighbour interpolation, by setting order=0
+          - You are using linear interpolation, by setting order=1, or
+          - You have already prefiltered the input array,
+          using the spline_filter or spline_filter1d functions.
 
     Comments:
       - in case of a small coordinate array, it might make sense to rechunk

From 513ba17973079d7d33711f4e044e3735e7f5b27f Mon Sep 17 00:00:00 2001
From: GenevieveBuckley <30920819+GenevieveBuckley@users.noreply.github.com>
Date: Tue, 11 Nov 2025 16:27:18 +1100
Subject: [PATCH 22/24] ndinterp.map_overlap make clear that GPU cupy arrays
 are not supported by this implementation

---
 dask_image/ndinterp/__init__.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 72ba4895..6fab0e36 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -22,6 +22,7 @@
     dispatch_spline_filter,
     dispatch_spline_filter1d,
 )
+from ..dispatch._utils import get_type
 from ..ndfilters._utils import _get_depth_boundary
 
 __all__ = [
@@ -687,6 +688,7 @@ def map_coordinates(input, coordinates, order=3,
     Wraps ndimage.map_coordinates.
 
     Both the input and coordinate arrays can be dask arrays.
+    GPU arrays are not supported.
 
     For each chunk in the coordinates array, the coordinates are computed
     and mapped onto the required slices of the input array. One task is
@@ -731,6 +733,8 @@ def map_coordinates(input, coordinates, order=3,
       - note the different default for `prefilter` compared to
         `scipy.ndimage.map_coordinates`, which is True by default.
     """
+    if "cupy" in str(get_type(input)) or "cupy" in str(get_type(coordinates)):
+        raise NotImplementedError("GPU cupy arrays are not supported by dask_image.ndinterp.map_overlap")
 
     # if coordinate array is not a dask array, convert it into one
     if type(coordinates) is not da.Array:

From 59310061cfe172b6841e812e73216daec14649a8 Mon Sep 17 00:00:00 2001
From: GenevieveBuckley <30920819+GenevieveBuckley@users.noreply.github.com>
Date: Tue, 11 Nov 2025 16:39:09 +1100
Subject: [PATCH 23/24] Fix flake8 line too long errors

---
 dask_image/ndinterp/__init__.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/dask_image/ndinterp/__init__.py b/dask_image/ndinterp/__init__.py
index 6fab0e36..0175b7c4 100644
--- a/dask_image/ndinterp/__init__.py
+++ b/dask_image/ndinterp/__init__.py
@@ -734,7 +734,9 @@ def map_coordinates(input, coordinates, order=3,
         `scipy.ndimage.map_coordinates`, which is True by default.
     """
     if "cupy" in str(get_type(input)) or "cupy" in str(get_type(coordinates)):
-        raise NotImplementedError("GPU cupy arrays are not supported by dask_image.ndinterp.map_overlap")
+        raise NotImplementedError(
+            "GPU cupy arrays are not supported by "
+            "dask_image.ndinterp.map_overlap")
 
     # if coordinate array is not a dask array, convert it into one
     if type(coordinates) is not da.Array:
@@ -758,7 +760,7 @@ def map_coordinates(input, coordinates, order=3,
     # executed. Therefore two `compute` calls are required to obtain the
     # final result, one of which is peformed by
     # `_map_single_coordinates_array_chunk`
-    # Discussion https://dask.discourse.group/t/passing-dask-objects-to-delayed-computations-without-triggering-compute/1441
+    # Discussion https://dask.discourse.group/t/passing-dask-objects-to-delayed-computations-without-triggering-compute/1441 # noqa: E501
     output = da.map_blocks(
         _map_single_coordinates_array_chunk,
         delayed(da.Array)(input.dask, input.name, input.chunks, input.dtype),

From 6fbaeb7284c8c400f56c3b2bb1d519b0309925d5 Mon Sep 17 00:00:00 2001
From: GenevieveBuckley <30920819+GenevieveBuckley@users.noreply.github.com>
Date: Tue, 11 Nov 2025 16:59:43 +1100
Subject: [PATCH 24/24] Fix coverage.rst table syntax

---
 docs/coverage.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/coverage.rst b/docs/coverage.rst
index a43465c4..6b4b52fd 100644
--- a/docs/coverage.rst
+++ b/docs/coverage.rst
@@ -190,6 +190,7 @@ This table shows which SciPy ndimage functions are supported by dask-image.
    * - ``map_coordinates``
      - ✓
      - ✓
+     -
    * - ``maximum``
      - ✓
      - ✓