From bae4a5050f372ba530bbbe920b899ac1f489a9e8 Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Tue, 30 Jun 2026 23:41:32 +0200 Subject: [PATCH] Implement map_coordinates functionality --- .../bioimage_cpp/transformation/affine.hxx | 85 +++++------- .../transformation/coordinate.hxx | 129 ++++++++++++++++++ src/bindings/transformation.cxx | 112 +++++++++++++++ src/bioimage_cpp/transformation/__init__.py | 2 + .../transformation/_transformation.py | 106 ++++++++++++++ tests/test_transformation.py | 121 ++++++++++++++++ 6 files changed, 507 insertions(+), 48 deletions(-) create mode 100644 include/bioimage_cpp/transformation/coordinate.hxx diff --git a/include/bioimage_cpp/transformation/affine.hxx b/include/bioimage_cpp/transformation/affine.hxx index bd8070f..c7391f8 100644 --- a/include/bioimage_cpp/transformation/affine.hxx +++ b/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 +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 +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(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(value); cz += m02; cy += m12; diff --git a/include/bioimage_cpp/transformation/coordinate.hxx b/include/bioimage_cpp/transformation/coordinate.hxx new file mode 100644 index 0000000..7b5eec9 --- /dev/null +++ b/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 +void require_coordinates(const ConstArrayView &coordinates, const ArrayView &output) { + if (coordinates.ndim() != static_cast(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(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 +void map_coordinates_2d( + const ConstArrayView &input, + ArrayView &output, + const ConstArrayView &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(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(value); + } +} + +// ----- 3D entry point ------------------------------------------------------- + +template +void map_coordinates_3d( + const ConstArrayView &input, + ArrayView &output, + const ConstArrayView &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(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(value); + } +} + +} // namespace bioimage_cpp::transformation diff --git a/src/bindings/transformation.cxx b/src/bindings/transformation.cxx index 8b74f63..716fecb 100644 --- a/src/bindings/transformation.cxx +++ b/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 @@ -113,6 +114,86 @@ void bind_affine_for_dtype(nb::module_ &m, const char *name_2d, const char *name ); } +template +OutputArray map_coordinates_t( + ConstArray input, + MatrixArray coordinates, + OutputArray 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 input_view{input.data(), input_shape, input_strides}; + ArrayView output_view{output.data(), output_shape, output_strides}; + ConstArrayView coordinates_view{ + coordinates.data(), coordinates_shape, coordinates_strides + }; + + { + nb::gil_scoped_release release; + if constexpr (D == 2) { + transformation::map_coordinates_2d( + input_view, output_view, coordinates_view, order, fill_value + ); + } else { + transformation::map_coordinates_3d( + input_view, output_view, coordinates_view, order, fill_value + ); + } + } + + return output; +} + +template +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( m, "_affine_transform_2d_float64", "_affine_transform_3d_float64" ); + + bind_map_coordinates_for_dtype( + m, "_map_coordinates_2d_uint8", "_map_coordinates_3d_uint8" + ); + bind_map_coordinates_for_dtype( + m, "_map_coordinates_2d_uint16", "_map_coordinates_3d_uint16" + ); + bind_map_coordinates_for_dtype( + m, "_map_coordinates_2d_uint32", "_map_coordinates_3d_uint32" + ); + bind_map_coordinates_for_dtype( + m, "_map_coordinates_2d_uint64", "_map_coordinates_3d_uint64" + ); + bind_map_coordinates_for_dtype( + m, "_map_coordinates_2d_int8", "_map_coordinates_3d_int8" + ); + bind_map_coordinates_for_dtype( + m, "_map_coordinates_2d_int16", "_map_coordinates_3d_int16" + ); + bind_map_coordinates_for_dtype( + m, "_map_coordinates_2d_int32", "_map_coordinates_3d_int32" + ); + bind_map_coordinates_for_dtype( + m, "_map_coordinates_2d_int64", "_map_coordinates_3d_int64" + ); + bind_map_coordinates_for_dtype( + m, "_map_coordinates_2d_float32", "_map_coordinates_3d_float32" + ); + bind_map_coordinates_for_dtype( + m, "_map_coordinates_2d_float64", "_map_coordinates_3d_float64" + ); } } // namespace bioimage_cpp::bindings diff --git a/src/bioimage_cpp/transformation/__init__.py b/src/bioimage_cpp/transformation/__init__.py index b63a096..31c79af 100644 --- a/src/bioimage_cpp/transformation/__init__.py +++ b/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", ] diff --git a/src/bioimage_cpp/transformation/_transformation.py b/src/bioimage_cpp/transformation/_transformation.py index 934acd2..205f77c 100644 --- a/src/bioimage_cpp/transformation/_transformation.py +++ b/src/bioimage_cpp/transformation/_transformation.py @@ -36,6 +36,32 @@ np.dtype("float64"): _core._affine_transform_3d_float64, } +_MAP_COORDINATES_2D_BY_DTYPE = { + np.dtype("uint8"): _core._map_coordinates_2d_uint8, + np.dtype("uint16"): _core._map_coordinates_2d_uint16, + np.dtype("uint32"): _core._map_coordinates_2d_uint32, + np.dtype("uint64"): _core._map_coordinates_2d_uint64, + np.dtype("int8"): _core._map_coordinates_2d_int8, + np.dtype("int16"): _core._map_coordinates_2d_int16, + np.dtype("int32"): _core._map_coordinates_2d_int32, + np.dtype("int64"): _core._map_coordinates_2d_int64, + np.dtype("float32"): _core._map_coordinates_2d_float32, + np.dtype("float64"): _core._map_coordinates_2d_float64, +} + +_MAP_COORDINATES_3D_BY_DTYPE = { + np.dtype("uint8"): _core._map_coordinates_3d_uint8, + np.dtype("uint16"): _core._map_coordinates_3d_uint16, + np.dtype("uint32"): _core._map_coordinates_3d_uint32, + np.dtype("uint64"): _core._map_coordinates_3d_uint64, + np.dtype("int8"): _core._map_coordinates_3d_int8, + np.dtype("int16"): _core._map_coordinates_3d_int16, + np.dtype("int32"): _core._map_coordinates_3d_int32, + np.dtype("int64"): _core._map_coordinates_3d_int64, + np.dtype("float32"): _core._map_coordinates_3d_float32, + np.dtype("float64"): _core._map_coordinates_3d_float64, +} + def _normalize_matrix(matrix, ndim: int) -> np.ndarray: array = np.asarray(matrix, dtype=np.float64) @@ -340,3 +366,83 @@ def affine_transform( output = _validate_or_allocate_out(out, output_shape, dtype) run(contiguous, normalized_matrix, starts, output, order, typed_fill) return output + + +def _normalize_coordinates(coordinates, ndim: int) -> tuple[np.ndarray, tuple[int, ...]]: + array = np.ascontiguousarray(coordinates, dtype=np.float64) + if array.ndim != ndim + 1: + raise ValueError( + f"coordinates must have ndim={ndim + 1} (a leading axis of length {ndim} plus the " + f"output shape), got ndim={array.ndim}" + ) + if array.shape[0] != ndim: + raise ValueError( + f"coordinates.shape[0] must equal the data dimension {ndim}, got {array.shape[0]}" + ) + output_shape = tuple(int(s) for s in array.shape[1:]) + return array, output_shape + + +def map_coordinates( + data: np.ndarray, + coordinates, + *, + order: int = 1, + fill_value: Number = 0, + out: np.ndarray | None = None, +) -> np.ndarray: + """Resample a 2D or 3D array at explicit, per-output-voxel source coordinates. + + Like ``scipy.ndimage.map_coordinates``, but specialized to a same-rank resampling: for each + output voxel the source coordinate to sample is read from ``coordinates`` and the input is + interpolated there. ``coordinates`` has shape ``(ndim, *output_shape)`` in NumPy axis order, so + ``coordinates[d]`` holds the axis-``d`` source coordinate of every output voxel. The output has + the same number of dimensions as ``data`` (a 2D or 3D resampling), so ``output_shape`` = + ``coordinates.shape[1:]`` has ``ndim`` entries. This is the deformation-field counterpart of + :func:`affine_transform` (which derives the same per-voxel source coordinate from an affine + matrix) and shares its interpolation backend. + + Supported interpolation orders match :func:`affine_transform`: + + - ``0`` — nearest neighbour. + - ``1`` — linear (bi-/tri-linear). + - ``2`` — quadratic B-spline (3 taps per axis). + - ``3`` — local Keys cubic convolution (4 taps; Catmull-Rom, ``a=-0.5``); reproduces input + exactly at integer coordinates. + - ``4`` — quartic B-spline (5 taps per axis). + - ``5`` — quintic B-spline (6 taps per axis). + + Orders ``2``, ``4``, ``5`` evaluate the cardinal B-spline kernel directly on the input samples + (no prefilter); they are low-pass smoothing, not interpolating. See :func:`affine_transform` + for the full discussion. + + The output preserves the input dtype; integer outputs round to nearest and clamp to the dtype + range. Coordinates that map outside the input contribute ``fill_value``. + + Pass a pre-allocated, C-contiguous, writable NumPy array as ``out`` to write the result in + place; it must have shape ``coordinates.shape[1:]`` and dtype equal to ``data.dtype``. + """ + array = np.asarray(data) + if array.ndim not in (2, 3): + raise ValueError(f"data must be 2D or 3D, got ndim={array.ndim}") + + order = int(order) + if order not in (0, 1, 2, 3, 4, 5): + raise ValueError(f"order must be in 0..5, got {order}") + + table = _MAP_COORDINATES_2D_BY_DTYPE if array.ndim == 2 else _MAP_COORDINATES_3D_BY_DTYPE + dtype = array.dtype + try: + run = table[dtype] + except KeyError as error: + supported = ", ".join(str(name) for name in table) + raise TypeError( + f"data must have one of dtypes ({supported}), got dtype={dtype}" + ) from error + + coords, output_shape = _normalize_coordinates(coordinates, array.ndim) + contiguous = np.ascontiguousarray(array) + typed_fill = _normalize_fill_value(fill_value, dtype) + output = _validate_or_allocate_out(out, output_shape, dtype) + run(contiguous, coords, output, order, typed_fill) + return output diff --git a/tests/test_transformation.py b/tests/test_transformation.py index 36853c1..5300220 100644 --- a/tests/test_transformation.py +++ b/tests/test_transformation.py @@ -408,3 +408,124 @@ def test_invalid_inputs_raise(): bic.transformation.affine_transform(data, _matrix(2), bounding_box=(slice(None, None, 2), slice(None))) with pytest.raises(TypeError, match="dtype"): bic.transformation.affine_transform(np.zeros((4, 4), dtype=np.bool_), _matrix(2)) + + +# --------------------------------------------------------------------------- +# map_coordinates +# --------------------------------------------------------------------------- + +_ALL_DTYPES = [ + np.uint8, np.uint16, np.uint32, np.uint64, + np.int8, np.int16, np.int32, np.int64, + np.float32, np.float64, +] + + +def _affine_coords(matrix, shape): + """The (ndim, *shape) coordinate field that reproduces `matrix` as an explicit deformation.""" + ndim = len(shape) + grid = np.indices(shape, dtype=np.float64) + coords = np.empty_like(grid) + for d in range(ndim): + acc = np.full(shape, matrix[d, ndim], dtype=np.float64) + for k in range(ndim): + acc = acc + matrix[d, k] * grid[k] + coords[d] = acc + return coords + + +@pytest.mark.parametrize("order", [0, 1, 2, 3, 4, 5]) +@pytest.mark.parametrize("dtype", _ALL_DTYPES) +def test_map_coordinates_matches_affine_2d(order, dtype): + # A coordinate field built from an affine matrix must reproduce affine_transform exactly: same + # interpolation backend, so the results are bit-identical for every order and dtype. + data = (np.arange(5 * 7) % 17).astype(dtype).reshape(5, 7) + matrix = _matrix(2, translation=[0.5, -1.25]) + coords = _affine_coords(matrix, data.shape) + got = bic.transformation.map_coordinates(data, coords, order=order, fill_value=0) + ref = bic.transformation.affine_transform(data, matrix, order=order, fill_value=0) + np.testing.assert_array_equal(got, ref) + + +@pytest.mark.parametrize("order", [0, 1, 2, 3, 4, 5]) +@pytest.mark.parametrize("dtype", [np.uint8, np.uint16, np.int32, np.float32, np.float64]) +def test_map_coordinates_matches_affine_3d(order, dtype): + data = (np.arange(4 * 5 * 6) % 13).astype(dtype).reshape(4, 5, 6) + matrix = _matrix(3, translation=[0.25, -0.5, 1.0]) + coords = _affine_coords(matrix, data.shape) + got = bic.transformation.map_coordinates(data, coords, order=order, fill_value=0) + ref = bic.transformation.affine_transform(data, matrix, order=order, fill_value=0) + np.testing.assert_array_equal(got, ref) + + +@pytest.mark.parametrize("order", [0, 1]) +@pytest.mark.parametrize("ndim", [2, 3]) +def test_map_coordinates_matches_scipy(order, ndim): + sp = pytest.importorskip("scipy.ndimage") + rng = np.random.default_rng(0) + shape = (9, 11) if ndim == 2 else (6, 7, 8) + data = rng.random(shape).astype(np.float64) + # random source coordinates strictly inside the volume, so boundary handling is not involved. + coords = np.stack([rng.uniform(1.0, s - 2.0, size=shape) for s in shape]) + got = bic.transformation.map_coordinates(data, coords, order=order) + ref = sp.map_coordinates(data, coords, order=order, mode="nearest", prefilter=False) + np.testing.assert_allclose(got, ref, atol=1e-6) + + +@pytest.mark.parametrize("order", [0, 1, 3]) +@pytest.mark.parametrize("shape", [(5, 7), (4, 5, 6)]) +def test_map_coordinates_identity_round_trip(order, shape): + # Sampling at the integer grid reproduces the input exactly for the interpolating orders. + data = np.arange(int(np.prod(shape)), dtype=np.float32).reshape(shape) + coords = np.indices(shape, dtype=np.float64) + got = bic.transformation.map_coordinates(data, coords, order=order, fill_value=-1) + np.testing.assert_array_equal(got, data) + + +@pytest.mark.parametrize("order", [0, 1, 3]) +def test_map_coordinates_out_of_bounds_uses_fill(order): + data = np.arange(16, dtype=np.float32).reshape(4, 4) + coords = np.full((2, 3, 3), -50.0) # every output voxel maps far outside the input + got = bic.transformation.map_coordinates(data, coords, order=order, fill_value=7.0) + np.testing.assert_array_equal(got, np.full((3, 3), 7.0, dtype=np.float32)) + + +def test_map_coordinates_writes_into_out(): + data = np.arange(20, dtype=np.float64).reshape(4, 5) + coords = np.indices(data.shape, dtype=np.float64) + out = np.empty(data.shape, dtype=np.float64) + returned = bic.transformation.map_coordinates(data, coords, order=1, out=out) + assert returned is out + np.testing.assert_array_equal(out, data) + + +def test_map_coordinates_accepts_noncontiguous_coordinates(): + data = np.arange(20, dtype=np.float64).reshape(4, 5) + matrix = _matrix(2, translation=[0.5, 0.5]) + coords = _affine_coords(matrix, data.shape) + noncontig = np.asfortranarray(coords) # coerced to contiguous float64 internally + assert not noncontig.flags.c_contiguous + got = bic.transformation.map_coordinates(data, noncontig, order=1) + ref = bic.transformation.affine_transform(data, matrix, order=1) + np.testing.assert_array_equal(got, ref) + + +def test_map_coordinates_invalid_inputs_raise(): + data = np.zeros((4, 5), dtype=np.float32) + good = np.indices(data.shape).astype(np.float64) + with pytest.raises(ValueError, match="2D or 3D"): + bic.transformation.map_coordinates(np.zeros((4,), dtype=np.float32), np.zeros((1, 4))) + with pytest.raises(ValueError, match="order"): + bic.transformation.map_coordinates(data, good, order=6) + # coordinates must have ndim == data.ndim + 1 + with pytest.raises(ValueError, match="ndim"): + bic.transformation.map_coordinates(data, np.zeros((2, 4), dtype=np.float64)) + # coordinates leading axis must equal data ndim + with pytest.raises(ValueError, match=r"shape\[0\]"): + bic.transformation.map_coordinates(data, np.zeros((3, 4, 5), dtype=np.float64)) + with pytest.raises(TypeError, match="dtype"): + bic.transformation.map_coordinates(np.zeros((4, 5), dtype=np.bool_), good) + with pytest.raises(ValueError, match="shape"): + bic.transformation.map_coordinates(data, good, out=np.empty((3, 3), dtype=np.float32)) + with pytest.raises(TypeError, match="dtype"): + bic.transformation.map_coordinates(data, good, out=np.empty((4, 5), dtype=np.float64))