Approximation von Daten mit einer kubischen Bezierkurve aus mehreren Segmenten und einer Distanz sowie einer Krümmungsbeschränkung

Ich habe einige Geodaten (das Bild unten zeigt den Verlauf eines Flusses als rote Punkte), die ich mit einer kubischen Bezierkurve mit mehreren Segmenten approximieren möchte. Durch andere Fragen zu Stackoverflow hier und hier bin ich auf den Algorithmus von Philip J. Schneider von “Graphics Gems” gestoßen. Ich habe es erfolgreich implementiert und kann berichten, dass es selbst mit Tausenden von Punkten sehr schnell ist. Leider hat diese Geschwindigkeit einige Nachteile, nämlich dass die Anpassung ziemlich schlampig durchgeführt wird. Betrachten Sie die folgende Grafik:

Die roten Punkte sind meine Originaldaten und die blaue Linie ist das Multi-Segment-Bezier, das durch den Algorithmus von Schneider erstellt wurde. Wie Sie sehen können, war die Eingabe in den Algorithmus eine Toleranz, die mindestens so hoch ist, wie die grüne Linie anzeigt. Trotzdem erzeugt der Algorithmus eine Bezierkurve, die zu viele scharfe Kurven hat. Sie sehen auch diese unnötig scharfen Kurven im Bild. Es ist leicht, sich eine Bezierkurve mit weniger scharfen Kurven für die angezeigten Daten vorzustellen, während die maximale Toleranzbedingung beibehalten wird (schieben Sie die Bezierkurve einfach ein wenig in Richtung der magentafarbenen Pfeile). Das Problem scheint zu sein, dass der Algorithmus Datenpunkte aus meinen Originaldaten als Endpunkte der einzelnen Bezierkurven auswählt (die magentafarbenen Pfeile weisen auf einige Verdächtige hin). Wenn die Endpunkte der Bezierkurven so eingeschränkt sind, ist es klar, dass der Algorithmus manchmal ziemlich scharfe Krümmungen erzeugt.

Was ich suche, ist ein Algorithmus, der meine Daten mit einer Bezierkurve mit mehreren Segmenten mit zwei Einschränkungen approximiert:

• die mehrsegmentige Bezierkurve darf nie mehr als einen bestimmten Abstand von den Datenpunkten entfernt sein (dies wird durch den Algorithmus von Schneider bereitgestellt)
• Die Bezierkurve mit mehreren Segmenten darf niemals zu scharfe Krümmungen erzeugen. Eine Möglichkeit, dieses Kriterium zu überprüfen, besteht darin, einen Kreis mit dem minimalen Krümmungsradius entlang der Bezierkurve mit mehreren Segmenten zu rollen und zu prüfen, ob er alle Teile der Kurve entlang seines Pfads berührt. Obwohl es scheint, dass es eine bessere Methode gibt, die das Kreuzprodukt der ersten und zweiten Ableitung verwendet

Die Lösungen, die ich gefunden habe und die bessere Passungen erzeugen, funktionieren leider entweder nur für einzelne Bezierkurven (und lassen die Frage aus, wie man gute Start- und Endpunkte für jede Bezierkurve in der Bezierkurve mit mehreren Segmenten findet) oder erlauben keine minimale Krümmungsbeschränkung. Ich denke, dass die minimale Krümmungsbeschränkung hier die schwierige Bedingung ist.

Hier ein weiteres Beispiel (dieses ist handgezeichnet und nicht 100% genau):

Nehmen wir an, dass Abbildung 1 sowohl die Krümmungsbeschränkung (der Kreis muss entlang der gesamten Kurve passen) als auch den maximalen Abstand eines beliebigen Datenpunkts von der Kurve (der zufällig der Radius des Kreises in Grün ist) zeigt. Eine erfolgreiche Annäherung an den roten Pfad in Abbildung 2 ist blau dargestellt. Diese Annäherung berücksichtigt die Krümmungsbedingung (der Kreis kann innerhalb der gesamten Kurve rollen und berührt sie überall) sowie die Abstandsbedingung (in grün dargestellt). Abbildung 3 zeigt eine andere Annäherung an den Pfad. Es erfüllt zwar die Abstandsbedingung, aber es ist klar, dass der Kreis nicht mehr in die Krümmung passt. Abbildung 4 zeigt einen Weg, der mit den gegebenen Einschränkungen nicht angenähert werden kann, weil er zu spitz ist. Dieses Beispiel soll veranschaulichen, dass es zur korrekten Annäherung einiger spitzer Kurven im Pfad erforderlich ist, dass der Algorithmus Kontrollpunkte auswählt, die nicht Teil des Pfads sind. Abbildung 3 zeigt, dass die Krümmungsbeschränkung nicht mehr erfüllt werden kann, wenn Kontrollpunkte entlang des Pfads ausgewählt wurden. Dieses Beispiel zeigt auch, dass der Algorithmus bei einigen Eingaben beendet werden muss, da es nicht möglich ist, ihn mit den gegebenen Einschränkungen zu approximieren.

Gibt es eine Lösung für dieses Problem? Die Lösung muss nicht schnell sein. Wenn es einen Tag dauert, 1000 Punkte zu verarbeiten, dann ist das in Ordnung. Die Lösung muss auch nicht in dem Sinne optimal sein, dass sie zu einer Anpassung nach der Methode der kleinsten Quadrate führen muss.

Am Ende werde ich das in C und Python implementieren, aber ich kann auch die meisten anderen Sprachen lesen.

Ich habe die Lösung gefunden, die meine Kriterien erfüllt. Die Lösung besteht darin, zuerst einen B-Spline zu finden, der die Punkte im Sinne der kleinsten Quadrate annähert, und diesen Spline dann in eine Bezierkurve mit mehreren Segmenten umzuwandeln. B-Splines haben den Vorteil, dass sie im Gegensatz zu Bezier-Kurven nicht durch die Kontrollpunkte verlaufen und eine Möglichkeit bieten, eine gewünschte “Glätte” der Annäherungskurve zu spezifizieren. Die erforderliche Funktionalität zum Generieren eines solchen Splines ist in der FITPACK-Bibliothek implementiert, für die scipy eine Python-Bindung anbietet. Nehmen wir an, ich lese meine Daten in die Listen ein x und ydann kann ich tun:

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate
tck,u = interpolate.splprep([x,y],s=3)
unew = np.arange(0,1.01,0.01)
out = interpolate.splev(unew,tck)
plt.figure()
plt.plot(x,y,out[0],out[1])
plt.show()

Das Ergebnis sieht dann so aus:

Wenn ich die Kurve glatter haben möchte, dann kann ich die erhöhen s Parameter zu splprep. Wenn ich die Annäherung näher an die Daten möchte, kann ich die verringern s Parameter für weniger Glätte. Indem man mehrere durchläuft s Parameter kann ich programmgesteuert einen guten Parameter finden, der den gegebenen Anforderungen entspricht.

Die Frage ist jedoch, wie man dieses Ergebnis in eine Bezier-Kurve umwandelt. Die Antwort drin diese E-Mail von Zachary Pincus. Ich werde seine Lösung hier replizieren, um meine Frage vollständig zu beantworten:

def b_spline_to_bezier_series(tck, per = False):
  """Convert a parametric b-spline into a sequence of Bezier curves of the same degree.

  Inputs:
    tck : (t,c,k) tuple of b-spline knots, coefficients, and degree returned by splprep.
    per : if tck was created as a periodic spline, per *must* be true, else per *must* be false.

  Output:
    A list of Bezier curves of degree k that is equivalent to the input spline. 
    Each Bezier curve is an array of shape (k+1,d) where d is the dimension of the
    space; thus the curve includes the starting point, the k-1 internal control 
    points, and the endpoint, where each point is of d dimensions.
  """
  from fitpack import insert
  from numpy import asarray, unique, split, sum
  t,c,k = tck
  t = asarray
  try:
    c[0][0]
  except:
    # I can't figure out a simple way to convert nonparametric splines to 
    # parametric splines. Oh well.
    raise TypeError("Only parametric b-splines are supported.")
  new_tck = tck
  if per:
    # ignore the leading and trailing k knots that exist to enforce periodicity 
    knots_to_consider = unique(t[k:-k])
  else:
    # the first and last k+1 knots are identical in the non-periodic case, so
    # no need to consider them when increasing the knot multiplicities below
    knots_to_consider = unique(t[k+1:-k-1])
  # For each unique knot, bring it's multiplicity up to the next multiple of k+1
  # This removes all continuity constraints between each of the original knots, 
  # creating a set of independent Bezier curves.
  desired_multiplicity = k+1
  for x in knots_to_consider:
    current_multiplicity = sum(t == x)
    remainder = current_multiplicity%desired_multiplicity
    if remainder != 0:
      # add enough knots to bring the current multiplicity up to the desired multiplicity
      number_to_insert = desired_multiplicity - remainder
      new_tck = insert(x, new_tck, number_to_insert, per)
  tt,cc,kk = new_tck
  # strip off the last k+1 knots, as they are redundant after knot insertion
  bezier_points = numpy.transpose(cc)[:-desired_multiplicity]
  if per:
    # again, ignore the leading and trailing k knots
    bezier_points = bezier_points[k:-k]
  # group the points into the desired bezier curves
  return split(bezier_points, len(bezier_points) / desired_multiplicity, axis = 0)

Also B-Splines, FITPACK, numpy und scipy haben mir den Tag gerettet 🙂

1. Daten polygonisieren

   Finden Sie die Reihenfolge der Punkte, so dass Sie einfach die nächsten Punkte finden und versuchen, sie “durch Linien” zu verbinden. Vermeiden Sie es, zum Ursprungspunkt zurückzukehren

2. Berechnung der Ableitung entlang des Pfades

   Es ist die Richtungsänderung der ‘Linien’, wo Sie das lokale Minimum oder Maximum treffen, dort ist Ihr Kontrollpunkt … Tun Sie dies, um Ihre Eingabedaten zu reduzieren (lassen Sie nur Kontrollpunkte übrig).

3. Kurve

   Verwenden Sie diese Punkte nun als Kontrollpunkte. Ich empfehle dringend das Interpolationspolynom für beide x und y separat zum Beispiel so etwas:

   x=a0+a1*t+a2*t*t+a3*t*t*t
   y=b0+b1*t+b2*t*t+b3*t*t*t
   

   wo a0..a3 werden so berechnet:

   d1=0.5*(p2.x-p0.x);
   d2=0.5*(p3.x-p1.x);
   a0=p1.x;
   a1=d1;
   a2=(3.0*(p2.x-p1.x))-(2.0*d1)-d2;
   a3=d1+d2+(2.0*(-p2.x+p1.x));
   

   • b0 .. b3 werden auf die gleiche Weise berechnet, verwenden aber natürlich y-Koordinaten
   • p0..p3 sind Kontrollpunkte für die kubische Interpolationskurve
   • t =<0.0,1.0> ist Kurvenparameter von p1 zu p2

   Dies stellt sicher, dass Position und erste Ableitung stetig sind (c1) und Sie können auch BEZIER verwenden, aber es wird nicht so gut übereinstimmen.

[edit1] Zu scharfe Kanten sind ein großes Problem

Um es zu lösen, können Sie Punkte aus Ihrem Datensatz entfernen, bevor Sie die Kontrollpunkte erhalten. Ich kann mir jetzt zwei Möglichkeiten vorstellen, es zu tun … wähle, was für dich besser ist

1. Punkte mit zu hoher erster Ableitung aus Datensatz entfernen

   dx/dl oder dy/dl wo x,y sind Koordinaten und l ist die Kurvenlänge (entlang ihres Pfades). Die exakte Berechnung des Krümmungsradius aus der Kurvenableitung ist knifflig

2. Punkte aus dem Datensatz entfernen, die zu einem zu kleinen Krümmungsradius führen

   Berechnen Sie den Schnittpunkt benachbarter Liniensegmente (schwarze Linien) als Mittelpunkt. Senkrechte Achsen wie auf dem Bild (rote Linien), der Abstand davon und der Verbindungspunkt (blaue Linie) ist Ihr Krümmungsradius. Wenn der Krümmungsradius kleiner ist als Ihr Limit, entfernen Sie diesen Punkt …

   

   Wenn Sie jetzt wirklich nur BEZIER-Kubik brauchen, können Sie meine Interpolationskubik wie folgt in BEZIER-Kubik umwandeln:

//  ---------------------------------------------------------------------------
//  x=cx[0]+(t*cx[1])+(tt*cx[2])+(ttt*cx[3]); // cubic x=f
//  ---------------------------------------------------------------------------
//  cubic matrix                           bz4 = it4
//  ---------------------------------------------------------------------------
//  cx[0]=                            (    x0) =                    (    X1)
//  cx[1]=                   (3.0*x1)-(3.0*x0) =           (0.5*X2)         -(0.5*X0)
//  cx[2]=          (3.0*x2)-(6.0*x1)+(3.0*x0) = -(0.5*X3)+(2.0*X2)-(2.5*X1)+(    X0)
//  cx[3]= (    x3)-(3.0*x2)+(3.0*x1)-(    x0) =  (0.5*X3)-(1.5*X2)+(1.5*X1)-(0.5*X0)
//  ---------------------------------------------------------------------------
    const double m=1.0/6.0;
    double x0,y0,x1,y1,x2,y2,x3,y3;
    x0 = X1;           y0 = Y1;
    x1 = X1-(X0-X2)*m; y1 = Y1-(Y0-Y2)*m;
    x2 = X2+(X1-X3)*m; y2 = Y2+(Y1-Y3)*m;
    x3 = X2;           y3 = Y2;

Falls Sie die umgekehrte Konvertierung benötigen, siehe:

• Bezier-Kurve mit Kontrollpunkten innerhalb der Kurve

Die Frage wurde vor langer Zeit gepostet, aber hier ist eine einfache Lösung, die auf splprep basiert und den minimalen Wert von s findet, der es ermöglicht, ein Kriterium für einen minimalen Krümmungsradius zu erfüllen.

Route ist die Menge der Eingabepunkte, wobei die erste Dimension die Anzahl der Punkte ist.

import numpy as np
from scipy.interpolate import splprep, splev

#The minimum curvature radius we want to enforce
minCurvatureConstraint = 2000
#Relative tolerance on the radius
relTol = 1.e-6

#Initial values for bisection search, should bound the solution
s_0 = 0
minCurvature_0 = 0
s_1 = 100000000 #Should be high enough to produce curvature radius larger than constraint
s_1 *= 2
minCurvature_1 = np.float('inf')

while np.abs(minCurvature_0 - minCurvature_1)>minCurvatureConstraint*relTol:
    s = 0.5 * (s_0 + s_1)
    tck, u  = splprep(np.transpose(route), s=s)
    smoothed_route = splev(u, tck)
    
    #Compute radius of curvature
    derivative1 = splev(u, tck, der=1)
    derivative2 = splev(u, tck, der=2) 
    xprim = derivative1[0]
    xprimprim = derivative2[0]
    yprim = derivative1[1]
    yprimprim = derivative2[1]
    curvature = 1.0 / np.abs((xprim*yprimprim - yprim* xprimprim) / np.power(xprim*xprim + yprim*yprim, 3 / 2))
    minCurvature = np.min(curvature)
    print("s is %g => Minimum curvature radius is %g"%(s,np.min(curvature)))

    #Perform bisection
    if minCurvature > minCurvatureConstraint:
        s_1 = s
        minCurvature_1 = minCurvature
    else:
        s_0 = s
        minCurvature_0 = minCurvature

Es kann einige Verfeinerungen wie Iterationen erfordern, um ein geeignetes s_1 zu finden, aber es funktioniert.