Implement map_coordinates functionality

include/bioimage_cpp/transformation/affine.hxx

@@ -495,6 +495,40 @@ inline double cubic_3d(
     return value;
 }
 
+// Dispatch to the per-voxel sampler for the requested interpolation order. Shared by the affine and
+// map_coordinates kernels so the interpolation backend lives in exactly one place. `order` must be in
+// 0..5 (validated by the public entry points before the sampling loop).
+template <class T>
+inline double sample_2d(
+    const T *data, std::ptrdiff_t in_h, std::ptrdiff_t in_w,
+    double cy, double cx, const int order, double fill
+) {
+    switch (order) {
+        case 0: return nearest_2d(data, in_h, in_w, cy, cx, fill);
+        case 1: return linear_2d(data, in_h, in_w, cy, cx, fill);
+        case 2: return bspline_2d<2>(data, in_h, in_w, cy, cx, fill);
+        case 3: return cubic_2d(data, in_h, in_w, cy, cx, fill);
+        case 4: return bspline_2d<4>(data, in_h, in_w, cy, cx, fill);
+        default: return bspline_2d<5>(data, in_h, in_w, cy, cx, fill);  // 5
+    }
+}
+
+template <class T>
+inline double sample_3d(
+    const T *data,
+    std::ptrdiff_t in_d, std::ptrdiff_t in_h, std::ptrdiff_t in_w,
+    double cz, double cy, double cx, const int order, double fill
+) {
+    switch (order) {
+        case 0: return nearest_3d(data, in_d, in_h, in_w, cz, cy, cx, fill);
+        case 1: return linear_3d(data, in_d, in_h, in_w, cz, cy, cx, fill);
+        case 2: return bspline_3d<2>(data, in_d, in_h, in_w, cz, cy, cx, fill);
+        case 3: return cubic_3d(data, in_d, in_h, in_w, cz, cy, cx, fill);
+        case 4: return bspline_3d<4>(data, in_d, in_h, in_w, cz, cy, cx, fill);
+        default: return bspline_3d<5>(data, in_d, in_h, in_w, cz, cy, cx, fill);  // 5
+    }
+}
+
 } // namespace detail
 
 // ----- 2D entry point -------------------------------------------------------
@@ -554,27 +588,7 @@ void affine_transform_2d(
         double cy = row_y;
         double cx = row_x;
         for (std::ptrdiff_t j = 0; j < out_w; ++j) {
-            double value;
-            switch (order) {
-                case 0:
-                    value = detail::nearest_2d(in_data, in_h, in_w, cy, cx, fill);
-                    break;
-                case 1:
-                    value = detail::linear_2d(in_data, in_h, in_w, cy, cx, fill);
-                    break;
-                case 2:
-                    value = detail::bspline_2d<2>(in_data, in_h, in_w, cy, cx, fill);
-                    break;
-                case 3:
-                    value = detail::cubic_2d(in_data, in_h, in_w, cy, cx, fill);
-                    break;
-                case 4:
-                    value = detail::bspline_2d<4>(in_data, in_h, in_w, cy, cx, fill);
-                    break;
-                default:  // 5
-                    value = detail::bspline_2d<5>(in_data, in_h, in_w, cy, cx, fill);
-                    break;
-            }
+            const double value = detail::sample_2d(in_data, in_h, in_w, cy, cx, order, fill);
             *out_ptr++ = detail::to_output<T>(value);
             cy += m01;
             cx += m11;
@@ -655,33 +669,8 @@ void affine_transform_3d(
             double cy = row_y;
             double cx = row_x;
             for (std::ptrdiff_t j = 0; j < out_w; ++j) {
-                double value;
-                switch (order) {
-                    case 0:
-                        value = detail::nearest_3d(in_data, in_d, in_h, in_w,
-                                                   cz, cy, cx, fill);
-                        break;
-                    case 1:
-                        value = detail::linear_3d(in_data, in_d, in_h, in_w,
-                                                  cz, cy, cx, fill);
-                        break;
-                    case 2:
-                        value = detail::bspline_3d<2>(in_data, in_d, in_h, in_w,
-                                                      cz, cy, cx, fill);
-                        break;
-                    case 3:
-                        value = detail::cubic_3d(in_data, in_d, in_h, in_w,
-                                                 cz, cy, cx, fill);
-                        break;
-                    case 4:
-                        value = detail::bspline_3d<4>(in_data, in_d, in_h, in_w,
-                                                      cz, cy, cx, fill);
-                        break;
-                    default:  // 5
-                        value = detail::bspline_3d<5>(in_data, in_d, in_h, in_w,
-                                                      cz, cy, cx, fill);
-                        break;
-                }
+                const double value = detail::sample_3d(in_data, in_d, in_h, in_w,
+                                                       cz, cy, cx, order, fill);
                 *out_ptr++ = detail::to_output<T>(value);
                 cz += m02;
                 cy += m12;

include/bioimage_cpp/transformation/coordinate.hxx

@@ -0,0 +1,129 @@
+#pragma once
+
+// Coordinate-based resampling (the analogue of ``scipy.ndimage.map_coordinates``). For every output
+// voxel the source coordinate to sample is read from an explicit ``coordinates`` array, instead of
+// being computed from an affine matrix. The per-voxel interpolation reuses the affine samplers
+// (``detail::sample_2d`` / ``detail::sample_3d``) so the interpolation backend lives in one place.
+//
+// This kernel is pure and in-memory (NumPy in, NumPy out); reading the source data and producing the
+// coordinate (deformation) field are the caller's responsibility.
+
+#include "bioimage_cpp/transformation/affine.hxx"
+
+namespace bioimage_cpp::transformation {
+
+namespace detail {
+
+// Validate that ``coordinates`` is a ``(D, *output_shape)`` field, i.e. it carries one source
+// coordinate per output voxel along its leading axis.
+template <std::size_t D, class T>
+void require_coordinates(const ConstArrayView<double> &coordinates, const ArrayView<T> &output) {
+    if (coordinates.ndim() != static_cast<std::ptrdiff_t>(D + 1)) {
+        throw std::invalid_argument(
+            "coordinates must have ndim=" + std::to_string(D + 1) +
+            ", got ndim=" + std::to_string(coordinates.ndim())
+        );
+    }
+    if (coordinates.shape[0] != static_cast<std::ptrdiff_t>(D)) {
+        throw std::invalid_argument(
+            "coordinates.shape[0] must equal the data dimension " + std::to_string(D) +
+            ", got " + std::to_string(coordinates.shape[0])
+        );
+    }
+    for (std::size_t axis = 0; axis < D; ++axis) {
+        if (coordinates.shape[axis + 1] != output.shape[axis]) {
+            throw std::invalid_argument(
+                "coordinates spatial shape (coordinates.shape[1:]) must match the output shape"
+            );
+        }
+    }
+}
+
+} // namespace detail
+
+// ----- 2D entry point -------------------------------------------------------
+
+template <class T>
+void map_coordinates_2d(
+    const ConstArrayView<T> &input,
+    ArrayView<T> &output,
+    const ConstArrayView<double> &coordinates,
+    const int order,
+    const T fill_value
+) {
+    detail::require_views<2, T>(input, output);
+    detail::require_coordinates<2, T>(coordinates, output);
+    if (order < 0 || order > 5) {
+        throw std::invalid_argument(
+            "order must be in 0..5, got " + std::to_string(order)
+        );
+    }
+
+    const auto out_h = output.shape[0];
+    const auto out_w = output.shape[1];
+    if (out_h == 0 || out_w == 0) return;
+
+    const auto in_h = input.shape[0];
+    const auto in_w = input.shape[1];
+
+    // coordinates is C-contiguous with shape (2, out_h, out_w), so axis-d coordinates form a
+    // contiguous block of n_out doubles starting at d * n_out.
+    const std::ptrdiff_t n_out = out_h * out_w;
+    const double *cy = coordinates.data;
+    const double *cx = coordinates.data + n_out;
+
+    const double fill = static_cast<double>(fill_value);
+    const T *in_data = input.data;
+    T *out_ptr = output.data;
+
+    for (std::ptrdiff_t p = 0; p < n_out; ++p) {
+        const double value = detail::sample_2d(in_data, in_h, in_w, cy[p], cx[p], order, fill);
+        out_ptr[p] = detail::to_output<T>(value);
+    }
+}
+
+// ----- 3D entry point -------------------------------------------------------
+
+template <class T>
+void map_coordinates_3d(
+    const ConstArrayView<T> &input,
+    ArrayView<T> &output,
+    const ConstArrayView<double> &coordinates,
+    const int order,
+    const T fill_value
+) {
+    detail::require_views<3, T>(input, output);
+    detail::require_coordinates<3, T>(coordinates, output);
+    if (order < 0 || order > 5) {
+        throw std::invalid_argument(
+            "order must be in 0..5, got " + std::to_string(order)
+        );
+    }
+
+    const auto out_d = output.shape[0];
+    const auto out_h = output.shape[1];
+    const auto out_w = output.shape[2];
+    if (out_d == 0 || out_h == 0 || out_w == 0) return;
+
+    const auto in_d = input.shape[0];
+    const auto in_h = input.shape[1];
+    const auto in_w = input.shape[2];
+
+    // coordinates is C-contiguous with shape (3, out_d, out_h, out_w).
+    const std::ptrdiff_t n_out = out_d * out_h * out_w;
+    const double *cz = coordinates.data;
+    const double *cy = coordinates.data + n_out;
+    const double *cx = coordinates.data + 2 * n_out;
+
+    const double fill = static_cast<double>(fill_value);
+    const T *in_data = input.data;
+    T *out_ptr = output.data;
+
+    for (std::ptrdiff_t p = 0; p < n_out; ++p) {
+        const double value = detail::sample_3d(in_data, in_d, in_h, in_w,
+                                               cz[p], cy[p], cx[p], order, fill);
+        out_ptr[p] = detail::to_output<T>(value);
+    }
+}
+
+} // namespace bioimage_cpp::transformation

src/bindings/transformation.cxx

@@ -3,6 +3,7 @@
 #include "bioimage_cpp/array_view.hxx"
 #include "bioimage_cpp/detail/grid.hxx"
 #include "bioimage_cpp/transformation/affine.hxx"
+#include "bioimage_cpp/transformation/coordinate.hxx"
 
 #include <nanobind/ndarray.h>
 
@@ -113,6 +114,86 @@ void bind_affine_for_dtype(nb::module_ &m, const char *name_2d, const char *name
     );
 }
 
+template <std::size_t D, class T>
+OutputArray<T> map_coordinates_t(
+    ConstArray<T> input,
+    MatrixArray coordinates,
+    OutputArray<T> output,
+    const int order,
+    const T fill_value
+) {
+    if (input.ndim() != D) {
+        throw std::invalid_argument(
+            "input must have ndim=" + std::to_string(D) +
+            ", got ndim=" + std::to_string(input.ndim())
+        );
+    }
+    if (output.ndim() != D) {
+        throw std::invalid_argument(
+            "output must have ndim=" + std::to_string(D) +
+            ", got ndim=" + std::to_string(output.ndim())
+        );
+    }
+    if (coordinates.ndim() != D + 1) {
+        throw std::invalid_argument(
+            "coordinates must have ndim=" + std::to_string(D + 1) +
+            ", got ndim=" + std::to_string(coordinates.ndim())
+        );
+    }
+
+    const auto input_shape = shape_of(input);
+    const auto input_strides = detail::c_order_strides(input_shape);
+    const auto coordinates_shape = shape_of(coordinates);
+    const auto coordinates_strides = detail::c_order_strides(coordinates_shape);
+    const auto output_shape = shape_of(output);
+    const auto output_strides = detail::c_order_strides(output_shape);
+
+    ConstArrayView<T> input_view{input.data(), input_shape, input_strides};
+    ArrayView<T> output_view{output.data(), output_shape, output_strides};
+    ConstArrayView<double> coordinates_view{
+        coordinates.data(), coordinates_shape, coordinates_strides
+    };
+
+    {
+        nb::gil_scoped_release release;
+        if constexpr (D == 2) {
+            transformation::map_coordinates_2d<T>(
+                input_view, output_view, coordinates_view, order, fill_value
+            );
+        } else {
+            transformation::map_coordinates_3d<T>(
+                input_view, output_view, coordinates_view, order, fill_value
+            );
+        }
+    }
+
+    return output;
+}
+
+template <class T>
+void bind_map_coordinates_for_dtype(nb::module_ &m, const char *name_2d, const char *name_3d) {
+    m.def(
+        name_2d,
+        &map_coordinates_t<2, T>,
+        nb::arg("input"),
+        nb::arg("coordinates"),
+        nb::arg("output"),
+        nb::arg("order"),
+        nb::arg("fill_value"),
+        "Apply a 2D coordinate-based resampling into a pre-allocated NumPy array."
+    );
+    m.def(
+        name_3d,
+        &map_coordinates_t<3, T>,
+        nb::arg("input"),
+        nb::arg("coordinates"),
+        nb::arg("output"),
+        nb::arg("order"),
+        nb::arg("fill_value"),
+        "Apply a 3D coordinate-based resampling into a pre-allocated NumPy array."
+    );
+}
+
 } // namespace
 
 void bind_transformation(nb::module_ &m) {
@@ -146,6 +227,37 @@ void bind_transformation(nb::module_ &m) {
     bind_affine_for_dtype<double>(
         m, "_affine_transform_2d_float64", "_affine_transform_3d_float64"
     );
+
+    bind_map_coordinates_for_dtype<std::uint8_t>(
+        m, "_map_coordinates_2d_uint8", "_map_coordinates_3d_uint8"
+    );
+    bind_map_coordinates_for_dtype<std::uint16_t>(
+        m, "_map_coordinates_2d_uint16", "_map_coordinates_3d_uint16"
+    );
+    bind_map_coordinates_for_dtype<std::uint32_t>(
+        m, "_map_coordinates_2d_uint32", "_map_coordinates_3d_uint32"
+    );
+    bind_map_coordinates_for_dtype<std::uint64_t>(
+        m, "_map_coordinates_2d_uint64", "_map_coordinates_3d_uint64"
+    );
+    bind_map_coordinates_for_dtype<std::int8_t>(
+        m, "_map_coordinates_2d_int8", "_map_coordinates_3d_int8"
+    );
+    bind_map_coordinates_for_dtype<std::int16_t>(
+        m, "_map_coordinates_2d_int16", "_map_coordinates_3d_int16"
+    );
+    bind_map_coordinates_for_dtype<std::int32_t>(
+        m, "_map_coordinates_2d_int32", "_map_coordinates_3d_int32"
+    );
+    bind_map_coordinates_for_dtype<std::int64_t>(
+        m, "_map_coordinates_2d_int64", "_map_coordinates_3d_int64"
+    );
+    bind_map_coordinates_for_dtype<float>(
+        m, "_map_coordinates_2d_float32", "_map_coordinates_3d_float32"
+    );
+    bind_map_coordinates_for_dtype<double>(
+        m, "_map_coordinates_2d_float64", "_map_coordinates_3d_float64"
+    );
 }
 
 } // namespace bioimage_cpp::bindings

src/bioimage_cpp/transformation/__init__.py

@@ -3,11 +3,13 @@
 from ._transformation import (
     affine_transform,
     compute_anti_aliasing_sigma,
+    map_coordinates,
     resample,
 )
 
 __all__ = [
     "affine_transform",
     "compute_anti_aliasing_sigma",
+    "map_coordinates",
     "resample",
 ]