127 lines
5.1 KiB
Python
127 lines
5.1 KiB
Python
# import csv
|
|
import math
|
|
import numpy as np
|
|
|
|
|
|
class bilinearInterpolator:
|
|
"""
|
|
This class takes a collection of 3-dimensional points from a .csv file.
|
|
It contains a bilinear interpolator to find unknown points within the grid.
|
|
"""
|
|
@property
|
|
def probedGrid(self):
|
|
return self._probedGrid
|
|
|
|
"""
|
|
Constructor takes a file with a .csv extension and creates an evenly-spaced 'ideal' grid from the data points.
|
|
This is done to get around any floating point errors that may exist in the data
|
|
"""
|
|
def __init__(self, pointsFile):
|
|
|
|
self.pointsFile = pointsFile
|
|
self.points = np.loadtxt(self.pointsFile, delimiter=',')
|
|
|
|
self.xMin, self.xMax, self.xSpacing, self.xCount = self._axisParams(0)
|
|
self.yMin, self.yMax, self.ySpacing, self.yCount = self._axisParams(1)
|
|
|
|
# generate ideal grid to match actually probed points -- this is due to floating-point error issues
|
|
idealGrid = ([
|
|
[(x, y) for x in np.linspace(self.xMin, self.xMax, self.xCount, True)]
|
|
for y in np.linspace(self.yMin, self.yMax, self.yCount, True)
|
|
])
|
|
|
|
self._probedGrid = [[0] * self.yCount for i in range(0, self.xCount)]
|
|
|
|
# align ideal grid indices with probed data points
|
|
for rowIndex, row in enumerate(idealGrid):
|
|
for colIndex, idealPoint in enumerate(row):
|
|
minSqDist = math.inf
|
|
for probed in self.points:
|
|
# find closest point in ideal grid that corresponds to actual tested point
|
|
# put z value in correct index
|
|
sqDist = pow(probed[0] - idealPoint[0], 2) + pow(probed[1] - idealPoint[1], 2)
|
|
if sqDist <= minSqDist:
|
|
minSqDist = sqDist
|
|
indexX = rowIndex
|
|
indexY = colIndex
|
|
closestProbed = probed
|
|
self.probedGrid[indexY][indexX] = closestProbed
|
|
|
|
def Interpolate(self, point):
|
|
"""
|
|
Bilinear interpolation method to determine unknown z-values within grid of known z-values.
|
|
|
|
NOTE: If one axis is outside the grid, linear interpolation is used instead.
|
|
If both axes are outside of the grid, the z-value of the closest corner of the grid is returned.
|
|
"""
|
|
lin = False
|
|
|
|
if point[0] < self.xMin:
|
|
ix1 = 0
|
|
lin = True
|
|
elif point[0] > self.xMax:
|
|
ix1 = self.xCount-1
|
|
lin = True
|
|
else:
|
|
ix1 = math.floor((point[0] - self.xMin)/self.xSpacing)
|
|
ix2 = math.ceil((point[0] - self.xMin)/self.xSpacing)
|
|
|
|
def interpolatePoint(p1, p2, pt, axis):
|
|
return (p2[2]*(pt[axis] - p1[axis]) + p1[2]*(p2[axis] - pt[axis]))/(p2[axis] - p1[axis])
|
|
|
|
if point[1] < self.yMin:
|
|
if lin:
|
|
return self.probedGrid[ix1][0][2]
|
|
return interpolatePoint(self.probedGrid[ix1][0], self.probedGrid[ix2][0], point, 0)
|
|
elif point[1] > self.yMax:
|
|
if lin:
|
|
return self.probedGrid[ix1][self.yCount - 1][2]
|
|
return interpolatePoint(
|
|
self.probedGrid[ix1][self.yCount - 1], self.probedGrid[ix2][self.yCount - 1], point, 0)
|
|
else:
|
|
iy1 = math.floor((point[1] - self.yMin)/self.ySpacing)
|
|
iy2 = math.ceil((point[1] - self.yMin)/self.ySpacing)
|
|
# if x was at an extrema, but y was not, perform linear interpolation on x axis
|
|
if lin:
|
|
return interpolatePoint(self.probedGrid[ix1][iy1], self.probedGrid[ix1][iy2], point, 1)
|
|
|
|
def specialDiv(a, b):
|
|
if b == 0:
|
|
return 0.5
|
|
else:
|
|
return a/b
|
|
|
|
x1 = self.probedGrid[ix1][iy1][0]
|
|
x2 = self.probedGrid[ix2][iy1][0]
|
|
y1 = self.probedGrid[ix2][iy1][1]
|
|
y2 = self.probedGrid[ix2][iy2][1]
|
|
|
|
Q11 = self.probedGrid[ix1][iy1][2]
|
|
Q12 = self.probedGrid[ix1][iy2][2]
|
|
Q21 = self.probedGrid[ix2][iy1][2]
|
|
Q22 = self.probedGrid[ix2][iy2][2]
|
|
|
|
r1 = specialDiv(point[0]-x1, x2-x1)*Q21 + specialDiv(x2-point[0], x2-x1)*Q11
|
|
r2 = specialDiv(point[0]-x1, x2-x1)*Q22 + specialDiv(x2-point[0], x2-x1)*Q12
|
|
p = specialDiv(point[1]-y1, y2-y1)*r2 + specialDiv(y2-point[1], y2-y1)*r1
|
|
|
|
return p
|
|
|
|
# Returns the min, max, spacing and size of one axis of the 2D grid
|
|
def _axisParams(self, sortAxis):
|
|
# sort the set and eliminate the previous, unsorted set
|
|
srtSet = sorted(self.points, key=lambda x: x[sortAxis])
|
|
|
|
dists = []
|
|
for item0, item1 in zip(srtSet[:(len(srtSet)-2)], srtSet[1:]):
|
|
dists.append(float(item1[sortAxis]) - float(item0[sortAxis]))
|
|
axisSpacing = max(dists)
|
|
|
|
# add an extra one for axisCount to account for the starting point
|
|
axisMin = float(min(srtSet, key=lambda x: x[sortAxis])[sortAxis])
|
|
axisMax = float(max(srtSet, key=lambda x: x[sortAxis])[sortAxis])
|
|
axisRange = axisMax - axisMin
|
|
axisCount = round((axisRange/axisSpacing) + 1)
|
|
|
|
return axisMin, axisMax, axisSpacing, axisCount
|