Senkrechtester Pixel in Punktwolke

[Tobi (Gast)]

Hallo,

ich habe eine Punktwolke von einer Kamera. Nun guckt diese Kamera auf 
eine Wand und ich möchte bestimmen, welcher Pixel der senkrechteste zur 
Wand ist.

Ich denke der "beste" Ansatz wäre, eine ebene in eine gemittelte 
Punktwolke zu fitten und dann für jeden Pixel den Schnittwinkel mit 
dieser Ebene zu berechnen. Für jeden Pixel habe ich die (x y 
z)-Koordinaten (bezogen auf die Kamera) oder den radialen Abstand zur 
Verfügung.

Der Schnittwinkel sollte ja nach

$$\alpha = sin^{-1}(\frac{r \cdot n}{|r| \cdot |n|}) \ mit \ r = (x \ y \ z) \ und \ n = Normale \ der \ Ebene$$

zu berechnen sein.

Nun habe ich aber keine Ahnung, wie ich am besten eine Ebene in die 
Punktwolke fitten kann. Im Dreidimensionalen habe ich so etwas noch nie 
gemacht.
Das Ganze soll nachher in C# und später eventuell auch in MATLAB 
implementiert werden.

Hat jemand da einen Ansatz?

[Karl H. (kbuchegg)]

Tobi schrieb:
> Hallo,
>
> ich habe eine Punktwolke von einer Kamera. Nun guckt diese Kamera auf
> eine Wand und ich möchte bestimmen, welcher Pixel der senkrechteste zur
> Wand ist.

Und von der Wand weißt du nichts - richtig?

> Ich denke der "beste" Ansatz wäre, eine ebene in eine gemittelte
> Punktwolke zu fitten

Kann man machen.
Das Stichwort dazu lautet zb "Orthogonal Distance Regression Planes"

http://mathforum.org/library/drmath/view/63765.html

> und dann für jeden Pixel den Schnittwinkel mit
> dieser Ebene zu berechnen.

Zu kompliziert.
Was genau verstehst du unter 'dem Senkrechtesten'?
Jeder Punkt (sofern) er nicht gerade in der Ebene liegt, steht auf eine 
Ebene senkrecht.
Ich denke du meinst: du suchst dann den Punkt, der zu dieser Ebene am 
nächsten liegt, also den, der den kleinsten Abstand zur Ebene hat. Das 
ist eber einfach: Einfach den Punkt  in die Ebenengleichung einsetzen 
und heraus kommt der Normalabstand des Punktes von dieser Ebene.

Eine Ebene ist ja definiert als die Menge aller Punkte, die Abstand 0 
von dieser Ebene haben

    a*x + b*y + c*z + d = 0

a, b, c, d sind die Parameter der Ebenengleichung. Punkt x/y/z einsetzen 
und wenn 0 rauskommt, dann liegt der Punkt auf der Ebene. Kommt da nicht 
0 raus, dann ergibt sich so die Normaldistanz des Punktes von der Ebene, 
wobei dir das Vorzeichen auch noch verrät, auf welcher Seite der Ebene 
sich der Punkt befindet.

[Tobi (Gast)]

Danke für den Link, werde ich mir gleich direkt angucken.

Unter "am senkrechtesten" verstehe ich folgendes:
Hat man eine gerade Wand und richtet eine Kamera perfekt parallel zur 
Wand aus, sodass sie zur Wand guckt, dann steht der mittlere Pixel genau 
senkrecht auf der Wand. Alle anderen Pixel werden ja weiter nach außen 
emittiert. Daher treffen diese dann nicht mehr senkrecht auf die Wand.
Ist die Kamera jetzt nicht perfekt ausgerichtet, so ist es nicht mehr 
der mittlere Pixel der senkrecht zur Wand steht, sondern ein anderer.

Genau diesem Pixel wollte ich bestimmen. Aber das mit dem Winkel scheint 
doch etwas zu verkompliziert zu sein. Der kleinste Abstand sollte ja 
genau dieser Pixel sein.

Also folgendes vorgehen:
- n Bilder aufnehmen
- Mittelwert über alle n Aufnahme für alle 3 Koordinaten (x y z)
- Ebene bestimmen, welche am besten in die Punktwolke passt (hier muss 
ich noch gucken wie ich das am besten mache)
- Über alle Punkte iterieren und in die Ebenengleichung einsetzen
- Punkt mit kleinstem Abstand

$$min(a \cdot x_i + b \cdot y_i + c \cdot z_i + d)$$

 bestimmen

[zoggl (Gast)]

du kennst den nodalpunkt der kamera und hast deine ebene bestimmt.

jetzt projeziert du alle punkte (auch den nodalpunkt) entlang der 
Nodalpunktnormalen auf die ebene und suchst den punkt mit dem kleinsten 
abstand zum nodalpunkt.

ist unter umständen schneller, da du besser suchen kannst.

[Tobi (Gast)]

Ich habe mir das Verfahren hinter dem Link mal angeguckt. So ganz 
verstehe ich das noch nicht.

Also ich brauche den Punkt des Flächenschwerpunkts (x0 y0 z0). 
Allerdings habe ich keine Ahnung wie ich diese bestimme.

Anschließend könnte ich die Matrix M aufstellen mit
                 --                             --
                 | x1 - x0    y1 - y0    z1 - z0 |
                 | x2 - x0    y2 - y0    z2 - z0 |
            M =  |    .          .          .    |
                 |    .          .          .    |
                 |    .          .          .    |
                 | xn - x0    yn - y0    zn - z0 |
                 --                             --

Für M muss dann eine Singulärwertzerlegung gemacht werden. Davon habe 
ich keine Ahnung. Daher würde ich mir eine fertige Bibliothek suchen, 
die das kann. Spontan habe ich diese gefunden:
http://www.codeproject.com/Articles/51470/Advanced-Matrix-Library-in-C-NET
Wenn jemand eine andere empfehlen kann, dann wär ich für Tipps dankbar.

Anschließend wird A = M^T * M mit Hilfe der Singulärwertzerlegung 
bestimmt und ergibt sich zu A = V * S^2 * V^T.

S enthält die Singulärwerte (bzw. das Quadrat). Der Singulärvektor (in 
den Spalten von V enthalten), welcher zum kleinesten Singulärwert gehört 
ist dann der Normalvektor (n) der Ebene.

Somit hätte ich dann die Ebene in der Normalenform:

$$(\vec{r} - (x_0\ y_0\ z_0)^T) \cdot \vec{n} = 0$$

[Karl H. (kbuchegg)]

Tobi schrieb:
> Ich habe mir das Verfahren hinter dem Link mal angeguckt. So ganz
> verstehe ich das noch nicht.
>
> Also ich brauche den Punkt des Flächenschwerpunkts (x0 y0 z0).

Nicht ganz.
Centroid ist nicht der Flächenschwerpunkt.

Für x0 addierst du einfach alle X Werte und dividierst durch die Anzahl.
Für y0     -"-                  Y
Für z0     -"-                  Z

Also einfach immer nur der Durchschnitt aller Punkte.

[Tobi (Gast)]

Also einfach nur:

List<List<ScanPoint>> points; // (dynamische) Matrix von Scanpunkten

[...]

double x0 = 0, y0 = 0, z0 = 0;
int n = points.Count * points[0].Count;

foreach (var point in points)
{
    x0 += point.X;
    y0 += point.Y;
    z0 += point.Z;
}

x0 /= n;
y0 /= n;
z0 /= n;

OK, danke für die Hilfe. Dann sollte ich soweit die Ebene fitten können.

[Tobi (Gast)]

Mir fällt gerade auf, dass das komplette rumgerechne mit A = M^T * M 
doch eigentlich überflüssig sein müsste.
Es interessiert ja letztendlich nur der Singulärvektor der zum kleinsten 
Singulärwert gehört. Diesen bekomm ich doch schon direkt aus der 
Singulärwertzerlegung von M. Also müsste ich doch eigentlich nach der 
Singulärwertzerlegung von M die Spalte in S mit dem kleinsten Wert 
suchen, und dann die Zeile (mit dem bestimmten Spaltenindex) auslesen 
(Zeile weil man V^T hat und nicht V) und dann hab ich doch schon direkt 
(x0 y0 z0).

Oder übersehe ich da was?