Бикубическая интерполяция поверхности

Столкнулся с проблемой интерполяции поверхности кубическим сплайном, выкопал весь интернет, но кроме банальной интерполяции в квадрате [0,1] x [[0,1] ничего не нашел, но мне нужно интерполировать значение внутри произвольной области (будь то 16 баллов или 4 с указанием всех производных и т. д.), я приведу пример:
Мы берем точки, которые удовлетворяют уравнению эллиптического параболоида: x = {100..102}, y = {100..102}

х = 100, у = 100, f (х, у) = 1805,56

х = 100, у = 101, ф (х, у) = 1816,72

х = 100, у = 102, ф (х, у) = 1828

х = 101, у = 100, ф (х, у) = 1830,68

х = 101, у = 101, f (х, у) = 1841,85

х = 101, у = 102, ф (х, у) = 1853,12

х = 102, у = 100, ф (х, у) = 1856,06

х = 102, у = 101, ф (х, у) = 1867,22

х = 102, у = 102, ф (х, у) = 1878,5

Как я могу вычислить значение в точке {101.5, 101.5}?
Я прочитал статью в Википедии, но ничего не сказано о том, как она работает в произвольных точках, за пределами квадрата [0,1] x [0,1] Хорошая реализация, которая работает на одном квадрате и использует последовательную интерполяцию здесь

double cubicInterpolate (double p[4], double x) {
return p[1] + 0.5 * x*(p[2] - p[0] + x*(2.0*p[0] - 5.0*p[1] + 4.0*p[2] -
p[3] + x*(3.0*(p[1] - p[2]) + p[3] - p[0])));
}

double bicubicInterpolate (double p[4][4], double x, double y) {
double arr[4];
arr[0] = cubicInterpolate(p[0], y);
arr[1] = cubicInterpolate(p[1], y);
arr[2] = cubicInterpolate(p[2], y);
arr[3] = cubicInterpolate(p[3], y);
return cubicInterpolate(arr, x);
}