Es muß gelten
fi(ti + 1) | = | fi + 1(ti + 1) | für | i = 1,...,n - 2 | stetig |
fi'(ti + 1) | = | fi + 1'(ti + 1) | für | i = 1,...,n - 2 | differenzierbar |
fi''(ti + 1) | = | fi + 1''(ti + 1) | für | i = 1,...,n - 2 | keine Steigungsänderung |
Nach Vorgabe von a1 und an - 1 entsteht ein Gleichungssystem
mit n - 3 Gleichungen und n - 3 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)) .
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];
}