加速曲线网格到矩形网格插值的点所属单元查找

我正尝试搭建一套从曲线网格到规则网格的点插值方法，核心需求是确定规则网格的每个点所属的曲线网格单元（曲线网格单元为随机四边形）。

网格创建示例（含不利场景） #

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微秒，2000×2000网格的单点查找耗时7毫秒，若查找2000×2000个点，预估总耗时达28000秒，3D场景耗时更甚。

我认为可通过跳过曲线网格外的点来节省时间，但暂无简单判断指标；同时推测存在更优的全局查找方法，但尚不了解。

尝试的二分法优化及问题 #

我尝试通过针对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):
            (...)

咨询问题 #

如何借助Cython/Numpy提升该点所属单元查找的效率？

内容的提问来源于stack exchange，提问作者Ipse Lium