bilinear interpolation for angles

Question

I have a 2d array of directional data. I need to interpolate over a higher resolution grid however the ready made functions like scipy interp2d, etc don't account for the discontinuity between 0 and 360.

I have code for doing this for a single grid of 4 points (thanks How to perform bilinear interpolation in Python and Rotation Interpolation) however I would like it to accept big data sets at once - just like the interp2d function. How can I encorporate this into the code below in a way which doesn't just loop over all of the data?

Thanks!

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

I'm going to reuse some personal Point and Point3D simplified classes for this example:

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) 

It represents a 2D point. It has a simple constructor, x and y properties that ensure they always store floats, magic methods for string representation as (x,y) and comparison to make them sortable (sorts by x, then by y). My original class has additional features such as addition and substraction (vector behaviour) magic methods both but they are not needed for this example.

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)

Same as Point but for 3D points. It also includes an additional constructor classmethod that takes a Point and its z value as arguments.

Linear interpolation

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

It takes a first position argument that will determine the x whose y value we want to calculate, and an infinite amount of Point instances from where we want to interpolate. A keyword argument (extrapolate) allows to turn on extrapolation. A Point instance is returned with the requested x and the calculated y values.

Bilinear interpolation

I offer two alternatives, both of them have a similar signature to the previous interpolation function. A Point whose z value we want to calculate, a keyword argument (extrapolate) that turns on extrapolation and return a Point3D instance with the requested and calculated data. The difference between these two approaches are how the values that will be used to interpolate are provided:

Approach 1

The first approach takes a two-levels-deep nested dict. The first level keys represent the x values, the second level keys the y values and the second level values the z values.

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

Approach 2

The second approach takes an infinite number of Point3D instances.

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

You can see from both approaches that they use the previously defined linear_interpoaltion function, and that they always set extrapolation to True as they already raised an exception if it was False and the requested point was outside the provided interval.