Skip to content

Commit

Permalink
Allow float32 in gradient search
Browse files Browse the repository at this point in the history
  • Loading branch information
@mraspaud
mraspaud committed Sep 3, 2024
1 parent 27e04a4 commit 03cec67
Show file tree
Hide file tree
Showing 2 changed files with 35 additions and 23 deletions.
4 changes: 2 additions & 2 deletions pyresample/gradient/__init__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -414,14 +414,14 @@ def parallel_gradient_search(data, src_x, src_y, dst_x, dst_y,
else:
is_pad = False
res = dask.delayed(_gradient_resample_data)(
arr.astype(np.float64),
arr,
src_x[i], src_y[i],
src_gradient_xl[i], src_gradient_xp[i],
src_gradient_yl[i], src_gradient_yp[i],
dst_x[i], dst_y[i],
method=method)
res = da.from_delayed(res, (num_bands, ) + dst_x[i].shape,
dtype=np.float64).astype(arr.dtype)
dtype=arr.dtype)
if dst_mosaic_locations[i] in chunks:
if not is_pad:
chunks[dst_mosaic_locations[i]].append(res)
Expand Down
54 changes: 33 additions & 21 deletions pyresample/gradient/_gradient_search.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,18 +23,26 @@

import numpy as np

cimport cython
cimport numpy as np
from libc.math cimport fabs, isinf

DTYPE = np.double
ctypedef np.double_t DTYPE_t
cimport cython
from libc.math cimport fabs, isinf

ctypedef fused data_type:
double
float

ctypedef double float_index
f_index_dtype = float
ctypedef int i_index_type


np.import_array()

@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline void nn(const DTYPE_t[:, :, :] data, int l0, int p0, double dl, double dp, int lmax, int pmax, DTYPE_t[:] res) noexcept nogil:
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 int nnl, nnp
cdef size_t z_size = res.shape[0]
cdef size_t i
Expand All @@ -54,7 +62,7 @@ cdef inline void nn(const DTYPE_t[:, :, :] data, int l0, int p0, double dl, doub

@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline void bil(const DTYPE_t[:, :, :] data, int l0, int p0, double dl, double dp, int lmax, int pmax, DTYPE_t[:] res) noexcept nogil:
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 int l_a, l_b, p_a, p_b
cdef double w_l, w_p
cdef size_t z_size = res.shape[0]
Expand Down Expand Up @@ -84,18 +92,18 @@ cdef inline void bil(const DTYPE_t[:, :, :] data, int l0, int p0, double dl, dou

@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline void indices_xy(const DTYPE_t[:, :, :] data, int l0, int p0, double dl, double dp, int lmax, int pmax, DTYPE_t[:] res) noexcept nogil:
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 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 DTYPE_t[:, :, :] data, int l0, int p0, double dl, double dp, int lmax, int pmax, DTYPE_t[:] res) noexcept nogil
ctypedef void (*FN)(const data_type[:, :, :] data, int l0, int p0, double dl, double dp, int lmax, int pmax, data_type[:] res) noexcept nogil

@cython.boundscheck(False)
@cython.wraparound(False)
cpdef one_step_gradient_search(const DTYPE_t [:, :, :] data,
cpdef one_step_gradient_search(const data_type[:, :, :] data,
DTYPE_t [:, :] src_x,
DTYPE_t [:, :] src_y,
DTYPE_t [:, :] xl,
Expand All @@ -118,10 +126,14 @@ cpdef one_step_gradient_search(const DTYPE_t [:, :, :] data,
cdef size_t y_size = dst_y.shape[0]
cdef size_t x_size = dst_x.shape[1]

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

# output image array --> needs to be (lines, pixels) --> y,x
image = np.full([z_size, y_size, x_size], np.nan, dtype=DTYPE)
cdef DTYPE_t [:, :, :] image_view = image
image = np.full([z_size, y_size, x_size], np.nan, dtype=dtype)
cdef data_type[:, :, :] image_view = image
with nogil:
one_step_gradient_search_no_gil(data,
src_x, src_y,
Expand All @@ -135,7 +147,7 @@ cpdef one_step_gradient_search(const DTYPE_t [:, :, :] data,
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef void one_step_gradient_search_no_gil(const DTYPE_t[:, :, :] data,
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,
Expand All @@ -147,7 +159,7 @@ cdef void one_step_gradient_search_no_gil(const DTYPE_t[:, :, :] data,
const size_t x_size,
const size_t y_size,
FN fun,
DTYPE_t[:, :, :] result_array) noexcept nogil:
data_type[:, :, :] result_array) noexcept nogil:

# pixel max ---> data is expected in [lines, pixels]
cdef int pmax = src_x.shape[1] - 1
Expand Down Expand Up @@ -239,17 +251,17 @@ cpdef one_step_gradient_indices(DTYPE_t [:, :] src_x,
cdef size_t x_size = dst_x.shape[1]

# output indices arrays --> needs to be (lines, pixels) --> y,x
indices = np.full([2, y_size, x_size], np.nan, dtype=DTYPE)
cdef DTYPE_t [:, :, :] indices_view = indices
indices = np.full([2, y_size, x_size], np.nan, dtype=f_index_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 DTYPE_t [:, :, :] fake_data = np.full([1, 1, 1], np.nan, dtype=DTYPE)
cdef float_index [:, :, :] fake_data = np.full([1, 1, 1], np.nan, dtype=f_index_dtype)

with nogil:
one_step_gradient_search_no_gil(fake_data,
src_x, src_y,
xl, xp, yl, yp,
dst_x, dst_y,
x_size, y_size,
indices_xy, indices_view)
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)
return indices

0 comments on commit 03cec67

Please sign in to comment.