Formel für kubische Bezier-Kurve nach t umformen
Roman Bobik
Ich schreibe gerade an einem Programm in dem es notwendig ist, dass man
Bezier-Kurven erstellen kann.
Nun bin ich auf ein Problem gestoßen.
Ich gehe von der folgenden Formel für eine kubische Bezier-Kurve aus:
x = x1*b1(t) + x2*b2(t) + x3*b3(t) + x4*b4(t)
wobei:
b1(t) := t ^ 3
b2(t) := 3 * t ^ 2 * (1 - t)
b3(t) := 3 * t * (1 - t) ^ 2
b4(t) := (1 - t) ^ 3
also:
x = x1 * t ^ 3 + x2 * 3 * t ^ 2 * (1 - t) + x3 * 3 * t * (1 - t) ^ 2 + x4 *
(1 - t) ^ 3
Soweit so gut. Eine Funktion im Programm benötigt allerdings diese nach t
umgeformte Formel.
Hab mir gedacht dass das nicht so kompliziert sein kann und hab mir Derive 5
Trial heruntergeladen und damit diese Formel umgeformt.
Nach 8 Minuten (!) Rechenzeit auf meinem K6-2 400 bekam ich folgende kurze
und übersichtliche Formel (bzw. 3 Lösungen):
t=(SQRT(3)*SQRT(x1*(x4-x3)+x2^2-x2*(x4+x3)+x3^2)*SIGN(x1-3*x2-x4+3*x3)*COS(A
TAN((x1*(x4-x3)+x2^2-x2*(x4+x3)+x3^2)^(3/2)*(x*(x1-3*x2-x4+3*x3)^2-x1^2*x4+x
1*(3*x2*(x4+x3)-x4^2+3*x3*(x4-2*x3))-2*x2^3+3*x2^2*(x3-2*x4)+x3*(3*x2*(x4+x3
)-2*x3^2))*SQRT(-(x^2*(x1-3*x2-x4+3*x3)^2-2*x*(x1^2*x4-x1*(3*x2*(x4+x3)-x4^2
+3*x3*(x4-2*x3))+2*x2^3+3*x2^2*(2*x4-x3)-3*x2*x3*(x4+x3)+2*x3^3)+x1^2*x4^2+2
*x1*x3*(2*x3^2-3*x2*x4)+x2^2*(4*x2*x4-3*x3^2))/(x1*(x4-x3)+x2^2-x2*(x4+x3)+x
3^2)^3)/((x1-3*x2-x4+3*x3)*(x^2*(x1-3*x2-x4+3*x3)^2-2*x*(x1^2*x4-x1*(3*x2*(x
4+x3)-x4^2+3*x3*(x4-2*x3))+2*x2^3+3*x2^2*(2*x4-x3)-3*x2*x3*(x4+x3)+2*x3^3)+x
1^2*x4^2+2*x1*x3*(2*x3^2-3*x2*x4)+x2^2*(4*x2*x4-3*x3^2))))/3)-SQRT(x1*(x4-x3
)+x2^2-x2*(x4+x3)+x3^2)*SIGN(x1-3*x2-x4+3*x3)*SIN(ATAN((x1*(x4-x3)+x2^2-x2*(
x4+x3)+x3^2)^(3/2)*(x*(x1-3*x2-x4+3*x3)^2-x1^2*x4+x1*(3*x2*(x4+x3)-x4^2+3*x3
*(x4-2*x3))-2*x2^3+3*x2^2*(x3-2*x4)+x3*(3*x2*(x4+x3)-2*x3^2))*SQRT(-(x^2*(x1
-3*x2-x4+3*x3)^2-2*x*(x1^2*x4-x1*(3*x2*(x4+x3)-x4^2+3*x3*(x4-2*x3))+2*x2^3+3
*x2^2*(2*x4-x3)-3*x2*x3*(x4+x3)+2*x3^3)+x1^2*x4^2+2*x1*x3*(2*x3^2-3*x2*x4)+x
2^2*(4*x2*x4-3*x3^2))/(x1*(x4-x3)+x2^2-x2*(x4+x3)+x3^2)^3)/((x1-3*x2-x4+3*x3
)*(x^2*(x1-3*x2-x4+3*x3)^2-2*x*(x1^2*x4-x1*(3*x2*(x4+x3)-x4^2+3*x3*(x4-2*x3)
)+2*x2^3+3*x2^2*(2*x4-x3)-3*x2*x3*(x4+x3)+2*x3^3)+x1^2*x4^2+2*x1*x3*(2*x3^2-
3*x2*x4)+x2^2*(4*x2*x4-3*x3^2))))/3)-x2-x4+2*x3)/(x1-3*x2-x4+3*x3) OR
t=-(SQRT(3)*SQRT(x1*(x4-x3)+x2^2-x2*(x4+x3)+x3^2)*SIGN(x1-3*x2-x4+3*x3)*COS(
ATAN((x1*(x4-x3)+x2^2-x2*(x4+x3)+x3^2)^(3/2)*(x*(x1-3*x2-x4+3*x3)^2-x1^2*x4+
x1*(3*x2*(x4+x3)-x4^2+3*x3*(x4-2*x3))-2*x2^3+3*x2^2*(x3-2*x4)+x3*(3*x2*(x4+x
3)-2*x3^2))*SQRT(-(x^2*(x1-3*x2-x4+3*x3)^2-2*x*(x1^2*x4-x1*(3*x2*(x4+x3)-x4^
2+3*x3*(x4-2*x3))+2*x2^3+3*x2^2*(2*x4-x3)-3*x2*x3*(x4+x3)+2*x3^3)+x1^2*x4^2+
2*x1*x3*(2*x3^2-3*x2*x4)+x2^2*(4*x2*x4-3*x3^2))/(x1*(x4-x3)+x2^2-x2*(x4+x3)+
x3^2)^3)/((x1-3*x2-x4+3*x3)*(x^2*(x1-3*x2-x4+3*x3)^2-2*x*(x1^2*x4-x1*(3*x2*(
x4+x3)-x4^2+3*x3*(x4-2*x3))+2*x2^3+3*x2^2*(2*x4-x3)-3*x2*x3*(x4+x3)+2*x3^3)+
x1^2*x4^2+2*x1*x3*(2*x3^2-3*x2*x4)+x2^2*(4*x2*x4-3*x3^2))))/3)+SQRT(x1*(x4-x
3)+x2^2-x2*(x4+x3)+x3^2)*SIGN(x1-3*x2-x4+3*x3)*SIN(ATAN((x1*(x4-x3)+x2^2-x2*
(x4+x3)+x3^2)^(3/2)*(x*(x1-3*x2-x4+3*x3)^2-x1^2*x4+x1*(3*x2*(x4+x3)-x4^2+3*x
3*(x4-2*x3))-2*x2^3+3*x2^2*(x3-2*x4)+x3*(3*x2*(x4+x3)-2*x3^2))*SQRT(-(x^2*(x
1-3*x2-x4+3*x3)^2-2*x*(x1^2*x4-x1*(3*x2*(x4+x3)-x4^2+3*x3*(x4-2*x3))+2*x2^3+
3*x2^2*(2*x4-x3)-3*x2*x3*(x4+x3)+2*x3^3)+x1^2*x4^2+2*x1*x3*(2*x3^2-3*x2*x4)+
x2^2*(4*x2*x4-3*x3^2))/(x1*(x4-x3)+x2^2-x2*(x4+x3)+x3^2)^3)/((x1-3*x2-x4+3*x
3)*(x^2*(x1-3*x2-x4+3*x3)^2-2*x*(x1^2*x4-x1*(3*x2*(x4+x3)-x4^2+3*x3*(x4-2*x3
))+2*x2^3+3*x2^2*(2*x4-x3)-3*x2*x3*(x4+x3)+2*x3^3)+x1^2*x4^2+2*x1*x3*(2*x3^2
-3*x2*x4)+x2^2*(4*x2*x4-3*x3^2))))/3)+x2+x4-2*x3)/(x1-3*x2-x4+3*x3) OR
t=(2*SQRT(x1*(x4-x3)+x2^2-x2*(x4+x3)+x3^2)*SIGN(x1-3*x2-x4+3*x3)*SIN(ATAN((x
1*(x4-x3)+x2^2-x2*(x4+x3)+x3^2)^(3/2)*(x*(x1-3*x2-x4+3*x3)^2-x1^2*x4+x1*(3*x
2*(x4+x3)-x4^2+3*x3*(x4-2*x3))-2*x2^3+3*x2^2*(x3-2*x4)+x3*(3*x2*(x4+x3)-2*x3
^2))*SQRT(-(x^2*(x1-3*x2-x4+3*x3)^2-2*x*(x1^2*x4-x1*(3*x2*(x4+x3)-x4^2+3*x3*
(x4-2*x3))+2*x2^3+3*x2^2*(2*x4-x3)-3*x2*x3*(x4+x3)+2*x3^3)+x1^2*x4^2+2*x1*x3
*(2*x3^2-3*x2*x4)+x2^2*(4*x2*x4-3*x3^2))/(x1*(x4-x3)+x2^2-x2*(x4+x3)+x3^2)^3
)/((x1-3*x2-x4+3*x3)*(x^2*(x1-3*x2-x4+3*x3)^2-2*x*(x1^2*x4-x1*(3*x2*(x4+x3)-
x4^2+3*x3*(x4-2*x3))+2*x2^3+3*x2^2*(2*x4-x3)-3*x2*x3*(x4+x3)+2*x3^3)+x1^2*x4
^2+2*x1*x3*(2*x3^2-3*x2*x4)+x2^2*(4*x2*x4-3*x3^2))))/3)-x2-x4+2*x3)/(x1-3*x2
-x4+3*x3)
so. Mein erster Eindruck war natürlich das das nicht stimmen kann...
Jedoch versuchte ich diese Formeln zu verwenden und es ist tatsächlich so
das min. eine dieser Lösungen (fast) immer das richtige Ergebnis liefert.
Jedoch gibt es Grenzsituationen bei der keine der Lösungen richtig ist (Div
durch 0, Wurzel aus neg. Zahl, oder t nicht zwischen 0 und 1, wies
eigentlich sein sollte). Das ist vor allem dann der Fall, wenn der Wert x2
"nahe" am Wert x4 ist oder x3 "nahe" an x1. "nahe" bedeutet, dass z.B.
(x4-x2) < (x2-x1) (x2 ist näher an x4 als an x1...). Wann das genau auftritt
und was der Grund ist weiß ich nicht genau.
Jetzt würde ich gerne wissen ob ich etwas falsch mache, ob ich an die
Grenzen von Derive gestoßen bin oder ob es da vielleicht eine einfachere
Lösung gibt... "Näherungsverfahren" hab ich versucht ist allerdings nicht
sehr performant, da diese Berechnung min. 50 mal/sec durchgeführt werden
soll. Vielleicht hat jemand eine Idee?
Vielen Dank!
Mit freundlichen Grüßen!
Roman Bobik
Hermann Kremer
<3bcf1bfc$0$33462$5039...@newsreader01.highway.telekom.at>...
>Hi Mathe-Profis!
Hallo Roman,
>Ich schreibe gerade an einem Programm in dem es notwendig ist, dass man
>Bezier-Kurven erstellen kann.
>Nun bin ich auf ein Problem gestoßen.
>Ich gehe von der folgenden Formel für eine kubische Bezier-Kurve aus:
>
>x = x1*b1(t) + x2*b2(t) + x3*b3(t) + x4*b4(t)
>
>wobei:
>b1(t) := t ^ 3
>b2(t) := 3 * t ^ 2 * (1 - t)
>b3(t) := 3 * t * (1 - t) ^ 2
>b4(t) := (1 - t) ^ 3
>
>also:
>
> x = x1 * t ^ 3 + x2 * 3 * t ^ 2 * (1 - t) + x3 * 3 * t * (1 - t) ^ 2 +
> + x4 * (1 - t) ^ 3
OK, das übliche Bernstein-Polynom ...
>Soweit so gut. Eine Funktion im Programm benötigt allerdings diese nach t
>umgeformte Formel.
>Hab mir gedacht dass das nicht so kompliziert sein kann ...
Hmmm ...
>dass mind. eine dieser Lösungen (fast) immer das richtige Ergebnis liefert.
OK, Du hast eine kubische Gleichung, in allgemeiner Form
a_3 * t^3 + a_2 * t^2 + a_1 * t + a_0 = 0 ,
und die hat genau 3 Lösungen. Eine dieser Lösungen ist immer reell -
m.E. müßte das die dritte oben sein, und die beiden anderen Lösungen
sind entweder beide reell, oder sie sind konjugiert komplex. Da Du vermutlich
aber nur die reelle Lösung im Intervall 0 <= t <= 1 brauchen kannst, mußt
Du dafür sorgen, daß die komplexen Lösungen gar nicht erst berechnet werden.
Damit Du eine Vorstellung erhältst, zeichne Dir mal den Graphen für irgend
eine Bézier-Kurve und zeichne dann Parallelen zur t-Achse bei verschiedenen
x-Werten. Du siehst jetzt, wie oft und wo eine solche Parallele den Bézier-
Graphen schneidet - und überall dort ist eine Lösung für das jeweilige x .
>Jedoch gibt es Grenzsituationen bei der keine der Lösungen richtig ist (Div
>durch 0, Wurzel aus neg. Zahl, oder t nicht zwischen 0 und 1, wies
>eigentlich sein sollte). Das ist vor allem dann der Fall, wenn der Wert x2
>"nahe" am Wert x4 ist oder x3 "nahe" an x1. "nahe" bedeutet, dass z.B.
>(x4-x2) < (x2-x1) (x2 ist näher an x4 als an x1...). Wann das genau auftritt
>und was der Grund ist weiß ich nicht genau.
OK, es gibt algebraische Verfahren, mit denen Du vorher bestimmen
kannst, ob für ein gegebenes x im Intervall 0 <= t <= 1 drei, zwei, eine
oder gar keine Lösungen vorhanden sind (Sturm'sche Kette u.a); ein solcher
Test ist jedoch ziemlich aufwendig.
>Jetzt würde ich gerne wissen ob ich etwas falsch mache, ob ich an die
>Grenzen von Derive gestoßen bin oder ob es da vielleicht eine einfachere
>Lösung gibt... "Näherungsverfahren" hab ich versucht ist allerdings nicht
>sehr performant, da diese Berechnung min. 50 mal/sec durchgeführt werden
>soll. Vielleicht hat jemand eine Idee?
Du hast nix falsch gemacht, und Dein Derive war auch brav - Du hast das
Problem nur sehr ungeschickt angepackt. Die Derive-Formel ist zwar
mathematisch völlig OK, aber numerisch denkbar ungeeignet, da
dort jede Menge Zeug mehrfach berechnet wird und dadurch ganz üble
Stellenauslöschungsfehler auftreten können.
Ich würde folgendes vorschlagen:
1) Du rechnest Dein Bézier-Polynom in die Form
x = a_3*t^3 + a_2*t^2 + a_1*t + a_0
um, wobei die a_i natürlich von den Bézier-Knoten x1 ... x4 abhängen;
für jeden Satz von Bézier-Knoten mußt Du das nur einmal machen.
2) Für ein gegebenes x berechnest Du die reelle(n) Wurzeln(n) der
Polynomgleichung
a_3*t^3 + a_2*t^2 + a_1*t + a_0 - x = 0
Wie das geht, kannst Du in
http://hem.passagen.se/ceem/cubic/index.htm
nachlesen; und wenn Du in C (oder C++) programmierst, kann ich Dir dafür
auch eine fertige function posten.
Grüße
Hermann
--
Marcus Roeckrath
Roman Bobik wrote:
> Ich schreibe gerade an einem Programm in dem es notwendig ist, dass man
> Bezier-Kurven erstellen kann.
Was willst Du genau? Bezierkurven zeichnen?
Da gibt es sehr schöne Alghorithmen, die nur über Teilung des durch die
Gewichtspunkte gebildeten Polygons arbeiten.
--
Gruss Marcus
Marcus Roeckrath -- Vikarsbusch 8 -- D-48308 Senden -- Germany
Phone : +49-2536-9944 -- Mailer/BBS/Fax : +49-2536-9943 (V34, X75)
FidoNet: 2:2449/523
E-Mail : marcus.r...@gmx.de
WWW : http://home.foni.net/~marcusroeckrath/
Jan C. Hoffmann
"Roman Bobik" <roman...@aon.at> schrieb im Newsbeitrag news:3bcf1bfc$0$33462$5039...@newsreader01.highway.telekom.at...
> Hi Mathe-Profis!
>
> Ich schreibe gerade an einem Programm in dem es notwendig ist, dass man
> Bezier-Kurven erstellen kann.
> Nun bin ich auf ein Problem gestoßen.
> Ich gehe von der folgenden Formel für eine kubische Bezier-Kurve aus:
>
> x = x1*b1(t) + x2*b2(t) + x3*b3(t) + x4*b4(t)
>
> wobei:
> b1(t) := t ^ 3
> b2(t) := 3 * t ^ 2 * (1 - t)
> b3(t) := 3 * t * (1 - t) ^ 2
> b4(t) := (1 - t) ^ 3
>
> also:
>
> x = x1 * t ^ 3 + x2 * 3 * t ^ 2 * (1 - t) + x3 * 3 * t * (1 - t) ^ 2 + x4 *
> (1 - t) ^ 3
>
> Soweit so gut. Eine Funktion im Programm benötigt allerdings diese nach t
> umgeformte Formel.
Hallo Roman,
das ist einfach mit Newton Iteration.
Beispiel
0 = x_1 * t ^ 3 + x_2 * 3 * t ^ 2 * (1 - t) + x_3 * 3 * t * (1 - t) ^ 2 + x_4 * (1 - t) ^ 3 - x
t = 1,166667
für
x_1 = 1,000000
x_2 = 2,000000
x_3 = 3,000000
x_4 = 4,000000
x = 0,500000
--
Gruss Jan C. Hoffmann (c/o MTEC)
http://ourworld.compuserve.com/homepages/MTEC/
Roman Bobik
> Roman Bobik wrote:
>
> > Ich schreibe gerade an einem Programm in dem es notwendig ist, dass man
> > Bezier-Kurven erstellen kann.
>
> Was willst Du genau? Bezierkurven zeichnen?
Danke erst mal für die Antwort. Nein, Bezier kurven zeichnen ist kein
Problem, da gibt es die Windows-Api-Funktion PolyBezier, die x1 .. xn als
Parameter bekommt. Das funktioniert ja ohne Problem. Ich weiß auch wie ich
"händisch" auf die Pixelkoordinaten komme (nämlich mit der oberen Formel:
einmal mit x, einmal mit y).
Angenommen ich zeichne mir jetzt eine Kurve, die von x=0 nach x=100 geht.
Die Kurve befindet sich immer zwischen diesen Koordinaten, diese müssen
immer Start und Endpunkt sein. Eine Funktion im Programm soll jetzt zu einem
x-Wert (0..100) den zugehörigen y-Punkt berechnen. Ich habe die
Kurven-Erstellung im Programm so eingeschränkt, dass keine "Loopings", also
keine 2 oder 3 Werte pro x-Punkt entstehen können. Hintergrund ist der, dass
ich per "Zeitschleife" den y-Punkt ausgeben will, wobei die Zeit durch x
gegeben ist. Wie ich bemerkt habe kann ich nicht einfach den y-Punkt für t =
Zeit ausrechnen, da t ja nicht linear wie x-Achse ist, sondern genau die
Position auf der Kurven-Linie angibt. Ich hoffe das ist zu verstehen... Habe
mit der oben genannten Formel ja auch das richtige Ergenis erziehlt, jedoch
in Grenzsituationen nicht.
> Da gibt es sehr schöne Alghorithmen, die nur über Teilung des durch die
> Gewichtspunkte gebildeten Polygons arbeiten.
ähh, bitte? Könntest du mir das genauer erklären?
Danke!
Roman Bobik
Roman Bobik
> OK, Du hast eine kubische Gleichung, in allgemeiner Form
>
> a_3 * t^3 + a_2 * t^2 + a_1 * t + a_0 = 0 ,
>
> und die hat genau 3 Lösungen. Eine dieser Lösungen ist immer reell -
> m.E. müßte das die dritte oben sein, und die beiden anderen Lösungen
> sind entweder beide reell, oder sie sind konjugiert komplex. Da Du
vermutlich
> aber nur die reelle Lösung im Intervall 0 <= t <= 1 brauchen kannst,
mußt
> Du dafür sorgen, daß die komplexen Lösungen gar nicht erst berechnet
werden.
ja stimmt, müsste die dritte sein. Die ersten zwei unterscheiden sich nur
durch das Vorzeichen. Hab versucht mittels VC++ das ganze mit komplexen
Zahlen zu rechnen, hab aber die gleich (falschen) Ergebnisse bekommen.
OK, das mit der Formel hab ich glaub ich geschnallt.
Laut meinen Berechnungen müsste:
a_3 = x1 - 3*x2 + 3*x3 + x4
a_2 = 3 * x2 - 6 * x3 + 3 * x4
a_1 = 3 * x3 - 3 * x4 * t
a_0 = x4
sein. (sollte stimmen wenn ich mich verrechnet hab werd ichs ja noch
merken...)
> Damit Du eine Vorstellung erhältst, zeichne Dir mal den Graphen für irgend
> eine Bézier-Kurve und zeichne dann Parallelen zur t-Achse bei
verschiedenen
> x-Werten. Du siehst jetzt, wie oft und wo eine solche Parallele den
Bézier-
> Graphen schneidet - und überall dort ist eine Lösung für das jeweilige
x .
ja hab ich schon mal gemacht.
> >Jedoch gibt es Grenzsituationen bei der keine der Lösungen richtig ist
(Div
> >durch 0, Wurzel aus neg. Zahl, oder t nicht zwischen 0 und 1, wies
> >eigentlich sein sollte). Das ist vor allem dann der Fall, wenn der Wert
x2
> >"nahe" am Wert x4 ist oder x3 "nahe" an x1. "nahe" bedeutet, dass z.B.
> >(x4-x2) < (x2-x1) (x2 ist näher an x4 als an x1...). Wann das genau
auftritt
> >und was der Grund ist weiß ich nicht genau.
>
> OK, es gibt algebraische Verfahren, mit denen Du vorher bestimmen
> kannst, ob für ein gegebenes x im Intervall 0 <= t <= 1 drei,
zwei, eine
> oder gar keine Lösungen vorhanden sind (Sturm'sche Kette u.a); ein
solcher
> Test ist jedoch ziemlich aufwendig.
Benötige ich das? Es sollte ja sowieso Lösung 3 richtig sein, die anderen
fallen ja deshal weg, weil sie entweder im Reelen nicht lösbar sind oder
nicht im bereich 0 .. 1 liegen.
> >Jetzt würde ich gerne wissen ob ich etwas falsch mache, ob ich an die
> >Grenzen von Derive gestoßen bin oder ob es da vielleicht eine einfachere
> >Lösung gibt... "Näherungsverfahren" hab ich versucht ist allerdings nicht
> >sehr performant, da diese Berechnung min. 50 mal/sec durchgeführt werden
> >soll. Vielleicht hat jemand eine Idee?
>
> Du hast nix falsch gemacht, und Dein Derive war auch brav - Du hast das
> Problem nur sehr ungeschickt angepackt. Die Derive-Formel ist zwar
> mathematisch völlig OK, aber numerisch denkbar ungeeignet, da
> dort jede Menge Zeug mehrfach berechnet wird und dadurch ganz üble
> Stellenauslöschungsfehler auftreten können.
Ja das mit der Mehrfachberechnung ist mir auch aufgefallen und das
möglicherweise Stellenauslöschungsfehler zu falschen Ergebnissen führen
könnten, war mir auch klar, allerdings wenn der Wertebereich/Genauigkeit von
double nicht ausreicht, was soll man da als Programmierer machen... Aber ok
wenn die obige Lösung einfacher ist dann lässt sich das umgehen.
> Ich würde folgendes vorschlagen:
>
> 1) Du rechnest Dein Bézier-Polynom in die Form
>
> x = a_3*t^3 + a_2*t^2 + a_1*t + a_0
>
> um, wobei die a_i natürlich von den Bézier-Knoten x1 ... x4
abhängen;
> für jeden Satz von Bézier-Knoten mußt Du das nur einmal machen.
hab ich gemacht.
> 2) Für ein gegebenes x berechnest Du die reelle(n) Wurzeln(n) der
> Polynomgleichung
>
> a_3*t^3 + a_2*t^2 + a_1*t + a_0 - x = 0
>
> Wie das geht, kannst Du in
> http://hem.passagen.se/ceem/cubic/index.htm
> nachlesen; und wenn Du in C (oder C++) programmierst, kann ich Dir
dafür
> auch eine fertige function posten.
Das Programm wird zwar nicht in C entwickelt, allerdings kann ich es auch.
Wäre nett wenn du mir diese Funktion posten könntest!
Danke für die ausführliche Antwort!
Roman
Roman Bobik
> ... da t ja nicht linear wie x-Achse ist, sondern genau die
> Position auf der Kurven-Linie angibt ...
stimmt natürlich auch nicht, sorry. Außerdem war das ganze jetzt für eine
XY-Bezier-Kurve gedacht. Die Idee, wie ich auf den y-Wert komme ist, ich
berechne mit t aus der x-Formel und setze das Ergebnis in die y-Formel ein.
mfg Roman
Roman Bobik
> das ist einfach mit Newton Iteration.
du meinst das Newton'sche Näherungsverfahren? Hab ich mir auch schon gedacht
(wie oben erwähnt), allerdings benötige ich max. 50 * 512 solche
Berechnungen pro Sekunde, und das ist ja nicht alles was das Programm
gleichzeitig macht... Es ist also rechnerisch ziemlich aufwendig, vor allem
weil pro Sekunde 50 Berechnugen durch geführt werden, bei denen sich nur x
ändert und sonst alles gleich bleibt.
> Beispiel
>
> 0 = x_1 * t ^ 3 + x_2 * 3 * t ^ 2 * (1 - t) + x_3 * 3 * t * (1 - t) ^ 2 +
x_4 * (1 - t) ^ 3 - x
>
> t = 1,166667
>
> für
>
> x_1 = 1,000000
> x_2 = 2,000000
> x_3 = 3,000000
> x_4 = 4,000000
>
> x = 0,500000
das verstehe ich nicht ganz. x muss ja zwischen x_1 und x_4 liegen, sonst
ist es nicht auf der Kurve, d.h. es gibt dafür auch "keinen" t-Wert (bzw.
keinen gültigen: 0 <= t <= 1).
danke.
Roman
Marcus Roeckrath
Roman Bobik wrote:
>> Da gibt es sehr schöne Alghorithmen, die nur über Teilung des durch die
>> Gewichtspunkte gebildeten Polygons arbeiten.
>
> ähh, bitte? Könntest du mir das genauer erklären?
Da kram ich jetzt mal in meinem alterschwachen Gehirn rum, ist nämlich
schon über 17 Jahre her, als ich mich im Rahmen meines Mathematikstudiums
in einer CAD-Forschungsgruppe auch mit B-Splines und Bezierkurven
auseinandergesetzt habe.
Satz: Algorithmus von Casteljau (vgl. auch Algorithmus von Cox und de Boor
für B-Splines)
Sei P(m) ein Polygon mit den Bezierpunkten b(0), ...,b(m) und 0<=t<=m fest
aber beliebig gewählt, dann gilt der Kurvenpunkt X(t) ergibt sich als
x(t) = b(m,m)
mit
b(i,j)=(1-t)b(i-1,j-1)+t*b(i,j-1)
mit j=1,...,m
und i=j,...,m
Anmerkung: Auf P, b und X gehören Vekorpfeile.
Was nun vielleicht noch etwas kompliziiert aussieht, ist ein schneller
Algorithmus, der praktisch so funktioniert.
Man verbinde die Ecken des Polygons und markiere auf jeder Polygonseite den
Punkt der dem Parameter t (z. B. 0.4) entspricht, wobei der Anfangspunkt
des Polygons bei t=0 der Endpunkt bei T=1 liegt.
Das macht man mit jeder Polygonseite und erhält so einen neuen Polygonzug,
der nun ein Segment weniger hat.
Diesen Teilungsschritt führt man nun so lange durch, bis nur noch ein Punkt
übrigbleibt, daß ist der gesuchte Kurvenpunkt zu t. Den gesamten Kurvenzug
erhält man durch Variation von t.
Eigentlich haben wir nun erst ein Beziersegment. Eine Bezierkurve erhält
man durch Aneinanderhängen mehrerer Segmente z. B. Segmente mit 4
Bezierpunkten (Grad 3), die jedoch in den Trennstellen gewissen Regeln
genügen müssen, damit in den Trennstellen Steigkeit und auch
Differenzierbarkeit mehreren Grades gegeben ist.
Bei 4 Punkten pro Segment kann man maximal 2-mal stetige
Differenzierbarkeit voraussetzen.
Dabei verläßt man zur Konstruktion die Bezierpunkte und definiert
sogenannte Gewichtspunkte, aus denen dann, bis auf den ersten und
allerletzten Bezierpunkt der Bezierkurve sämtliche Bezierpunkte berechnet
und dann durch obigen Algorithmus die Kurve konstruiert werden kann:
Algorithmus zur Bestimmung der Bezierpolygone aus dem Gewichtspolygon für
kubische Beziersegmente und C²-Stetigkeit in den Übergangsstellen der
Segmente.
Seinen k+1 Gewichtspunkte d(0), ..., d(k) gegeben, so errechnen sich die
3k-1 inneren Bezierpunkte einer kubischen Bezierkurve wie folgt:
b(3l-2)=1/3(2d(l-1) + d(l) l=1,...,k
b(3l)=1/6(d(l-1 + 4d(l) + d(l+1) l=1,...,k-1
b(3l+2)=1/3(d(l) + 2d(l+1) l=0,...,k-1
Anmerkung: d und b sind Vektoren.
b(0) und b(3k) sind frei wählbar.
Hermann Kremer
>"Hermann Kremer" <hermann...@online.de> schrieb:
>> OK, Du hast eine kubische Gleichung, in allgemeiner Form
>>
>> a_3 * t^3 + a_2 * t^2 + a_1 * t + a_0 = 0 ,
[ ... ]
>OK, das mit der Formel hab ich glaub ich geschnallt.
>Laut meinen Berechnungen müsste:
>
> a_3 = x1 - 3*x2 + 3*x3 - x4
> a_2 = 3*x2 - 6*x3 + 3*x4
> a_1 = 3*x3 - 3*x4
> a_0 = x4
>
>sein. (sollte stimmen wenn ich mich verrechnet hab werd ichs ja noch
>merken...)
Stimmt bis auf einen kleinen Vorzeichenfehler in a_3 (oben schon korrigiert).
[ ... ]
>>.... ein solcher Test ist jedoch ziemlich aufwendig.
>
>Benötige ich das? Es sollte ja sowieso Lösung 3 richtig sein, die anderen
>fallen ja deshal weg, weil sie entweder im Reelen nicht lösbar sind oder
>nicht im bereich 0 .. 1 liegen.
Ja, in abgespeckter Form schon - Du mußt ja vorher wissen, was Du
rechnen mußt und was nicht ...
>Ja das mit der Mehrfachberechnung ist mir auch aufgefallen und das
>möglicherweise Stellenauslöschungsfehler zu falschen Ergebnissen führen
>könnten, war mir auch klar, allerdings wenn der Wertebereich/Genauigkeit von
>double nicht ausreicht, was soll man da als Programmierer machen...
Hmm, nach meiner Erfahrung reicht die Genauigkeit von "double" *fast*
immer, wenn man den Rechenverlauf bezüglich numerischer Stabilität
optimiert und entsprechend programmiert ...
>> Ich würde folgendes vorschlagen:
>>
>> 1) Du rechnest Dein Bézier-Polynom in die Form
>>
>> x = a_3*t^3 + a_2*t^2 + a_1*t + a_0
>>
>> um, wobei die a_i natürlich von den Bézier-Knoten x1 ... x4
>> abhängen; für jeden Satz von Bézier-Knoten mußt Du das
>> nur einmal machen.
>
>hab ich gemacht.
OK, fein
>> 2) Für ein gegebenes x berechnest Du die reelle(n) Wurzeln(n) der
>> Polynomgleichung
>>
>> a_3*t^3 + a_2*t^2 + a_1*t + a_0 - x = 0
>>
>> Wie das geht, kannst Du in
>> http://hem.passagen.se/ceem/cubic/index.htm
>> nachlesen; und wenn Du in C (oder C++) programmierst, kann ich Dir
>> dafür auch eine fertige function posten.
>
>Das Programm wird zwar nicht in C entwickelt, allerdings kann ich es auch.
In welcher Sprache denn? In Java dürfte es kein Problem geben, das hat
beinahe das gleiche Gleitkomma-Modell wie C und C++.
Pascal / Delphi / Modula dürfte auch keine Probleme machen.
Bei VB kenne ich mich nicht aus - I love Linux ...
>Wäre nett wenn du mir diese Funktion posten könntest!
OK, mach ich - ich muß sie nur erst wieder ausbuddeln ;-) Ich nehme an,
ich kann sie Dir morgen nachmittag posten ...
Grüße
Hermann
--
Roman Bobik
> >möglicherweise Stellenauslöschungsfehler zu falschen Ergebnissen führen
> >könnten, war mir auch klar, allerdings wenn der Wertebereich/Genauigkeit
von
> >double nicht ausreicht, was soll man da als Programmierer machen...
>
> Hmm, nach meiner Erfahrung reicht die Genauigkeit von "double" *fast*
> immer, wenn man den Rechenverlauf bezüglich numerischer Stabilität
> optimiert und entsprechend programmiert ...
tja genau das Optimieren wollte ich mir ja ersparen... (hab mich auf Derive
verlassen..., aber das macht halt nur seine Aufgabe, dass es eine Formel für
das Ergebnis berechnet und mir das nicht in einzelne Teile aufteilt...) aber
anscheinend gehts nicht anders.
> >Das Programm wird zwar nicht in C entwickelt, allerdings kann ich es
auch.
>
>
> In welcher Sprache denn? In Java dürfte es kein Problem geben, das hat
> beinahe das gleiche Gleitkomma-Modell wie C und C++.
> Pascal / Delphi / Modula dürfte auch keine Probleme machen.
> Bei VB kenne ich mich nicht aus - I love Linux ...
Bingo, VB bzw .NET. Aber das hat auch das gleiche Gleitkomma-Modell. double
nach IEEE, 8 Byte, ein paar Funktionen heißen anders. Ist aber kein Problem
wenn du mir die C Funktion postest. Sollte nach VB portierbar sein ansonsten
kann ich das ja auch in eine DLL per VC++ packen... I love Windows... hoffe
aber dass du mir deswegen trotzdem die Funktion gibst *g*
> >Wäre nett wenn du mir diese Funktion posten könntest!
>
> OK, mach ich - ich muß sie nur erst wieder ausbuddeln ;-) Ich nehme
an,
> ich kann sie Dir morgen nachmittag posten ...
ok, ich komme vor Montag eh nicht zum ausprobieren.
Danke mal bis jetzt!
Roman
Roman Bobik
[...]
sorry, das ganze ist ein bisschen zu hoch für mich. Habe nicht den Vorteil
das ich auf ein Mathematikstudium zurückblicken darf... Brauchte schon eine
ganze weile bis ich im Internet (für mich) brauchbare Informationen zum
"erstellen" der kubischen Bezier-Formel fand... Es gab nur irgendwo die
quadratische und in irgendeiner Diplomarbeit stand wie man sich eine solche
Formel herleitet.
Danke für deine Mühe, ich glaube Hermann hat mir schon helfen können, vor
allem was das umsetzen in ein Programm betrifft.
Roman
Marcus Roeckrath
Roman Bobik wrote:
> sorry, das ganze ist ein bisschen zu hoch für mich. Habe nicht den Vorteil
> das ich auf ein Mathematikstudium zurückblicken darf.
Der Algorithmus ist halb so wild, nur ein wenig grundlegende
Vekorarithmetik. Hier mal ein Beispiel für ein kubisches Beziersegment:
Gegeben 4 Bezierpunkte
b0=(0; 0)
b1=(2; 4,5)
b2=(8,5; 6,5)
b3=(11; 2)
die wir als Vektoren in der normalen Ebene auffassen, und berechnen mal den
Bezierkurvenpunkt zum Parameterwert t=0,6.
Im ersten Schritt reduzieren wir die vier Punkte auf 3 folgendermassen:
b1' = (1-0,6)*b0 + 0,6*b1 = 0,4*b0 + 0,6*b1 = 0,4*(0; 0) + 0,6*(2; 4,5) =
(1,2; 2,7)
b2' = 0,4*(2; 4,5) + 0,6*(8,5; 6,5) = (5,9; 5,7)
b3' = 0,4*(8,5; 6,5) + 0,6*(11; 2) = (10; 3,8)
Nun das gleiche mit den 3 neuen Punkten:
b2'' = 0,4*(1,2; 2,7) + 0,6*(5,9; 5,7) = (4,02; 4,5)
b3'' = 0,4*(5,9; 5,7) + 0,6*(10; 3,8) = (8,36; 4,56)
Aus diesen 2 Punkten wird nun im dritten Schritt der Kurvenpunkt selbst:
b3''' = 0,4*(4,02; 4,5) + 0,6*(8,36; 4,56) = (6,63; 4,54)
Die im weiteren Verlauf meiner vorigen Mail angegebenen Formel zur
Umrechnung der Gewichtspunkte in Bezierpunkte sind auch nichts weiter als
einfache Vektoroperationen.
Der Übergang von den Bezierpunkten zu den Gewichtspunkten ist allerdings
notwendig, um beim Anschluß mehrerer Beziersegmente gewisse Eigenschaften
zu erhalten, die sogenannte C²-Stetigkeit, also stetig, und 2mal stetig
differenzierbar, damit die Kurve keine Spitzen hat bzw. die Beschleunigung
keine Sprünge aufweist.
Über die Herleitung braucht man sich aber als Anwender keine Gedanken zu
machen, der Algorithmis dieser Umrechnung gewährleistet das.
Hermann Kremer
>"Hermann Kremer" <hermann...@online.de> schrieb:
>> Roman Bobik schrieb in Nachricht
>> <3bd0420b$0$33864$5039...@newsreader01.highway.telekom.at>...
Hallo Roman,
[ ... ]
>> Bei VB kenne ich mich nicht aus - I love Linux ...
>
>Bingo, VB bzw .NET. Aber das hat auch das gleiche Gleitkomma-Modell. double
>nach IEEE, 8 Byte, ein paar Funktionen heißen anders. Ist aber kein Problem
>wenn du mir die C Funktion postest. Sollte nach VB portierbar sein ansonsten
>kann ich das ja auch in eine DLL per VC++ packen
Ja, Du brauchst aber nicht nur double, sondern auch float ...
>... I love Windows... hoffe aber dass du mir deswegen trotzdem die
> Funktion gibst *g*
Hmmm, das muß ich mir aber noch sehr überlegen ;-))
OK, hier ist sie
~~~~~~~~~~~~~~~~~~~~~~ SNIP ~~~~~~~~~~~~~~~~~~~~~~~~~~
#ifdef DOCUMENTATION
Die Function InvB3P(a0, a1, a2, a3, x) erwartet die Koeffizienten und
den Funktionswert des Bezier-Polynoms
x = a3*t^3 + a2*t^2 + a1*t + a0
und liefert den Wert von t zurueck. Dabei ist vorausgesetzt, dass das
Polynom streng monoton ist, d.h. dass
x2 > x1 fuer alle 1 >= t2 > t1 >= 0
oder
x2 < x1 fuer alle 1 >= t2 > t1 >= 0
gilt. Ist dies nicht der Fall, so ist der zurueckgelieferte Wert u.U.
falsch.
Sowohl die Argumente als auch der Funktionswerte sind einfach genaue
Gleitkommazahlen (float); intern wird doppelt genau gerechnet (double).
Dies ist Absicht, um richtig gerundete Ein- und Ausgabewerte zu
garantieren.
In TEST1 werden die 4 Knoten (x0,y0) ... (x3,y3) eingegeben, daraus
werden die Koeffizienten der Bezierpolynome
x = a3*t^3 + a2*t^2 + a1*t + a0
y = b3*t^3 + b2*t^2 + b1*t + b0
berechnet und damit die Bezierkurve (x(t), y(t)), 0 <= t <= 1.
Fuer jedes t werden dann aus den entsprechenden Polynomkoeffizienten
und Funktionswerten mittels InvB3P(...) die t-Werte
tau(x) = InvB3P(a3, a2, a1, a0, x)
und
tau(y) = InvB3P(b3, b2, b1, b0, y)
zurueckgerechnet und zusammen mit den Abweichungen bzgl. t ausgedruckt.
Beispiel 1: x(t), y(t) streng monoton
======================================
P0=(1.0,1.0), P1=(100.0,100.0), P2=(150.0,101.0), P3=(200.0,300.0)
a3 = 49.0000, a2 = -147.0000, a1 = 297.0000, a0 = 1.0000
b3 = 296.0000, b2 = -294.0000, b1 = 297.0000, b0 = 1.0000
------+----------------------++------------------++---------------------
t | x y || tau(x) tau(y) || t-tau(x) t-tau(y)
------+----------------------++------------------++---------------------
0.00 | 1.00000 1.00000 || 0.00000 0.00000 ||-2.25e-016 -1.94e-017
0.05 | 15.48863 15.15200 || 0.05000 0.05000 || 0.00e+000 0.00e+000
0.10 | 29.27900 28.05600 || 0.10000 0.10000 || 0.00e+000 0.00e+000
0.15 | 42.40788 39.93400 || 0.15000 0.15000 || 0.00e+000 0.00e+000
0.20 | 54.91200 51.00800 || 0.20000 0.20000 || 0.00e+000 0.00e+000
0.25 | 66.82813 61.50000 || 0.25000 0.25000 || 0.00e+000 0.00e+000
0.30 | 78.19300 71.63200 || 0.30000 0.30000 || 0.00e+000 0.00e+000
0.35 | 89.04338 81.62600 || 0.35000 0.35000 || 0.00e+000 0.00e+000
0.40 | 99.41601 91.70401 || 0.40000 0.40000 || 0.00e+000 0.00e+000
0.45 | 109.34763 102.08801 || 0.45000 0.45000 || 0.00e+000 0.00e+000
0.50 | 118.87501 113.00001 || 0.50000 0.50000 || 0.00e+000 0.00e+000
0.55 | 128.03489 124.66202 || 0.55000 0.55000 || 0.00e+000 0.00e+000
0.60 | 136.86401 137.29602 || 0.60000 0.60000 || 0.00e+000 0.00e+000
0.65 | 145.39914 151.12403 || 0.65000 0.65000 || 0.00e+000 0.00e+000
0.70 | 153.67702 166.36803 || 0.70000 0.70000 || 0.00e+000 0.00e+000
0.75 | 161.73439 183.25004 || 0.75000 0.75000 || 0.00e+000 0.00e+000
0.80 | 169.60802 201.99205 || 0.80000 0.80000 || 0.00e+000 0.00e+000
0.85 | 177.33465 222.81606 || 0.85000 0.85000 || 0.00e+000 0.00e+000
0.90 | 184.95102 245.94408 || 0.90000 0.90000 || 0.00e+000 0.00e+000
0.95 | 192.49390 271.59809 || 0.95000 0.95000 || 0.00e+000 0.00e+000
1.00 | 200.00002 300.00007 || 1.00000 1.00000 || 0.00e+000 0.00e+000
------+----------------------++------------------++---------------------
Beispiel 2: x(t) streng monoton, y(t) nicht monoton
====================================================
P0=(1.0,1.0), P1=(100.0,101.0), P2=(150.0,100.0), P3=(200.0,1.0)
a3 = 49.0000, a2 = -147.0000, a1 = 297.0000, a0 = 1.0000
b3 = 3.0000, b2 = -303.0000, b1 = 300.0000, b0 = 1.0000
------+----------------------++------------------++---------------------
t | x y || tau(x) tau(y) || t-tau(x) t-tau(y)
------+----------------------++------------------++---------------------
0.00 | 1.00000 1.00000 || 0.00000 1.00000 ||-2.25e-016 -1.00e+000
0.05 | 15.48863 15.24288 || 0.05000 0.94952 || 0.00e+000 -9.00e-001
0.10 | 29.27900 27.97300 || 0.10000 0.89909 || 0.00e+000 -7.99e-001
0.15 | 42.40788 39.19263 || 0.15000 0.84871 || 0.00e+000 -6.99e-001
0.20 | 54.91200 48.90400 || 0.20000 0.79839 || 0.00e+000 -5.98e-001
0.25 | 66.82813 57.10938 || 0.25000 0.74811 || 0.00e+000 -4.98e-001
0.30 | 78.19300 63.81100 || 0.30000 0.69789 || 0.00e+000 -3.98e-001
0.35 | 89.04338 69.01113 || 0.35000 0.64771 || 0.00e+000 -2.98e-001
0.40 | 99.41601 72.71200 || 0.40000 0.59759 || 0.00e+000 -1.98e-001
0.45 | 109.34763 74.91588 || 0.45000 0.54751 || 0.00e+000 -9.75e-002
0.50 | 118.87501 75.62500 || 0.50000 0.50000 || 0.00e+000 0.00e+000
0.55 | 128.03489 74.84162 || 0.55000 0.55000 || 0.00e+000 0.00e+000
0.60 | 136.86401 72.56799 || 0.60000 0.60000 || 0.00e+000 0.00e+000
0.65 | 145.39914 68.80637 || 0.65000 0.65000 || 0.00e+000 0.00e+000
0.70 | 153.67702 63.55899 || 0.70000 0.70000 || 0.00e+000 0.00e+000
0.75 | 161.73439 56.82811 || 0.75000 0.75000 || 0.00e+000 0.00e+000
0.80 | 169.60802 48.61598 || 0.80000 0.80000 || 0.00e+000 0.00e+000
0.85 | 177.33465 38.92485 || 0.85000 0.85000 || 0.00e+000 0.00e+000
0.90 | 184.95102 27.75696 || 0.90000 0.90000 || 0.00e+000 0.00e+000
0.95 | 192.49390 15.11458 || 0.95000 0.95000 || 0.00e+000 0.00e+000
1.00 | 200.00001 0.99999 || 1.00000 1.00000 || 0.00e+000 0.00e+000
------+----------------------++------------------++---------------------
Beispiel 3: x(t) nicht monoton, y(t) streng monoton
====================================================
P0=(1.0,1.0), P1=(101.0,100.0), P2=(100.0,150.0), P3=(1.0,200.0)
a3 = 3.0000, a2 = -303.0000, a1 = 300.0000, a0 = 1.0000
b3 = 49.0000, b2 = -147.0000, b1 = 297.0000, b0 = 1.0000
------+----------------------++------------------++---------------------
t | x y || tau(x) tau(y) || t-tau(x) t-tau(y)
------+----------------------++------------------++---------------------
0.00 | 1.00000 1.00000 || 1.00000 0.00000 ||-1.00e+000 -2.25e-016
0.05 | 15.24288 15.48863 || 0.94952 0.05000 ||-9.00e-001 0.00e+000
0.10 | 27.97300 29.27900 || 0.89909 0.10000 ||-7.99e-001 0.00e+000
0.15 | 39.19263 42.40788 || 0.84871 0.15000 ||-6.99e-001 0.00e+000
0.20 | 48.90400 54.91200 || 0.79839 0.20000 ||-5.98e-001 0.00e+000
0.25 | 57.10938 66.82813 || 0.74811 0.25000 ||-4.98e-001 0.00e+000
0.30 | 63.81100 78.19300 || 0.69789 0.30000 ||-3.98e-001 0.00e+000
0.35 | 69.01113 89.04338 || 0.64771 0.35000 ||-2.98e-001 0.00e+000
0.40 | 72.71200 99.41601 || 0.59759 0.40000 ||-1.98e-001 0.00e+000
0.45 | 74.91588 109.34763 || 0.54751 0.45000 ||-9.75e-002 0.00e+000
0.50 | 75.62500 118.87501 || 0.50000 0.50000 || 0.00e+000 0.00e+000
0.55 | 74.84162 128.03489 || 0.55000 0.55000 || 0.00e+000 0.00e+000
0.60 | 72.56799 136.86401 || 0.60000 0.60000 || 0.00e+000 0.00e+000
0.65 | 68.80637 145.39914 || 0.65000 0.65000 || 0.00e+000 0.00e+000
0.70 | 63.55899 153.67702 || 0.70000 0.70000 || 0.00e+000 0.00e+000
0.75 | 56.82811 161.73439 || 0.75000 0.75000 || 0.00e+000 0.00e+000
0.80 | 48.61598 169.60802 || 0.80000 0.80000 || 0.00e+000 0.00e+000
0.85 | 38.92485 177.33465 || 0.85000 0.85000 || 0.00e+000 0.00e+000
0.90 | 27.75696 184.95102 || 0.90000 0.90000 || 0.00e+000 0.00e+000
0.95 | 15.11458 192.49390 || 0.95000 0.95000 || 0.00e+000 0.00e+000
1.00 | 0.99999 200.00001 || 1.00000 1.00000 || 0.00e+000 0.00e+000
------+----------------------++------------------++---------------------
Beispiel 4: x(t) nicht monoton, y(t) nicht monoton
===================================================
P0=(1.0,1.0), P1=(200.0,300.0), P2=(150.0,100.0), P3=(100.0,101.0)
a3 = 249.0000, a2 = -747.0000, a1 = 597.0000, a0 = 1.0000
b3 = 700.0000, b2 = -1497.0000, b1 = 897.0000, b0 = 1.0000
------+----------------------++------------------++---------------------
t | x y || tau(x) tau(y) || t-tau(x) t-tau(y)
------+----------------------++------------------++---------------------
0.00 | 1.00000 1.00000 || 0.00000 0.00000 ||-6.04e-017 -4.35e-017
0.05 | 29.01363 42.19500 || 0.05000 0.05000 || 0.00e+000 0.00e+000
0.10 | 53.47900 76.43000 || 0.10000 0.10000 || 0.00e+000 0.00e+000
0.15 | 74.58288 104.23000 || 0.15000 0.92039 || 0.00e+000 -7.70e-001
0.20 | 92.51200 126.12000 || 0.20000 0.75528 || 0.00e+000 -5.55e-001
0.25 | 107.45313 142.62500 || 0.95011 0.65725 ||-7.00e-001 -4.07e-001
0.30 | 119.59300 154.27000 || 0.86533 0.57983 ||-5.65e-001 -2.80e-001
0.35 | 129.11838 161.58000 || 0.79065 0.51439 ||-4.41e-001 -1.64e-001
0.40 | 136.21600 165.08000 || 0.72345 0.45739 ||-3.23e-001 -5.74e-002
0.45 | 141.07263 165.29500 || 0.66219 0.45000 ||-2.12e-001 0.00e+000
0.50 | 143.87500 162.75000 || 0.60587 0.50000 ||-1.06e-001 0.00e+000
0.55 | 144.80988 157.96999 || 0.55378 0.55000 ||-3.78e-003 0.00e+000
0.60 | 144.06400 151.47999 || 0.60000 0.60000 || 0.00e+000 0.00e+000
0.65 | 141.82412 143.80498 || 0.65000 0.65000 || 0.00e+000 0.00e+000
0.70 | 138.27699 135.46998 || 0.70000 0.70000 || 0.00e+000 0.00e+000
0.75 | 133.60936 126.99998 || 0.75000 0.75000 || 0.00e+000 0.00e+000
0.80 | 128.00798 118.91998 || 0.80000 0.80000 || 0.00e+000 0.00e+000
0.85 | 121.65961 111.75498 || 0.85000 0.85000 || 0.00e+000 0.00e+000
0.90 | 114.75098 106.02999 || 0.90000 0.90000 || 0.00e+000 0.00e+000
0.95 | 107.46885 102.26999 || 0.95000 0.95000 || 0.00e+000 0.00e+000
1.00 | 99.99999 101.00000 || 1.00000 0.99500 || 0.00e+000 5.00e-003
------+----------------------++------------------++---------------------
In TEST2 wird der x-Bereich x0 <= x <= x3 durchlaufen, fuer jeden
x-Wert der zugehoerige t-Wert mittels InvB3P(...) bestimmt und damit
der entsprechende y-Wert berechnet.
Beispiel 1: x(t), y(t) streng monoton
======================================
P0=(1.0,1.0), P1=(100.0,100.0), P2=(150.0,101.0), P3=(200.0,300.0)
a3 = 49.0000, a2 = -147.0000, a1 = 297.0000, a0 = 1.0000
b3 = 296.0000, b2 = -294.0000, b1 = 297.0000, b0 = 1.0000
----------------------------------------------------------
x = 1.00000 t = 0.00000 y = 1.00000
x = 10.95000 t = 0.03407 y = 10.78914
x = 20.90000 t = 0.06933 y = 20.27578
x = 30.85000 t = 0.10586 y = 29.49579
x = 40.80000 t = 0.14374 y = 38.49626
x = 50.75000 t = 0.18309 y = 47.33833
x = 60.70000 t = 0.22399 y = 56.10061
x = 70.65000 t = 0.26655 y = 64.88342
x = 80.60000 t = 0.31090 y = 73.81389
x = 90.55000 t = 0.35713 y = 83.05199
x = 100.50000 t = 0.40535 y = 92.79745
x = 110.45000 t = 0.45568 y = 103.29730
x = 120.40000 t = 0.50819 y = 114.85347
x = 130.35000 t = 0.56294 y = 127.82937
x = 140.30000 t = 0.61993 y = 142.65344
x = 150.25000 t = 0.67913 y = 159.81734
x = 160.20000 t = 0.74038 y = 179.86520
x = 170.15000 t = 0.80348 y = 203.37106
x = 180.09999 t = 0.86808 y = 230.90292
x = 190.04999 t = 0.93376 y = 262.97575
x = 199.99999 t = 1.00000 y = 299.99989
Beispiel 2: x(t) streng monoton, y(t) nicht monoton
====================================================
P0=(1.0,1.0), P1=(100.0,101.0), P2=(150.0,100.0), P3=(200.0,1.0)
a3 = 49.0000, a2 = -147.0000, a1 = 297.0000, a0 = 1.0000
b3 = 3.0000, b2 = -303.0000, b1 = 300.0000, b0 = 1.0000
----------------------------------------------------------
x = 1.00000 t = 0.00000 y = 1.00000
x = 10.95000 t = 0.03407 y = 10.86931
x = 20.90000 t = 0.06933 y = 20.34288
x = 30.85000 t = 0.10586 y = 29.36496
x = 40.80000 t = 0.14374 y = 37.87131
x = 50.75000 t = 0.18309 y = 45.78770
x = 60.70000 t = 0.22399 y = 53.02842
x = 70.65000 t = 0.26655 y = 59.49453
x = 80.60000 t = 0.31090 y = 65.07204
x = 90.55000 t = 0.35713 y = 69.63012
x = 100.50000 t = 0.40535 y = 73.01954
x = 110.45000 t = 0.45568 y = 75.07170
x = 120.40000 t = 0.50819 y = 75.59882
x = 130.35000 t = 0.56294 y = 74.39607
x = 140.30000 t = 0.61993 y = 71.24655
x = 150.25000 t = 0.67913 y = 65.93025
x = 160.20000 t = 0.74038 y = 58.23790
x = 170.15000 t = 0.80348 y = 47.98950
x = 180.09999 t = 0.86808 y = 35.05623
x = 190.04999 t = 0.93376 y = 19.38205
x = 199.99999 t = 1.00000 y = 1.00005
Beispiel 3: x(t) nicht monoton, y(t) streng monoton
====================================================
P0=(1.0,1.0), P1=(101.0,100.0), P2=(100.0,150.0), P3=(1.0,200.0)
a3 = 3.0000, a2 = -303.0000, a1 = 300.0000, a0 = 1.0000
b3 = 49.0000, b2 = -147.0000, b1 = 297.0000, b0 = 1.0000
----------------------------------------------------------
x = 1.00000 t = 1.00000 y = 200.00000
x = 6.00000 t = 0.98287 y = 197.43095
x = 11.00000 t = 0.96513 y = 194.76676
x = 16.00000 t = 0.94668 y = 191.99442
x = 21.00000 t = 0.92744 y = 189.09803
x = 26.00000 t = 0.90731 y = 186.05789
x = 31.00000 t = 0.88614 y = 182.84903
x = 36.00000 t = 0.86375 y = 179.43914
x = 41.00000 t = 0.83991 y = 175.78511
x = 46.00000 t = 0.81427 y = 171.82729
x = 51.00000 t = 0.78638 y = 167.47901
x = 56.00000 t = 0.75548 y = 162.60595
x = 61.00000 t = 0.72034 y = 156.97896
x = 66.00000 t = 0.67847 y = 150.14253
x = 71.00000 t = 0.62330 y = 140.87591
x = 76.00000 t = 100.00253 y =47563339.83264
x = 81.00000 t = 100.00269 y =47563581.69516
x = 86.00000 t = 100.00286 y =47563823.55851
x = 91.00000 t = 100.00303 y =47564065.42267
x = 96.00000 t = 100.00320 y =47564307.28766
x = 101.00000 t = 100.00336 y =47564549.15346
Beispiel 4: x(t) nicht monoton, y(t) nicht monoton
===================================================
P0=(1.0,1.0), P1=(200.0,300.0), P2=(150.0,100.0), P3=(100.0,101.0)
a3 = 249.0000, a2 = -747.0000, a1 = 597.0000, a0 = 1.0000
b3 = 700.0000, b2 = -1497.0000, b1 = 897.0000, b0 = 1.0000
----------------------------------------------------------
x = 1.00000 t = 0.00000 y = 1.00000
x = 10.95000 t = 0.01703 y = 15.84299
x = 20.90000 t = 0.03483 y = 30.45920
x = 30.85000 t = 0.05352 y = 44.82688
x = 40.80000 t = 0.07321 y = 58.92003
x = 50.75000 t = 0.09406 y = 72.70707
x = 60.70000 t = 0.11626 y = 86.14883
x = 70.65000 t = 0.14007 y = 99.19563
x = 80.60000 t = 0.16585 y = 111.78246
x = 90.55000 t = 0.19408 y = 123.82093
x = 100.50000 t = 0.99667 y = 100.99667
x = 110.45000 t = 0.92976 y = 103.52184
x = 120.40000 t = 0.85938 y = 110.55483
x = 130.35000 t = 0.77999 y = 122.07349
x = 140.30000 t = 0.67362 y = 139.91778
x = 150.25000 t = 1.90810 y =1125.18001
x = 160.20000 t = 1.92882 y =1184.90482
x = 170.15000 t = 1.94840 y =1243.35892
x = 180.09999 t = 1.96699 y =1300.70127
x = 190.04999 t = 1.98472 y =1357.05935
x = 199.99999 t = 2.00167 y =1412.53872
Bei den Beispielen 3 und 4 sieht man an den berechneten t-Werten
deutlich den Effekt der nicht strengen Monotonizitaet von x(t).
#endif /* DOCUMENTATION */
/*
* 17. 5. 95 /Kr.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifndef PI
#define PI 3.141592653589
#endif
float InvB3P(float a0, float a1, float a2, float a3, float x)
/*
****************************************************************************
*
* Bestimmung der (einer) reellen Wurzel des kubischen Polynoms
*
* a3*t^3 + a2*t^2 + a1*t + a0 = x
*
* mittels einer trigonometrischen Version der Cardano'schen Formel.
*
****************************************************************************
*/
{
float c;
double h, p, q, D, R, S, F, t;
double w1 = 2.0*PI/3.0;
double w2 = 4.0*PI/3.0;
c = (float)(1.0 + a3); /* Quick & dirty hack to cover the case */
if (c == (float)1.0) /* that the 3-rd order polynomial dege- */
a3 = (float)1.0e-6; /* nerates to a 2-rd or 1-st order one */
h = a2/3.0/a3;
p = (3.0*a1*a3 - a2*a2)/3.0/a3/a3;
q = (2.0*a2*a2*a2 - 9.0*a1*a2*a3 - 27.0*a3*a3*(x - a0))/27.0/a3/a3/a3;
c = (float)(1.0 + p); /* Check for p being too near to zero */
if (c == (float)1.0)
{
c = (float)(1.0 + q); /* Check for q being too near to zero */
if (c == (float)1.0)
return( (float)(-h) );
t = -exp(log(fabs(q))/3.0);
if (q < 0.0)
t = -t;
t -= h;
return( (float)t );
}
R = sqrt(fabs(p)/3.0);
S = fabs(q)/2.0/R/R/R;
R = -2.0*R;
if (q < 0.0)
R = -R;
if (p < 0.0)
{
D = p*p*p/27.0 + q*q/4.0;
if (D <= 0.0)
{
F = acos(S)/3.0;
t = R*cos(F + w2) - h;
if ((t < -0.00005) || (t > 1.00005))
{
t = R*cos(F + w1) - h;
if ((t < -0.00005) || (t > 1.00005))
{
t = R*cos(F) - h;
if (t < -0.00005)
t = 0.0;
if (t > 1.00005)
t = 1.0;
}
}
}
else
t = R*cosh(log(S + sqrt((S + 1.0)*(S - 1.0)))/3.0) - h; /* arcosh */
}
else
t = R*sinh(log(S + sqrt(S*S + 1.0))/3.0) - h; /* arsinh */
return( (float)t );
}
main()
{
float x0, y0, x1, y1, x2, y2, x3, y3;
float a0, a1, a2, a3;
float b0, b1, b2, b3;
float t, x, y, taux, tauy, xmin, xmax;
AGAIN:
printf("x0/y0: ");
if (scanf("%f %f", &x0, &y0) != 2)
exit(1);
printf("x1/y1: ");
if (scanf("%f %f", &x1, &y1) != 2)
exit(1);
printf("x2/y2: ");
if (scanf("%f %f", &x2, &y2) != 2)
exit(1);
printf("x3/y3: ");
if (scanf("%f %f", &x3, &y3) != 2)
exit(1);
printf("\n");
a0 = x0; b0 = y0;
a1 = 3.0*(x1 - x0); b1 = 3.0*(y1 - y0);
a2 = 3.0*(x2 - 2.0*x1 + x0); b2 = 3.0*(y2 - 2.0*y1 + y0);
a3 = x3 - 3.0*x2 + 3.0*x1 - x0; b3 = y3 - 3.0*y2 + 3.0*y1 - y0;
printf("P0=(%.1f,%.1f), P1=(%.1f,%.1f), P2=(%.1f,%.1f), P3=(%.1f,%.1f)\n",
x0,y0, x1,y1, x2,y2, x3,y3);
printf("a3 = %10.4f, a2 = %10.4f, a1 = %10.4f, a0 = %10.4f\n",
a3, a2, a1, a0);
printf("b3 = %10.4f, b2 = %10.4f, b1 = %10.4f, b0 = %10.4f\n",
b3, b2, b1, b0);
#ifdef TEST1
printf("------+----------------------++------------------++---------------------\n");
printf(" t | x y || tau(x) tau(y) || t-tau(x) t-tau(y)
\n");
printf("------+----------------------++------------------++---------------------\n");
for (t = 0.0; t <= 1.01; t += 0.05)
{
x = ((a3*t + a2)*t + a1)*t + a0;
y = ((b3*t + b2)*t + b1)*t + b0;
taux = InvB3P(a0, a1, a2, a3, x);
tauy = InvB3P(b0, b1, b2, b3, y);
printf("%5.2f |%10.5f %10.5f ||%8.5f %8.5f ||%10.2e %10.2e\n",
(float)t,
(float)x, (float)y,
(float)taux, (float)tauy, (float)(t-taux), (float)(t-tauy));
}
#endif /* TEST1 */
#ifdef TEST2
printf("\n----------------------------------------------------------\n");
xmin = min(x0, min(x1, min(x2, x3)));
xmax = max(x0, max(x1, max(x2, x3)));
for (x = xmin; x <= xmax+0.001; x += (xmax-xmin)/20.0)
{
t = InvB3P(a0, a1, a2, a3, x);
y = ((b3*t + b2)*t + b1)*t + b0;
printf("x =%10.5f t = %8.5f y =%10.5f\n",
(float)x, (float)t, (float)y);
}
#endif /* TEST2 */
goto AGAIN;
}
~~~~~~~~~~~~~~~~~~~~~~ SNIP ~~~~~~~~~~~~~~~~~~~~~~~~~~
Grüße
Hermann
--
Marcus Roeckrath
Marcus Roeckrath wrote:
>> sorry, das ganze ist ein bisschen zu hoch für mich. Habe nicht den
>> Vorteil das ich auf ein Mathematikstudium zurückblicken darf.
>
> Der Algorithmus ist halb so wild, nur ein wenig grundlegende
> Vekorarithmetik. Hier mal ein Beispiel für ein kubisches Beziersegment:
Wer das Programm DynaGeo hat, kann die beschriebene Konstruktion der
kubischen Bezierkurve dort als Beispieldatei finden.
--
Gruss Marcus
Roman Bobik
> Roman Bobik schrieb in Nachricht
> <3bd16897$0$9880$6e36...@newsreader02.highway.telekom.at>...
> >"Hermann Kremer" <hermann...@online.de> schrieb:
> >> Roman Bobik schrieb in Nachricht
> >> <3bd0420b$0$33864$5039...@newsreader01.highway.telekom.at>...
>
double
> >nach IEEE, 8 Byte, ein paar Funktionen heißen anders. Ist aber kein
Problem
> >wenn du mir die C Funktion postest. Sollte nach VB portierbar sein
ansonsten
> >kann ich das ja auch in eine DLL per VC++ packen
>
> Ja, Du brauchst aber nicht nur double, sondern auch float ...
Natürlich hat VB auch float (= Single)...
> >... I love Windows... hoffe aber dass du mir deswegen trotzdem die
> > Funktion gibst *g*
>
> Hmmm, das muß ich mir aber noch sehr überlegen ;-))
>
> OK, hier ist sie
>
> [...]
Danke für den Code! Hab ihn nach VB portiert und funktioniert perfekt. Ich
konnte nun damit die selben Ergebnisse erziehlen wie mit der 5kb-Formel
oben...., jedoch ohne die Nebeneffekte. Kurz: So gibt es keine Fehler mehr
in Grenzbereichen.
Nochmals Danke, weil alleine wäre ich vermutlich nie drauf gekommen...
Hatte anfangs nur Probleme, da ich meine x umgekehrt Indiziert hatte, aber
Dank deiner Beispiele bin ich schnell auf den Fehler gekommen...
Wer Interesse am nach Basic/Visual Basic portierten Code hat, hier ist er:
'Hyperbolic Cosine
Function CosH(ByVal x As Double) As Double
CosH = (Exp(x) + Exp(-x)) / 2
End Function
'Hyperbolic Sine
Function SinH(ByVal x As Double) As Double
SinH = (Exp(x) - Exp(-x)) / 2
End Function
'Inverse Hyperbolic Cos
Function ArcCosH(ByVal x As Double) As Double
ArcCosH = Log(x + Sqr(x * x - 1))
End Function
'Inverse Hyperbolic Sine
Function ArcSinH(ByVal x As Double) As Double
ArcSinH = Log(x + Sqr(x * x + 1))
End Function
'Inverse Cosine
Function ArcCos(ByVal x As Double) As Double
On Error Resume Next
ArcCos = Atn(-x / Sqr(-x * x + 1)) + 2 * Atn(1)
End Function
Function PI() As Double
PI = 4 * Atn(1)
End Function
Public Function SolveSymbolicCubicRoot(ByVal a0 As Single, ByVal a1 As
Single, ByVal a2 As Single, ByVal a3 As Single, ByVal x As Single) As Single
'***************************************************************************
*
'*
'* Bestimmung der (einer) reellen Wurzel des kubischen Polynoms
'*
'* a3*t^3 + a2*t^2 + a1*t + a0 = x
'*
'* mittels einer trigonometrischen Version der Cardano'schen Formel.
'*
'* [--- Original-Code (in C) von Hermann Kremer
hermann...@online.de) ---]
'***************************************************************************
*
Dim c As Single
Dim h As Double, p As Double, q As Double, D As Double, R As Double, _
S As Double, F As Double, t As Double, w1 As Double, w2 As Double
w1 = 2# * PI / 3#
w2 = 4# * PI / 3#
c = CSng(1# + a3) 'Quick & dirty hack to cover the case
If c = CSng(1#) Then _
a3 = CSng(0.000001) 'that the 3-rd order polynomial dege-
'nerates to a 2-rd or 1-st order one
h = a2 / 3# / a3
p = (3# * a1 * a3 - a2 * a2) / 3# / a3 / a3
q = (2# * a2 * a2 * a2 - 9# * a1 * a2 * a3 - 27# * a3 * a3 * (x - a0)) /
27# / a3 / a3 / a3
c = CSng(1# + p) 'Check for p being too near to zero
If c = CSng(1#) Then
c = CSng(1# + q) 'Check for q being too near to zero
If c = CSng(1#) Then
SolveSymbolicCubicRoot = CSng(-h)
Exit Function
End If
t = -Exp(Log(Abs(q)) / 3#)
If q < 0# Then t = -t
t = t - h
SolveSymbolicCubicRoot = CSng(t)
Exit Function
End If
R = Sqr(Abs(p) / 3#)
S = Abs(q) / 2# / R / R / R
R = -2# * R
If q < 0# Then R = -R
If p < 0# Then
D = p * p * p / 27# + q * q / 4#
If D <= 0# Then
F = ArcCos(S) / 3#
t = R * Cos(F + w2) - h
If (t < -0.00005) Or (t > 1.00005) Then
t = R * Cos(F + w1) - h
If (t < -0.00005) Or (t > 1.00005) Then
t = R * Cos(F) - h
If t < -0.00005 Then t = 0#
If t > 1.00005 Then t = 1#
End If
End If
Else
t = R * CosH(ArcCosH(S) / 3#) - h
End If
Else
t = R * SinH(ArcSinH(S) / 3#) - h 'arsinh
End If
SolveSymbolicCubicRoot = CSng(t)
End Function
Private Function GetTFromBezierEquation(ByVal X As Double, x4 As Double, x3
As Double, x2 As Double, x1 As Double) As Double
Dim a(3) As Double
a(0) = x1 - 3 * x2 + 3 * x3 - x4
a(1) = 3 * x2 - 6 * x3 + 3 * x4
a(2) = 3 * x3 - 3 * x4
a(3) = x4
GetTFromBezierEquation = SolveSymbolicCubicRoot(a(3), a(2), a(1), a(0),
X)
End Function
Die Lösung ist zusätzlich noch um einiges schneller als meine fehlerhafte
alte.
Danke!
Roman Bobik