Skip to content

Commit

Permalink
Refactor gradient search cython
Browse files Browse the repository at this point in the history
  • Loading branch information
@mraspaud
mraspaud committed Sep 3, 2024
1 parent 03cec67 commit fdea688
Showing 1 changed file with 51 additions and 44 deletions.
95 changes: 51 additions & 44 deletions pyresample/gradient/_gradient_search.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,22 +27,19 @@ cimport cython
cimport numpy as np
from libc.math cimport fabs, isinf

ctypedef np.double_t DTYPE_t

ctypedef fused data_type:
double
float

ctypedef double float_index
f_index_dtype = float
ctypedef int i_index_type

ctypedef fused float_index:
double
float

np.import_array()

@cython.boundscheck(False)
@cython.wraparound(False)
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:
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:
cdef int nnl, nnp
cdef size_t z_size = res.shape[0]
cdef size_t i
Expand All @@ -62,9 +59,9 @@ cdef inline void nn(const data_type[:, :, :] data, i_index_type l0, int p0, doub

@cython.boundscheck(False)
@cython.wraparound(False)
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:
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:
cdef int l_a, l_b, p_a, p_b
cdef double w_l, w_p
cdef float_index w_l, w_p
cdef size_t z_size = res.shape[0]
cdef size_t i
if dl < 0:
Expand Down Expand Up @@ -92,26 +89,28 @@ cdef inline void bil(const data_type[:, :, :] data, int l0, int p0, double dl, d

@cython.boundscheck(False)
@cython.wraparound(False)
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:
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:
cdef int nnl, nnp
cdef size_t z_size = res.shape[0]
cdef size_t i
res[1] = dl + l0
res[0] = dp + p0

ctypedef void (*FN)(const data_type[:, :, :] data, int l0, int p0, double dl, double dp, int lmax, int pmax, data_type[:] res) noexcept nogil

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


@cython.boundscheck(False)
@cython.wraparound(False)
cpdef one_step_gradient_search(const data_type[:, :, :] data,
DTYPE_t [:, :] src_x,
DTYPE_t [:, :] src_y,
DTYPE_t [:, :] xl,
DTYPE_t [:, :] xp,
DTYPE_t [:, :] yl,
DTYPE_t [:, :] yp,
DTYPE_t [:, :] dst_x,
DTYPE_t [:, :] dst_y,
float_index [:, :] src_x,
float_index [:, :] src_y,
float_index [:, :] xl,
float_index [:, :] xp,
float_index [:, :] yl,
float_index [:, :] yp,
float_index [:, :] dst_x,
float_index [:, :] dst_y,
str method='bilinear'):
"""Gradient search, simple case variant."""
cdef FN fun
Expand Down Expand Up @@ -144,18 +143,19 @@ cpdef one_step_gradient_search(const data_type[:, :, :] data,
# return the output image
return image


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef void one_step_gradient_search_no_gil(const data_type[:, :, :] data,
const DTYPE_t[:, :] src_x,
const DTYPE_t[:, :] src_y,
const DTYPE_t[:, :] xl,
const DTYPE_t[:, :] xp,
const DTYPE_t[:, :] yl,
const DTYPE_t[:, :] yp,
const DTYPE_t[:, :] dst_x,
const DTYPE_t[:, :] dst_y,
const float_index[:, :] src_x,
const float_index[:, :] src_y,
const float_index[:, :] xl,
const float_index[:, :] xp,
const float_index[:, :] yl,
const float_index[:, :] yp,
const float_index[:, :] dst_x,
const float_index[:, :] dst_y,
const size_t x_size,
const size_t y_size,
FN fun,
Expand All @@ -173,7 +173,7 @@ cdef void one_step_gradient_search_no_gil(const data_type[:, :, :] data,
# intermediate variables:
cdef int l_a, l_b, p_a, p_b
cdef size_t i, j, elt
cdef double dx, dy, d, dl, dp
cdef float_index dx, dy, d, dl, dp
cdef int col_step = -1
# number of iterations
cdef int cnt = 0
Expand Down Expand Up @@ -230,16 +230,17 @@ cdef void one_step_gradient_search_no_gil(const data_type[:, :, :] data,
p0 = int(p0 + dp)
j += col_step


@cython.boundscheck(False)
@cython.wraparound(False)
cpdef one_step_gradient_indices(DTYPE_t [:, :] src_x,
DTYPE_t [:, :] src_y,
DTYPE_t [:, :] xl,
DTYPE_t [:, :] xp,
DTYPE_t [:, :] yl,
DTYPE_t [:, :] yp,
DTYPE_t [:, :] dst_x,
DTYPE_t [:, :] dst_y):
cpdef one_step_gradient_indices(float_index [:, :] src_x,
float_index [:, :] src_y,
float_index [:, :] xl,
float_index [:, :] xp,
float_index [:, :] yl,
float_index [:, :] yp,
float_index [:, :] dst_x,
float_index [:, :] dst_y):
"""Gradient search, simple case variant, returning float indices.
This is appropriate for monotonous gradients only, i.e. not modis or viirs in satellite projection.
Expand All @@ -250,18 +251,24 @@ cpdef one_step_gradient_indices(DTYPE_t [:, :] src_x,
cdef size_t y_size = dst_y.shape[0]
cdef size_t x_size = dst_x.shape[1]


if float_index is double:
dtype = np.float64
else:
dtype = np.float32

# output indices arrays --> needs to be (lines, pixels) --> y,x
indices = np.full([2, y_size, x_size], np.nan, dtype=f_index_dtype)
indices = np.full([2, y_size, x_size], np.nan, dtype=dtype)
cdef float_index [:, :, :] indices_view_result = indices

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

with nogil:
one_step_gradient_search_no_gil[float_index](fake_data,
src_x, src_y,
xl, xp, yl, yp,
dst_x, dst_y,
x_size, y_size,
indices_xy, indices_view_result)
one_step_gradient_search_no_gil[float_index, float_index](fake_data,
src_x, src_y,
xl, xp, yl, yp,
dst_x, dst_y,
x_size, y_size,
indices_xy, indices_view_result)
return indices

0 comments on commit fdea688

Please sign in to comment.