角度的双线性插值
bilinear interpolation for angles
我有一个二维定向数据数组。我需要在更高分辨率的网格上进行插值，但是像 scipy interp2d 等现成的函数没有考虑 0 到 360 之间的不连续性。
我有针对单个 4 点网格执行此操作的代码（感谢 How to perform bilinear interpolation in Python and Rotation Interpolation），但我希望它能立即接受大数据集 - 就像 interp2d 函数一样。我怎样才能以一种不只是循环遍历所有数据的方式将其合并到下面的代码中？
谢谢！
def shortest_angle(beg,end,amount):
shortest_angle=((((end - beg) % 360) + 540) % 360) - 180
return shortest_angle*amount
def bilinear_interpolation_rotation(x, y, points):
'''Interpolate (x,y) from values associated with four points.
The four points are a list of four triplets: (x, y, value).
The four points can be in any order. They should form a rectangle.
'''
points = sorted(points) # order points by x, then by y
(x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = points
if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2:
raise ValueError('points do not form a rectangle')
if not x1 <= x <= x2 or not y1 <= y <= y2:
raise ValueError('(x, y) not within the rectangle')
# interpolate over the x value at each y point
fxy1 = q11 + shortest_angle(q11,q21,((x-x1)/(x2-x1)))
fxy2 = q12 + shortest_angle(q12,q22,((x-x1)/(x2-x1)))
# interpolate over the y values
fxy = fxy1 + shortest_angle(fxy1,fxy2,((y-y1)/(y2-y1)))
return fxy
我将在此示例中重用一些个人 Point 和 Point3D 简化的 classes：
Point
class Point:
#Constructors
def __init__(self, x, y):
self.x = x
self.y = y
# Properties
@property
def x(self):
return self._x
@x.setter
def x(self, value):
self._x = float(value)
@property
def y(self):
return self._y
@y.setter
def y(self, value):
self._y = float(value)
# Printing magic methods
def __repr__(self):
return "({p.x},{p.y})".format(p=self)
# Comparison magic methods
def __is_compatible(self, other):
return hasattr(other, 'x') and hasattr(other, 'y')
def __eq__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x == other.x) and (self.y == other.y)
def __ne__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x != other.x) or (self.y != other.y)
def __lt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) < (other.x, other.y)
def __le__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) <= (other.x, other.y)
def __gt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) > (other.x, other.y)
def __ge__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) >= (other.x, other.y)
代表一个二维点。它有一个简单的构造函数、确保它们始终存储 float 的 x 和 y 属性、将字符串表示为 (x,y) 的魔术方法以及使它们可排序的比较（排序依据x，然后 y)。我原来的 class 具有附加功能，例如加法和减法（矢量行为）魔法方法，但本示例不需要它们。
Point3D
class Point3D(Point):
# Constructors
def __init__(self, x, y, z):
super().__init__(x, y)
self.z = z
@classmethod
def from2D(cls, p, z):
return cls(p.x, p.y, z)
# Properties
@property
def z(self):
return self._z
@z.setter
def z(self, value):
self._z = (value + 180.0) % 360 - 180
# Printing magic methods
def __repr__(self):
return "({p.x},{p.y},{p.z})".format(p=self)
# Comparison magic methods
def __is_compatible(self, other):
return hasattr(other, 'x') and hasattr(other, 'y') and hasattr(other, 'z')
def __eq__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x == other.x) and (self.y == other.y) and (self.z == other.z)
def __ne__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x != other.x) or (self.y != other.y) or (self.z != other.z)
def __lt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) < (other.x, other.y, other.z)
def __le__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) <= (other.x, other.y, other.z)
def __gt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) > (other.x, other.y, other.z)
def __ge__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) >= (other.x, other.y, other.z)
与 Point 相同，但用于 3D 点。它还包括一个额外的构造函数 class 方法，该方法采用 Point 及其 z 值作为参数。
线性插值
def linear_interpolation(x, *points, extrapolate=False):
# Check there are a minimum of two points
if len(points) < 2:
raise ValueError("Not enought points given for interpolation.")
# Sort the points
points = sorted(points)
# Check that x is the valid interpolation interval
if not extrapolate and (x < points[0].x or x > points[-1].x):
raise ValueError("{} is not in the interpolation interval.".format(x))
# Determine which are the two surrounding interpolation points
if x < points[0].x:
i = 0
elif x > points[-1].x:
i = len(points)-2
else:
i = 0
while points[i+1].x < x:
i += 1
p1, p2 = points[i:i+2]
# Interpolate
return Point(x, p1.y + (p2.y-p1.y) * (x-p1.x) / (p2.x-p1.x))
它采用第一个位置参数来确定我们要计算其 y 值的 x，以及我们要从中插值的无限数量的 Point 个实例。关键字参数 (extrapolate) 允许打开外推。 Point 实例是 return 使用请求的 x 和计算的 y 值编辑的。
双线性插值
我提供了两个替代方案，它们都具有与之前的插值函数相似的签名。我们要计算其 z 值的 Point，一个打开外推的关键字参数 (extrapolate) 和 return 一个包含请求和计算数据的 Point3D 实例。这两种方法之间的区别在于如何提供将用于插值的值：
方法一
第一种方法采用两层深度嵌套 dict。第一级键表示 x 值，第二级键表示 y 值，第二级键表示 z 值。
def bilinear_interpolation(p, points, extrapolate=False):
x_values = sorted(points.keys())
# Check there are a minimum of two x values
if len(x_values) < 2:
raise ValueError("Not enought points given for interpolation.")
y_values = set()
for value in points.values():
y_values.update(value.keys())
y_values = sorted(y_values)
# Check there are a minimum of two y values
if len(y_values) < 2:
raise ValueError("Not enought points given for interpolation.")
# Check that p is in the valid interval
if not extrapolate and (p.x < x_values[0] or p.x > x_values[-1] or p.y < y_values[0] or p.y > y_values[-1]):
raise ValueError("{} is not in the interpolation interval.".format(p))
# Determine which are the four surrounding interpolation points
if p.x < x_values[0]:
i = 0
elif p.x > x_values[-1]:
i = len(x_values) - 2
else:
i = 0
while x_values[i+1] < p.x:
i += 1
if p.y < y_values[0]:
j = 0
elif p.y > y_values[-1]:
j = len(y_values) - 2
else:
j = 0
while y_values[j+1] < p.y:
j += 1
surroundings = [
Point(x_values[i ], y_values[j ]),
Point(x_values[i ], y_values[j+1]),
Point(x_values[i+1], y_values[j ]),
Point(x_values[i+1], y_values[j+1]),
]
for i, surrounding in enumerate(surroundings):
try:
surroundings[i] = Point3D.from2D(surrounding, points[surrounding.x][surrounding.y])
except KeyError:
raise ValueError("{} is missing in the interpolation grid.".format(surrounding))
p1, p2, p3, p4 = surroundings
# Interpolate
p12 = Point3D(p1.x, p.y, linear_interpolation(p.y, Point(p1.y,p1.z), Point(p2.y,p2.z), extrapolate=True).y)
p34 = Point3D(p3.x, p.y, linear_interpolation(p.y, Point(p3.y,p3.z), Point(p4.y,p4.z), extrapolate=True).y)
return Point3D(p.x, p12.y, linear_interpolation(p.x, Point(p12.x,p12.z), Point(p34.x,p34.z), extrapolate=True).y)
print(bilinear_interpolation(Point(2,3), {1: {2: 5, 4: 6}, 3: {2: 3, 4: 9}}))
方法二
第二种方法需要无限数量的 Point3D 个实例。
def bilinear_interpolation(p, *points, extrapolate=False):
# Check there are a minimum of four points
if len(points) < 4:
raise ValueError("Not enought points given for interpolation.")
# Sort the points into a grid
x_values = set()
y_values = set()
for point in sorted(points):
x_values.add(point.x)
y_values.add(point.y)
x_values = sorted(x_values)
y_values = sorted(y_values)
# Check that p is in the valid interval
if not extrapolate and (p.x < x_values[0] or p.x > x_values[-1] or p.y < y_values[0] or p.y > y_values[-1]):
raise ValueError("{} is not in the interpolation interval.".format(p))
# Determine which are the four surrounding interpolation points
if p.x < x_values[0]:
i = 0
elif p.x > x_values[-1]:
i = len(x_values) - 2
else:
i = 0
while x_values[i+1] < p.x:
i += 1
if p.y < y_values[0]:
j = 0
elif p.y > y_values[-1]:
j = len(y_values) - 2
else:
j = 0
while y_values[j+1] < p.y:
j += 1
surroundings = [
Point(x_values[i ], y_values[j ]),
Point(x_values[i ], y_values[j+1]),
Point(x_values[i+1], y_values[j ]),
Point(x_values[i+1], y_values[j+1]),
]
for point in points:
for i, surrounding in enumerate(surroundings):
if point.x == surrounding.x and point.y == surrounding.y:
surroundings[i] = point
for surrounding in surroundings:
if not isinstance(surrounding, Point3D):
raise ValueError("{} is missing in the interpolation grid.".format(surrounding))
p1, p2, p3, p4 = surroundings
# Interpolate
p12 = Point3D(p1.x, p.y, linear_interpolation(p.y, Point(p1.y,p1.z), Point(p2.y,p2.z), extrapolate=True).y)
p34 = Point3D(p3.x, p.y, linear_interpolation(p.y, Point(p3.y,p3.z), Point(p4.y,p4.z), extrapolate=True).y)
return Point3D(p.x, p12.y, linear_interpolation(p.x, Point(p12.x,p12.z), Point(p34.x,p34.z), extrapolate=True).y)
print(bilinear_interpolation(Point(2,3), Point3D(3,2,3), Point3D(1,4,6), Point3D(3,4,9), Point3D(1,2,5)))
你可以从这两种方法中看到它们使用了先前定义的 linear_interpoaltion 函数，并且它们总是将 extrapolation 设置为 True 因为如果它是 False 并且请求的点在提供的间隔之外。
我有一个二维定向数据数组。我需要在更高分辨率的网格上进行插值，但是像 scipy interp2d 等现成的函数没有考虑 0 到 360 之间的不连续性。
我有针对单个 4 点网格执行此操作的代码（感谢 How to perform bilinear interpolation in Python and Rotation Interpolation），但我希望它能立即接受大数据集 - 就像 interp2d 函数一样。我怎样才能以一种不只是循环遍历所有数据的方式将其合并到下面的代码中？
谢谢！
def shortest_angle(beg,end,amount):
shortest_angle=((((end - beg) % 360) + 540) % 360) - 180
return shortest_angle*amount
def bilinear_interpolation_rotation(x, y, points):
'''Interpolate (x,y) from values associated with four points.
The four points are a list of four triplets: (x, y, value).
The four points can be in any order. They should form a rectangle.
'''
points = sorted(points) # order points by x, then by y
(x1, y1, q11), (_x1, y2, q12), (x2, _y1, q21), (_x2, _y2, q22) = points
if x1 != _x1 or x2 != _x2 or y1 != _y1 or y2 != _y2:
raise ValueError('points do not form a rectangle')
if not x1 <= x <= x2 or not y1 <= y <= y2:
raise ValueError('(x, y) not within the rectangle')
# interpolate over the x value at each y point
fxy1 = q11 + shortest_angle(q11,q21,((x-x1)/(x2-x1)))
fxy2 = q12 + shortest_angle(q12,q22,((x-x1)/(x2-x1)))
# interpolate over the y values
fxy = fxy1 + shortest_angle(fxy1,fxy2,((y-y1)/(y2-y1)))
return fxy
我将在此示例中重用一些个人 Point 和 Point3D 简化的 classes：
Point
class Point:
#Constructors
def __init__(self, x, y):
self.x = x
self.y = y
# Properties
@property
def x(self):
return self._x
@x.setter
def x(self, value):
self._x = float(value)
@property
def y(self):
return self._y
@y.setter
def y(self, value):
self._y = float(value)
# Printing magic methods
def __repr__(self):
return "({p.x},{p.y})".format(p=self)
# Comparison magic methods
def __is_compatible(self, other):
return hasattr(other, 'x') and hasattr(other, 'y')
def __eq__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x == other.x) and (self.y == other.y)
def __ne__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x != other.x) or (self.y != other.y)
def __lt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) < (other.x, other.y)
def __le__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) <= (other.x, other.y)
def __gt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) > (other.x, other.y)
def __ge__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y) >= (other.x, other.y)
代表一个二维点。它有一个简单的构造函数、确保它们始终存储 float 的 x 和 y 属性、将字符串表示为 (x,y) 的魔术方法以及使它们可排序的比较（排序依据x，然后 y)。我原来的 class 具有附加功能，例如加法和减法（矢量行为）魔法方法，但本示例不需要它们。
Point3D
class Point3D(Point):
# Constructors
def __init__(self, x, y, z):
super().__init__(x, y)
self.z = z
@classmethod
def from2D(cls, p, z):
return cls(p.x, p.y, z)
# Properties
@property
def z(self):
return self._z
@z.setter
def z(self, value):
self._z = (value + 180.0) % 360 - 180
# Printing magic methods
def __repr__(self):
return "({p.x},{p.y},{p.z})".format(p=self)
# Comparison magic methods
def __is_compatible(self, other):
return hasattr(other, 'x') and hasattr(other, 'y') and hasattr(other, 'z')
def __eq__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x == other.x) and (self.y == other.y) and (self.z == other.z)
def __ne__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x != other.x) or (self.y != other.y) or (self.z != other.z)
def __lt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) < (other.x, other.y, other.z)
def __le__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) <= (other.x, other.y, other.z)
def __gt__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) > (other.x, other.y, other.z)
def __ge__(self, other):
if not self.__is_compatible(other):
return NotImplemented
return (self.x, self.y, self.z) >= (other.x, other.y, other.z)
与 Point 相同，但用于 3D 点。它还包括一个额外的构造函数 class 方法，该方法采用 Point 及其 z 值作为参数。
线性插值
def linear_interpolation(x, *points, extrapolate=False):
# Check there are a minimum of two points
if len(points) < 2:
raise ValueError("Not enought points given for interpolation.")
# Sort the points
points = sorted(points)
# Check that x is the valid interpolation interval
if not extrapolate and (x < points[0].x or x > points[-1].x):
raise ValueError("{} is not in the interpolation interval.".format(x))
# Determine which are the two surrounding interpolation points
if x < points[0].x:
i = 0
elif x > points[-1].x:
i = len(points)-2
else:
i = 0
while points[i+1].x < x:
i += 1
p1, p2 = points[i:i+2]
# Interpolate
return Point(x, p1.y + (p2.y-p1.y) * (x-p1.x) / (p2.x-p1.x))
它采用第一个位置参数来确定我们要计算其 y 值的 x，以及我们要从中插值的无限数量的 Point 个实例。关键字参数 (extrapolate) 允许打开外推。 Point 实例是 return 使用请求的 x 和计算的 y 值编辑的。
双线性插值
我提供了两个替代方案，它们都具有与之前的插值函数相似的签名。我们要计算其 z 值的 Point，一个打开外推的关键字参数 (extrapolate) 和 return 一个包含请求和计算数据的 Point3D 实例。这两种方法之间的区别在于如何提供将用于插值的值：
方法一
第一种方法采用两层深度嵌套 dict。第一级键表示 x 值，第二级键表示 y 值，第二级键表示 z 值。
def bilinear_interpolation(p, points, extrapolate=False):
x_values = sorted(points.keys())
# Check there are a minimum of two x values
if len(x_values) < 2:
raise ValueError("Not enought points given for interpolation.")
y_values = set()
for value in points.values():
y_values.update(value.keys())
y_values = sorted(y_values)
# Check there are a minimum of two y values
if len(y_values) < 2:
raise ValueError("Not enought points given for interpolation.")
# Check that p is in the valid interval
if not extrapolate and (p.x < x_values[0] or p.x > x_values[-1] or p.y < y_values[0] or p.y > y_values[-1]):
raise ValueError("{} is not in the interpolation interval.".format(p))
# Determine which are the four surrounding interpolation points
if p.x < x_values[0]:
i = 0
elif p.x > x_values[-1]:
i = len(x_values) - 2
else:
i = 0
while x_values[i+1] < p.x:
i += 1
if p.y < y_values[0]:
j = 0
elif p.y > y_values[-1]:
j = len(y_values) - 2
else:
j = 0
while y_values[j+1] < p.y:
j += 1
surroundings = [
Point(x_values[i ], y_values[j ]),
Point(x_values[i ], y_values[j+1]),
Point(x_values[i+1], y_values[j ]),
Point(x_values[i+1], y_values[j+1]),
]
for i, surrounding in enumerate(surroundings):
try:
surroundings[i] = Point3D.from2D(surrounding, points[surrounding.x][surrounding.y])
except KeyError:
raise ValueError("{} is missing in the interpolation grid.".format(surrounding))
p1, p2, p3, p4 = surroundings
# Interpolate
p12 = Point3D(p1.x, p.y, linear_interpolation(p.y, Point(p1.y,p1.z), Point(p2.y,p2.z), extrapolate=True).y)
p34 = Point3D(p3.x, p.y, linear_interpolation(p.y, Point(p3.y,p3.z), Point(p4.y,p4.z), extrapolate=True).y)
return Point3D(p.x, p12.y, linear_interpolation(p.x, Point(p12.x,p12.z), Point(p34.x,p34.z), extrapolate=True).y)
print(bilinear_interpolation(Point(2,3), {1: {2: 5, 4: 6}, 3: {2: 3, 4: 9}}))
方法二
第二种方法需要无限数量的 Point3D 个实例。
def bilinear_interpolation(p, *points, extrapolate=False):
# Check there are a minimum of four points
if len(points) < 4:
raise ValueError("Not enought points given for interpolation.")
# Sort the points into a grid
x_values = set()
y_values = set()
for point in sorted(points):
x_values.add(point.x)
y_values.add(point.y)
x_values = sorted(x_values)
y_values = sorted(y_values)
# Check that p is in the valid interval
if not extrapolate and (p.x < x_values[0] or p.x > x_values[-1] or p.y < y_values[0] or p.y > y_values[-1]):
raise ValueError("{} is not in the interpolation interval.".format(p))
# Determine which are the four surrounding interpolation points
if p.x < x_values[0]:
i = 0
elif p.x > x_values[-1]:
i = len(x_values) - 2
else:
i = 0
while x_values[i+1] < p.x:
i += 1
if p.y < y_values[0]:
j = 0
elif p.y > y_values[-1]:
j = len(y_values) - 2
else:
j = 0
while y_values[j+1] < p.y:
j += 1
surroundings = [
Point(x_values[i ], y_values[j ]),
Point(x_values[i ], y_values[j+1]),
Point(x_values[i+1], y_values[j ]),
Point(x_values[i+1], y_values[j+1]),
]
for point in points:
for i, surrounding in enumerate(surroundings):
if point.x == surrounding.x and point.y == surrounding.y:
surroundings[i] = point
for surrounding in surroundings:
if not isinstance(surrounding, Point3D):
raise ValueError("{} is missing in the interpolation grid.".format(surrounding))
p1, p2, p3, p4 = surroundings
# Interpolate
p12 = Point3D(p1.x, p.y, linear_interpolation(p.y, Point(p1.y,p1.z), Point(p2.y,p2.z), extrapolate=True).y)
p34 = Point3D(p3.x, p.y, linear_interpolation(p.y, Point(p3.y,p3.z), Point(p4.y,p4.z), extrapolate=True).y)
return Point3D(p.x, p12.y, linear_interpolation(p.x, Point(p12.x,p12.z), Point(p34.x,p34.z), extrapolate=True).y)
print(bilinear_interpolation(Point(2,3), Point3D(3,2,3), Point3D(1,4,6), Point3D(3,4,9), Point3D(1,2,5)))
你可以从这两种方法中看到它们使用了先前定义的 linear_interpoaltion 函数，并且它们总是将 extrapolation 设置为 True 因为如果它是 False 并且请求的点在提供的间隔之外。