加快在插值中查找点所属单元的速度，从曲线网格到矩形网格

Question

我正在尝试设置一个方法，将曲线网格上的点插值到规则网格上。为此，我需要找出规则网格上的每个点属于曲线网格的哪个单元。曲线网格的单元是随机的四边形。

网格的创建方式如下，举个例子，针对一个不太理想的情况：

import numpy as np
import matplotlib.pyplot as plt

def circle(xn, yn):
    """Polar grid."""
    R = (xn[-1, 0] - xn[0, 0]) / (2*np.pi)
    return (yn + R)*np.sin(xn/R), (yn + R)*np.cos(xn/R)


_x= np.linspace(0, 0.1, 150)
_y = np.linspace(0, 0.1, 100)

# Curvilinear grid
up, vp = circle(np.meshgrid(_x, _y, indexing='ij'))

# Cartesian grid
ui = np.linspace(up.min(), up.max(), 300)
vi = np.linspace(vp.min(), vp.max(), 200)

我能够使用经典的四边形点查找算法，找到每个网格之间的对应关系：

%%cython -f -c=-O2 -I./ --compile-args=-fopenmp --link-args=-fopenmp

cimport cython
from cython.parallel cimport prange
import numpy as np
cimport numpy as np


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
cpdef neighbors(double[:, ::1] xp, double[:, ::1] yp, double xi, double yi):
    """Clockwise winding Polygon :

        [[xp[ix, iy], yp[ix, iy]],
         [xp[ix, iy+1], yp[ix, iy+1]],
         [xp[ix+1, iy+1], yp[ix+1, iy+1]],
         [xp[ix+1, iy], yp[ix+1, iy]]])
    """
    cdef Py_ssize_t nx = xp.shape[0]
    cdef Py_ssize_t ny = xp.shape[1]
    cdef Py_ssize_t ix, iy = 0
    cdef Py_ssize_t tag = 0
    cdef double xmin, xmax, ymin, ymax
    cdef double x1, x2, x3, x4, y1, y2, y3, y4

    for ix in range(nx-1):

        for iy in range(ny-1):

            x1 = xp[ix, iy]
            x2 = xp[ix, iy+1]
            x3 = xp[ix+1, iy+1]
            x4 = xp[ix+1, iy]
            y1 = yp[ix, iy]
            y2 = yp[ix, iy+1]
            y3 = yp[ix+1, iy+1]
            y4 = yp[ix+1, iy]

            if xi - min(x1, x2, x3, x4) < 0.:
                continue

            if yi - min(y1, y2, y3, y4) < 0.:
                continue
            
            # Check that point (xi, yi) is on the same side of each edge of the quadrilateral
            if not (xi - x1) * (y2 - y1) - (x2 - x1) * (yi - y1) > 0:
                continue
            if not (xi - x2) * (y3 - y2) - (x3 - x2) * (yi - y2) > 0:
                continue
            if not (xi - x3) * (y4 - y3) - (x4 - x3) * (yi - y3) > 0 :
                continue
            if not (xi - x4) * (y1 - y4) - (x1 - x4) * (yi - y4) > 0:
                continue

            tag = 1
            break
        if tag == 1:
            break
    if tag == 1:
        return (ix, iy), (ix, iy+1), (ix+1, iy+1), (ix+1, iy)
    return None

而且这个方法是有效的：

fig, axes = plt.subplots(1, 1, figsize=(6, 6))

for i in np.arange(0, ui.shape[0]):
    axes.axvline(ui[i], color='k', linewidth=0.1)
for i in np.arange(0, vi.shape[0]):
    axes.axhline(vi[i], color='k', linewidth=0.1)

for i in np.arange(0, vp.shape[1]):
    axes.plot(up[:, i], vp[:, i], 'k', linewidth=0.5)
for i in np.arange(0, up.shape[0]):
    axes.plot(up[i, :], vp[i, :], 'k', linewidth=0.5)

for iu, iv in ((180, 120), (100, 90), (80, 40), (10, 10)):
    axes.plot(ui[iu], vi[iv], 'ro')
    idx = neighbors(up, vp, ui[iu], vi[iv])
    if idx:
        for _idx in idx:
            axes.plot(up[_idx], vp[_idx], 'go')
plt.show()

问题是，在我的系统上，每个点的查找平均需要15微秒，对于一个2000x2000的网格来说，总共需要7毫秒，也就是说，查找2000x2000个点的估计时间是28000秒（更别提三维的情况了 :））。

我觉得可以通过避免查找那些在曲线网格中不存在的点来节省时间，但我找不到一个简单的方法来跳过这些点的查找。

我认为对于这种研究，可能有更好的全局方法，但我还不太了解。

我该如何利用Cython/numpy来加快这个过程呢？

通过使用简单的二分法来定位x值，可以节省相当多的计算时间。不过，下面的代码如果x值不是单调的，可能会导致错误：

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
cpdef bint is_monotonic(double[:, ::1] v):
    cdef Py_ssize_t i

    for i in range(v.shape[0] - 1):
        if (v[i, 1] > v[i, 0]) != (v[i+1, 1] > v[i+1, 0]):
            return False
    return True


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
cpdef (double, double) minmax(double[:] arr):
    cdef double min = arr[0]
    cdef double max = arr[0]
    cdef Py_ssize_t i
    cdef Py_ssize_t n = arr.shape[0]

    for i in range(1, n):
        if arr[i] < min:
            min = arr[i]
        elif arr[i] > max:
            max = arr[i]
    return min, max


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
cpdef (Py_ssize_t, Py_ssize_t) dichotomy(double[:, ::1] xp, double xi, Py_ssize_t offset=1):

    cdef Py_ssize_t nx = xp.shape[0]
    cdef Py_ssize_t istart = 0
    cdef Py_ssize_t istop = nx - 1
    cdef Py_ssize_t middle = (istart + istop) // 2
    cdef double xmin1, xmin2, xmax1, xmax2
    cdef double xmin, xmax


    while (istop - istart) > 3:

        middle = (istart + istop) // 2

        xmin1, xmax1 = minmax(xp[istart, :])
        xmin2, xmax2 = minmax(xp[middle, :])
        xmin = min(xmin1, xmax1, xmin2, xmax2)
        xmax = max(xmin1, xmax1, xmin2, xmax2)

        if xmin <= xi <= xmax:
            if middle - istart < 3:
                return istart - offset, middle + offset
            istop = middle
        else:
            istart = middle

    xmin1, xmax1 = minmax(xp[middle, :])
    xmin2, xmax2 = minmax(xp[istop, :])
    xmin = min(xmin1, xmax1, xmin2, xmax2)
    xmax = max(xmin1, xmax1, xmin2, xmax2)

    if xmin <= xi <= xmax:
        return middle - offset, istop + offset
    return -1, -1

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
cpdef neighbor(double[:, ::1] xp, double[:, ::1] yp, double xi, double yi):

    (...)

    cdef Py_ssize_t istart = 0
    cdef Py_ssize_t istop = nx - 2

    if is_monotonic(xp):
        istart, istop = dichotomy(xp, xi)
        if istart == -1:
            return None

    for ix in range(istart, istop + 1):
        for iy in range(ny-1):
            (...)