Бикубическая сплайн-интерполяция гауссианы ⇐ Python

[Anonymous]

Я пытаюсь выполнить 2D-интерполяцию бикубическим сплайном согласно этой статье (стр. 15 сек. 6.2 Бикубическая интерполяция)
Мой тест --
У меня есть функция Гаусса, определенная для интервал $ x,y \in [-1,1] $, и я хочу использовать метод бикубической интерполяции, чтобы получить разумную интерполяцию функции.
Это моя реализация пока -

Код: Выделить всё

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)

Когда я визуализирую Z_interpolated (зеленый), он сильно отличается от Z_actual (красный) — производные не только имеют неправильную величину в углах, но и имеют неверный знак.
Я хочу, чтобы схема интерполяции лучше фиксировала производные в четырех указанных мной точках.
Это ограничение бикубического сплайна или есть ли ошибка в реализация?
Должен ли я включить больше точек сетки? Мне кажется, поскольку это схема кусочной интерполяции, должна ли она действительно улучшить ситуацию?

Подробнее здесь: https://stackoverflow.com/questions/793 ... a-gaussian