CODEWITHSUNDEEP
python-3.xinterpolation
ThaNoob
ThaNoob

Reputation: 622

Python implementation of bilinear quadrilateral interpolation

I'm trying to perform bilinear quadrilateral interpolation. So I have four nodes with known values and I want to find a value that lies in between those four nodes by interpolation, but the four nodes do not form a rectangle. 4-node sketch

I found several ways to solve this, but none of them is implemented in Python already. Does there exist somewhere an already finished python implementation? If not which of the two solutions below would you recommend? Or would you recommend another approach?

**************Different solutions*******************

Solution 1:

I found here, https://www.colorado.edu/engineering/CAS/courses.d/IFEM.d/IFEM.Ch16.d/IFEM.Ch16.pdf, that I should solve the following set of equations: set of equations with Ni being: N definition.

Finally this results in solving a set of equations of the form:

a*x+b*y+c*xy=z1
d*x+e*y+f*xy=z2

with x and y being the unknowns. This could be solved numerically using fsolve.

Solution 2:

This one is completely explained here: https://math.stackexchange.com/questions/828392/spatial-interpolation-for-irregular-grid

but it's quite complex and I think it will take me longer to code it.

Upvotes: 2

Views: 1070

Answers (1)

ThaNoob
ThaNoob

Reputation: 622

Due to a lack of answers I went for the first option. You can find the code below. Recommendations to improve this code are always welcome.

import numpy as np
from scipy.optimize import fsolve

def interpolate_quatrilateral(pt1,pt2,pt3,pt4,pt):
    '''Interpolates a value in a quatrilateral figure defined by 4 points. 
    Each point is a tuple with 3 elements, x-coo,y-coo and value.
    point1 is the lower left corner, point 2 the lower right corner,
    point 3 the upper right corner and point 4 the upper left corner.
    args is a list of coordinates in the following order:
     x1,x2,x3,x4 and x (x-coo of point to be interpolated) and y1,y2...
     code based on the theory found here:
     https://www.colorado.edu/engineering/CAS/courses.d/IFEM.d/IFEM.Ch16.d/IFEM.Ch16.pdf'''

    coos = (pt1[0],pt2[0],pt3[0],pt4[0],pt[0],
            pt1[1],pt2[1],pt3[1],pt4[1],pt[1]) #coordinates of the points merged in tuple
    guess = np.array([0,0]) #The center of the quadrilateral seem like a good place to start
    [eta, mu] = fsolve(func=find_local_coo_equations, x0=guess, args=coos)

    densities = (pt1[2], pt2[2], pt3[2], pt4[2])
    density = find_density(eta,mu,densities)

    return density

def find_local_coo_equations(guess, *args):
    '''This function creates the transformed coordinate equations of the quatrilateral.'''

    eta = guess[0]
    mu = guess[1]

    eq=[0,0]#Initialize eq
    eq[0] = 1 / 4 * (args[0] + args[1] + args[2] + args[3]) - args[4] + \
            1 / 4 * (-args[0] - args[1] + args[2] + args[3]) * mu + \
            1 / 4 * (-args[0] + args[1] + args[2] - args[3]) * eta + \
            1 / 4 * (args[0] - args[1] + args[2] - args[3]) * mu * eta
    eq[1] = 1 / 4 * (args[5] + args[6] + args[7] + args[8]) - args[9] + \
            1 / 4 * (-args[5] - args[6] + args[7] + args[8]) * mu + \
            1 / 4 * (-args[5] + args[6] + args[7] - args[8]) * eta + \
            1 / 4 * (args[5] - args[6] + args[7] - args[8]) * mu * eta
    return eq

def find_density(eta,mu,densities):
    '''Finds the final density based on the eta and mu local coordinates calculated
    earlier and the densities of the 4 points'''
    N1 = 1/4*(1-eta)*(1-mu)
    N2 = 1/4*(1+eta)*(1-mu)
    N3 = 1/4*(1+eta)*(1+mu)
    N4 = 1/4*(1-eta)*(1+mu)
    density = densities[0]*N1+densities[1]*N2+densities[2]*N3+densities[3]*N4
    return density

pt1= (0,0,1)
pt2= (1,0,1)
pt3= (1,1,2)
pt4= (0,1,2)
pt= (0.5,0.5)
print(interpolate_quatrilateral(pt1,pt2,pt3,pt4,pt))

Upvotes: 3

Related Questions