From a4ce5b09d85b1aede9ff1c82e08010af96ef84c6 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Tue, 14 Feb 2023 18:35:07 +0100 Subject: _interp_c_template: allow interpolating close to the boundary --- _interp_c_template.c | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/_interp_c_template.c b/_interp_c_template.c index e1c5f8a..aa2796d 100644 --- a/_interp_c_template.c +++ b/_interp_c_template.c @@ -3,6 +3,8 @@ #define FUNC2(name, stencil) JOIN2(name, stencil) #define FUNC(name) FUNC2(name, STENCIL) +#define FUZZ 1e-3 + void FUNC(interp2d)(const double *src_start, const double *src_step, const double *src_val, const int *src_stride, const int *src_len, @@ -13,18 +15,27 @@ void FUNC(interp2d)(const double *src_start, const double *src_step, const double x = x_vals[pt]; const double y = y_vals[pt]; - const int idx_x = (x - src_start[0]) / src_step[0] - (STENCIL / 2 - 1); - const int idx_y = (y - src_start[1]) / src_step[1] - (STENCIL / 2 - 1); - - const double src_x = idx_x * src_step[0] + src_start[0]; - const double src_y = idx_y * src_step[1] + src_start[1]; - double fact[2][STENCIL]; + double src_x, src_y, ret; + int idx_x, idx_y; - double ret = 0.0; + idx_x = (x - src_start[0]) / src_step[0] - (STENCIL / 2 - 1); + if (idx_x < 0) + idx_x = 0; + else if (idx_x > src_len[0] - STENCIL) + idx_x = src_len[0] - STENCIL; - if (idx_x < 0 || idx_x > src_len[0] - STENCIL || - idx_y < 0 || idx_y > src_len[1] - STENCIL) { + idx_y = (y - src_start[1]) / src_step[1] - (STENCIL / 2 - 1); + if (idx_y < 0) + idx_y = 0; + else if (idx_y > src_len[1] - STENCIL) + idx_y = src_len[1] - STENCIL; + + src_x = idx_x * src_step[0] + src_start[0]; + src_y = idx_y * src_step[1] + src_start[1]; + if (x < src_x - FUZZ || y < src_y - FUZZ || + x > src_x + STENCIL * src_step[0] + FUZZ || + y > src_y + STENCIL * src_step[1] + FUZZ) { ret = NAN; goto output; } @@ -45,6 +56,7 @@ void FUNC(interp2d)(const double *src_start, const double *src_step, fact[1][i] = fact_y; } + ret = 0.0; for (int j = 0; j < STENCIL; j++) { double val = 0.0; @@ -57,3 +69,8 @@ output: dst[pt] = ret; } } + +#undef FUZZ +#undef FUNC +#undef FUNC2 +#undef JOIN2 -- cgit v1.2.3