diff options
Diffstat (limited to 'interp.py')
| -rw-r--r-- | interp.py | 64 |
1 files changed, 64 insertions, 0 deletions
@@ -1,3 +1,5 @@ +import ctypes + import numpy as np def interp1d(src_start, src_step, src_val, dst_coords, stencil): @@ -50,3 +52,65 @@ def interp2d(src_start, src_step, src_val, dst_coords, stencil): ret += fact[0][i] * val return ret + +_doubleptr = ctypes.POINTER(ctypes.c_double) + +class Interp2D_C: + _lib = None + _interp_func = None + + _src_start_c = None + _src_step_c = None + _src_val_c = None + _src_stride_c = None + _src_len_c = None + + _x = None + _y = None + _ret = None + _x_c = None + _y_c = None + _ret_c = None + + def __init__(self, src_start, src_step, src_val, stencil): + self._lib = ctypes.CDLL('lib_interp_c.so') + + + func = getattr(self._lib, 'interp2d_%d' % stencil) + func.argtypes = [_doubleptr, _doubleptr, _doubleptr, + ctypes.POINTER(ctypes.c_int), ctypes.POINTER(ctypes.c_int), + ctypes.c_int, _doubleptr, _doubleptr, _doubleptr] + + self._interp_func = func + + self._src_start_c = (ctypes.c_double * 2)(src_start[0], src_start[1]) + self._src_step_c = (ctypes.c_double * 2)(src_step[0], src_step[1]) + self._src_val_c = ctypes.cast(np.ctypeslib.as_ctypes(src_val), _doubleptr) + self._src_stride_c = (ctypes.c_int * 2)(src_val.strides[0] // 8, src_val.strides[1] // 8) + self._src_len_c = (ctypes.c_int * 2)(src_val.shape[0], src_val.shape[1]) + + self._setup_arrays((1,)) + + def _setup_arrays(self, shape): + self._x = np.empty(shape) + self._y = np.empty(shape) + self._ret = np.empty(shape) + + self._x_c = ctypes.cast(np.ctypeslib.as_ctypes(self._x), _doubleptr) + self._y_c = ctypes.cast(np.ctypeslib.as_ctypes(self._y), _doubleptr) + self._ret_c = ctypes.cast(np.ctypeslib.as_ctypes(self._ret), _doubleptr) + + def interp(self, x, y): + if x.shape != self._x.shape: + self._setup_arrays(x.shape) + + np.copyto(self._x, x) + np.copyto(self._y, y) + self._interp_func(self._src_start_c, self._src_step_c, self._src_val_c, + self._src_stride_c, self._src_len_c, + self._x.shape[0], + self._x_c, self._y_c, self._ret_c) + return self._ret.copy() + + def __call__(self, x, y): + return self.interp(x, y) |