Skip to content

Commit fdea688

Browse files
mraspaudmraspaud
committed
Refactor gradient search cython
1 parent 03cec67 commit fdea688

File tree

1 file changed

+51
-44
lines changed

1 file changed

+51
-44
lines changed

‎pyresample/gradient/_gradient_search.pyx

+51-44
Original file line numberDiff line numberDiff line change
@@ -27,22 +27,19 @@ cimport cython
2727
cimport numpy as np
2828
from libc.math cimport fabs, isinf
2929

30-
ctypedef np.double_t DTYPE_t
31-
3230
ctypedef fused data_type:
3331
double
3432
float
3533

36-
ctypedef double float_index
37-
f_index_dtype = float
38-
ctypedef int i_index_type
39-
34+
ctypedef fused float_index:
35+
double
36+
float
4037

4138
np.import_array()
4239

4340
@cython.boundscheck(False)
4441
@cython.wraparound(False)
45-
cdef inline void nn(const data_type[:, :, :] data, i_index_type l0, int p0, double dl, double dp, int lmax, int pmax, data_type[:] res) noexcept nogil:
42+
cdef inline void nn(const data_type[:, :, :] data, int l0, int p0, float_index dl, float_index dp, int lmax, int pmax, data_type[:] res) noexcept nogil:
4643
cdef int nnl, nnp
4744
cdef size_t z_size = res.shape[0]
4845
cdef size_t i
@@ -62,9 +59,9 @@ cdef inline void nn(const data_type[:, :, :] data, i_index_type l0, int p0, doub
6259

6360
@cython.boundscheck(False)
6461
@cython.wraparound(False)
65-
cdef inline void bil(const data_type[:, :, :] data, int l0, int p0, double dl, double dp, int lmax, int pmax, data_type[:] res) noexcept nogil:
62+
cdef inline void bil(const data_type[:, :, :] data, int l0, int p0, float_index dl, float_index dp, int lmax, int pmax, data_type[:] res) noexcept nogil:
6663
cdef int l_a, l_b, p_a, p_b
67-
cdef double w_l, w_p
64+
cdef float_index w_l, w_p
6865
cdef size_t z_size = res.shape[0]
6966
cdef size_t i
7067
if dl < 0:
@@ -92,26 +89,28 @@ cdef inline void bil(const data_type[:, :, :] data, int l0, int p0, double dl, d
9289

9390
@cython.boundscheck(False)
9491
@cython.wraparound(False)
95-
cdef inline void indices_xy(const data_type[:, :, :] data, int l0, int p0, double dl, double dp, int lmax, int pmax, data_type[:] res) noexcept nogil:
92+
cdef inline void indices_xy(const data_type[:, :, :] data, int l0, int p0, float_index dl, float_index dp, int lmax, int pmax, data_type[:] res) noexcept nogil:
9693
cdef int nnl, nnp
9794
cdef size_t z_size = res.shape[0]
9895
cdef size_t i
9996
res[1] = dl + l0
10097
res[0] = dp + p0
10198

102-
ctypedef void (*FN)(const data_type[:, :, :] data, int l0, int p0, double dl, double dp, int lmax, int pmax, data_type[:] res) noexcept nogil
99+
100+
ctypedef void (*FN)(const data_type[:, :, :] data, int l0, int p0, float_index dl, float_index dp, int lmax, int pmax, data_type[:] res) noexcept nogil
101+
103102

104103
@cython.boundscheck(False)
105104
@cython.wraparound(False)
106105
cpdef one_step_gradient_search(const data_type[:, :, :] data,
107-
DTYPE_t [:, :] src_x,
108-
DTYPE_t [:, :] src_y,
109-
DTYPE_t [:, :] xl,
110-
DTYPE_t [:, :] xp,
111-
DTYPE_t [:, :] yl,
112-
DTYPE_t [:, :] yp,
113-
DTYPE_t [:, :] dst_x,
114-
DTYPE_t [:, :] dst_y,
106+
float_index [:, :] src_x,
107+
float_index [:, :] src_y,
108+
float_index [:, :] xl,
109+
float_index [:, :] xp,
110+
float_index [:, :] yl,
111+
float_index [:, :] yp,
112+
float_index [:, :] dst_x,
113+
float_index [:, :] dst_y,
115114
str method='bilinear'):
116115
"""Gradient search, simple case variant."""
117116
cdef FN fun
@@ -144,18 +143,19 @@ cpdef one_step_gradient_search(const data_type[:, :, :] data,
144143
# return the output image
145144
return image
146145

146+
147147
@cython.boundscheck(False)
148148
@cython.wraparound(False)
149149
@cython.cdivision(True)
150150
cdef void one_step_gradient_search_no_gil(const data_type[:, :, :] data,
151-
const DTYPE_t[:, :] src_x,
152-
const DTYPE_t[:, :] src_y,
153-
const DTYPE_t[:, :] xl,
154-
const DTYPE_t[:, :] xp,
155-
const DTYPE_t[:, :] yl,
156-
const DTYPE_t[:, :] yp,
157-
const DTYPE_t[:, :] dst_x,
158-
const DTYPE_t[:, :] dst_y,
151+
const float_index[:, :] src_x,
152+
const float_index[:, :] src_y,
153+
const float_index[:, :] xl,
154+
const float_index[:, :] xp,
155+
const float_index[:, :] yl,
156+
const float_index[:, :] yp,
157+
const float_index[:, :] dst_x,
158+
const float_index[:, :] dst_y,
159159
const size_t x_size,
160160
const size_t y_size,
161161
FN fun,
@@ -173,7 +173,7 @@ cdef void one_step_gradient_search_no_gil(const data_type[:, :, :] data,
173173
# intermediate variables:
174174
cdef int l_a, l_b, p_a, p_b
175175
cdef size_t i, j, elt
176-
cdef double dx, dy, d, dl, dp
176+
cdef float_index dx, dy, d, dl, dp
177177
cdef int col_step = -1
178178
# number of iterations
179179
cdef int cnt = 0
@@ -230,16 +230,17 @@ cdef void one_step_gradient_search_no_gil(const data_type[:, :, :] data,
230230
p0 = int(p0 + dp)
231231
j += col_step
232232

233+
233234
@cython.boundscheck(False)
234235
@cython.wraparound(False)
235-
cpdef one_step_gradient_indices(DTYPE_t [:, :] src_x,
236-
DTYPE_t [:, :] src_y,
237-
DTYPE_t [:, :] xl,
238-
DTYPE_t [:, :] xp,
239-
DTYPE_t [:, :] yl,
240-
DTYPE_t [:, :] yp,
241-
DTYPE_t [:, :] dst_x,
242-
DTYPE_t [:, :] dst_y):
236+
cpdef one_step_gradient_indices(float_index [:, :] src_x,
237+
float_index [:, :] src_y,
238+
float_index [:, :] xl,
239+
float_index [:, :] xp,
240+
float_index [:, :] yl,
241+
float_index [:, :] yp,
242+
float_index [:, :] dst_x,
243+
float_index [:, :] dst_y):
243244
"""Gradient search, simple case variant, returning float indices.
244245
245246
This is appropriate for monotonous gradients only, i.e. not modis or viirs in satellite projection.
@@ -250,18 +251,24 @@ cpdef one_step_gradient_indices(DTYPE_t [:, :] src_x,
250251
cdef size_t y_size = dst_y.shape[0]
251252
cdef size_t x_size = dst_x.shape[1]
252253

254+
255+
if float_index is double:
256+
dtype = np.float64
257+
else:
258+
dtype = np.float32
259+
253260
# output indices arrays --> needs to be (lines, pixels) --> y,x
254-
indices = np.full([2, y_size, x_size], np.nan, dtype=f_index_dtype)
261+
indices = np.full([2, y_size, x_size], np.nan, dtype=dtype)
255262
cdef float_index [:, :, :] indices_view_result = indices
256263

257264
# fake_data is not going to be used anyway as we just fill in the indices
258-
cdef float_index [:, :, :] fake_data = np.full([1, 1, 1], np.nan, dtype=f_index_dtype)
265+
cdef float_index [:, :, :] fake_data = np.full([1, 1, 1], np.nan, dtype=dtype)
259266

260267
with nogil:
261-
one_step_gradient_search_no_gil[float_index](fake_data,
262-
src_x, src_y,
263-
xl, xp, yl, yp,
264-
dst_x, dst_y,
265-
x_size, y_size,
266-
indices_xy, indices_view_result)
268+
one_step_gradient_search_no_gil[float_index, float_index](fake_data,
269+
src_x, src_y,
270+
xl, xp, yl, yp,
271+
dst_x, dst_y,
272+
x_size, y_size,
273+
indices_xy, indices_view_result)
267274
return indices

0 commit comments

Comments
 (0)