Verfasst am: 24.02.2013, 19:47
Titel: Schnelle Berechnung des Flächeninhaltes einer Punktewolke
Moin,
ich hoffe ihr könnt mir weiterhelfen. Ich möchte gerne die Fläche eines Objektes im Raum berechnen, dargestellt aus X-Y-Z Punkten (Siehe Bild a.jpg). Auf der Linken Seite ist so ein Objekt dargestellt. Die Punkte liegen als 3 Vektoren vor. Es ist quasi ein "Z-Map-Model", wo keine Hinterschnitte vorkommen
Was ganz angenehm ist: die einzelnen Punkte in der X bzw. Y Richtung sind äquidestant zueinander (siehe Bild rechte Seite). Bisher war mein Ansatz immer 4 Punkte zu nehmen, Vektoren zu bilden und das Kreuzprodukt zu berechnen - aber dabei verzettele ich mich in den Schleifen
Die Funktion die sowas berechnet habe ich leider noch nicht gefunden. Hättet ihr eine Lösung, die möglichst schnell durchgeführt wird?
Verfasst am: 26.02.2013, 13:04
Titel: Re: Schnelle Berechnung des Flächeninhaltes einer Punktewol
Hallo dotwinX,
Wie ist denn die Fläche genau definiert? Am Rand gibt es ja mehrere Methoden um die Punkte zu verbinden. Eine exakte Bestimmung der Umrandung ist aber notwendig, um die Fläche zu erhalten.
Jup, es geht wirklich um die Fläche, nicht um das Volumen.
Die Fläche hat keine genaue Definition. Also eine konvexe Hülle kann ich dafür leider nicht nehmen, da sie von "oben betrachtet" auch wie ein L z.B. aussehen kann. Hinterschneidungen sind aber nicht möglich.
Ist das mit dem Kreuzprodukt denn kein guter Ansatz?
Dann gibt es eine einfache Strategie: Eine Fläche ohne genaue Definition hat keinen genauen Inhalt. Es könnte schlicht jeder beliebeige Wert sein.
Du benötist also als aller erstes eine Definition, wie denn aus den Punkten eine Fläche konstruiert werden darf und kann. Eine lokale konvexe Hülle ist denkbar, bei der nur eine gewisse Umgebung mitgezählt wird.
Wenn aber die Fläche von oben wie ein L aussehen kann, wie sieht es denn dann am Knick aus? Welche Dicke hat das Objekt oder anders gefragt: hat das Objekt eine endliche Dicke, so dass der Ranb auch eine Fläche hat?
Zitat:
4 Punkte zu nehmen, Vektoren zu bilden und das Kreuzprodukt zu berechnen
Was heißt das genau? Welche 4 Punkte nimmst Du, wie bildest Du daraus Vektoren und wieso erwartest Du, dass das Kreuzprodukt die Fläche ergibt? Bei einem Quadrat könnte das die Fläche 0 ergeben, wenn Du benachbarte Ecken verwendest, bei gegenüberliegenden Ecken ist das schon anders.
Für 4 beliebige Punkte im Raum ist aber ein Flächeninhalt nicht einmal definiert, solange sie nicht in einer Ebene liegen. Wenn Du mit den 4 Punkten aber in 2 Dreiecke definierst, hättest Du wieder eine (also zwei) Flächen. Aber dies ist wieder nicht eindeutig definiert, da es ja zwei Dreiecksaufteilungen gibt.
Die Frage ist also wirklich nicht trivial. Gerade die Fehlende Definition des Randes ist tückisch.
ich glaub bzgl. der Definition der Fläche hatte ich dich falsch verstanden
Mal schaun ob ich es klarer beschreiben kann.
Also die Fläche ist unendlich dünn, ergo gibt es keine Fläche am Rand, welche orthogonal zur eigentlichen Fläche ist. Bei der Berechnung sind numerische Ungenauigkeiten von +-5% kein Problem. Daher auch mein flapsiger Hinweis, dass es keine "Definition" gäbe.
Wie ich mir die Fläche vorstelle habe ich mal an einem Bild veranschaulicht. Blau sind die Punkte und rot umrandet die Fläche die ich berechnen möchte. Stellt euch nur jetzt vor, dass die Punkte in X-Y Richtung auch unterschiedliche Z-Werte haben können.
Mit den 4 Punkten hast du natürlich Recht. Aber da ich ja recht viele Punkte habe und die nebeneinanderliegenden Punkte keine großen Sprünge in Z-Richtung machen, wäre bspw. eine Approximation durch 2 Dreiecke zulässig.
Meine Frage geht eher in die Richtung wie ich sowas programmiertechnisch klug umsetzen kann.
ich glaube ich habe eine Methodik gefunden die wirklich gut für meinen Anwendungsfall ist und näherungsweise korrekt. Ich nutze einfach die Viereckformel!
Wenn das alles fertig ist, dann kann ich das gerne posten. Nur leider hat mein Code im 1. Schritt bisher ein großen Nachteil an einer Stelle: Die Geschwindigkeit...
Ich muss aus der Punktewolke aus 3 separaten Vektoren (pktwolke) ein Matrix generieren (mesh). Und meine bisherige Lösung ist sehr langsam:
Code:
%Entferne alle doppelten Einträge
x = unique(pktwolke_X_werte);
y = unique(pktwolke_Y_werte);
[X,Y] = meshgrid(x,y);
%Fülle die Z-Matrix durch zeilen/spaltenweisen Vergleich
Z=zeros(size(X)); %Pre-Allo
Ist "pktwolke_Z_werte" ein Skalar? Andernfalls würde "Z(index_x, index_y)=pktwolke_Z_werte;" aber einen Fehler bewirken. Aber wenn dies kein Skalar ist, hat die Schleife nur eine Iteration. Deshalb verstehe ich nicht, was er Code bewirkt und rate mal ins Blaue hinein:
Code:
aX = X(1, :);
aY = Y(:, 1);
for k = 1:length(pktwolke_Z_werte)
index_x = find(aX == pktwolke_X_werte(k));
index_y = find(aY == pktwolke_Y_werte(k));
Z(index_x, index_y) = pktwolke_Z_werte(k);
end toc
Bitte erkläre "langsam" immer in absoluten Einheiten. Es ist wichtig, ob Du um Microsekunden oder Minuten kämpfst.
INTERSECT ist ineffizient, wenn nach einem Skalar gesucht wird. Kannst Du eventuell gleich nach "intersect(aX, pktwolke_Y_werte)" suchen?
Vielen Dank für deine Antwort. pktwolke_Z_werte ist ein recht langer Vektor (mehr als tausend Einträge). Dein Vorschlag ist schon mal sehr gut!
Vorher:
Elapsed time is 4.323416 seconds.
Jetzt:
Elapsed time is 0.343593 seconds.
Und es geht leider wirklich um Millisekunden, nicht um Minuten :/
Daher ist das schon ein sehr guter Anfang, aber noch ein bisschen schneller wäre natürlich besser.
Ich bin gerade am überlegen, ob dies auch ohne Schleife möglich ist die Werte zu finden.
Durch unique gehen ja schon Infos verloren...
Wenn ich die X-Y pktwolke als 2xn Matrix aufziehe und dann die X und Y-Einträge direkt in den Matrizen X und Y suche?...
Sowas wie ...
Code:
clear('i'); % It should be the imaginary unit
v1 = X + i*Y;
v2 = pktwolke_X_werte + i*pktwolke_Y_werte;
[~, indices] = ismember(v1, v2);
Ich habe immer noch nicht verstanden, was das Programm überhaupt bewirken soll. Wenn Du einen relevanten lauffähigen Code postest, sind Verbesserungsvorschläge möglich.
ursprünglich war meine Frage wie ich einen Flächeninhalt schnell bererechne. Dies hab ich mit der Vierecksformel gelöst. Ich wollte bzgl. meiner 2. Frage keinen neuen Thread aufmachen.
Bei der 2. Frage geht es darum, wie ich einen Teil des Codes verschnellere. Hierbei hab ich ursprünglich 3 Vektoren, die X, Y und Z Einträge von Punkten im Raum darstellen (pktwolke). Diese möchte ich in einem Array anordnen für die Berechnung mit der Viereckformel. Oder anders ausgedrückt: Statt plot3 würde ich gerne mesh/surf verwenden.
Nur ist aufgrund der Datenmengen (Diese Berechnung wird sehr häufig ausgeführt) diese Lösung noch nicht schnell genug.
In meinem Post davor hatte ich nur etwas gebrainstormed - bin aber zu keinem besseren Ergebnis gekommen :/
Und da war meine Kernfrage, ob der Schritt mit unique überhaupt notwendig ist um aus plot3 --> mesh zu machen.
Hallo, ich hab jetzt deinen Code noch nicht durchgesehen (bin bis zum Bild gekommen), leider schreibst du auch nix über die Äquidistanz und ob die Fläche wirklich planar liegt (also nicht im Raum gekrümmt ist) - aber wie wäre es, zu jedem Punkt die im Rasterradius befindlichen Punkte (oder mit Skalierung auch noch größer als Rasterradius) zu zählen und diese Summe als Summe kleiner Quadratflächen zu sehen? pi*daumen käme das hin.
[edit] Ahhh, ich hab das erste Bild "a" übersehen und nur das untitled angeguckt...
Also hast du eine Rasterung in -sagen wir mal- x&y.
Je nach dem wie deine Daten vorliegen klingt Kreuzprodukt doch nach der besten Lösung...
[edit2] okay, in dem Code geht es eher um die Umformung von Koordinatenlisten zu X,Y,Z. Warum denkst du zwingend, dass es so rum leichter geht?
Du könntest bspw. erstmal eine in Spalten ("nach hinten") offene Matrix bauen und die zu einem jeden Punkt nächst gelegenen Punkte aus der Koordinatenliste suchen (und deren Indizes in diese neue Matrix schreiben). Dann gehst du alle nächsten Punkte ab und berechnest die dadruch definierten Dreiecke. Oder?
_________________
Ich hasse es wenn die Leute Fragen stellen, man dann versucht sich Mühe zu geben, und diejenigen ihren Thread nie wieder besuchen...
vielen Dank für deine Antwort. Ich weiß nicht ob ich alles verstanden habe was du geschrieben hast
Aber so wie ich es verstanden habe, müsste ich somit jede Punktewolke einzelnd durchsuchen und jeder Berechnungsschritt eines Rasterpunktes wäre für sich alleine. Oder verpeile ich jetzt was? Mein großes Problem ist ja, dass die Berechnung schnell ablaufen muss. Daher würde ich eine einzelne Abfrage der Rasterpunkte gern vermeiden.
Ich hatte auch länger den Thread nicht besucht, weil ich eine Lösung gefunden habe die recht fix ist (könnte aber noch besser sein ). Ich möchte die erstmal kurz erläutern und zur Diskussion stellen.
Also vorweg: Ich nutze jetzt zur Berechnung der Fläche keine Dreiecke mehr, sondern Vierecke. Das ist der Einfachkeit geschuldet, weil ich bei Dreiecken ganz häufig Fallunterscheidungen machen muss. Wenn du untitled anschaust, dann kann solch ein Dreieck mal auf der einen, mal auf der anderne Seite, mal auf beiden Seiten von einem Rasterpunkt existieren.
Aus dem Grund hab ich einfach angenommen, dass wenn ich einen Raster mit 4 Punkten betrachte, die Diagonalen eines solchen Vierecks sich einfach schneiden. Dies ist eine Vereinfachung, aber aufgrund der vorliegenden Daten ist das eine gute Näherung mit der ich leben kann. Mit den Diagonalvektoren e und f (vgl. http://de.wikipedia.org/wiki/Viereck unter "Formeln") kann ich im dreidimensonalen Fall die Fläche berechnen.
Dies ist für jeden Rasterpunkt möglich. Wenn nun die Oberflächenpunkte nicht als Punktewolke (3xn Vektor mit X,Y,Z für Plot3) vorliegt, sondern als Matrix (Mx,My,Mz für Surf/Mesh usw), dann ist die Gleichung bei Wiki auch für Matrizen verwendbar.
Also 1. Schritt: Wandle Punktewolke in Matrizen um (Das dauert noch am längsten im Code).
Code:
%Schritt 1: Erstelle ein Mesh der realen Einträge mittels unique
x = unique(Pktwolke(:,1));
y = unique(Pktwolke(:,2));
[X,Y] = meshgrid(x,y);
%Schritt 2: Fülle die Z-Matrix durch zeilen/spaltenweisen Vergleich
Z=zeros(size(X)); %Pre-Allo
%http://www.gomatlab.de/viewtopic,p,110463.html#110463
aX = X(1, :); %X(1,:) reicht, weil die Spalteneinträge alle gleich sind;
aY = Y(:, 1); %Y(:,1) reicht, weil die Zeileneinträge alle gleich sind;
%Schneide die zu großgeratene Matrix Z (Ist durch die Schleife oben immer eine Square-Matrix, dies ist aber nicht nötig)
Z(~any(Z,2),:)=[];
Z(:,~any(Z))=[]; Z=Z';
Schritt 2: Wende die Gleichung mit e und f an und zwar mit folgendem "Trick". Statt den Punkten A, B, C, D (vgl. http://de.wikipedia.org/wiki/Viereck) hab ich die Matrix der Oberfläche Z. Diese ist um je einen Eintrag "verschoben" und zwar so:
Nun berechne ich die Diagonallängen e und f (Mit der Berücksichtigung, dass dies Matrizen A,B,C,D sind), das gleiche mit dem Skalarprodukt (auch das ist dann eine Matrix). Final:
vielen Dank für deine Antwort. Ich weiß nicht ob ich alles verstanden habe was du geschrieben hast
Was gibts denn da nicht zu verstehen?
Sorry, war da echt zu doof deinen ersten Post zu verarbeiten.
dotwinX hat Folgendes geschrieben:
Aber so wie ich es verstanden habe, müsste ich somit jede Punktewolke einzelnd durchsuchen und jeder Berechnungsschritt eines Rasterpunktes wäre für sich alleine. Oder verpeile ich jetzt was?
Was ist für dich jede Punktewolke? Die Nächste-Nachbar-Tabelle muss für jeden Punkt einen Eintrag haben, ja. Sofern die Punkte nahe genug gerastert wurden, kannst du die auch ohne Radius einspeichern, das zu rechnen würde ja nochmal Zeit kosten...
zu deiner Lösung:
* Passen denn die Punkte halbwegs ins Mx,My,Mz?
* hab grad kein Matlab hier, aber meshgrid sollte auch mit X, Y UND Z funktionieren...
* irgendwie verstehe ich nicht, warum du die Meshgrid-Variablen wieder zu Zeilenvektoren machst - das läuft ja dann wieder entgegen der Vektorisierung - Wie genau berechnest du dann e und f?
_________________
Ich hasse es wenn die Leute Fragen stellen, man dann versucht sich Mühe zu geben, und diejenigen ihren Thread nie wieder besuchen...
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS