|
1
|
|
|
import numpy as np |
|
|
|
|
|
|
2
|
|
|
|
|
3
|
|
|
from numba import jit |
|
4
|
|
|
|
|
5
|
|
|
|
|
6
|
|
|
class FastBilinearInterpolation(object): |
|
7
|
|
|
""" |
|
8
|
|
|
A super fast bilinar interpolation implementation which exploits the fact that we are always interpolating in the |
|
9
|
|
|
same grid. For example, if we always go from the same flat sky projection to the same Healpix map, we can precompute |
|
10
|
|
|
the weights for the interpolation and then apply them to any new data instead of recomputing them every time. |
|
11
|
|
|
""" |
|
12
|
|
|
|
|
13
|
|
|
def __init__(self, input_shape, new_points): |
|
14
|
|
|
|
|
15
|
|
|
self._gridx = np.arange(input_shape[0]) |
|
16
|
|
|
self._gridy = np.arange(input_shape[1]) |
|
17
|
|
|
|
|
18
|
|
|
self._x_bounds = [self._gridx.min(), self._gridx.max()] |
|
19
|
|
|
self._y_bounds = [self._gridy.min(), self._gridy.max()] |
|
20
|
|
|
|
|
21
|
|
|
self._data_shape = (self._gridx.shape[0], self._gridy.shape[0]) |
|
22
|
|
|
|
|
23
|
|
|
self._bs, self._flat_points = self.compute_coefficients(new_points) |
|
24
|
|
|
|
|
25
|
|
|
@staticmethod |
|
26
|
|
|
def _find_bounding_box(xaxis, yaxis, xs, ys): |
|
|
|
|
|
|
27
|
|
|
# Find lower left corner of bounding box |
|
28
|
|
|
xidx = np.searchsorted(xaxis, xs, 'left') - 1 |
|
29
|
|
|
yidx = np.searchsorted(yaxis, ys, 'left') - 1 |
|
30
|
|
|
|
|
31
|
|
|
lower_left_x = xaxis[xidx] |
|
32
|
|
|
lower_left_y = yaxis[yidx] |
|
33
|
|
|
|
|
34
|
|
|
upper_left_x = xaxis[xidx] |
|
35
|
|
|
upper_left_y = yaxis[yidx + 1] |
|
36
|
|
|
|
|
37
|
|
|
upper_right_x = xaxis[xidx + 1] |
|
38
|
|
|
upper_right_y = yaxis[yidx + 1] |
|
39
|
|
|
|
|
40
|
|
|
lower_right_x = xaxis[xidx + 1] |
|
41
|
|
|
lower_right_y = yaxis[yidx] |
|
42
|
|
|
|
|
43
|
|
|
return (lower_left_x, lower_left_y, |
|
44
|
|
|
upper_left_x, upper_left_y, |
|
45
|
|
|
upper_right_x, upper_right_y, |
|
46
|
|
|
lower_right_x, lower_right_y) |
|
47
|
|
|
|
|
48
|
|
|
def compute_coefficients(self, p): |
|
|
|
|
|
|
49
|
|
|
|
|
50
|
|
|
xx = p[0] |
|
|
|
|
|
|
51
|
|
|
yy = p[1] |
|
|
|
|
|
|
52
|
|
|
|
|
53
|
|
|
# Find bounding boxes |
|
54
|
|
|
# We need a situation like this |
|
55
|
|
|
# x1 x2 |
|
56
|
|
|
# y1 q11 q21 |
|
57
|
|
|
# |
|
58
|
|
|
# y2 q12 q22 |
|
59
|
|
|
|
|
60
|
|
|
x1, y2, xx1, y1, x2, yy1, xx2, yy2 = self._find_bounding_box(self._gridx, self._gridy, xx, yy) |
|
|
|
|
|
|
61
|
|
|
|
|
62
|
|
|
bs = np.zeros((xx.shape[0], 4), np.float64) |
|
|
|
|
|
|
63
|
|
|
|
|
64
|
|
|
bs[:, 0] = (x2 - xx) * (y2 - yy) / (x2 - x1) * (y2 - y1) |
|
65
|
|
|
bs[:, 1] = (xx - x1) * (y2 - yy) / (x2 - x1) * (y2 - y1) |
|
66
|
|
|
bs[:, 2] = (x2 - xx) * (yy - y1) / (x2 - x1) * (y2 - y1) |
|
67
|
|
|
bs[:, 3] = (xx - x1) * (yy - y1) / (x2 - x1) * (y2 - y1) |
|
68
|
|
|
|
|
69
|
|
|
# Get the flat indexing for all the corners of the bounding boxes |
|
70
|
|
|
flat_upper_left = np.ravel_multi_index((x1, y1), self._data_shape) |
|
71
|
|
|
flat_upper_right = np.ravel_multi_index((x2, y1), self._data_shape) |
|
72
|
|
|
flat_lower_left = np.ravel_multi_index((x1, y2), self._data_shape) |
|
73
|
|
|
flat_lower_right = np.ravel_multi_index((x2, y2), self._data_shape) |
|
74
|
|
|
|
|
75
|
|
|
# Stack them so that they are in the right order, i.e.: |
|
76
|
|
|
# ul1, ur1, ll1, lr1, ul2, ur2, ll2, lr2 ... uln, urn, lln, lrn |
|
77
|
|
|
|
|
78
|
|
|
flat_points = np.vstack([flat_upper_left, |
|
79
|
|
|
flat_upper_right, |
|
80
|
|
|
flat_lower_left, |
|
81
|
|
|
flat_lower_right]).T.flatten() |
|
82
|
|
|
|
|
83
|
|
|
return bs, flat_points.astype(np.int64) |
|
84
|
|
|
|
|
85
|
|
|
def __call__(self, data): |
|
86
|
|
|
|
|
87
|
|
|
res = _apply_bilinar_interpolation(self._bs, self._flat_points, data) |
|
88
|
|
|
|
|
89
|
|
|
return res |
|
90
|
|
|
|
|
91
|
|
|
|
|
92
|
|
|
@jit("float64[:](float64[:, :], int64[:], float64[:, :])", nopython=True) |
|
93
|
|
|
def _apply_bilinar_interpolation(bs, points, data): # pragma: no cover |
|
|
|
|
|
|
94
|
|
|
|
|
95
|
|
|
vs = data.ravel()[points] |
|
|
|
|
|
|
96
|
|
|
|
|
97
|
|
|
return np.sum(bs * vs.reshape(bs.shape), axis=1).flatten() |
|
|
|
|
|
giacomov /
hawc_hal
Loading history...
The coding style of this project requires that you add a docstring to this code element. Below, you find an example for methods:
If you would like to know more about docstrings, we recommend to read PEP-257: Docstring Conventions.