_interp_c_template: allow interpolating close to the boundary

1 files 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