Problembeschrieb
Gegeben sind $n$ Punkte $p_0$ bis $p_{n-1}$, durch welche eine geschlossene Kurve gelegt werden soll. Diese Kurve soll möglichst “angenehm” zu fahren sein, d.h. die Beschleunigungen sind klein und vor allem soll die Beschleunigung stetig sein.
Die Kurve soll stückweise kubisch sein, d.h. durch Bezierkurven darstellbar.
Hermitsche Basis
Wir beschränken uns erst einmal auf den eindimensionalen Fall auf dem Intervall $t \in [0,1]$. Gesucht ist eine kubische Funktion $f(x)$ so, dass $$ f(0)=x_0,\quad f(1)=x_1, \quad f'(0) = m_0 \text{ und} \quad f'(1)=m_1. $$
Bestimmen Sie die kubischen Funktionen für folgende Parameter:
| Funktion | $x_0$ | $x_1$ | $m_0$ | $m_1$ | Funktionsterm |
|---|---|---|---|---|---|
| $h_{0,0}$ | 1 | 0 | 0 | 0 | |
| $h_{1,0}$ | 0 | 1 | 0 | 0 | |
| $h_{0,1}$ | 0 | 0 | 1 | 0 | |
| $h_{1,1}$ | 0 | 0 | 0 | 1 |
Zeigen Sie, dass $f(t) = x_0 h_{0,0}(t) + x_1 h_{1,0}(t)+ m_0 h{0,1}(t) + m_1 h_{1,1}(t)$ eine Lösung für das oben gestellte Problem ist.
Bestimmen Sie die zweiten Ableitungen der obigen Funktionen.
Zusammenfassung
| Funktion | Normalform | Faktorisiert | 2. Ableitung | 2. Ableitung für $t=0$ | 2. Ableitung für $t=1$ |
|---|---|---|---|---|---|
| $h_{0,0}$ | $2t^3-3t^2+1$ | $(1+2t)(1-t)^2$ | $12t-6 $ | -6 | 6 |
| $h_{1,0}$ | $-2t^3+3t^2$ | $t^2(3-2t)$ | $-12t+6$ | 6 | -6 |
| $h_{0,1}$ | $t^3-2t^2+t$ | $t(1-t)^2$ | $6t-4$ | -4 | 2 |
| $h_{1,1}$ | $t^3-t^2$ | $t^2(t-1)$ | $6t-2$ | -2 | 4 |
Für allgemeine Funktionen ergibt sich die zweite Ableitung für $t=0$ und $t=1$:
$f''(0) = -6x_0+6x_1-4m_0-2 m_1$ und $f''(1) = 6x_0-6x_1+2m_0+4m_1$.
Siehe auch: https://de.wikipedia.org/wiki/Kubisch_Hermitescher_Spline
Optimale Steigungen
Sind $n$ Werte $x_0$ bis $x_{n-1}$ gegeben sind die Steigungen $m_0$ bis $m_{n-1}$ so zu bestimmen, dass $$ f''_{i}(1) = f''_{i+1}(0) $$ (Die Indizies sind modulo $n$ aufzufassen, da die Kurve geschlossen werden soll, d.h. $f_{n}=f_{0}$ etc.)
Damit haben wir $n$ Unbekannte und $n$ Gleichungen der Form: $$ \begin{align*} f''_{i}(1) & = & f''_{i+1}(0) \\ 6x_i-6x_{i+1}+2m_i+4m_{i+1} & = & -6x_{i+1}+6x_{i+2}-4m_{i+1}-2m_{i+2}\\ 2m_i+4m_{i+1} & = & -6x_i+6x_{i+2}-4m_{i+1}-2m_{i+2}\\ 2m_i+8m_{i+1}+2m_{i+2} & = & -6x_i+6x_{i+2}\\ m_i+4m_{i+1}+m_{i+2} & = & -3x_i+3x_{i+2}\\ \end{align*} $$ In matrizieller Form geschrieben: \[ \begin{pmatrix} 1 & 4 & 1 & 0 & 0 & \cdots & 0 \\\\ 0 & 1 & 4 & 1 & 0 & \cdots & 0 \\\\ 0 & 0 & 1 & 4 & 1 & \cdots & 0 \\\\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\\\ 4 & 1 & 0 & 0 & 0 & \cdots & 1 \\\\ \end{pmatrix} \cdot \begin{pmatrix} m_0 \\\\ m_1 \\\\ m_2 \\\\ \vdots \\\\ m_{n-1} \end{pmatrix} = \begin{pmatrix} -3 & 0 & 3 & 0 & 0 & \cdots & 0 \\\\ 0 & -3 & 0 & 3 & 0 & \cdots & 0 \\\\ 0 & 0 & -3 & 0 & 3 & \cdots & 0 \\\\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\\\ 0 & 3 & 0 & 0 & 0 & \cdots & -3 \\\\ \end{pmatrix} \cdot \begin{pmatrix} x_0 \\\\ x_1 \\\\ x_2 \\\\ \vdots \\\\ x_{n-1} \\\\ \end{pmatrix} \]
Lineare Abschnitte mit bekannter Steigung
Es ist wünschenswert, gerade Segmente definieren zu können. Ein Segment trägt den Index der ersten Steigung. D.h. Wenn das Segment $i$ linear sein soll, dann sind $m_i$ und $m_{i+1}$ gegeben und zwar als $(x_{i+1}-x_{i})$. Die Einträge in der Matrix ändern sich dann für die Zeilen $i$ und $i+1$ wie folgt:
\[ \left(0 \ldots 0 \quad 1 \quad 0\ldots 0\right) \cdot \vec m = (0 \ldots 0 \quad -1 \quad 1 \quad 0 \ldots 0) \cdot \vec x \] mit $M_{i,i} = 1$ und $X_{i,i}=-1$ und $X_{i,i+1}=1$.
In Python (hier Blender Code) sieht das dann z.B. wie folgt aus:
import bpy import numpy as np from mathutils import Vector # Pfad heisst StartPath, gegebenenfalls anpassen oder auf aktive Auswahl anwenden obj = bpy.data.objects['StartPath'] # Array für Kurvenpunkte mypoints=[] if obj.type == 'CURVE': for subcurve in obj.data.splines: curvetype = subcurve.type if curvetype == 'POLY': for p in subcurve.points: mypoints.append(p.co) print("MyPoints") print(mypoints) n=len(mypoints) # Matrizen M*m = X*x M = np.zeros((n,n)) X = np.zeros((n,n)) x = np.zeros((n,3)) #Matrizen zeilenweise befüllen for i in range(n): M[i][i] = 1 M[i][(i+1)%n] = 4 M[i][(i+2)%n] = 1 # etc. #print(M) #print(X) m = np.linalg.solve(M,np.dot(X,x)) print(m) print(x) print(mypoints) korrektur = 1.0 #0.5519150244/0.5 handle_right = m/3*korrektur+x handle_left = x-m/3*korrektur print(Vector(x[0])) curvedata = bpy.data.curves.new(name="mySpline", type='CURVE') curvedata.dimensions = '3D' objectdata = bpy.data.objects.new("mySpline", curvedata) objectdata.location = (0,0,0) objectdata.data.bevel_depth = 0.03 bpy.context.scene.collection.objects.link(objectdata) #bpy.context.scene.objects.link(objectdata) # railsCol.objects.link(objectdata) polyline = curvedata.splines.new('BEZIER') polyline.bezier_points.add(n-1) polyline.use_cyclic_u = True for i in range(n): point = polyline.bezier_points[i] point.co = Vector(x[i]) point.handle_left = Vector(handle_left[i]) point.handle_right = Vector(handle_right[i]) # point.handle_left = Vector(x[i]) # point.handle_right = Vector(x[i]) point.handle_left_type = 'ALIGNED' point.handle_right_type = 'ALIGNED'
Umwandlung in Bezierkurven
Da wir wissen, dass für eine Bezierkurve $p(t)$ folgendes gilt: $p'(0) = 3(P_1-P_0)$. Und im eindimensionalen Fall gilt $p(0) = x_0$ und $p'(0)=m_0$, also
\begin{align*} p'(0) & = & m_0 \\ 3(P_1- x_0) & = & m_0 \\ P_1 & = & \frac{1}{3}m_0+x_0 \end{align*}
Analog für $P_2$:
\begin{align*} p'(1) & = & m_1 \\ 3(x_1- P_2) & = & m_1 \\ P_2 & = & -\frac{1}{3}m_1+x_1 \end{align*}
Die komplette Berechnung kann für alle Dimensionen unabhängig durchgeführt werden. Damit erhalten wir die Kontrollpunkte der Bezierkurven für eine ruckelfreie Bahn durch $n$ Punkte.
Ausblick
Wie man es noch viel viel besser machen könnte: https://www2.eecs.berkeley.edu/Pubs/TechRpts/1993/CSD-93-732.pdf Siehe z.B. S. 96 (was hier vorgeschlagen ist) und S. 97 (was möglich wäre, aber nicht wirklich in unserem Rahmen hier)…
Variante mit variablen Interallen $[0,a_i]$ anstatt [0,1]
I.A. sind die Punkte $p_i$ nicht äquidistant, daher macht es Sinn die kubischen Funktionen auf den Intervallen $[0,a_i]$ zu bestimmen, wobei $a_i$ z.B. die Distanz der Punkte (oder die Länge einer vorhergehenden Bezierkurve ist).
Gegeben: $a$, $x_0$, $x_1$, $m_0$, $m_1$. Gesucht: Kubische Funktion $g(t)$ so, dass $g(0)=x_0$, $g'(0)=m_0$, $g(a)=x_1$, $g'(a)=m_1$.
Ansatz: $g(t) = f\left(\frac{t}{a}\right)$ wobei $f$ auf dem Intervall $[0,1]$ definiert ist. Damit gilt $g'(t) = f'\left(\frac{t}{a}\right) \cdot \frac{1}{a}$ d.h. $a \cdot g'(t) = f'\left(\frac{t}{a}\right)$.
Wir suchen also $f$ kubisch mit $f'(0)=am_0$ und $f'(1)=am_1$.
Für die zweite Ableitung gilt:
$g'(t) = f''\left(\frac{t}{a}\right) \cdot \frac{1}{a^2}$, d.h. $a^2 \cdot g'(t) = f''\left(\frac{t}{a}\right)$.
Wir suchen also $g_i$ mit $g_i(0)=x_i$, $g_i(a_i) = x_{i+1}$, $g'_i(a_i) = g'_{i+1}(0)$ und $g''_i(a_i) = g''_{i+1}(0)$.
Seien $m_i$ die unbekannten Steigungen von $g_i$ bei $t=0$ (und damit von $g_{i-1}$ bei $t=a_{i-1}$.)
Damit ist $$ g_i(t) = x_i \cdot h_{0,0}\left(\frac{t}{a_i}\right) + x_{i+1} \cdot h_{1,0}\left(\frac{t}{a_i}\right) + a_i \cdot m_i \cdot h_{0,1}\left(\frac{t}{a_i}\right) + a_i \cdot m_{i+1} \cdot h_{1,1}\left(\frac{t}{a_i}\right) $$ und $$ g''_i(t) = \frac{1}{a_i^2} x_i \cdot h''_{0,0}\left(\frac{t}{a_i}\right) + \frac{1}{a_i^2} x_{i+1} \cdot h''_{1,0}\left(\frac{t}{a_i}\right) + \frac{1}{a_i} \cdot m_i \cdot h_{0,1}\left(\frac{t}{a}\right) + \frac{1}{a_i} \cdot m_{i+1} \cdot h_{1,1}\left(\frac{t}{a}\right) $$
Zur Erinnerung: $f''(0) = -6x_0+6x_1-4m_0-2 m_1$ und $f''(1) = 6x_0-6x_1+2m_0+4m_1$, wobei hier $m_0$ und $m_1$ durch $a_i\cdot m_i$ und $a_i\cdot m_{i+1}$ zu ersetzen sind und es gilt $g''_i(t) = \frac{1}{a_i^2} f''(t)$.
Bei $t=0$ und $t=a_i$ ausgewertet: $$ g''_i(0) = \frac{1}{a_i^2}(-6x_i+6x_{i+1}) + \frac{1}{a_i}(-4m_i-2m_{i+1}) $$ und $$ g''_i(a_i) = \frac{1}{a_i^2}(6x_i-6x_{i+1}) + \frac{1}{a_i}(2m_i+4m_{i+1}) $$
Schliesslich die $n$ Gleichungen für $m_i$: $$ g''_i(a_i) = g''_{i+1}(0) $$ Ausgeschrieben: $$ \frac{1}{a_i^2}(6x_i-6x_{i+1}) + \frac{1}{a_i}(2m_i+4m_{i+1}) = g''_{i+1}(0) = \frac{1}{a_{i+1}^2}(-6x_{i+1}+6x_{i+2}) + \frac{1}{a_{i+1}}(-4m_{i+1}-2m_{i+2}) $$ Variablen nach links, Konstanten nach rechts: $$ \frac{1}{a_i}(2m_i+4m_{i+1}) - \frac{1}{a_{i+1}}(-4m_{i+1}-2m_{i+2}) = \frac{1}{a_{i+1}^2}(-6x_{i+1}+6x_{i+2}) - \frac{1}{a_i^2}(6x_i-6x_{i+1}) $$ Zusammenfassen: $$ \frac{2}{a_i}m_i +m_{i+1} \cdot \left(\frac{4}{a_i}+\frac{4}{a_{i+1}}\right) +\frac{2}{a_{i+1}}\cdot m_{i+2} = \frac{1}{a_{i+1}^2}(-6x_{i+1}+6x_{i+2}) - \frac{1}{a_i^2}(6x_i-6x_{i+1}) $$
Simulation der Achterbahn
Zu jedem Zeitpunkt $t$ ist die Position $\vec x(t)$ und die Geschwindigkeit $\vec v(t)$ gesucht auf der Achterbahn gesucht. Der Betrag der Geschwindigkeit ergibt sich im einfachsten Fall mit dem Energieerhaltungssatz direkt aus der Höhe der Position (wenn man Reibung vernachlässigt).
Effektiv werden wir eine “framerate” festlegen (wohl 30 fps) und für jedes frame $i$ möchten wir die Position kennen (und die Geschwindigkeit). Daraus kann z.B. die entsprechende Strecke auf der Kurve vorwärts gegangen werden. Wir müssen also auf einer Bezierkurve eine gegebene Länge vorwärts gehen können.
Ansatz 1: Kleine Abstände von Punkten auf der Kurve aufsummieren: $L = \sum_{k=0}^{n} |\vec x(t_{k+1})-\vec x(t_{k})|$
Ansatz 2: Betrag der Geschwindigkeit integrieren: $L = \int_{a}^b |\vec v(t)|dt$
Die Komponenten von $\vec v(t)$ sind Polynome 2. Grades, d.h. man integriert eine Wurzel aus einem Polynom 4. Grades. Keine Ahnung, ob man das noch explizit gebacken kriegt (maxima streikt schon in ganz einfachen Fällen).
Wir werden das darum nummerisch machen und eine nette Klasse Bezier programmieren, die die nötigen Methoden zur Verfügung stellt.
Um die Schienen zu legen, müssen wir die effektive Beschleunigung kennen, damit die Schienen normal dazu gelegt werden können. Dazu brauchen wir die Normale zur Bahn. Die erhalten wir als Ableitung der normalisierten Geschwindigkeit.
Nötige Methoden für die Bezier Klasse
- x(t): Liefert die Position zum Parameter $t$.
- v(t): Liefert die Geschwindigkeit zum Parameter $t$ (exakt).
- len(tmin, tmax): Liefert die Länge der Kurve zwischen den Parametern tmin, tmax (angenähert).
- vnorm(t): Geschwindkeitsvektor mit Länge 1 zum Parameter $t$.
- adl(t): Die Ableitung von vnorm nach der Länge der Kurve
- bezier2.py
# -*- coding: utf-8 -*- # import importlib # importlib.import_module("./vector") from vector import Vector class Bezier: def __init__(self,pts): if len(pts)==4: self.pts = [p.copy() for p in pts] if len(pts)==3: self.pts = [pts[0].copy(), 1/3*pts[0]+2/3*pts[1], 1/3*pts[2]*2/3*pts[1], pts[2].copy()] if len(pts)==2: self.pts = [pts[0].copy(), 2/3*pts[0]+1/3*pts[1], 1/3*pts[0]+2/3*pts[1], pts[1].copy()] # Punkt für Parameter t (in [0,1]) def x(self,t): return (1.0-t)**3*self.pts[0]+\ 3*(1.0-t)**2*t*self.pts[1]+\ 3*(1.0-t)*t**2*self.pts[2]+\ t**3*self.pts[3] # Geschwindigkeit (dx/dt) in t-Parametrierung def v(self,t): return -3*(1.0-t)**2*self.pts[0]+\ 3*(1.0-4*t+3*t**2)*self.pts[1]+\ 3*(2*t-3*t**2)*self.pts[2]+\ 3*t**2*self.pts[3] # Zeichnet den Spline # Achtung: Funktioniert nur, wenn diese Datei ausgeführt wird... def draw(self, steps=1000): move(self.pts[0][0], self.pts[0][1]); for i in range(1,steps+1): p = self.x(1.0*i/steps) lineTo((p[0],p[1])) if __name__== "__main__": from gpanel import * makeGPanel(Size(600,600)) window(0, 599, 0, 599) b = Bezier([Vector([100,500,0]), Vector([200,0,0]), Vector([300,300,0]), Vector([500,300,0])]) b.draw()