7.2 Splines

Gegeben n Stützpunkte.
Wähle n - 1 Kurven kleiner Ordnung für den Verlauf innerhalb der n - 1 Intervalle.
Kurven erster Ordnung (Geraden) ergeben nur den Linienzug.
Kurven zweiter Ordnung (Parabeln) sind an den Intervallgrenzen nur einmal differenzierbar.
Kurven dritter Ordnung mit stetigen ersten und zweiten Ableitungen an den Intervallgrenzen (kubische Splines) können wie folgt bestimmt werden:

Gesucht:

Approximation einer Kurve f (t) , die durch n vorgegebene Punkte laufen soll.

Bewertung:

f ist in parametrischer Form gegeben (d.h. x = g(t),y = h(t) ), da die explizite Form y = f (x) pro x -Wert nur einen y -Wert erlaubt.

Die n Stützpunkte seien f (t_i) , die resultierenden Intervalle seien [t_i,t_{i + 1}] , die für das Intervall i zuständige Funktion sei f_i(t) , 1 i n - 1 . Jedes f_i hat die Form

Es muß gelten

f_i(t_{i + 1}) = f_{i + 1}(t_{i + 1}) für i = 1,...,n - 2 stetig

f_i'(t_{i + 1}) = f_{i + 1}'(t_{i + 1}) für i = 1,...,n - 2 differenzierbar

f_i''(t_{i + 1}) = f_{i + 1}''(t_{i + 1}) für i = 1,...,n - 2 keine Steigungsänderung

Nach Vorgabe von a₁ und a_{n - 1} entsteht ein Gleichungssystem mit n - 3 Gleichungen und n - 3 Unbekannten a₂,a₃,...,a_{n - 2} .
Die Werte von a₁ und a_{n - 1} können z.B. ermittelt werden durch Festlegung f''(t₁) = f''(t_n) = 0 (natürliche Spline-Interpolation).
Für die Folge t_i reicht jede monoton wachsende Folge. Geeignet ist z.B. die Folge der euklidischen Abstände zwischen den Stützpunkten. Statt den Abstand d = zu berechnen, verwendet man die Näherungsformel · (dx + dy + 2 · max (dx,dy)) .

Vom Spline-Algorithmus erzeugte Kurveninterpolation mit 13 Stützpunkten und 65 Interpolationspunkten (pro Intervall)

/****************************************************************************************/
/*               Spline-Interpolation durch Polynome 3. Grades                          */
/****************************************************************************************/

private static double[] calcfx = new double[MAX_POINTS];       	// x-Funktionswerte
private static double[] calcfy = new double[MAX_POINTS];       	// y-Funktionswerte
private static double[] x = new double[MAX_POINTS];    	// x-Koordinaten der Stuetzpunkte
private static double[] y = new double[MAX_POINTS];     // y-Koordinaten der Stuetzpunkte
private static double[] t = new double[MAX_POINTS];     // t-Werte der Stuetzpunkte

void besetze_arrays(                          // berechnet Intervallgrenzen
              int point_cnt,                  // Anzahl der Polygonpunkte
              Point punkt[])                  // Punktliste
{
      int i;
      double ax,ay,dd;
 
      for (i=0; i<point_cnt; i++) {           // besetze x- und y-Koordinaten
              x[i] = (double) punkt[i+1].x;
              y[i] = (double) punkt[i+1].y;}
      t[0] = 0.0;                             // besetze t-Werte mit dem
      for (i=1; i<point_cnt; i++) {           // Euklidischen Abstand der Stuetzpunke
              ax=Math.abs(x[i]-x[i-1]);
              ay=Math.abs(y[i]-y[i-1]);
              dd = ax + ay;
              if (ax > ay) dd = dd+2*ax;      // Naeherungsformel fuer Euklid:
                      else dd = dd+2*ay;      // Abstand = 1/3*(dx+dy+2*max{dx,dy})
              t[i] = t[i-1]+dd/3;
      }
}
 
void curve_fitting(                           // berechnet die Approximation der Kurve
      int point_cnt,                          // fuer point_cnt Stuetzpunkte
      Point punkt[],                          // uebergeben im Array punkt
      int spline_size)                        // fuer spline_size Interpolationswerte
{
      int i,n=0;
      besetze_arrays(point_cnt, punkt);                    // erzeuge Intervallgrenzen und
                                                           // Gleitkommaversion von punkt
      cubic_spline(point_cnt, t, x, spline_size, calcfx);  // bestimme die Interpol.werte
      cubic_spline(point_cnt, t, y, spline_size, calcfy);  // bestimme die Interpol.werte
      for (i=0; i<spline_size; i++)                        // und uebertrage sie
        spline_punkt[i+1] = new Point((int)calcfx[i],(int)calcfy[i]);
}
void cubic_spline(                    // loest Gleichungssystem zur Bestimmung
                                      // der Koeffizienten a,b,c
              int n,                  // Anzahl der Stuetzpunkte
              double t[], double f[], // Stuetzpunkte in der Form (t[i],f[t[i]])
              int spline_size,        // Anzahl der Interpolationswerte
              double calcf[])         // Ausgabearray: Interpolationswerte
{ 
      int i,j;
      double a[] = new double[MAX_POINTS];
      double b[] = new double[MAX_POINTS];
      double c[] = new double[MAX_POINTS];
      double delta_t[] = new double[MAX_POINTS];
      double D[] = new double[MAX_POINTS];
      double m[] = new double[MAX_POINTS];
      double k[] = new double[MAX_POINTS];
      double bh,dh,e,h,wt,dt;
      for (i=1; i<n; i++) {
              delta_t[i] = t[i]-t[i-1];
              D[i] = (f[i]-f[i-1])/delta_t[i];
      }
      m[0] = delta_t[2];
      delta_t[0] = t[2]-t[0];
      h = delta_t[1];
      k[0] = ((h + 2*delta_t[0])*D[1]*delta_t[2]+h*h*D[2])/delta_t[0];
      for (i=1; i<(n-1); i++) {
              h = -delta_t[i+1]/m[i-1];
              k[i] = h*k[i-1]+3*(delta_t[i]*D[i+1]+delta_t[i+1]*D[i]);
              m[i] = h *delta_t[i-1] + 2* (delta_t[i] + delta_t[i+1]);
      }
      h = t[n-1]-t[n-3];
      dh = delta_t[n-1];
      k[n-1] = ((dh+h+h)*D[n-1]*delta_t[n-2] +  dh*dh*D[n-2])/h;
      h = -h/m[n-2];
      m[n-1] = delta_t[n-2];
      m[n-1] = h*delta_t[n-2] + m[n-1];
      a[n-1] = (h*k[n-2]+k[n-1])/m[n-1];
      for (i=n-2; i>=0; i--) a[i] = (k[i]-delta_t[i]*a[i+1])/m[i];
      for (i=1; i<n; i++) {
              dh = D[i]; bh = delta_t[i];
              e = a[i-1]+a[i]-dh-dh;
              b[i-1] = 2*(dh-a[i-1]-e)/bh;
              c[i-1] = 6*e/(bh*bh);
      }
      wt = 0; j=0;
      dt = t[n-1]/(double)(spline_size-1);
      for (i=0; i<spline_size; i++) {
              while ((t[j+1]<wt)&&(j<spline_size)) j++;
              h = wt - t[j];
              calcf[i] = f[j]+h*(a[j]+h*(b[j]+h*c[j]/3)/2);
              wt = wt+dt;
      }
      calcf[spline_size-1]= f[n-1];
}