Korrekte Implementierung der kubischen Spline-Interpolation

Ich habe einen der vorgeschlagenen Algorithmen benutzt, aber die Ergebnisse sind sehr schlecht.

Ich habe den Wiki-Algorithmus implementiert

in Java (Code unten). x(0) ist points.get(0) , y(0) ist values[points.get(0)] , α ist alfa und μ ist mi . Der Rest ist der gleiche wie im Wiki-Pseudocode.

  public void createSpline(double[] values, ArrayList points){ a = new double[points.size()+1]; for (int i=0; i <points.size();i++) { a[i] = values[points.get(i)]; } b = new double[points.size()]; d = new double[points.size()]; h = new double[points.size()]; for (int i=0; i<points.size()-1; i++){ h[i] = points.get(i+1) - points.get(i); } alfa = new double[points.size()]; for (int i=1; i <points.size()-1; i++){ alfa[i] = (double)3 / h[i] * (a[i+1] - a[i]) - (double)3 / h[i-1] *(a[i+1] - a[i]); } c = new double[points.size()+1]; l = new double[points.size()+1]; mi = new double[points.size()+1]; z = new double[points.size()+1]; l[0] = 1; mi[0] = z[0] = 0; for (int i =1; i0; j--) { c[j] = z[j] - mi[j]*c[j-1]; b[j] = (a[j+1]-a[j]) - (h[j] * (c[j+1] + 2*c[j])/(double)3) ; d[j] = (c[j+1]-c[j])/((double)3*h[j]); } for (int i=0; i<points.size()-1;i++){ for (int j = points.get(i); j<points.get(i+1);j++){ // fk[j] = values[points.get(i)]; functionResult[j] = a[i] + b[i] * (j - points.get(i)) + c[i] * Math.pow((j - points.get(i)),2) + d[i] * Math.pow((j - points.get(i)),3); } } } 

Das Ergebnis, das ich bekomme, ist das Folgende:

Bildbeschreibung hier eingeben

aber es sollte ähnlich sein:

Bildbeschreibung hier eingeben


Ich versuche auch, den Algorithmus auf andere Weise zu implementieren: http://www.geos.ed.ac.uk/~yliu23/docs/level_spline.pdf

Zuerst zeigen sie, wie man linearen Spline macht, und es ist ziemlich einfach. Ich erstelle functionen, die A und B Koeffizienten berechnen. Dann erweitern sie den linearen Spline, indem sie eine zweite Ableitung hinzufügen. C und D Koeffizienten sind ebenfalls einfach zu berechnen.

Aber die Probleme beginnen, wenn ich versuche, die zweite Ableitung zu berechnen. Ich verstehe nicht, wie sie sie berechnen.

Also habe ich nur lineare Interpolation implementiert. Das Ergebnis ist:

Bildbeschreibung hier eingeben

Kann jemand den ersten Algorithmus korrigieren oder mir erklären, wie man die zweite Ableitung im zweiten Algorithmus berechnet?

Cubic B-Spline wurde kürzlich in einer Reihe von Arbeiten von Unser, Thévenaz et al. , siehe unter anderem

M. Unser, A. Aldroubi, M. Eden, “Schnelle B-Spline-Transformation für kontinuierliche Bilddarstellung und Interpolation”, IEEE Trans. Muster Anal. Maschinenintell. , vol. 13, n. 3, S. 277-285, März 1991.

M. Unser, “Splines, eine perfekte Lösung für Signal- und Bildverarbeitung”, IEEE Signal Proc. Mag. , S. 22-38, Nov. 1999.

und

P. Thévenaz, T. Blu, M. Unser, “Interpolation Revisited”, IEEE Trans. über medizinische Bildgebung , vol. 19, nein. 7, S. 739-758, Juli 2000.

Hier sind einige Richtlinien.

Was sind Splines?

Splines sind stückweise Polynome, die glatt miteinander verbunden sind. Für einen Spline vom Grad n ist jedes Segment ein Polynom vom Grad n . Die Teile sind so verbunden, dass der Spline bis zu seiner Ableitung vom Grad n-1 an den Knoten , nämlich den Verbindungspunkten der Polynomstücke, stetig ist.

Wie können Splines konstruiert werden?

Der Spline nullter Ordnung ist der folgende

Bildbeschreibung hier eingeben

Alle anderen Splines können als konstruiert werden

Bildbeschreibung hier eingeben

wo die Faltung n-1 Mal genommen wird.

Kubische Splines

Die beliebtesten Splines sind kubische Splines, deren Ausdruck ist

Bildbeschreibung hier eingeben

Spline-Interpolationsproblem

Bei einer gegebenen function f(x) die an den diskreten ganzzahligen Punkten k abgetastet wird, dient das Spline-Interpolationsproblem dazu, eine Approximation s(x) bis f(x) zu bestimmen, die in der folgenden Weise ausgedrückt wird

Bildbeschreibung hier eingeben

wo die ck Interpolationskoeffizienten sind und s(k) = f(k) .

Spline-Vorfilterung

Leider, ab n=3 ,

Bildbeschreibung hier eingeben

so dass die ck ‘s nicht die Interpolationskoeffizienten sind. Sie können bestimmt werden, indem das lineare Gleichungssystem getriggers wird, das durch Erzwingen von s(k) = f(k) , nämlich

Bildbeschreibung hier eingeben

Eine solche Gleichung kann in eine Faltungsform umformuliert und in dem transformierten z Raum als getriggers werden

Bildbeschreibung hier eingeben

woher

Bildbeschreibung hier eingeben

Entsprechend,

Bildbeschreibung hier eingeben

Das Vorgehen auf diesem Weg ist immer vorzuziehen, als die Lösung eines linearen Gleichungssystems, beispielsweise durch LU Zerlegung.

Die Lösung für die obige Gleichung kann bestimmt werden, indem dies bemerkt wird

Bildbeschreibung hier eingeben

woher

Bildbeschreibung hier eingeben

Die erste Fraktion ist repräsentativ für einen kausalen Filter , während die zweite für ein antikausales Filter repräsentativ ist. Beide sind in den folgenden Abbildungen dargestellt.

Kausalfilter

Bildbeschreibung hier eingeben

Antikausaler Filter

Bildbeschreibung hier eingeben

In der letzten Figur,

Bildbeschreibung hier eingeben

Die Ausgabe der Filter kann durch die folgenden rekursiven Gleichungen ausgedrückt werden

Bildbeschreibung hier eingeben

Die obigen Gleichungen können getriggers werden, indem zuerst “Anfangsbedingungen” für c- und c+ . Bei Annahme einer periodischen, gespiegelten Eingabefolge fk so, dass

Bildbeschreibung hier eingeben

dann kann gezeigt werden, dass die Anfangsbedingung für c+ wie folgt ausgedrückt werden kann

Bildbeschreibung hier eingeben

während die Anfangsbedingung für c+ ausgedrückt werden kann als

Bildbeschreibung hier eingeben

Tut mir leid, aber Ihr Quellcode ist wirklich ein unlesbares Durcheinander, also bleibe ich bei der Theorie. Hier sind einige Tipps:

  1. SPLINE-Würfel

    SPLINE ist keine Interpolation, sondern eine Approximation , um sie zu verwenden. Sie benötigen keine Ableitung. Wenn Sie zehn Punkte haben: p0,p1,p2,p3,p4,p5,p6,p7,p8,p9 dann beginnt / endet der kubische Spline mit dem Tripelpunkt. Wenn Sie eine function erstellen, um den kubischen SPLINE- Kurven-Patch zu “zeichnen”, dann wird die Call-Sequenz folgendermaßen aussehen, um Kontinuität zu gewährleisten:

      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); 

    Vergessen Sie nicht, dass die SPLINE-Kurve für p0,p1,p2,p3 nur die Kurve ‘zwischen’ p1,p2 !!!

  2. BEZIER Kubik

    4-Punkt- BEZIER- Koeffizienten können wie folgt berechnet werden:

      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. Interpolation

    Um die Interpolation zu verwenden, müssen Interpolationspolynome verwendet werden. Es gibt viele da draußen, aber ich bevorzuge kubische … ähnlich wie BEZIER / SPLINE / NURBS … so

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

    Das einzige, was noch zu tun bleibt, ist a0,a1,a2,a3 berechnen. Sie haben 2 Gleichungen ( p(t) und seine Ableitung durch t ) und 4 Punkte aus dem Datensatz. Sie müssen auch die Stetigkeit sicherstellen … Also muss die erste Ableitung für Join-Punkte für benachbarte Kurven gleich sein ( t=0,t=1 ). Dies führt zu einem System von 4 linearen Gleichungen (verwenden Sie t=0 und t=1 ). Wenn Sie es ableiten, erzeugt es eine einfache Gleichung, die nur von den Koordinaten des Eingabepunkts abhängt:

      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)); 

    Die Aufrufreihenfolge ist dieselbe wie für SPLINE

[Anmerkungen]

  1. Der Unterschied zwischen Interpolation und Approximation ist folgender:

    Die Interpolationskurve geht immer durch die Kontrollpunkte, aber die Polynome höherer Ordnung tendieren dazu, zu oszillieren, und die Approximation nähert sich nur den Kontrollpunkten (in einigen Fällen können sie sich kreuzen, aber normalerweise nicht).

  2. Variablen:

    • a0,... p0,... sind Vektoren (die Anzahl der Dimensionen muss mit den Eingabepunkten übereinstimmen)
    • t ist skalar
  3. Würfel aus den Koeffizienten a0,..a3 zu zeichnen

    mach einfach so etwas:

      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(px,py); } 

Siehe Spline-Interpolation , obwohl sie nur ein brauchbares 3×3-Beispiel ergeben. Für mehr Beispielpunkte, sagen wir N + 1, x[0..N] mit Werten y[0..N] muss man das folgende System für das unbekannte 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]; 

woher

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

Dies kann mit der Gauß-Seidel-Iteration (wurde diese nicht genau für dieses System erfunden?) Oder mit Ihrem bevorzugten Krylov-Raumlöser getriggers werden.