角度的双线性插值

Question

我有一个二维方向数据数组。我需要在更高分辨率的网格上进行插值，但是像 scipy interp2d 等现成的函数并不能解释 0 到 360 之间的不连续性。

我有代码为 4 点的单个网格执行此操作（感谢 How to perform bilinear interpolation in Python 和 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

Answer 1

我将在此示例中重用一些个人 Point 和 Point3D 简化类：

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) 

它代表一个二维点。它有一个简单的构造函数，x 和y 属性确保它们始终存储floats，用于字符串表示的魔术方法(x,y) 和比较以使它们可排序（按x 排序，然后按@987654330 @)。我的原始类具有附加功能，例如加法和减法（向量行为）魔术方法，但本示例不需要它们。

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 点。它还包括一个额外的构造函数类方法，该方法接受 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 实例与请求的 x 和计算的 y 值一起返回。

双线性插值

我提供了两个替代方案，它们都具有与之前的插值函数相似的签名。一个 Point 我们要计算其 z 值，一个关键字参数 (extrapolate) 打开外推并返回一个带有请求和计算数据的 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}}))

方法2

第二种方法采用无限数量的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 并且请求的点在外部，它们已经引发了异常提供的间隔。