Bikubische Spline-Interpolation einer Gaußschen Funktion ⇐ Python

[Guest]

Ich versuche, im Anschluss an diesen Artikel (S. 15 Sek. 6.2 Bikubische Interpolation) eine zweidimensionale bikubische Spline-Interpolation durchzuführen.
Mein Test –
Ich habe eine Gaußsche Funktion auf dem definiert Intervall $ x,y \in [-1,1] $ und ich möchte die bikubische Interpolationsmethode verwenden, um eine vernünftige Interpolation der Funktion zu erhalten.
Dies ist meine Implementierung bisher -

Code: Select all

import numpy as np
from numdifftools import Hessian, Gradient

def gaussian(pt):
x, y = pt
return np.exp(-x**2 - y**2)

A_inv = np.array([
[1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0],
[-3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0],
[2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0],
[0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0],
[0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0],
[0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0],
[-3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0],
[0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0],
[9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1],
[-6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1],
[2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0],
[0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0],
[-6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1],
[4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1]
])

def get_coefficient_tensor(func_array):
alpha_array = A_inv @ func_array
return np.array([
[alpha_array[0], alpha_array[1], alpha_array[2], alpha_array[3]],
[alpha_array[4], alpha_array[5], alpha_array[6], alpha_array[7]],
[alpha_array[8], alpha_array[9], alpha_array[10], alpha_array[11]],
[alpha_array[12], alpha_array[13], alpha_array[14], alpha_array[15]]
]).T

def constructed_function(x, y, x_range, y_range, func_array):
'''
Returns the value of the constructed interpolation function at the point (x,y)

Parameters:
x : float
x coordinate of the point
y : float
y coordinate of the point
x_range : list / touple of 2 floats
range of the x coordinate of the current square grid being considered
y_range : list / touple of 2 floats
range of the y coordinate of the current square grid being considered
func_array :  numpy array of size 16
array containing the user defined quantities - 16x1 dimensional array
- first 4 elements are the values of the function at the 4 corners of the square
- next 8 elements are the values of the partial derivatives in the two directions at the 4 corners of the square (susceptibilities)
- last 4 elements are the values of the mixed partial derivatives at the 4 corners of the square (I don't know how to interpret this yet)
'''
t = (x - x_range[0]) / (x_range[1] - x_range[0])
u = (y - y_range[0]) / (y_range[1] - y_range[0])
coefficient_tensor = get_coefficient_tensor(func_array)
return np.sum([
coefficient_tensor[i,j] * t**i * u**j for i in range(4) for j in range(4)
])

def get_func_array(func, x_range, y_range):
# function values at the 4 corners
f_00 = func([x_range[0], y_range[0]])
f_01 = func([x_range[0], y_range[1]])
f_10 = func([x_range[1], y_range[0]])
f_11 = func([x_range[1], y_range[1]])

# partial derivatives at the 4 corners
f_x_00, f_y_00 = Gradient(func)([x_range[0], y_range[0]])
f_x_01, f_y_01 = Gradient(func)([x_range[0], y_range[1]])
f_x_10, f_y_10 = Gradient(func)([x_range[1], y_range[0]])
f_x_11, f_y_11 = Gradient(func)([x_range[1], y_range[1]])

# cross partial derivatives at the 4 corners
f_xy_00 = Hessian(func)([x_range[0],y_range[0]])[0,1]
f_xy_01 = Hessian(func)([x_range[0],y_range[1]])[0,1]
f_xy_10 = Hessian(func)([x_range[1],y_range[0]])[0,1]
f_xy_11 = Hessian(func)([x_range[1],y_range[1]])[0,1]
return np.array([f_00,f_01,f_10,f_11,f_x_00,f_x_01,f_x_10,f_x_11,f_y_00,f_y_01,f_y_10,f_y_11,f_xy_00,f_xy_01,f_xy_10,f_xy_11]).T

x_range = [-1,1]
y_range = [-1,1]

xs = np.linspace(x_range[0], x_range[1], 100)
ys = np.linspace(y_range[0], y_range[1], 100)

X, Y = np.meshgrid(xs, ys)

func_array = get_func_array(gaussian, x_range, y_range)

Z = np.array([
constructed_function(x, y, x_range, y_range, func_array) for x in xs for y in ys
]).reshape(X.shape)

Z_actual = np.array([
gaussian([x,y]) for x in xs for y in ys
]).reshape(X.shape)

Wenn ich mir Z_interpolated (grün) vorstelle, unterscheidet es sich stark von Z_actual (rot) – nicht nur haben die Ableitungen an den Ecken die falsche Größe, sie haben auch das falsche Vorzeichen.
Ich möchte, dass das Interpolationsschema die Ableitungen an den vier von mir angegebenen Punkten besser erfasst.
Ist dies eine Einschränkung des bikubischen Splines oder? Gibt es einen Fehler in der Implementierung?
Soll ich mehr Rasterpunkte einschließen? Ich denke, da es sich um ein stückweises Interpolationsschema handelt, sollte es die Dinge wirklich verbessern?