python: интерполяция: поиск прямоугольника, минимально включающего точку

Question

Я реализую билинейную интерполяцию, как в Как выполнить билинейную интерполяцию в Python

У меня есть отсортированный список точек, которые являются вершинами моей регулярной сетки.

[[x1,y1,z1],[x2,y2,z2],[x3,y3,z3],[x4,y4,z4],[x5,y5,z5],...]

Я хочу линейно интерполировать в точку (x,y).Я написал следующий код

def f(x, y, points):
    for i in range(len(points)-1, -1, -1):
        if (x>points[i][0])and(y>points[i][1]):
            break
    try:
        pp = [points[i], points[i+1]]
    except IndexError:
        pp = [points[i], points[i-1]]

    for j in range(len(points)):
        if (x<points[j][0])and(y<points[j][1]):
            break
    pp.append(points[j-1])
    pp.append(points[j])

    (x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = pp
    return (q11 * (x2 - x) * (y2 - y) +
            q21 * (x - x1) * (y2 - y) +
            q12 * (x2 - x) * (y - y1) +
            q22 * (x - x1) * (y - y1)) / ((x2 - x1) * (y2 - y1))

, но этот код не работает на границах.Я думаю, что это обычная проблема в интерполяции, поэтому мне было интересно, как мне выбрать наименьший прямоугольник точек вокруг (x, y) из моей регулярной сетки.

Answer 1

Следуя предложениям в комментариях, я написал следующий код:

def f(x, y, points):
    points = sorted(points)

    xunique = np.unique([point[0] for point in points])
    yunique = np.unique([point[1] for point in points])
    xmax    = np.max(xunique)
    ymax    = np.max(yunique)
    deltax  = xunique[1] - xunique[0]
    deltay  = yunique[1] - yunique[0]
    x0      = xunique[0]
    y0      = yunique[0]
    ni      = len(xunique)
    nj      = len(yunique)

    i1 = int(np.floor((x-x0)/deltax))
    if i1 == ni:
        i1 = i1 - 1
    i2 = int(np.ceil((x-x0)/deltax))
    if i2 == ni:
        i2 = i2 - 1
    j1 = int(np.floor((y-y0)/deltay))
    if j1 == nj:
        j1 = j1 - 1
    j2 = int(np.ceil((y-y0)/deltay))
    if j2 == ni:
        j2 = j2 - 1

    pp=[]
    if (i1==i2):
        if i1>0:
            i1=i1-1
        else:
            i2=i2+1
    if (j1==j2):
        if j1>0:
            j1=j1-1
        else:
            j2=j2+1

    pp=[points[i1 * nj + j1], points[i1 * nj + j2], 
            points[i2 * nj + j1], points[i2 * nj + j2]]

    (x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = pp
    return (q11 * (x2 - x) * (y2 - y) +
                q21 * (x - x1) * (y2 - y) +
                q12 * (x2 - x) * (y - y1) +
                q22 * (x - x1) * (y - y1)) / ((x2 - x1) * (y2 - y1))

Answer 2

Ваша сетка регулярна, поэтому вам не нужно проходить все точки, чтобы определить индексы ячеек.Просто разделите координаты на размер ячейки и округлите результат до меньшего целого числа.Пример 1D: если первая точка имеет координату 1 и размер ячейки равен 2, точка 6 находится на int (6-1)/2 = 2 -ом интервале

Ограничить индекс результата, чтобы убедиться, что он находится в пределах сетки - поэтому точки вне сетки будут использовать границуячейки

 i = int((x - points[i][0]) / xsize)  #not sure what is the best way in Python
 if (i < 0):
     i = 0
 if (i >=  XCount):
     i = XCount - 1
 // same for j and y-coordinate