![]() |
---|
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 Parameter bestimmen muss, um auf die Form ax^2 by^2 cz^2 exz fyz gx hy iz 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 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." |
![]() |
![]() |
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 Punkte um auf deine Gleichungen für die Variablen zu kommen, oder? Dann hättest du dann eine Quadrik im bestimmt, welche genau diese 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. |
![]() |
Ä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 . tausend Mal durchlaufen lasse - um so erst einmal zu schauen, ob sich eine Tendenz abzeichnet (zB. kriege ich dann 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. |
![]() |
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 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 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 oder 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) 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. |
![]() |
Meine Punktwolke besteht aus 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. |
![]() |
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 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 . Immerhin nur 6 Parameter, was die Rechenzeit stark verkürzen sollte. 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 |
![]() |
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. |
![]() |
Hallo, ich habe doch noch eine Rückfrage zum Thema: Mit dem . . best fit Algorithmus habe ich mir die zehn Quadrikkoeffizienten ax^2 by^2 cz^2 exz fyz gx hy iz bestimmt. Es handelt sich bei meiner Quadrik um ein Ellipsoid - wie komme ich mit diesen zehn Parametern denn nun auf die Form ? Viele Grüße |
![]() |
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 . |
![]() |
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 . 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) bzw. svd(A) heraus. Bei svd sogar mit einem geänderten Vorzeichen. Laut Beispiel sind die Eigenwerte jedoch also eine andere Reihenfolge 2. Die Hauptachsentrafo der Fläche im Bsp. lautet . Woher kommt die ? Wird die aus dem Beispiel übernommen? |
![]() |
ad ist definitiv kein Eigenwert der Matrix ist schon richtig und wird dir ja auch von Matlab geliefert. Der Vektor, den svd() . liefert besteht aus den singulären Werten von die positiven Quadratwurzeln der Eigenwerte von . Das sind nicht die Eigenwerte von 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 Ja. Beachte das "Kochrezept" auf Seite und da insbesondere den Summanden . |
![]() |
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 Parameter erstellt. Für diese Punktwolke versuche ich das gerade nachzurechnen, da ich dort ja die Sollwerte habe. Ich erhalte zuerst einmal die Koeffizienten Wenn ich jetzt die Spektralzerlegung der Matrix durchführe, erhalte ich nun mit Matlab die Eigenwerte und die Eigenvektoren Nach Kochrezept rechne ich ja nun EV' EV, was eine Matrix ergibt, die auf den Nebendiag und auf der Diag wiederum die Eigenwerte hat. Und im letzten Schritt rechne ich diag(D)/a10, was mich zum Ergebnis von führt, was überhaupt nicht mit meiner Solllösung übereinstimmt. Wo liegt denn jetzt mein Fehler? Hilflose Grüße |
![]() |
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 und 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 und bringt die Gleichung in Normalform, erhält man . 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 und unterwerfe es der einen oder anderen Drehung und dann schau, ob du wieder zurück findest. |
![]() |
"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. |
![]() |
Genau das habe ich ja eigentlich probiert - meine Sollparameter sind und . Und welche Dreh- und Schiebeoperationen hast du nun auf die Quadrik 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 bis stellt das, wie schon geschrieben, kein Ellipsoid mehr dar. |
![]() |
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)]; |
![]() |
" meshgrid(-100:100,-100:100); sqrt(abs(-a^2 " 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 und liegen. Und dann erzeust du einen yz-Raster jeweils von bis . Du erzeugst also Punkte, von denen aber nur höchstens auf deinem Ellipsoid liegen. Wenn man es sich genauer ansieht, sind es sogar nur genau . 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 und 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 und 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 mit und oder, wenn du nur das vordere Halbellipsoid haben möchtest, denn eben . Natürlich kannst du, da bei dir ja die x-Richtung ausgezeichnet ist, die Gleichungen für und 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. |
![]() |
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.
|