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 zu einer Gesamtkurve mit weichen Übergängen an den Nahtstellen zusammengesetzt werden.
Die n Stützpunkte seien f (ti) , die resultierenden Intervalle
seien
[ti,ti + 1] , die für das Intervall i
zuständige Funktion sei fi(t) ,
1 i
n - 1 .
Jedes fi hat die Form
Es muß gelten
fi(ti + 1) | = | fi + 1(ti + 1) | für | i = 1,...,n - 2 | gleicher Wert |
fi'(ti + 1) | = | fi + 1'(ti + 1) | für | i = 1,...,n - 2 | gleiche Steigung |
fi''(ti + 1) | = | fi + 1''(ti + 1) | für | i = 1,...,n - 2 | gleiche Steigungsänderung |
Nach Vorgabe von a1 und an - 1 entsteht ein Gleichungssystem mit n - 2 Gleichungen und n - 2 Unbekannten a2,a3,...,an - 2 .
Die Werte von a1 und an - 1 können z.B. ermittelt werden durch Festlegung f''(t1) = f''(tn) = 0 (natürliche Spline-Interpolation).
Für die Folge ti 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)) .
/*************************************************************************************/ /* 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 Interpol.werte cubic_spline(point_cnt, t, y, spline_size, calcfy); // bestimme 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]; }