Правильная реализация кубической сплайн-интерполяции

Извините, но Ваш исходный код действительно нечитаемый беспорядок для меня, поэтому я придерживаюсь теории. Вот несколько советов:

1. SPLINE кубики

   SPLINE это не интерполяция, а приближение чтобы использовать их, вам не нужно никакого деривации. Если у вас есть десять очков: p0,p1,p2,p3,p4,p5,p6,p7,p8,p9 затем кубический сплайн начинается / заканчивается тройной точкой. Если вы создаете функцию «рисовать» СПЛАЙН Затем исправьте кубическую кривую, чтобы обеспечить непрерывность последовательности вызовов:

       spline(p0,p0,p0,p1);
   spline(p0,p0,p1,p2);
   spline(p0,p1,p2,p3);
   spline(p1,p2,p3,p4);
   spline(p2,p3,p4,p5);
   spline(p3,p4,p5,p6);
   spline(p4,p5,p6,p7);
   spline(p5,p6,p7,p8);
   spline(p6,p7,p8,p9);
   spline(p7,p8,p9,p9);
   spline(p8,p9,p9,p9);
   

   не забывайте, что кривая SPLINE для p0,p1,p2,p3 нарисовать только кривую «между» p1,p2 !!!

2. Безье кубики

   4-точечный Безье Коэффициенты могут быть вычислены следующим образом:

       a0=                              (    p0);
   a1=                    (3.0*p1)-(3.0*p0);
   a2=          (3.0*p2)-(6.0*p1)+(3.0*p0);
   a3=(    p3)-(3.0*p2)+(3.0*p1)-(    p0);
   

3. интерполирование

   чтобы использовать интерполяцию, вы должны использовать интерполяционные полиномы. Есть много, но я предпочитаю использовать кубики … похоже на Безье / СПЛАЙН / NURBS … так

   • p(t) = a0+a1*t+a2*t*t+a3*t*t*t где t = <0,1>

   Осталось только вычислить a0,a1,a2,a3, У вас есть 2 уравнения (p(t) и его вывод по t) и 4 балла из набора данных. Вы также должны обеспечить непрерывность … Таким образом, первый вывод для точек соединения должен быть одинаковым для соседних кривых (t=0,t=1). Это приводит к системе 4 линейных уравнений (использование t=0 а также t=1). Если вы получите его, он создаст простое уравнение, зависящее только от координат входной точки:

       double  d1,d2;
   d1=0.5*(p2-p0);
   d2=0.5*(p3-p1);
   a0=p1;
   a1=d1;
   a2=(3.0*(p2-p1))-(2.0*d1)-d2;
   a3=d1+d2+(2.0*(-p2+p1));
   

   последовательность вызовов такая же, как для СПЛАЙН

[Заметки]

1. Разница между интерполяцией и аппроксимацией заключается в том, что:

   Кривая интерполяции всегда проходит через контрольные точки, но многочлены высокого порядка имеют тенденцию колебаться, и приближение только приближается к контрольным точкам (в некоторых случаях может пересекать их, но обычно нет).

2. переменные:

   • a0,... p0,... являются векторами (число измерений должно соответствовать входным точкам)
   • t является скалярным

3. рисовать кубики из коэффициентов a0,..a3

   просто сделайте что-то вроде этого:

       MoveTo(a0.x,a0.y);
   for (t=0.0;t<=1.0;t+=0.01)
   {
   tt=t*t;
   ttt=tt*t;
   p=a0+(a1*t)+(a2*tt)+(a3*ttt);
   LineTo(p.x,p.y);
   }
   

Увидеть сплайн-интерполяция, хотя они дают только полезный пример 3х3. Для большего количества точек выборки, скажем, N + 1 перечислил x[0..N] со значениями y[0..N] нужно решить следующую систему для неизвестного k[0..N]

              2*    c[0]    * k[0] + c[0] * k[1] == 3*   d[0];
c[0] * k[0] + 2*(c[0]+c[1]) * k[1] + c[1] * k[2] == 3*(d[0]+d[1]);
c[1] * k[1] + 2*(c[1]+c[2]) * k[2] + c[2] * k[3] == 3*(d[1]+d[2]);
c[2] * k[2] + 2*(c[2]+c[3]) * k[3] + c[3] * k[4] == 3*(d[2]+d[3]);
...
c[N-2]*k[N-2]+2*(c[N-2]+c[N-1])*k[N-1]+c[N-1]*k[N] == 3*(d[N-2]+d[N-1]);
c[N-2]*k[N-1] + 2*c[N-1]        *  k[N]          == 3*   d[N-1];

где

c[k]=1/(x[k+1]-x[k]); d[k]=(y[k+1]-y[k])*c[k]*c[k];

Эту проблему можно решить с помощью итерации Гаусса-Зейделя (не было ли она изобретена именно для этой системы?) Или вашего любимого космического решателя Крылова.