Kubische Splines

Um die Punkte in einem Koordinatensystem miteinander möglichst gut zu verbinden ohne scharfe Kanten zu erhalten, können Polynome 3ten Grades verwendet werden.

Jede Verbindung zwischen 2 Punkten erhält ein Polynom 3ten Grades (Spline). Dieses Polynom stellt die Verbindung zwischen den beiden Punkten möglichst gut dar.

Es gibt hierzu noch viele mathematische Definitionen (2te Ableitung, stetig differenzierbar…) die ich bei Gelegenheit nachtragen werde ;)

Fürs erste zusammengefasst, müssen für die Ermittlung der kubischen Splines folgende Schritte durchgeführt werden:

  • Die 4 Koeffizienten für alle Polynome 3ten Grades berechnen. (ak, bk, ck, dk).

    • alle aks erhält man direkt aus den yk’s

    • die cks bestimmt man durch das Lösen eines linearen Gleichungssystems mit dem Thomas-Algorithmus (auch genannt TDMA/vereinfachter Gauss Algorithmus). Dieser Algorithmus ist speziell für das Lösen von Tridiagonal-Matrizen geeignet. Also Matrizen die nur auf 3 Diagonalen definiert sind.

    • Die bks und dks kann man nun über die cks berechnen.

  • Jeder Übergang zwischen zwei Punkten hat nun ein Polynom 3ten Grades zugewiesen bekommen. Dieses Polynom kann mittels dem Horner Schema für beliebige x die zwischen diesen beiden Punkten liegen gelöst werden.

  • Zeichnen der Polynome für alle x über einen Bereich des Koordinatensystems

Es folgt der Java Code:

import java.applet.Applet;  
import java.awt.BorderLayout;  
import java.awt.Color;  
import java.awt.Graphics;  
import java.awt.Graphics2D;  
import java.awt.event.MouseEvent;  
import java.awt.geom.Point2D;  
import java.awt.geom.Rectangle2D;  
import java.awt.geom.Point2D.Double;  
import java.util.ArrayList;  
import java.util.Collections;  
import java.util.LinkedList;  
import java.util.List;

import javax.swing.JButton;  
import javax.swing.JPanel;  
import javax.swing.event.MouseInputAdapter;

public class SplineInterpolation extends Applet {  
  public class Punkt extends Point2D.Double implements Comparable<Punkt> {

    public Punkt(double d, double e) {
      super(d, e);
    }

    @Override
    public int compareTo(Punkt arg0) {
      if (arg0 == null) {
        throw new IllegalArgumentException();
      }

      if (getX() < arg0.getX())
        return -1;
      else if (getX() > arg0.getX())
        return 1;
      else
        return 0;
    }

  }

  private static final long serialVersionUID = -4431496457566318110L;
  private List<Punkt> stuetzpunkte = new ArrayList<Punkt>();
  private List<Punkt> splinepunkte = new ArrayList<Punkt>();
  private List<String> polynoms = new ArrayList<String>();

  private double[] x; // Stützstellen (x-Komponente)
  private double[] y; // Stützstellen (y-Komponente)
  private double[] a; // Koeffizienten ai
  private double[] b; // Koeffizienten bi
  private double[] c; // Koeffizienten ci
  private double[] d; // Koeffizienten di
  private double[] h; // h = (x-xi)

  private boolean clear = false;

  public void init() {
    MouseInputAdapter mouseListener = new MouseInputAdapter() {

      public synchronized void mouseClicked(MouseEvent event) {

        if (event.getSource() instanceof JButton) {
          if (stuetzpunkte.size() >= 2) {
            clear = true;
            startSplineInterpolation();
          }
        } else {
          stuetzpunkte.add(new Punkt(event.getPoint().getX() / 20,
              (200 - event.getPoint().getY()) / 20));
        }
        repaint();
        return;
      }
    };

    setSize(750, 600);
    this.addMouseListener(mouseListener);
    this.setLayout(new BorderLayout());

    JPanel rightPanel = new JPanel();
    JButton button = new JButton();
    button.setText("Start Spline Interpolation");
    button.addMouseListener(mouseListener);
    rightPanel.add(button);
    this.add(rightPanel, BorderLayout.NORTH);

  }

  public void startSplineInterpolation() {
    Collections.sort(stuetzpunkte);
    createSplines((Double[]) stuetzpunkte.toArray(new Double[stuetzpunkte
                                                             .size()]));
    for (int i = 1; i <= 29; i++) {

      for (double j = 0; j <= 1; j = j + 0.1) {
        double value = apply(i + j);
        Punkt pt = new Punkt(i + j, value);
        splinepunkte.add(pt);
      }
    }
  }

  @Override
  public void paint(Graphics g) {
    Graphics2D g2d = (Graphics2D) g;
    super.paint(g2d);

    // Gitterlinien malen
    for (int i = 0; i < 31; i++) {
      g2d.setColor(Color.BLACK);
      g2d.drawString(String.valueOf(i), 0 + i * 20, 200);
      g2d.setColor(Color.GRAY);
      g2d.drawLine(0 + i * 20, 35, 0 + i * 20, 200);
    }

    for (int i = 0; i < 8; i++) {
      g2d.setColor(Color.BLACK);
      g2d.drawString(String.valueOf(i), 0, 200 - i * 20);
      g2d.setColor(Color.GRAY);
      g2d.drawLine(0, 200 - i * 20, 600, 200 - i * 20);
    }

    // farbe ändern
    g.setColor(Color.RED);

    // punkte malen
    for (Point2D.Double p : splinepunkte) {
      g2d.draw(new Rectangle2D.Double(p.x * 20, 200 - p.y * 20, 1, 1));
    }

    for (Point2D.Double p : stuetzpunkte) {
      g2d.setColor(Color.BLUE);
      g2d.draw(new Rectangle2D.Double(p.x * 20 - 2, 200 - p.y * 20 - 2,
          4, 4));
    }

    // String-Liste der Polynome malen
    for (int i = 0; i < polynoms.size(); i++) {
      String polynom = polynoms.get(i);
      g2d.setColor(Color.BLUE);
      g2d.drawString(polynom, 20, 240 + i * 20);
    }

    if (clear) {
      stuetzpunkte.clear();
      splinepunkte.clear();
      polynoms.clear();
      clear = false;
    }
  }

  public void createSplines(Point2D.Double[] knots) {
    // Übergabeparameter: knots = gegebene Stützstellen

    this.x = new double[knots.length]; // Stützstellen (x-Koordinate)
    this.y = new double[knots.length]; // Stützstellen (y-Koordinate)
    this.a = new double[knots.length]; // Koeffizienten a[0..n]
    this.b = new double[knots.length]; // Koeffizienten b[0..n]
    this.c = new double[knots.length]; // Koeffizienten c[0..n]
    this.d = new double[knots.length]; // Koeffizienten d[0..n]
    this.h = new double[knots.length]; // Koeffizienten h[0..n]

    for (int i = 0; i <= knots.length - 1; i++) {
      x[i] = knots[i].x; // Stützstellen (x-Komponenten)
      y[i] = knots[i].y; // Stützstellen (y-Komponenten)
    }

    qSplines();
  }

  public double apply(double arg) {
    // Index der passenden Splinefunktion si(x) finden
    int i = 0;
    for (i = 0; i < this.x.length - 1; i++) {
      if (arg < this.x[i + 1])
        break;
    }

    // si(arg) auswerten (Hornerschema)
    double h = (arg - x[i]);

    double koeff1 = d[i] * h + c[i];
    double koeff2 = koeff1 * h + b[i];
    double koeff3 = koeff2 * h + a[i];
    return koeff3;
  }

  public void qSplines() {

    // alle hks und alle aks bestimmen
    for (int i = 0; i < x.length - 1; i++) {
      h[i] = x[i + 1] - x[i];
      a[i] = y[i];
    }

    // die Matrix zur Bestimmung der cks aufbauen
    double[] m1 = new double[x.length];
    double[] m2 = new double[x.length];
    double[] m3 = new double[x.length];
    double[] l = new double[x.length];

    for (int i = 0; i < x.length - 2; i++) {
      m1[i + 1] = h[i + 1];
      m2[i] = 2 * (h[i + 1] + h[i + 2]);
      m3[i + 1] = h[i + 1];
      l[i] = 3 * (((1 / h[i + 1]) * (y[i + 2] - y[i + 1])) - ((1 / h[i]) * (y[i + 1] - y[i])));
    }

    m1[0] = 0;
    m3[0] = 0;
    m1[m1.length - 1] = 0;
    m3[m3.length - 1] = 0;

    // die cks durch Auflösen der Matrix nach dem Thomas Algorithmus
    // (anderer Name: TDMA / vereinfachter Gauss Algorithmus für Matrizen
    // mit 3 Diagonalen) bestimmen.
    TridiagonalSolve(m1, m2, m3, c, l);

    // die erhaltenen cks müssen noch so verschoben werden dass Platz für
    // die erste und letzte 0 ist.
    for (int i = x.length - 1; i > 0; i--)
      c[i] = c[i - 1];

    // erstes und letztes ck ist immer 0.
    c[0] = 0;
    c[x.length - 1] = 0;

    for (int i = 0; i < x.length - 1; i++) {
      b[i] = ((y[i + 1] - y[i]) / h[i])
          - ((h[i] * (c[i + 1] + 2 * c[i])) / 3);
      d[i] = (1 / (3 * h[i])) * (c[i + 1] - c[i]);
    }

    for (int i = 0; i < x.length; i++) {
      polynoms.add("S" + i + "(x): " + a[i] + " + " + b[i] + "*(x-"
          + x[i] + ") + " + c[i] + "*(x-" + x[i] + ")^2 +" + d[i]
              + "*(x-" + x[i] + ")^3)");
    }
  }

  /**
   * der hier implementierte Thomas Algorithmus stammt aus dem Wikipedia
   * Eintrag zum Thomas Algorithmus:
   * http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
   */
  static void TridiagonalSolve(double[] a, double[] b, double[] c,
      double[] x, double[] d) {
    /* Fills solution into x. Warning: will modify c and d! */
    int i;
    int n = d.length - 1;

    /* Modify the coefficients. */
    c[0] = c[0] / b[0];
    d[0] = c[0] / b[0];
    for (i = 1; i < n; i++) {
      double id = (b[i] - c[i - 1] * a[i]);
      c[i] = c[i] / id;
      d[i] = (d[i] - d[i - 1] * a[i]) / id;
    }

    /* Now back substitute. */
    x[n - 1] = d[n - 1];
    for (i = n - 2; i >= 0; i--)
      x[i] = d[i] - c[i] * x[i + 1];
  }
}
comments powered by Disqus