Länge einer kubischen Bezierkurve integrieren



  • Wie mir bekannt ist, lässt sich die Länge einer kubischen Bezierkurve bzw. allgemein die Länge einer Kurve der Art

    C(t) = A t^3 + B t^2 + C t + D mit 0 <= t <= 1

    nicht analytisch berechnen, da kein Integral gefunden werden kann.

    Nun bin ich auf der Suche nach einer schnellen Methode dies numerisch zu berechnen. Allerdings kommt man mit stupiden zerlegen von C(t) in viele kleine Geradenstücke nicht viel weiter, da der Rechenaufwand groß und die Genauigkeit gering ist.

    Nun habe ich schon probiert das Ganze nach de Casteljau (Divide and Conquer) zu rechnen, wobei ich gerade Stücken ohne besonders tiefe Zerlegung bestimmen kann. Dadurch wird das Verhalten schon ein gutes Stück optimaler. Mit double-Präzision erreiche ich hier eine Durchschnittliche Laufzeit von 2 μs.

    Jedoch frage ich mich, ob es da nicht noch eine elegantere Methode gibt. Ich habe zwar diverse Dokumente zu diesem Thema gefunden, aber kaum eines davon setzte sich mit einer Implementierung auseinander. Irgendeine Idee nach was ich dort noch suchen könnte, oder gar einen Vorschlag für einen passenden Algorithmus?

    PS: Weiterhin werde ich wohl auf t für eine bestimmte Länge und die Länge für ein bestimmtes t bestimmen müssen.



  • Ich hab mich jetzt nicht näher damit auseinandergesetzt, aber wäre es nicht viel naheliegender, dass Integral numerisch zu berechnen?



  • Genau das mache ich doch gerade. Die Frage ist nur, welches Verfahren dafür geeignet ist. De Casteljau eignet sich ganz gut, ich weiß nur nicht, ob ein einfacheres Verfahren (schneller zu berechnen) vergleichbare oder gar bessere Ergebnisse bringt. Insbesondere ärgert mich die Wurzel bei der numerischen Integration. Da ich ja L = int(0,1) | C'(t) | lösen muss, was ja der Distanz zwischen zwei Punkten im R² entspricht, also sqrt(x²+y²).



  • Niabot schrieb:

    Genau das mache ich doch gerade.

    Das sieht mir aber nicht so aus. Du versuchst, die Bogenlänge zu berechnen, indem du die Kurve stückweise linear approximierst und die einzelnen Längen zusammenrechnest.

    Ich meinte, dass du das Integral $$ab1+(C˙(t))2dt\int_a^b \sqrt{1+(\dot{C}(t))^2} dt$$ mit irgendeinem numerischen Verfahren integrierst.

    Wie gesagt, ich hab mir das nicht näher angesehen, vielleicht kommt es ja auf das gleiche raus. Vielleicht aber auch nicht. Meine Idee hat immerhin den Vorteil, dass es eine Standardaufgabe ist, wofür es viele Verfahren mit unterschiedlichen Tradeoffs zwischen Aufwand und Genauigkeit gibt.



  • Ist eigentlich die summierte Sehnentrapezformel (also n Trapeze unter die Funktion stellen) eigentlich ein gutes Verfahren zur num. Integration?



  • numeriker2 schrieb:

    Ist eigentlich die summierte Sehnentrapezformel (also n Trapeze unter die Funktion stellen) eigentlich ein gutes Verfahren zur num. Integration?

    Eher nicht, weil mit man mit subminimal mehr Aufwand gleich die Simpson-Regel nehmen kann, also statt stets eine kleine Gerade auf zwei Stützpunkte zu legen, stets eine kleine Parabel auf drei Stützpunkte legen.



  • volkard schrieb:

    numeriker2 schrieb:

    Ist eigentlich die summierte Sehnentrapezformel (also n Trapeze unter die Funktion stellen) eigentlich ein gutes Verfahren zur num. Integration?

    Eher nicht, weil mit man mit subminimal mehr Aufwand gleich die Simpson-Regel nehmen kann, also statt stets eine kleine Gerade auf zwei Stützpunkte zu legen, stets eine kleine Parabel auf drei Stützpunkte legen.

    Dein "besser" bezieht sich aber nur auf "Genauigkeit des Ergebnsisses", oder? Die Simpsonregel scheint ja offenbar teurer (in Laufzeitkosten) zu sein als die Trapezregel.
    Wenn Geschwindigkeit sehr wichtig ist, dann ist die Trapezregel doch ganz gut, oder? (Quasi ein guter Kompromiss aus Genauigkeit und Geschwindigkeit?)



  • numeriker2 schrieb:

    Wenn Geschwindigkeit sehr wichtig ist, dann ist die Trapezregel doch ganz gut, oder? (Quasi ein guter Kompromiss aus Genauigkeit und Geschwindigkeit?)

    Nein. Hattu Simpsonregel angeschaut?
    Simpson hat Gewichte 1 4 2 4 2 4 2...2 4 1
    Sehnentrapez hat Gewichte 1 2 2 2 2 2...2 2 1
    Sonst kein Unterschied, ich glaube nicht, daß *4 lahmer ist als *2. Die Schleife um zwei zu unwinden tut gut, da ist Simpson potentiell schneller, was aber ein unfairer Vergleich ist. Aber irgendwie nicht langsamer.Oder in Zweierschritten auf zwei unterschiedliche Variablen summieren und die erst zum Schluß *2 und *4 machen? Geht alles unter im Vergleich zur Datenbeschaffung, fürchte ich.



  • volkard schrieb:

    numeriker2 schrieb:

    Wenn Geschwindigkeit sehr wichtig ist, dann ist die Trapezregel doch ganz gut, oder? (Quasi ein guter Kompromiss aus Genauigkeit und Geschwindigkeit?)

    Nein. Hattu Simpsonregel angeschaut?
    Simpson hat Gewichte 1 4 2 4 2 4 2...2 4 1
    Sehnentrapez hat Gewichte 1 2 2 2 2 2...2 2 1
    Sonst kein Unterschied, ich glaube nicht, daß *4 lahmer ist als *2. Die Schleife um zwei zu unwinden tut gut, da ist Simpson potentiell schneller, was aber ein unfairer Vergleich ist. Aber irgendwie nicht langsamer.Oder in Zweierschritten auf zwei unterschiedliche Variablen summieren und die erst zum Schluß *2 und *4 machen? Geht alles unter im Vergleich zur Datenbeschaffung, fürchte ich.

    Was bedeuten diese Gewichte? Und wie kommst du auf die Simpson Gewichte von 1,4,2..? Laut Wiki: http://de.wikipedia.org/wiki/Newton-Cotes-Formeln
    hat das Ding Gewichte von 1/6, 1/4, 1/6



  • numeriker2 schrieb:

    Laut Wiki: http://de.wikipedia.org/wiki/Newton-Cotes-Formeln
    hat das Ding Gewichte von 1/6, 1/4, 1/6

    Nee. 1/6 4/6 1/6
    Häng von denen mal ein paar aneinander.
    1/6 4/6 2/6 4/6 2/6 ... 2/6 4/6 1/6
    Und dann können wir auch gleich 1 4 2 4 2 ... 2 4 1 nehmen und erst ganz am Ende durch 6 teilen.



  • Bashar schrieb:

    Das sieht mir aber nicht so aus. Du versuchst, die Bogenlänge zu berechnen, indem du die Kurve stückweise linear approximierst und die einzelnen Längen zusammenrechnest.

    Ich meinte, dass du das Integral $$ab1+(C˙(t))2dt\int_a^b \sqrt{1+(\dot{C}(t))^2} dt$$ mit irgendeinem numerischen Verfahren integrierst.

    Wie gesagt, ich hab mir das nicht näher angesehen, vielleicht kommt es ja auf das gleiche raus. Vielleicht aber auch nicht. Meine Idee hat immerhin den Vorteil, dass es eine Standardaufgabe ist, wofür es viele Verfahren mit unterschiedlichen Tradeoffs zwischen Aufwand und Genauigkeit gibt.

    Genau die Bogenlänge will ich ja berechnen. Nicht etwa die Fläche unter der Kurve. Deshalb gilt ja im zweidimensionalen Fall $$_abC˙(t)dt=_abC˙_x(t)2+C˙_y(t)2dt\int\_a^b |\dot{C}(t)| dt = \int\_a^b \sqrt{\dot{C}\_x(t)^2 + \dot{C}\_y(t)^2} dt$$. Ich habe es nun schon einmal nach Simpson probiert. Nur ist das im Vergleich sehr viel langsamer. Zumindest bei annähernd gleicher Präzision.



  • Aufbauend auf de Casteljau hatte ich folgenden Code geschrieben. Er ist recht schnell, wobei das natürlich von den Konstanten abhängt. Natürlich ist die Fehlerabschätzung dann schwierig, weshalb man diesen wohl nur anhand zufälliger Ergebnisse mit einem anderen Verfahren vergleichen kann. Welche Ideen hättet ihr denn wie man es noch berechnen könnte. Etwas Code wäre natürlich auch nicht schlecht, damit man sie in der Laufzeit miteinander vergleichen kann. 😋

    private static final int MAX_DEPTH = 10;
    	private static final double RELATIVE_TOLLERANCE = 0.001;
    
    	public static double approximateArcLength(double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3) {
    		return approximateArcLengthRecursive(x0, y0, x1, y1, x2, y2, x3, y3, 0);
    	}
    
    	private static double approximateArcLengthRecursive(double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, int depth) {
    		double dx = x3 - x0;
    		double dy = y3 - y0;
    		double d1x = x1 - x0;
    		double d1y = y1 - y0;
    		double d2x = x2 - x3;
    		double d2y = y2 - y3;
    
    		double dlenSq = dx*dx + dy*dy;
    		double dist1Sq, dist2Sq;
    		if (dlenSq > 0.0) {
    			dist1Sq = distanceFromVectorLineSegmentSquared(d1x, d1y, dx, dy, dlenSq);
    			dist2Sq = distanceFromVectorLineSegmentSquared(d2x, d2y, -dx, -dy, dlenSq);
    		} else {
    			dist1Sq = 0.0;
    			dist2Sq = 0.0;
    		}
    
    		if ((depth >= 1 && (dist1Sq + dist2Sq) < RELATIVE_TOLLERANCE * dlenSq) || (depth >= MAX_DEPTH)) {
    			return Math.sqrt(dlenSq);
    		}
    
    		double x01 = (x0 + x1) / 2.0;
    		double y01 = (y0 + y1) / 2.0;
    		double x12 = (x1 + x2) / 2.0;
    		double y12 = (y1 + y2) / 2.0;
    		double x23 = (x2 + x3) / 2.0;
    		double y23 = (y2 + y3) / 2.0;
    
    		double x012 = (x01 + x12) / 2.0;
    		double y012 = (y01 + y12) / 2.0;
    		double x123 = (x12 + x23) / 2.0;
    		double y123 = (y12 + y23) / 2.0;
    
    		double x0123 = (x012 + x123) / 2.0;
    		double y0123 = (y012 + y123) / 2.0;
    
    		return approximateArcLengthRecursive(x0, y0, x01, y01, x012, y012, x0123, y0123, depth+1)
    		     + approximateArcLengthRecursive(x0123, y0123, x123, y123, x23, y23, x3, y3, depth+1);
    	}
    
    	/**
    	 * Returns the squared distance from a Point <b>p</b> to a line segment
    	 * from the origin (0,0) to the endpoint <b>v</b>.
    	 * @param px the x coordinate of the point
    	 * @param py the y coordinate of the point
    	 * @param vx the x component of the direction vector from the line
    	 * @param vy the y component of the direction vector from the line
    	 * @param vlen2 the squared length of the vector v
    	 * @return a double value of the squared distance of p to the line
    	 */
    	private static double distanceFromVectorLineSegmentSquared(double px, double py, double vx, double vy, double vlen2) {
    		double dot = vx*px + vy*py;
    		if (dot > 0.0) {
    			px -= vx;
    			py -= vy;
    			dot = vx*px + vy*py;
    			if (dot > 0.0)
    				return px*px + py*py;
    			double cross = px*vy - py*vx;
    			return (cross * cross) / vlen2;
    		}
    		return px*px + py*py;
    	}
    

Log in to reply