1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
import ctypes
import numpy as np
def interp1d(src_start, src_step, src_val, dst_coords, stencil):
idx_src = ((dst_coords - src_start) / src_step - (stencil / 2 - 1)).astype(np.int)
src_coord = np.linspace(src_start, src_start + src_step * (src_val.shape[0] - 1),
src_val.shape[0], dtype = src_val.dtype)
fact = np.zeros((stencil, ) + idx_src.shape, dtype = src_val.dtype) + 1.0
for i in range(stencil):
for j in range(stencil):
if i == j:
continue
fact[i] *= (dst_coords - src_coord[idx_src + j]) / (src_coord[idx_src + i] - src_coord[idx_src + j])
ret = np.zeros_like(dst_coords, dtype = src_val.dtype)
for i in range(stencil):
ret += fact[i] * src_val[idx_src + i]
return ret
def interp2d(src_start, src_step, src_val, dst_coords, stencil):
ret = np.zeros_like(dst_coords[0], dtype = src_val.dtype)
idx_src = []
fact = []
for dim in range(len(src_start)):
src_coord = np.linspace(src_start[dim], src_start[dim] + src_step[dim] * (src_val.shape[dim] - 1),
src_val.shape[dim], dtype = src_val.dtype)
idx = ((dst_coords[dim] - src_start[dim]) / src_step[dim] - (stencil / 2 - 1)).astype(np.int)
f = np.zeros((stencil, ) + idx.shape, dtype = src_val.dtype) + 1.0
for i in range(stencil):
for j in range(stencil):
if i == j:
continue
f[i] *= (dst_coords[dim] - src_coord[idx + j]) / (src_coord[idx + i] - src_coord[idx + j])
idx_src.append(idx)
fact.append(f)
for i in range(stencil):
val = np.zeros_like(ret)
for j in range(stencil):
val += fact[1][j] * src_val[idx_src[0] + i, idx_src[1] + j]
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):
try:
x.shape[0]
except (AttributeError, IndexError):
x = np.array([x])
try:
y.shape[0]
except (AttributeError, IndexError):
y = np.array([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)
|