Berechnung einer Ellipse und deren Halbachsen aus 5 Punkten

22.04.2013, 18:48

Ultimade

Berechnung einer Ellipse und deren Halbachsen aus 5 Punkten

Meine Frage:
Hallo,

ich möchte eine Ellipse aus fünf Punkten bestimmen, und dann die Länge der Haupt und Nebenachse berechnen.

Die Ellipse liegt frei im Koordinatensystem. D.h. der Mittelpunkt der Ellipse liegt nicht im Zentrum des KOS und die Halbachsen sind nicht parallel mit der x bzw y Achse.

Meine Ideen:
Soweit bin ich bis jetzt.
Ich habe in der allgemeinen Gleichung:

ax²+bxy+cy²+dx+ey+f=0

a=1 gesetzt.

Anschließend habe ich meine fünf Punkte nach der Reihe eingesetzt und damit 5 Gleichungen erhalten. Diese dann gelöst, so dass ich für b,c,d,e und f meine Werte erhalten habe. Eingesetzt in die allg. Gleichung erhalte ich meine Ellipsengleichung.
Soweit so gut. Jetzt kann ich zwar jeden Punkt auf der Ellipse bestimmen kenne aber noch nicht die Länge der Hauptachse oder kann die Brennpunkte bestimmen.
Nun habe ich gelesen:
http://www.herder-oberschule.de/madincea/aufg1213/kurven.pdf
dass man die Ellipse in den Nullpunkt verschieben muss und anschließend drehen. (Richtig so?)

Auf Seite 4 des PDF`s (LINK) ist beschrieben,
dass man x`=x-u und y´=y-v einsetzt und danach die Formel löst, indem man die linearen Terme elliminiert.

Kann mir das jemand zeigen? Gerne mit den Werten die im Link verwendet werden.

Merci Ultimade

22.04.2013, 23:50

Lampe16

Auf diesen Beitrag antworten »

Hier kann die Hauptachsentransformation helfen.

23.04.2013, 08:12

Ultimade

Auf diesen Beitrag antworten »

Ja sowas meine ich. Weißt du evtl. wo es zu sowas ein ausführliches Beispiel gibt?

23.04.2013, 08:22

Lampe16

Auf diesen Beitrag antworten »

Hier ist ein Beispiel. Keine Gewähr für die Qualität!

23.04.2013, 08:32

Ultimade

Auf diesen Beitrag antworten »

Danke.

Ich kämpfe mich mal durch und schau wie weit ich komme.

23.04.2013, 08:49

Lampe16

Auf diesen Beitrag antworten »

Meyberg/vachenauer Band 1 kommt noch für die Theorie in Frage. Der schreibt präzise, aber verständlich.

23.04.2013, 15:42

Ultimade

Auf diesen Beitrag antworten »

So ich habe jetzt die Matrix diagonalisiert und aus den Eigenräumen die Eigenvektoren bestimmt.

Nun werden hier: http://www.stud.uni-hannover.de/~fmodler...nsformation.pdf
die Vetktoren normiert. Mit der Orthonormalbasis des R² im Beispiel 1/sprt(10)
Frage: Kann man diesen Wert immer für eine zweidimensionale Ellipse verwenden?

Weiter geht es mit der Substitution. Hier stecke ich leider fest.
Es wird der Vektor D aufgestellt aus transponierte S, S und A.
Eigentlich doch nur die drei Matrizen miteinander multiplizieren oder?
Ich komme auf jeden fall nicht auf die Lösung von D mit
-5/2 0
0 5/2

23.04.2013, 16:23

Lampe16

Auf diesen Beitrag antworten »

Ich müsste mich da erst wieder einarbeiten. Vielleicht kann hier ein anderer Helfer eingreifen.

23.04.2013, 17:15

Ultimade

Auf diesen Beitrag antworten »

Das mit Der Matrix D hat sich erledigt. (Matrizen sind hald meistens nicht kommutativ ;-) )

Bin jetzt in Zeile fünf im Unterpunkt Substitution.
Hier wurden die Matrizen eingesetzt.
Warum wurde die Matrix S vertikal gespiegelt?

23.04.2013, 18:35

Lampe16

Auf diesen Beitrag antworten »

RE: Berechnung einer Ellipse und deren Halbachsen aus 5 Punkten

Zitat:

Original von Ultimade
ich möchte eine Ellipse aus fünf Punkten bestimmen, und dann die Länge der Haupt und Nebenachse berechnen.

Zwischenfrage am Rande: Wenn die fünf Punkte z.B. die Ecken und der Mittelpunkt eines Rechtecks sind, gibt es sicher keine passende Ellipse. Woher weißt du, ob die fünf Punkte überhaupt mit einer Ellipse verträglich sind?

Bei Meyberg/vachenauer Band 1 habe ich noch einen alternativen, sehr kurzen Transformationsweg gefunden.

Der Scilab-Code demonstriert ihn an zwei dir bekannten Beispielen.

//Hyperbel-Beispiel aus http://www.stud.uni-hannover.de/~fmodler...nsformation.pdf
A=[-2 1.5; 1.5 2];
a=[-3;1];
Beta=-4;
m=inv(A)*(-a/2)//m=[-0.6;0.2]//Zentrum der Quadrik
[EV,EW]=spec(A)//Eigenwerte EW(1,1)=-2.5, EW(2,2)=2.5 (spec(A) berechnet Eigenvektoren und -werte)
gam=Beta+1/2*a'*m//gam=-3
//Die Normalform der Quadrik (bis auf die Normierung des konstanten Glieds) ist
// EW(1,1)*z1^2+EW(2,2)*z2^2+gam=0 in einem gedrehten und verschobenen z1-z2-Koordinatensystem

//Kreiskegel-Beispiel aus http://de.wikipedia.org/wiki/Hauptachsentransformation
A=[1 -1 -1; -1 1 -1; -1 -1 1];
a=[2;2;2];
Beta=-3;
m=inv(A)*(-a/2)//m=[1;1;1]//Zentrum der Quadrik
[EV,EW]=spec(A)//Eigenwerte EW(1,1)=-1, EW(2,2)=2, EW(3,3)=2 (spec(A) berechnet Eigenvektoren und -werte)
gam=Beta+1/2*a'*m//gam=0
//Die Normalform der Quadrik ist
// EW(1,1)*z1^2+EW(2,2)*z2^2++EW(3,3)*z3^2+gam=0 in einem gedrehten und verschobenen z1-z2-z3-Koordinatensystem

Der Weg sollte auch für deine Ellipse funktionieren.

24.04.2013, 18:20

Ultimade

Auf diesen Beitrag antworten »

Zitat:

Berechtigte Frage. Mathematisch müsste man das zuerst kontrollieren. Praktisch sind meine Punkte allerdings aus einer Messung, so dass ich davon ausgehe, dass nicht alle Punkte auf einer Geraden liegen oder ähnliches. Sondern eine Kreisähnliche Ellipse bilden.

Der Scilab code ist super. Danke schön. Bekomme so die Werte super raus.

Zu meiner Frage

Zitat:

Nun werden hier: http://www.stud.uni-hannover.de/~fmodler...nsformation.pdf die Vetktoren normiert. Mit der Orthonormalbasis des R² im Beispiel 1/sprt(10) Frage: Kann man diesen Wert immer für eine zweidimensionale Ellipse verwenden?

würde ich inzwischen sagen NEIN kann man nicht. Ich vermute, dass man ihn über den Phytagoras errechnet. Wie genau und weshalb????

Zu:

Zitat:

Warum wurde die Matrix S vertikal gespiegelt?

Ich nehme an, da man die Eigenvektoren in ihrer Richtung Ändern kann, und die Matrix S aus diesen Eigenvektoren besteht ist das zulässig.

Nun noch eine vorerst letzte Frage:
Im Scilab wird gam Berechnet aus:
Beta+1/2*a'*m
Wenn wir beim ersten Beispiel bleiben:
Beta=-4
m=[-0,6;0,2]
was ist a'

diese hat sich im Moment auch erledigt es ist die Transversale von a.

Habe jetzt einmal alles soweit rausgebracht. Warum man das so rechnet weiß ich allerdings noch nicht so ganz verwirrt

. Aber dafür das ich Begriffe wei Orthogonalbasis zum ersten mal gehört habe, liefs ganz gut.
Muss das noch sauber aufschreiben, evtl kommt dann da noch die ein oder andere Frage.

Vielen Dank nochmal Lampe16.

24.04.2013, 21:19

Lampe16

Auf diesen Beitrag antworten »

Zitat:

Original von Ultimade
Praktisch sind meine Punkte allerdings aus einer Messung, so dass ich davon ausgehe, dass nicht alle Punkte auf einer Geraden liegen oder ähnliches. Sondern eine Kreisähnliche Ellipse bilden.

Wenn du Daten in das Korsett einer Ellipse zwingen willst, bietet sich ein Ausgleichsverfahren an. Dann wird die Ellipse gefunden, welche die Messpunkte am besten annähert. Das kann evtl. die ganze Aufgabenstellung umkrempeln.

Zitat:

Original von Ultimade
Der Scilab code ist super. Danke schön. Bekomme so die Werte super raus.

Freut mich - und dass du zu Scilab gefunden hast. Die Beschäftigung damit lohnt sich.

Zitat:

Original von Ultimade
Zu meiner Frage... die Vetktoren normiert. Mit der Orthonormalbasis des R² im Beispiel 1/sqrt(10) Frage: Kann man diesen Wert immer für eine zweidimensionale Ellipse verwenden?

Nein, Normierung bedeutet hier, das der Vektor zu einem Einheitsvektor getrimmt wird.

Zitat:

Original von Ultimade
Warum wurde die Matrix S vertikal gespiegelt?

Das kann ich ohne die ganze Transformation aufzutischen nicht beantworten. Das Ganze ist relativ komplex. Das beste, was du machen kannst, ist ein Lehrbuch deines Vertrauens zu lesen und dann zu schauen, ob du mit einem Beispiel klar kommst. Ich bevorzuge den schon genannten Meyberg/Vachenauer.

Zitat:

Original von Ultimade
Warum man das so rechnet weiß ich allerdings noch nicht so ganz ?

Nochmal mein Hinweis auf ein Lehrbuch, wo man die Zusammenhänge sieht. Das "Schnellverfahren", das ich in Scilab kodiert habe, führt Meyberg nebenher unter "Bemerkung" an. Das wird nicht weiter erklärt.
Ein wesentliches Kapitel, das bei der Transformation von Quadriken auf Normalform gebraucht wird, ist das Matrizen-Eigenwertproblem. Alles keine leichte Kost!

25.04.2013, 10:24

Ultimade

Auf diesen Beitrag antworten »

Zitat:

Genau so etwas will ich am Ende machen.
Mein erster gedanke war:
Ich Messe 5 Punkte und dann gibt es genau eine Ellipse die da reinpasst. Aus der kann ich anschließend die Halbachsen ermitteln.
Praktisch kann ich auch mehr als 5 Punkte Messen, z.B. 10. Wenn man aus diesen dann eine Ellipse annähert, würde natürlich ein "ausreißer" in der Messung nicht so stark ins Gewicht fallen.

Wie müsste man den dann vorgehen?

P.S. Buch ist schon bestellt

25.04.2013, 12:44

Lampe16

Auf diesen Beitrag antworten »

Zitat:

Original von Ultimade
Genau so etwas will ich am Ende machen.
Mein erster gedanke war: Ich Messe 5 Punkte und dann gibt es genau eine Ellipse die da reinpasst. Aus der kann ich anschließend die Halbachsen ermitteln.
Praktisch kann ich auch mehr als 5 Punkte Messen, z.B. 10. Wenn man aus diesen dann eine Ellipse annähert, würde natürlich ein "ausreißer" in der Messung nicht so stark ins Gewicht fallen.
Wie müsste man den dann vorgehen?

Ich stelle es mir etwa so vor:
$\begin{eqnarray*} x \end{eqnarray*}$ und $\begin{eqnarray*} y \end{eqnarray*}$ sind die Messwerte, in Scilab jeweils ein Vektor mit mindestens 5 Elementen.
Jetzt unterstelle ich, dass es einen vektorwertigen Parameter $\begin{eqnarray*} p=(a,b,\alpha,m_x,m_y) \end{eqnarray*}$ gibt, der die $\begin{eqnarray*} (x,y) \end{eqnarray*}$ -Werte in eine Näherungsellipse in Normallage verwandelt, d. h. dreht und verschiebt.
Dabei bedeuten:
$\begin{eqnarray*} a,b: \end{eqnarray*}$ Halbachsen der Näherungsellipse
$\begin{eqnarray*} \alpha: \end{eqnarray*}$ Winkel zwischen den großen Halbachsen in Original- und Normallage
$\begin{eqnarray*} m_x,m_y: \end{eqnarray*}$ Verschiebung von Original- nach Normallage

Man schätzt $\begin{eqnarray*} p \end{eqnarray*}$ . Das tut später der Minimumsucher optim oder fminsearch in Scilab.
Durch Multiplikation der $\begin{eqnarray*} (x,y) \end{eqnarray*}$ mit der $\begin{eqnarray*} \alpha \end{eqnarray*}$ -Drehmatrix und anschließender Verschiebung um $\begin{eqnarray*} (m_x,m_y) \end{eqnarray*}$ gelangt man zu den Punkten $\begin{eqnarray*} (x_n,y_n) \end{eqnarray*}$ , die eine Ellipse in Normallage beschreiben, wenn der Parameter richtig war.
Der Defekt
$\begin{eqnarray*} d=b^2x_n^2+a^2y_n^2-a^2b^2 \end{eqnarray*}$
ist ein Maß für die Güte des Parameters. Wäre $\begin{eqnarray*} d \end{eqnarray*}$ gleich null, erfüllten die $\begin{eqnarray*} (x,y) \end{eqnarray*}$ und die $\begin{eqnarray*} (x_n,y_n) \end{eqnarray*}$ eine Ellipse genau.
Die klassische Hauptachsentransformation fällt hier gar nicht mehr an, weil $\begin{eqnarray*} a \end{eqnarray*}$ und $\begin{eqnarray*} b \end{eqnarray*}$ Elemente des Parameters sind. Das hatte ich mit dem "Umkrempeln" gemeint.

Aber Vorsicht: Dies ist nur eine nicht erprobte Idee. Gib mir mal 5 bis N<100 Beispielwerte für $\begin{eqnarray*} x \end{eqnarray*}$ und $\begin{eqnarray*} y \end{eqnarray*}$ an. Dann werde ich es probieren und soviel von der Lösung verraten, wie du anforderst.

Zitat:

Original von Ultimade
P.S. Buch ist schon bestellt

Gute Entscheidung!

25.04.2013, 14:24

Ultimade

Auf diesen Beitrag antworten »

Klinkt spannend.
Theoretisch kann ich dir folgen.

Eine echte Messung habe ich noch nicht ganz. Da fehlen mir noch ein paar Teile.
Ich habe aber im CAD mal eine Ellipse gezeichnet aus 5 Punkten siehe Anhang.
Die Punkte sind: Zn (xn/yn)

xn yn
-3 0,5
-9 1
-15 -1
-19 -5
-5 -2

=> Mittelpunkt (~-4,167/~-9,995)

Habe dann einfach mal Punkte in die näher der Ellipse gezeichnet, so wie sie bei einer Messung entstehen würden.

xn yn
-3,4 0,8
-8,3 0,9
-11,7 0,5
-14 -0,8
-18,1 -2,9
-18,4 -7,5
-14,6 -9,7
-9,5 -8,8
-5,9 -8,3
-3,4 -5,9
-0,6 -2,9

Ich denke mit dem müsstest du was anfangen können. Wenn nicht einfach melden.

Noch ne andere Frage:
Kann ich den Scilab code auch in einen C oder C++ code umwandeln?

25.04.2013, 19:16

Lampe16

Auf diesen Beitrag antworten »

Mein Ansatz funktioniert. Der aktive Teil des Codes hat ca. 15 Zeilen. Ich komme heute Abend später nochmal drauf zurück, wenn du nachfragst.

Wegen "C": Damit habe ich mich bisher nicht beschäftigt, aber vermutlich geht das.

25.04.2013, 21:23

Lampe16

Auf diesen Beitrag antworten »

So sieht das Ergebnis für deine Daten aus.
[attach]29757[/attach]

29.04.2013, 12:49

Lampe16

Auf diesen Beitrag antworten »

Hier noch der Code (sce-Datei) zum Bild:

//Ausgleichsverfahren für Ellipse
clear, clc(), mode(0)
function [xR,yR]=!Dz2(x,y,w)//Drehen um den Winkel w von Punkten in der x-y-Ebene um die z-Achse im rechtshändigen Koo.-System
//x,y: Koordinaten der Punkte, Matrizen gleichen Typs
//w: Winkel/rad, skalar
//xR,yR: Koordinaten der gedrehten Punkte, Typ wie x
//Meyberg-Vachenauer Bd. 1, S. 318, dort D_3
xR=cos(w)*x - sin(w)*y
yR=sin(w)*x + cos(w)*y
endfunction

function d=!circanull(p,x,y)//Defektberechnung
//p=[a; b; alpha; mx; my]
x=x-p(4)
y=y-p(5)
[x,y]=!Dz2(x,y,-p(3))
d=p(2)^2*x^2 + p(1)^2*y^2 - p(1)^2*p(2)^2//Vektor
d=d*d'//Eigentlicher Defekt, d. h. Summe der Abweichungsquadrate
endfunction

//Messdaten
if 0 then
Dat=[..
-3 0.5
-9 1
-15 -1
-19 -5
-5 -2];
else
Dat=[..
-3.4 0.8
-8.3 0.9
-11.7 0.5
-14 -0.8
-18.1 -2.9
-18.4 -7.5
-14.6 -9.7
-9.5 -8.8
-5.9 -8.3
-3.4 -5.9
-0.6 -2.9];
end

//Dat=Dat+3*(rand(Dat)-0.5) //weitere zufällige Fehler

x=Dat(:,1)'; y=Dat(:,2)';//Messwerte
//p=[a b alpha mx my]//Bedeutung der Parameter
ab=(max(x)-min(x) + max(y)-min(y))/4;//Schätzwert für große und kleine Halbachse
coefs=regress(x,y);//Regressionsanalyse wegen Drehwinkel-Schätzwert
p0=[ab; ab; atan(coefs(2)); mean(x); mean(y)]//Startparameter (Anfangsschätzung)
[dOpt,pOpt]=optim(list(NDcost,!circanull,x,y),p0) // Ermittlung der optimalen Parameter für minimalen Defekt
phi=linspace(-%pi,%pi,50); //Winkel zur Generierung von dichtliegenden Ellipsenpunkten
x_N=pOpt(1)*cos(phi); y_N=pOpt(2)*sin(phi);//Ausgleichs-Ellipse in Normallage
[x_,y_]=!Dz2(x_N,y_N,pOpt(3));//Drehschritt zur Ausgleichs-Ellipse in Originallage
x_=x_+pOpt(4); y_=y_+pOpt(5);//Verschiebungsschritt zur Ausgleichs-Ellipse in Originallage
xdel();
plot( x,y,"ro", x_,y_,"g", x_N,y_N,"b" ),
legend("Messpunkte","Ausgleichsellipse","Ausgleichsellipse in Normallage");
title("a = " + string(pOpt(1)) + " b = " + string(pOpt(2)))
cax=gca();
cax.isoview="on";
p0_pOpt=[p0 pOpt]
//Ausgleichsverfahren für Ellipse

Rückfragen willkommen.

Neue Frage »

Antworten »

Berechnung einer Ellipse und deren Halbachsen aus 5 Punkten

Verwandte Themen