Что не так с расчетами?

[WorldMaster]

Здравствуйте.
Никак не пойму то-ли глаза замылились то-ли еще что-то, но не пойму в чем ошибка при расчете интерполяции.

Задача:
Нужно визуализировать матрицу значений как показано на рисунке.
Если я правильно понял то это можно достигнуть бикубической интерполяцией.

Написал код. Для проверки взял простую матрицу 4на4 и пытаюсь интерполировать ее на матрицу 16 на 16 для проверки алгоритма.

public class Bicubic
    {
        double[][] M;           
        double[] X3, X2, X1, X0;
        double[][] C;   

        public Bicubic(double[][] G)
        {
            M = new double[4][];
            M[0] = new double[] { -0.5, 1.5, -1.5, 0.5 };
            M[1] = new double[] { 1, -2.5, 2, -0.5 };
            M[2] = new double[] { -0.5, 0, 0.5, 0 };
            M[3] = new double[] { 0, 1, 0, 0 };

            C = new double[4][];
            C[0] = new double[4];
            C[1] = new double[4];
            C[2] = new double[4];
            C[3] = new double[4];

            X3 = C[0];
            X2 = C[1];
            X1 = C[2];
            X0 = C[3];

            double[][] T = new double[4][];
            T[0] = new double[4];
            T[1] = new double[4];
            T[2] = new double[4];
            T[3] = new double[4];

            for (int i = 0; i < 4; i++)    // T = G * MT
                for (int j = 0; j < 4; j++)
                    for (int k = 0; k < 4; k++)
                        T[i][j] += G[i][k] * M[j][k];

            for (int i = 0; i < 4; i++)    // C = M * T
                for (int j = 0; j < 4; j++)
                    for (int k = 0; k < 4; k++)
                        C[i][j] += M[i][k] * T[k][j];
        }

        public double eval(double x, double y)
        {
            return x * (x * (x * (y * (y * (y * X3[0] + X3[1]) + X3[2]) + X3[3])
                         + (y * (y * (y * X2[0] + X2[1]) + X2[2]) + X2[3]))
                      + (y * (y * (y * X1[0] + X1[1]) + X1[2]) + X1[3]))
                   + (y * (y * (y * X0[0] + X0[1]) + X0[2]) + X0[3]);
        }
    }

 double[][] arr = new double[4][];  

            arr[0] = new double[] { 1, 1, 1, 1 };
            arr[1] = new double[] { 1, 2, 2, 1 };
            arr[2] = new double[] { 1, 2, 2, 1 };
            arr[3] = new double[] { 1, 1, 1, 1 };

            double[,] dst = new double[16, 16];

            Bicubic ipol = new Bicubic(arr);

            string res = "";
            double step = 0.25; // шаг для размерности 16 на 16
            for (int i = 0; i < 16; i++)
            {
                for (int j = 0; j < 16; j++)
                {
                    double ax = ((i*step) - 1);
                    double ay = ((j*step) - 1);
                    double val = ipol.eval(ax, ay); // вычисляем значение в квадрате [-1,0,1,2] X [-1,0,1,2]
                    dst[i, j] = val;
                    res += dst[i, j].ToString("f4") + "\t";
                }
                res += Environment.NewLine;
            }

в Итоге получаю какую то фигню:
Почему то интерполяция выполняется на меньшую размерность чем 16.
В чем может быть ошибка то??
Если шаг сделать 1 то я получаю исходную матрицу.

[Pavia]

Сообщение от WorldMaster

double ax = ((i*step) - 1);
double ay = ((j*step) - 1);

Неправильное масштабирование. Алгоритм правильного масштабирования переносим в центр умножаем на масштабируемый коэффициент возвращаем назад.
double ax=(i-(16-1)/2)*step+(4-1)/2;
double ay=(j-(16-1)/2)*step+(4-1)/2;

Советую найти в сети готовый код.

[WorldMaster]

Сообщение от Pavia

double ax=(i-(16-1)/2)*step+(4-1)/2;
double ay=(j-(16-1)/2)*step+(4-1)/2;

Советую найти в сети готовый код.

Ничего не изменилось. Может быть я плохо роюсь но примеров рабочих не нашел. Сами интерполяторы вроде попадаются но вот использование ... Хотя и делаю по аналогии с теми что находятся.

[Pavia]

Ничего если на Delphi?

Коэффициенты для InterpolationCubic найдены следующим способом.

f=c3*t^3+c2*t^2+c1*t+c0;
Составим систему и найдём коэффициенты
Дано:
При t=0.0, f==a0;
При t=1.0, f==a1;
При t=2.0, f==a2;
При t=3.0, f==a3;
Найти c0, c1, c2, c3.

Система:
A*X=Y

  [c3]
  [c4]
X=[c1]
  [c0]
A- Матрица Вандерморда
  [ 0 0 0 1]
A=[ 1 1 1 1]
  [ 8 4 2 1]
  [27 9 3 1]

  [a0]
  [a1]
Y=[a2]
  [a3]


решение по методу Гауса
1.
c0=a0;
c3+c2+c1+c0=a1;
c3*8+c2*4+c1*2+c0=a2;
c3*27+c2*9+c1*3+c0=a3;
2.
c0=a0;
3.
c3+c2+c1=a1-a0;
c3*8+c2*4+c1*2=a2-a0;
c3*27+c2*9+c1*3=a3-a0;
4.
c3+c2+c1=a1-a0;
  -c2*4-c1*6=a2-8*a1+7*a0;
  -c2*18-c1*24=a3-27*a1+26*a0;
5.
c3+c2+c1=a1-a0;
  -c2*4-c1*6=a2-8*a1+7*a0;
       +c1*3=a3-a2*4.5+9*a1-5.5*a0;

6.
c1=1/3*(a3-a2*4.5+9*a1-5.5*a0)
7.
  -c2*4-c1*6=a2-8*a1+7*a0;
       +c1*6=2*a3-a2*9+18*a1-11*a0;
8.
-c2*4=2*a3-8*a2+10*a1-4*a0
c2*4=-2*a3+8*a2-10*a1+4*a0
c2*2=-a3+4*a2-5*a1+2*a0
c2=-0.5*a3+2*a2-2.5*a1+a0
9.
c3+c2+c1=a1-a0;
c3=a1-a0+0.5*a3-2*a2+2.5*a1-1*a0-1/3*(a3-a2*4.5+9*a1-5.5*a0);
c3=1/6*a3-0.5*a2+0.5*a1-1/6*a0;

Ответ:
c0=a0;
c1=1/3*(a3-a2*4.5+9*a1-5.5*a0)
c2=-0.5*a3+2*a2-2.5*a1+a0
c3=1/6*a3-0.5*a2+0.5*a1-1/6*a0

function InterpolationCubic(t:Real; a0,a1,a2, a3:Real):Real; Overload;
begin
Result:=(1/6*a3-0.5*a2+0.5*a1-1/6*a0)*t*t*t+(-0.5*a3+2*a2-2.5*a1+a0)*t*t+1/3*(a3-a2*4.5+9*a1-5.5*a0)*t+a0;
end;

// Интерполируем входной массив в выходной.
// Размеры массивов должны быть заданы до вызова функции, !ОСОБЕННО для выходного.
procedure InterpolationCubic(var OutA:TArraySmallint;InA:TArraySmallint); Overload;
var OutR,InR:TArrayReal;
 dT, FactorSize:Real;
 OldLength,NewLength:Integer;
 a0,a1,a2,a3:Real;
 i,j:Integer;
 T,F:Real;
begin
InR:=ToReal(InA);
SetLength(OutR, Length(OutA));
  OldLength:=Length(InA);
  NewLength:=Length(OutA);
  if (OldLength<1) then exit;
  FactorSize:=(NewLength-1)/(OldLength-1);
  dT:=1/FactorSize;
  for j:=0 to NewLength-1 do
    begin
    T:=j*dt;
    i:=Trunc(T); // Берём целую часть числа Trunc(1.7)=1 
    F:=Frac(T);  // Берём дробную часть числа Frac(1.7)=0.7

    // Смещаем на -1, и проверяем выход за приделы массива. Для таких элементов берём конечный элемент.
    if i-1<=0 then a0:=InR[0]
       else  a0:=InR[i-1];
    a1:=InR[i];
    if i+1>=OldLength then a2:=InR[OldLength-1]
       else  a2:=InR[i+1];
    if i+2>=OldLength then a3:=InR[OldLength-1]
       else  a3:=InR[i+2];

    // Так как сместили на 1, то надо прибавить.
    OutR[j]:=InterpolationCubic(F+1, a0, a1, a2, a3);
    end;

OutA:=ToSmallint(OutR);
end;

Для двух-мерного случая вначале интеполируете по строчкам. Потом по столбцам.