python: interpolation: finding a rectangle that minimally include a point

Question

I am implementing a bilinear interpolation as in How to perform bilinear interpolation in Python

I have a sorted list of points that are the vertexes of my regular grid.

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

I want to interpolate linearly on the point (x,y). I have written the following code

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))

but this code doesn't work on the boundaries. I would think this is common problem in interpolation, so I was wondering how I should select the smallest rectangle of points around (x,y) from my regular grid.

Find the closest point to x, y and select the grid square surrounding it. — Mad Physicist, May 29 '18 at 13:37
@MadPhysicist yes this is the idea, I was wondering if you have a pseudocode for that — simona, May 29 '18 at 13:37
Define your grid in terms of an origin and unit vectors. It's easy to project onto that, and immediately gives you the cell you're interested in in O(1) time (no searching involved) — Mad Physicist, May 29 '18 at 13:39

score 0 · Answer 1 · answered May 30 '18 at 04:09

Your grid is regular, so you don't need to traverse all points to determine cell indexes. Just divide coordinates by cell size and round result to smaller integer. 1D example: if first point has coordinate 1 and cell size is 2, point 6 lies at int (6-1)/2 = 2-nd interval

Restrict result index to ensure that it is in grid limits - so points outside grid will use border cells

 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

score 0 · Accepted Answer · answered May 30 '18 at 08:10

Following the suggestions in the comments I have written the following code:

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))

python: interpolation: finding a rectangle that minimally include a point

2 Answers2