Mathematik online lernen im Mathe-Forum. Nachhilfe online
Startseite » Forum » Flächen zweiter Ordnung für Punktmenge finden

Flächen zweiter Ordnung für Punktmenge finden

Universität / Fachhochschule

Körper

Polynome

Tags: Fläche zweiter Ordnung, Punktmenge

 
Antworten Neue Frage stellen Im Forum suchen
Neue Frage
alissa19

alissa19 aktiv_icon

11:58 Uhr, 25.08.2015

Antworten
Hallo zusammen,

ich bin auf der Suche nach einer geeigneten Fläche 2. Ordnung, die ich durch meine Punktwolke legen kann. Was ich bisher so gelesen habe, ist, dass ich als erstes 10 Parameter bestimmen muss, um auf die Form ax^2 + by^2 + cz^2 +dxy+ exz + fyz + gx + hy + iz +j=0 zu kommen

Da ich zehn Parameter habe, bräuchte ich ja eigentlich vier Punkte à 3 Koordinaten, um das Gleichungssystem zu lösen, korrekt?

Gibt es dafür einen geschickteren Wert, als seitenlange Terme umzuformen, um auf die einzelnen Parameter zu kommen? Damit habe ich jetzt zwar begonnen, nur ist das ganze so dermaßen groß und unübersichtlich, dass mir da zwangsweise Fehler passieren müssen.
Als beispielhafte Punkte habe ich

3.858425.77679.1883
3.897526.11787.8360
3.943119.91685.4890
3.893925.97067.7307
3.366229.377813.9810

Für eure Antworten wäre ich sehr dankbar.
Liebe Grüße


Für alle, die mir helfen möchten (automatisch von OnlineMathe generiert):
"Ich möchte die Lösung in Zusammenarbeit mit anderen erstellen."
Online-Nachhilfe in Mathematik
Antwort
Roman-22

Roman-22

12:23 Uhr, 25.08.2015

Antworten
> Da ich zehn Parameter habe, bräuchte ich ja eigentlich vier Punkte à 3 Koordinaten, um das Gleichungssystem zu lösen, korrekt?
Du brauchst doch wohl eher 10 Punkte um auf deine 10 Gleichungen für die 10 Variablen zu kommen, oder?
Dann hättest du dann eine Quadrik im 3 bestimmt, welche genau diese 10 Punkte enthält.
Und ja, "zu Fuß" ist das mühsam und abschreckend. Da sollte man sich besser um maschinelle Hilfe umsehen.

Aber ich vermute mal, dass deine Punktwolke aus hunderten oder tausenden von Punkten besteht. Da ist es nicht sinnvoll, zehn davon auszuwählen und die anderen zu ignorieren.

Schätze also, dass ein best fit gesucht sein wird. Also eine Quadrik, bei der zB die Summe der Quadrate der Abstände der Punktes von ihrer Oberfläche minimal ist. Und auch diese Regression, diese Minimierung der Residuen, wird man wohl eher den numerischen Algorithmen eines Programms anvertrauen.

R

alissa19

alissa19 aktiv_icon

12:28 Uhr, 25.08.2015

Antworten
Ähm... ja, zehn Punkte...

Du hast Recht, meine Punktwolke besteht aus knapp 1 Mio Punkten, ich wollte jedoch so vorgehen, dass ich mir für zehn zufällige Punkte die Quadrik bestimme und das ganze dann maschinell z.B. tausend Mal durchlaufen lasse - um so erst einmal zu schauen, ob sich eine Tendenz abzeichnet (zB. kriege ich dann 900x das Ergebnis, dass ein elliptisches Paraboloid die richtige Form sei)

An einem best fit Algorithmus habe ich mich versucht, nur den Lösungsweg verstehe ich wirklich so gar nicht - ich hatte gehofft, dass dieser Weg etwas einfacher sei.

Antwort
Roman-22

Roman-22

12:42 Uhr, 25.08.2015

Antworten
Wie sind denn die Daten verteilt? Könnte es eine polynomiale Regression zweiten Grades vielleicht auch tun? Da wäre dann die Einschränkung, dass es zu jedem Paar (x;y) nur einen eindeutigen z-Wert gibt.

Ansonsten wüsste ich jetzt ad hoc kein Programm, mit dem das out of the box funktionieren könnte. Am ehesten kann man sich da sicher in Matlab, Mathcad oder womit auch immer was zusammenbasteln. Also Formel für den Abstand eines Punktes von der (von 10 Parametern abhängigen Quadrik) aufstellen, Funktion für die Summe der Quadrat der Abstände aller Punkte definieren und diese dann mit den, dem Programm eigenen numerischen Algorithmen, minimieren lassen.

Dein Ansatz, aus den Millionen Punkten gerade einmal 1000 oder 10000 Mal ein Zehnerpack rauszugreifen find ich weniger zielführend. Vor allem, was nützt dir die Info, dass es möglicherweise ein elliptisches Paraboloid ist? Welches von den (von dir beispielhaft angegebenen) 900 nimmst du dann denn? Die grundlegende Form kannst du doch vielleicht auch schon selbst optisch erkennen, wenn du den entsprechenden Scatterplot ansiehst. Bei 1 Million Punkten sollte man da schon eine Tendenz auch optisch erkennen sollen.

R


alissa19

alissa19 aktiv_icon

12:50 Uhr, 25.08.2015

Antworten
Meine Punktwolke besteht aus x,y,z Koordinaten einer Kirchenwand, die ich selber vermessen habe. Ich bin auf der Suche nach einem bestmöglichen geometrischen Prinzip, um diese zu beschreiben. Eine Ebene reicht aufgrund der Krümmung hier leider nicht aus.

Mit der Erkenntnis, dass es sich vermutlich um ein elliptisches Paraboloid handelt, wollte ich dann im nächsten Schritt einen RANSAC Algorithmus über die Punktwolke laufen lassen. Der berechnet dann x-mal für zufällig ausgewählte Punkte eben ein Paraboloid und guckt über die Abstände aller Punkte zu dem jeweilig gefundenen Paraboloid, welches das insgesamt beste Ergebnis ist und wählt dieses aus. So grob wollte ich ja eben dieses RANSAC Prinzip jetzt auch erst einmal zur Geometriefindung nutzen.
Antwort
Roman-22

Roman-22

13:48 Uhr, 25.08.2015

Antworten
Wenn ich das richtig in Erinnerung habe (und da bewege ich mich durchaus auf dünnem Eis) liefert RANSAC die Modellparameter jenes Modells (aus bei dir genau 10 Punkten), bei dem möglichst viele der restlichen Punkte innerhalb eines zu wählenden Toleranzanstands zur Modellfläche liegen. Üblicherweise wird daher nach RANSAC ja doch eine Regression durchgeführt (Methode der kleinsten Summe der Abstandsquadrate), die die Parameter dann idR wesentlich verbessert.

Für eine Ersteinschätzung des zu wählenden Modells sollte doch vielleicht die optische Begutachtung der Kirchenwand reichen und wenn die nicht gerade Wellen schlägt könnte ich mir gut vorstellen, dass eine Quadrik reicht.

Über der yz-Ebene wird die Kirchenwand ja vermutlich einer eindeutigen Funktion entsprechen. Vielleicht reicht daher schon der Ansatz x=ay2+bz2+cyz+dy+ez+f. Immerhin nur 6 Parameter, was die Rechenzeit stark verkürzen sollte.

R

Ohne das Tool näher angesehen zu haben: Könnte eventuell die freie formfittingtoolbox der kostenlosen Sodtware JAG3D etwas für dich sein? derletztekick.com/software/formanalyse
Siehe auch diegeodaeten.de/formanalyse.html

Frage beantwortet
alissa19

alissa19 aktiv_icon

19:31 Uhr, 28.08.2015

Antworten
Hallo roman,

vielen Dank für deine Mühen bisher. Ich habe mich jetzt doch wie von dir vorgeschlagen für einen best-fit Algorithmus entschieden, um mein Problem zu lösen.

alissa19

alissa19 aktiv_icon

10:25 Uhr, 02.09.2015

Antworten
Hallo,

ich habe doch noch eine Rückfrage zum Thema:

Mit dem o. g. best fit Algorithmus habe ich mir die zehn Quadrikkoeffizienten
ax^2 + by^2 + cz^2 +dxy+ exz + fyz + gx + hy + iz +j=0
bestimmt.

Es handelt sich bei meiner Quadrik um ein Ellipsoid - wie komme ich mit diesen zehn Parametern denn nun auf die Form
x2a2+y2b2+z2c2=1?

Viele Grüße

Antwort
Roman-22

Roman-22

13:37 Uhr, 02.09.2015

Antworten
Das zugehörige Stichwort ist "Hauptachsentransformation".

zB web.student.tuwien.ac.at/~e0325258/studium/linalg.pdf
de.wikipedia.org/wiki/Hauptachsentransformation
http//www.mathebibel.de/hauptachsentransformation
....

R
alissa19

alissa19 aktiv_icon

12:05 Uhr, 03.09.2015

Antworten
Danke für die Links, ich muss gestehen, dass ich die ganze Thematik nie in der Vorlesung/Schule hatte und mich dementsprechend einlesen muss.

Ich habe mich an diesem Beispiel mal versucht:
http//www.mathematik.uni-mainz.de/Members/froehli/skripte/ss2010/teil2kap12.pdf
(S. 100ff), dazu habe ich zwei Fragen:

1. Ich muss das ganze zwingend mit Matlab programmieren, die Eigenwerte der gegebenen Matrix A bekomme ich jedoch je nach Befehl mit eig(A) =-4,9,36 bzw. svd(A) =36,9,4 heraus. Bei svd sogar mit einem geänderten Vorzeichen. Laut Beispiel sind die Eigenwerte jedoch 9,36,-4- also eine andere Reihenfolge

2. Die Hauptachsentrafo der Fläche im Bsp. lautet (S. 102)9u2+36v2-4w2-46=0
Woher kommt die 36? Wird die 1:1 aus dem Beispiel (S.101) übernommen?
Antwort
Roman-22

Roman-22

13:04 Uhr, 03.09.2015

Antworten
ad 1)+4 ist definitiv kein Eigenwert der Matrix A,-4 ist schon richtig und wird dir ja auch von Matlab geliefert. Der Vektor, den svd() u.a. liefert besteht aus den singulären Werten von A, die positiven Quadratwurzeln der Eigenwerte von AT¯A. Das sind nicht die Eigenwerte von A, sondern nur deren Beträge.
Die Reihenfolge der Eigenwerte ist unerheblich und es ist nicht unvernünftig, dass Matlab sie der Größe nach ordnet. Eine Änderung der Reihen bei der Hauptachsentransformation ändert nichts Wesentliches. Die Drehmatrix ändert sich und die Achsen. Bei deinem Ellipsoid ist dann halt einmal die u-Achse jene mit der größten Ausdehnung, bei einer anderen Reihenfolge der Eigenwerte ist es dann vielleicht die w-Achse.

ad 2) Ja. Beachte das "Kochrezept" auf Seite 100 und da insbesondere den Summanden b.

R

alissa19

alissa19 aktiv_icon

13:09 Uhr, 04.09.2015

Antworten
Hi,

das Beispiel habe ich jetzt nachvollziehen können. Ich habe jedoch immer noch arge Probleme mit der Transferleistung auf meine Problematik. Da ich für meine Daten die ermittelten Ellipsoidparameter nicht auf den ersten Blick plausibilisieren kann, habe ich mir selber ein Ellipsoid mit frei gewählten a,b,c Parameter erstellt. Für diese Punktwolke versuche ich das gerade nachzurechnen, da ich dort ja die Sollwerte habe.

Ich erhalte zuerst einmal die 10 Koeffizienten [0.3618,-0.9278,-0.0907,0.0000,-0.0000,0.0000,0.0989,-0.0000,0.0000,3.3882]

Wenn ich jetzt die Spektralzerlegung der Matrix B=[a1,a42,a52;a42,a2,a62;a52,a62,a3] durchführe, erhalte ich nun mit Matlab die Eigenwerte [0.9278,-0.0907,0.3618] und die Eigenvektoren [0.0000,0.0000,-1.0000;
-1.0000,0.0000,-0.0000;0.0000,1.0000,0.0000]

Nach Kochrezept rechne ich ja nun D= EV' B EV, was eine 3x3 Matrix ergibt, die auf den Nebendiag =0 und auf der Diag wiederum die Eigenwerte hat.

Und im letzten Schritt rechne ich diag(D)/a10, was mich zum Ergebnis von [-0.2738,-0.0268,0.1068] führt, was überhaupt nicht mit meiner Solllösung übereinstimmt. Wo liegt denn jetzt mein Fehler? :(

Hilflose Grüße


Antwort
Roman-22

Roman-22

14:25 Uhr, 04.09.2015

Antworten
Hast du beachtet, dass eine Hauptachsentrafo aus zwei Teilen besteht? Einmal eine Translation, so dass der Mittelpunkt der Quadrik in den Ursprung gelangt und dann erst die Drehung. Beachte in dem von dir vorhin zitierten Beispiel, dass die Koeffizienten von x,y und z bereits Null sind.

Wenn ich die Reihenfolge deiner Koeffizienten richtig interpretiere (so wie in deinem ersten Post angegeben), dann handelt es sich ja nur um eine Quadrik, die ein wenig in x-Richtung verschoben wurde. Es wurde keine Drehung angewandt (sonst wären die Koeeffizienten von xy, xz bzw. yz nicht alle Null).
Macht man diese Verschiebung rückgängig (u=x+0,098920,3618) und bringt die Gleichung in Normalform, erhält man -u23,051762+y21,905712+z20,164072=1. Ein Ellipsoid ist das allerdings nicht, denn da müssten die Koeffizienten der quadratischen Glieder allesamt positiv sein.

Wenn du dir, um das Prinzip besser verstehen zu lernen, selbst ein Beispiel ausdenkst - warum dann eines mit solch "krummen" Werten? Nimm doch ein Ellipsoid mit "schönen" Werten für a,b und c, unterwerfe es der einen oder anderen Drehung und dann schau, ob du wieder zurück findest.

R
alissa19

alissa19 aktiv_icon

11:12 Uhr, 05.09.2015

Antworten
"Wenn du dir, um das Prinzip besser verstehen zu lernen, selbst ein Beispiel ausdenkst - warum dann eines mit solch "krummen" Werten? Nimm doch ein Ellipsoid mit "schönen" Werten für a,b und c, unterwerfe es der einen oder anderen Drehung und dann schau, ob du wieder zurück findest."

Genau das habe ich ja eigentlich probiert - meine Sollparameter sind a=4, b=2.5 und c=8.

Ich habe jetzt das Folgende in Matlab programmiert:

clc; clear all
s_dach = [0.3618,-0.9278,-0.0907,0.0000,-0.0000,0.0000,0.0989,-0.0000,0.0000,3.3882];
B = [s_dach(1) s_dach(4)/sqrt(2) s_dach(5)/sqrt(2); s_dach(4)/sqrt(2) s_dach(2) s_dach(6)/sqrt(2); s_dach(5)/sqrt(2) s_dach(6)/sqrt(2) s_dach(3)];
b = [s_dach(7) s_dach(8) s_dach(9)]';
b0 = s_dach(10);

%% Drehung
[EV,EW] = eig(B);

D=EV'*B*EV;
d = EV'*b;

% Terme d4-d6 sind jetzt eliminiert
d1 = D(1,1);
d2 = D(2,2);
d3 = D(3,3);
d7 = d(1,1);
d8 = d(2,1);
d9 = d(3,1);
d10 = s_dach(10);

%% Verschiebung

px2 = (((d7/d1)/2)^2)*d1;
py2 = (((d8/d2)/2)^2)*d2;
pz2 = (((d9/d3)/2)^2)*d3;
d = -px2 - py2 - pz2 + d10;

a = d1;
b = d2;
c = d3;

a/d
b/d
c/d

Kriege dort aber die Parameter a = -0,277, b = -0,0268, c = 0,1070 heraus.

Antwort
Roman-22

Roman-22

13:39 Uhr, 05.09.2015

Antworten
> Genau das habe ich ja eigentlich probiert - meine Sollparameter sind a=4,b=2.5 und c=8.
Und welche Dreh- und Schiebeoperationen hast du nun auf die Quadrik x2a2+y2b2+z2c2=1 angewandt - also wie bist du auf deinen Vektor s_dach gekommen?
Denn wenn ich die Werte darin so interpretiere, wie du das in deinem ersten Post angegeben hast (a bis j), stellt das, wie schon geschrieben, kein Ellipsoid mehr dar.

R

alissa19

alissa19 aktiv_icon

14:23 Uhr, 05.09.2015

Antworten
1) Konstruiert habe ich mir die Punktwolke so, dass ich Y und Z Werte von -100:100 laufen lassen habe und für dieses 2D Gitter mir dann den X-Wert berechnet habe mit meinen Sollparametern.

a = 4;
b = 2.5;
c = 8;
[Y,Z] = meshgrid(-100:100,-100:100);
X = sqrt(abs(-a^2 * (Y.^2/b^2 + Z.^2/c^2 - 1)));

Das Ergebnis habe ich auch mal als Bild aus CloudCompare angehängt.

2) s_dach bekomme ich aus meinen Algorithmus für die erstellte Punktwolke heraus, die Berechnung läuft mit Matlab so ab: (Der ALgorithmus basiert auf der Dissertation von Drixler, 1993: Analyse der Form und Lage von Objekten im Raum)

P = load('Ellipsoidwolke.txt');
% Reservierung von Speicherbedarf
A1=zeros(length(P),6);
A2=zeros(length(P),4);
U=sparse(zeros(length(P),length(P)*3));

% Aufstellen der Designmatrizen A1 und A2
for i=1:length(P)
A1(i,1:6) = [P(i,1)^2 P(i,2)^2 P(i,3)^2 sqrt(2)*P(i,1)*P(i,2) sqrt(2)*P(i,1)*P(i,3) sqrt(2)*P(i,2)*P(i,3)];
A2(i,1:4) = [P(i,1) P(i,2) P(i,3) 1];
end

% Vereinfachtes stochastisches Modell
P_ww = sparse(eye(length(P)));

% Überführung des Gesamtmodells in ein EWP
H=P_ww-P_ww*A2*inv(A2'*P_ww*A2)*A2'*P_ww;
Q_HH=A1'*H*A1;

% Eigenwertzerlegung
[M_QHH lambda_QHH]=svd(Q_HH);
lambda_QHH_help=[lambda_QHH(1,1) lambda_QHH(2,2) lambda_QHH(3,3) lambda_QHH(4,4) lambda_QHH(5,5) lambda_QHH(6,6)];
such=min(lambda_QHH_help);
[row col]=find(lambda_QHH_help==such);
row=col;

s_dach_1=[M_QHH(1,col); M_QHH(2,col); M_QHH(3,col); M_QHH(4,col); M_QHH(5,col); M_QHH(6,col)];
s_dach_2=-1*inv(A2'*P_ww*A2)*A2'*P_ww*A1*s_dach_1;

% Verzerte Genauigkeiten der ausgeglichenen Parameter
Q_sd1_sd1=inv(Q_HH)-inv(Q_HH)*s_dach_1*inv(s_dach_1'*inv(Q_HH)*s_dach_1)*s_dach_1'*inv(Q_HH);
Q_sd2_sd2=inv(A2'*P_ww*A2)*A2'*P_ww*A1*Q_sd1_sd1*A1'*P_ww*A2*inv(A2'*P_ww*A2);

%% Berechnung der Kovarianzmatrix der ausgeglichenen Parameter der Quadrik

% Reservierung von Speicherbedarf

for i=1:length(P)
U(i,i*3-2)=2*s_dach_1(1)*P(i,1)+sqrt(2)*s_dach_1(4)*P(i,2)+sqrt(2)*s_dach_1(5)*P(i,3)+s_dach_2(1);
U(i,i*3-1)=2*s_dach_1(2)*P(i,2)+sqrt(2)*s_dach_1(4)*P(i,1)+sqrt(2)*s_dach_1(6)*P(i,3)+s_dach_2(2);
U(i,i*3) =2*s_dach_1(3)*P(i,3)+sqrt(2)*s_dach_1(5)*P(i,1)+sqrt(2)*s_dach_1(6)*P(i,2)+s_dach_2(3);
end

% Berichtigtes stochastisches Modell
clear P_ww
P_ww=inv(U*U');

% Überführung des Gesamtmodells in ein EWP
clear H Q_HH
H=P_ww-P_ww*A2*inv(A2'*P_ww*A2)*A2'*P_ww;
Q_HH=A1'*H*A1;

% Eigenwertzerlegung
clear M_QHH lambda_QHH lambda_QHH_help such row col
[M_QHH lambda_QHH]=svd(Q_HH);
lambda_QHH_help=[lambda_QHH(1,1) lambda_QHH(2,2) lambda_QHH(3,3) lambda_QHH(4,4) lambda_QHH(5,5) lambda_QHH(6,6)];
such=min(lambda_QHH_help);
[row col]=find(lambda_QHH_help==such);
row=col;

clear s_dach_1 s_dach_2
s_dach_1=[M_QHH(1,col); M_QHH(2,col); M_QHH(3,col); M_QHH(4,col); M_QHH(5,col); M_QHH(6,col)];
s_dach_2=-1*inv(A2'*P_ww*A2)*A2'*P_ww*A1*s_dach_1;

s_dach=[s_dach_1(1); s_dach_1(2);s_dach_1(3); s_dach_1(4); s_dach_1(5); s_dach_1(6); s_dach_2(1); s_dach_2(2); s_dach_2(3); s_dach_2(4)];


Wolke
Antwort
Roman-22

Roman-22

19:41 Uhr, 05.09.2015

Antworten
"
a=4;
b=2.5;
c=8;
[Y,Z]= meshgrid(-100:100,-100:100);
X= sqrt(abs(-a^2 (Y.2b2+Z.2c2-1)));
"

Auweia! Alles klar. Das konnte so ja nicht gut gehen.

Du hast ein Ellipsoid, bei dem sich die Koordinaten der Punkte maximal in den Bereichen -4x4,-2,5y2,5 und -8z8 liegen.
Und dann erzeust du einen yz-Raster jeweils von -100 bis 100. Du erzeugst also 40401 Punkte, von denen aber nur höchstens 85(=517) auf deinem Ellipsoid liegen. Wenn man es sich genauer ansieht, sind es sogar nur genau 65.

Der große Rest der Punkte müsste eigentlich eine komplexe x-Koordinate haben, aber Dank des von dir eingefügten Absolutbetrags werden reelle Werte berechnet, die aber mit deinem Ellipsoid nichts zu tun haben. Kurz - die von dir erzeugte Punktwolke liegt nicht auf einem Ellipsoid. Du siehst ja auch an deinem Bild, dass das eher an einen hyperbolischen Zylinder erinnert, aber sicher an kein Ellipsoid.
Zur Verdeutlichung habe ich die Werte für a,b und c einmal verzehnfacht und jetzt siehst du in Bild_1 dein Ellipsoid und eben auch die unerwünschte Fläche außerhalb.

Um eine sinnvolle Punktwolke zu erhalten müsstest du bei Beibehaltung deiner Werte für a,b und c dein Netz kleiner und wesentlich enger wählen und musst jene Punkte, für die der Term unter der Wurzel negativ wird, aus deiner Wolke entfernen.
Dass du mit deinem Ansatz nur das vordere Halbellipsoid bekommst ist vermutlich in Hinblick auf die Anwendung (Kirchenwand) Absicht.

Ich denke dass es einfacher wäre, wenn du deine Punktwolke mithilfe einer Parameterdarstellung deines Ellipsoids erzeugen würdest, denn da treten dann von Haus aus keine ungültigen Punkte auf ( Bild2).
Die übliche Darstellung ist (x(u,v)=asin(u)cos(v)y(u,v)=bsin(u)sin(v)z(u,v)=ccos(u)) mit u[0;π) und v[0;2π) oder, wenn du nur das vordere Halbellipsoid haben möchtest, denn eben v[-π2;π2).
Natürlich kannst du, da bei dir ja die x-Richtung ausgezeichnet ist, die Gleichungen für x und z auch entsprechend vertauschen.

Eine Anmerkung noch. Das Ellipsoid, so wie du es angeben wolltest, befindet sich doch schon in Hauptachsenlage. Also der Mittelpunkt liegt im Ursprung und die Achsen sind die Koorinatenachsen. Da macht doch eine Hauptachsentransformation überhaupt keinen Sinn, denn das ist ja bereits genau die Situation, die man bei einer allgemein herumliegenden Quadrik mit der Trafo erst herzustellen versucht.

R



Bild1
Bild2
Antwort
Roman-22

Roman-22

19:20 Uhr, 09.09.2015

Antworten
Na, so schnell das Interesse wieder verloren und grußlos verschieden?

Diese Frage wurde automatisch geschlossen, da der Fragesteller kein Interesse mehr an der Frage gezeigt hat.