Skip to content

Commit 03179ee

Browse files
constantinpapeconstantinpape
authored
Implement map_coordinates functionality (#56)
1 parent 5c27072 commit 03179ee

6 files changed

Lines changed: 507 additions & 48 deletions

File tree

‎include/bioimage_cpp/transformation/affine.hxx‎

Lines changed: 37 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -495,6 +495,40 @@ inline double cubic_3d(
495495
return value;
496496
}
497497

498+
// Dispatch to the per-voxel sampler for the requested interpolation order. Shared by the affine and
499+
// map_coordinates kernels so the interpolation backend lives in exactly one place. `order` must be in
500+
// 0..5 (validated by the public entry points before the sampling loop).
501+
template <class T>
502+
inline double sample_2d(
503+
const T *data, std::ptrdiff_t in_h, std::ptrdiff_t in_w,
504+
double cy, double cx, const int order, double fill
505+
) {
506+
switch (order) {
507+
case 0: return nearest_2d(data, in_h, in_w, cy, cx, fill);
508+
case 1: return linear_2d(data, in_h, in_w, cy, cx, fill);
509+
case 2: return bspline_2d<2>(data, in_h, in_w, cy, cx, fill);
510+
case 3: return cubic_2d(data, in_h, in_w, cy, cx, fill);
511+
case 4: return bspline_2d<4>(data, in_h, in_w, cy, cx, fill);
512+
default: return bspline_2d<5>(data, in_h, in_w, cy, cx, fill); // 5
513+
}
514+
}
515+
516+
template <class T>
517+
inline double sample_3d(
518+
const T *data,
519+
std::ptrdiff_t in_d, std::ptrdiff_t in_h, std::ptrdiff_t in_w,
520+
double cz, double cy, double cx, const int order, double fill
521+
) {
522+
switch (order) {
523+
case 0: return nearest_3d(data, in_d, in_h, in_w, cz, cy, cx, fill);
524+
case 1: return linear_3d(data, in_d, in_h, in_w, cz, cy, cx, fill);
525+
case 2: return bspline_3d<2>(data, in_d, in_h, in_w, cz, cy, cx, fill);
526+
case 3: return cubic_3d(data, in_d, in_h, in_w, cz, cy, cx, fill);
527+
case 4: return bspline_3d<4>(data, in_d, in_h, in_w, cz, cy, cx, fill);
528+
default: return bspline_3d<5>(data, in_d, in_h, in_w, cz, cy, cx, fill); // 5
529+
}
530+
}
531+
498532
} // namespace detail
499533

500534
// ----- 2D entry point -------------------------------------------------------
@@ -554,27 +588,7 @@ void affine_transform_2d(
554588
double cy = row_y;
555589
double cx = row_x;
556590
for (std::ptrdiff_t j = 0; j < out_w; ++j) {
557-
double value;
558-
switch (order) {
559-
case 0:
560-
value = detail::nearest_2d(in_data, in_h, in_w, cy, cx, fill);
561-
break;
562-
case 1:
563-
value = detail::linear_2d(in_data, in_h, in_w, cy, cx, fill);
564-
break;
565-
case 2:
566-
value = detail::bspline_2d<2>(in_data, in_h, in_w, cy, cx, fill);
567-
break;
568-
case 3:
569-
value = detail::cubic_2d(in_data, in_h, in_w, cy, cx, fill);
570-
break;
571-
case 4:
572-
value = detail::bspline_2d<4>(in_data, in_h, in_w, cy, cx, fill);
573-
break;
574-
default: // 5
575-
value = detail::bspline_2d<5>(in_data, in_h, in_w, cy, cx, fill);
576-
break;
577-
}
591+
const double value = detail::sample_2d(in_data, in_h, in_w, cy, cx, order, fill);
578592
*out_ptr++ = detail::to_output<T>(value);
579593
cy += m01;
580594
cx += m11;
@@ -655,33 +669,8 @@ void affine_transform_3d(
655669
double cy = row_y;
656670
double cx = row_x;
657671
for (std::ptrdiff_t j = 0; j < out_w; ++j) {
658-
double value;
659-
switch (order) {
660-
case 0:
661-
value = detail::nearest_3d(in_data, in_d, in_h, in_w,
662-
cz, cy, cx, fill);
663-
break;
664-
case 1:
665-
value = detail::linear_3d(in_data, in_d, in_h, in_w,
666-
cz, cy, cx, fill);
667-
break;
668-
case 2:
669-
value = detail::bspline_3d<2>(in_data, in_d, in_h, in_w,
670-
cz, cy, cx, fill);
671-
break;
672-
case 3:
673-
value = detail::cubic_3d(in_data, in_d, in_h, in_w,
674-
cz, cy, cx, fill);
675-
break;
676-
case 4:
677-
value = detail::bspline_3d<4>(in_data, in_d, in_h, in_w,
678-
cz, cy, cx, fill);
679-
break;
680-
default: // 5
681-
value = detail::bspline_3d<5>(in_data, in_d, in_h, in_w,
682-
cz, cy, cx, fill);
683-
break;
684-
}
672+
const double value = detail::sample_3d(in_data, in_d, in_h, in_w,
673+
cz, cy, cx, order, fill);
685674
*out_ptr++ = detail::to_output<T>(value);
686675
cz += m02;
687676
cy += m12;

‎include/bioimage_cpp/transformation/coordinate.hxx‎

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
#pragma once
2+
3+
// Coordinate-based resampling (the analogue of ``scipy.ndimage.map_coordinates``). For every output
4+
// voxel the source coordinate to sample is read from an explicit ``coordinates`` array, instead of
5+
// being computed from an affine matrix. The per-voxel interpolation reuses the affine samplers
6+
// (``detail::sample_2d`` / ``detail::sample_3d``) so the interpolation backend lives in one place.
7+
//
8+
// This kernel is pure and in-memory (NumPy in, NumPy out); reading the source data and producing the
9+
// coordinate (deformation) field are the caller's responsibility.
10+
11+
#include "bioimage_cpp/transformation/affine.hxx"
12+
13+
namespace bioimage_cpp::transformation {
14+
15+
namespace detail {
16+
17+
// Validate that ``coordinates`` is a ``(D, *output_shape)`` field, i.e. it carries one source
18+
// coordinate per output voxel along its leading axis.
19+
template <std::size_t D, class T>
20+
void require_coordinates(const ConstArrayView<double> &coordinates, const ArrayView<T> &output) {
21+
if (coordinates.ndim() != static_cast<std::ptrdiff_t>(D + 1)) {
22+
throw std::invalid_argument(
23+
"coordinates must have ndim=" + std::to_string(D + 1) +
24+
", got ndim=" + std::to_string(coordinates.ndim())
25+
);
26+
}
27+
if (coordinates.shape[0] != static_cast<std::ptrdiff_t>(D)) {
28+
throw std::invalid_argument(
29+
"coordinates.shape[0] must equal the data dimension " + std::to_string(D) +
30+
", got " + std::to_string(coordinates.shape[0])
31+
);
32+
}
33+
for (std::size_t axis = 0; axis < D; ++axis) {
34+
if (coordinates.shape[axis + 1] != output.shape[axis]) {
35+
throw std::invalid_argument(
36+
"coordinates spatial shape (coordinates.shape[1:]) must match the output shape"
37+
);
38+
}
39+
}
40+
}
41+
42+
} // namespace detail
43+
44+
// ----- 2D entry point -------------------------------------------------------
45+
46+
template <class T>
47+
void map_coordinates_2d(
48+
const ConstArrayView<T> &input,
49+
ArrayView<T> &output,
50+
const ConstArrayView<double> &coordinates,
51+
const int order,
52+
const T fill_value
53+
) {
54+
detail::require_views<2, T>(input, output);
55+
detail::require_coordinates<2, T>(coordinates, output);
56+
if (order < 0 || order > 5) {
57+
throw std::invalid_argument(
58+
"order must be in 0..5, got " + std::to_string(order)
59+
);
60+
}
61+
62+
const auto out_h = output.shape[0];
63+
const auto out_w = output.shape[1];
64+
if (out_h == 0 || out_w == 0) return;
65+
66+
const auto in_h = input.shape[0];
67+
const auto in_w = input.shape[1];
68+
69+
// coordinates is C-contiguous with shape (2, out_h, out_w), so axis-d coordinates form a
70+
// contiguous block of n_out doubles starting at d * n_out.
71+
const std::ptrdiff_t n_out = out_h * out_w;
72+
const double *cy = coordinates.data;
73+
const double *cx = coordinates.data + n_out;
74+
75+
const double fill = static_cast<double>(fill_value);
76+
const T *in_data = input.data;
77+
T *out_ptr = output.data;
78+
79+
for (std::ptrdiff_t p = 0; p < n_out; ++p) {
80+
const double value = detail::sample_2d(in_data, in_h, in_w, cy[p], cx[p], order, fill);
81+
out_ptr[p] = detail::to_output<T>(value);
82+
}
83+
}
84+
85+
// ----- 3D entry point -------------------------------------------------------
86+
87+
template <class T>
88+
void map_coordinates_3d(
89+
const ConstArrayView<T> &input,
90+
ArrayView<T> &output,
91+
const ConstArrayView<double> &coordinates,
92+
const int order,
93+
const T fill_value
94+
) {
95+
detail::require_views<3, T>(input, output);
96+
detail::require_coordinates<3, T>(coordinates, output);
97+
if (order < 0 || order > 5) {
98+
throw std::invalid_argument(
99+
"order must be in 0..5, got " + std::to_string(order)
100+
);
101+
}
102+
103+
const auto out_d = output.shape[0];
104+
const auto out_h = output.shape[1];
105+
const auto out_w = output.shape[2];
106+
if (out_d == 0 || out_h == 0 || out_w == 0) return;
107+
108+
const auto in_d = input.shape[0];
109+
const auto in_h = input.shape[1];
110+
const auto in_w = input.shape[2];
111+
112+
// coordinates is C-contiguous with shape (3, out_d, out_h, out_w).
113+
const std::ptrdiff_t n_out = out_d * out_h * out_w;
114+
const double *cz = coordinates.data;
115+
const double *cy = coordinates.data + n_out;
116+
const double *cx = coordinates.data + 2 * n_out;
117+
118+
const double fill = static_cast<double>(fill_value);
119+
const T *in_data = input.data;
120+
T *out_ptr = output.data;
121+
122+
for (std::ptrdiff_t p = 0; p < n_out; ++p) {
123+
const double value = detail::sample_3d(in_data, in_d, in_h, in_w,
124+
cz[p], cy[p], cx[p], order, fill);
125+
out_ptr[p] = detail::to_output<T>(value);
126+
}
127+
}
128+
129+
} // namespace bioimage_cpp::transformation

‎src/bindings/transformation.cxx‎

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include "bioimage_cpp/array_view.hxx"
44
#include "bioimage_cpp/detail/grid.hxx"
55
#include "bioimage_cpp/transformation/affine.hxx"
6+
#include "bioimage_cpp/transformation/coordinate.hxx"
67

78
#include <nanobind/ndarray.h>
89

@@ -113,6 +114,86 @@ void bind_affine_for_dtype(nb::module_ &m, const char *name_2d, const char *name
113114
);
114115
}
115116

117+
template <std::size_t D, class T>
118+
OutputArray<T> map_coordinates_t(
119+
ConstArray<T> input,
120+
MatrixArray coordinates,
121+
OutputArray<T> output,
122+
const int order,
123+
const T fill_value
124+
) {
125+
if (input.ndim() != D) {
126+
throw std::invalid_argument(
127+
"input must have ndim=" + std::to_string(D) +
128+
", got ndim=" + std::to_string(input.ndim())
129+
);
130+
}
131+
if (output.ndim() != D) {
132+
throw std::invalid_argument(
133+
"output must have ndim=" + std::to_string(D) +
134+
", got ndim=" + std::to_string(output.ndim())
135+
);
136+
}
137+
if (coordinates.ndim() != D + 1) {
138+
throw std::invalid_argument(
139+
"coordinates must have ndim=" + std::to_string(D + 1) +
140+
", got ndim=" + std::to_string(coordinates.ndim())
141+
);
142+
}
143+
144+
const auto input_shape = shape_of(input);
145+
const auto input_strides = detail::c_order_strides(input_shape);
146+
const auto coordinates_shape = shape_of(coordinates);
147+
const auto coordinates_strides = detail::c_order_strides(coordinates_shape);
148+
const auto output_shape = shape_of(output);
149+
const auto output_strides = detail::c_order_strides(output_shape);
150+
151+
ConstArrayView<T> input_view{input.data(), input_shape, input_strides};
152+
ArrayView<T> output_view{output.data(), output_shape, output_strides};
153+
ConstArrayView<double> coordinates_view{
154+
coordinates.data(), coordinates_shape, coordinates_strides
155+
};
156+
157+
{
158+
nb::gil_scoped_release release;
159+
if constexpr (D == 2) {
160+
transformation::map_coordinates_2d<T>(
161+
input_view, output_view, coordinates_view, order, fill_value
162+
);
163+
} else {
164+
transformation::map_coordinates_3d<T>(
165+
input_view, output_view, coordinates_view, order, fill_value
166+
);
167+
}
168+
}
169+
170+
return output;
171+
}
172+
173+
template <class T>
174+
void bind_map_coordinates_for_dtype(nb::module_ &m, const char *name_2d, const char *name_3d) {
175+
m.def(
176+
name_2d,
177+
&map_coordinates_t<2, T>,
178+
nb::arg("input"),
179+
nb::arg("coordinates"),
180+
nb::arg("output"),
181+
nb::arg("order"),
182+
nb::arg("fill_value"),
183+
"Apply a 2D coordinate-based resampling into a pre-allocated NumPy array."
184+
);
185+
m.def(
186+
name_3d,
187+
&map_coordinates_t<3, T>,
188+
nb::arg("input"),
189+
nb::arg("coordinates"),
190+
nb::arg("output"),
191+
nb::arg("order"),
192+
nb::arg("fill_value"),
193+
"Apply a 3D coordinate-based resampling into a pre-allocated NumPy array."
194+
);
195+
}
196+
116197
} // namespace
117198

118199
void bind_transformation(nb::module_ &m) {
@@ -146,6 +227,37 @@ void bind_transformation(nb::module_ &m) {
146227
bind_affine_for_dtype<double>(
147228
m, "_affine_transform_2d_float64", "_affine_transform_3d_float64"
148229
);
230+
231+
bind_map_coordinates_for_dtype<std::uint8_t>(
232+
m, "_map_coordinates_2d_uint8", "_map_coordinates_3d_uint8"
233+
);
234+
bind_map_coordinates_for_dtype<std::uint16_t>(
235+
m, "_map_coordinates_2d_uint16", "_map_coordinates_3d_uint16"
236+
);
237+
bind_map_coordinates_for_dtype<std::uint32_t>(
238+
m, "_map_coordinates_2d_uint32", "_map_coordinates_3d_uint32"
239+
);
240+
bind_map_coordinates_for_dtype<std::uint64_t>(
241+
m, "_map_coordinates_2d_uint64", "_map_coordinates_3d_uint64"
242+
);
243+
bind_map_coordinates_for_dtype<std::int8_t>(
244+
m, "_map_coordinates_2d_int8", "_map_coordinates_3d_int8"
245+
);
246+
bind_map_coordinates_for_dtype<std::int16_t>(
247+
m, "_map_coordinates_2d_int16", "_map_coordinates_3d_int16"
248+
);
249+
bind_map_coordinates_for_dtype<std::int32_t>(
250+
m, "_map_coordinates_2d_int32", "_map_coordinates_3d_int32"
251+
);
252+
bind_map_coordinates_for_dtype<std::int64_t>(
253+
m, "_map_coordinates_2d_int64", "_map_coordinates_3d_int64"
254+
);
255+
bind_map_coordinates_for_dtype<float>(
256+
m, "_map_coordinates_2d_float32", "_map_coordinates_3d_float32"
257+
);
258+
bind_map_coordinates_for_dtype<double>(
259+
m, "_map_coordinates_2d_float64", "_map_coordinates_3d_float64"
260+
);
149261
}
150262

151263
} // namespace bioimage_cpp::bindings

‎src/bioimage_cpp/transformation/__init__.py‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,11 +3,13 @@
33
from ._transformation import (
44
affine_transform,
55
compute_anti_aliasing_sigma,
6+
map_coordinates,
67
resample,
78
)
89

910
__all__ = [
1011
"affine_transform",
1112
"compute_anti_aliasing_sigma",
13+
"map_coordinates",
1214
"resample",
1315
]

0 commit comments

Comments
 (0)