Oberflächennormale aus (x,y,z) Punktewolke

[Franziska]

Hallo,

ich bin relativ neu bei Python und suche eine einfache Möglichkeit, die Oberflächennormale einer Fläche aus 100.000 Punkten zu bestimmen, ähnlich des Befehls surfnorm(X,Y,Z) bei Matlab.

Die Daten liegen als X,Y,Z-Tabelle vor, z.B.:

X Y Z
1 1 4
1 2 1
1 3 3
2 1 3
2 2 5
2 3 3
3 1 2
3 2 5
3 3 5
...

Ich habe bisher nur die Newell-Methode gefunden und die scheint mir etwas umständlich. Habt ihr Anregungen, wie das vielleicht einfacher geht?

[ThomasL]

Hast du mal hier geschaut?
https://stackoverflow.com/questions/350 ... inalg-lsts

[Franziska]

Hallo,
vielen Dank für den Link. Ich habe die dritte Method auf dieser Seite verwendet: https://stackoverflow.com/questions/206 ... 3#20700063
Meine Testebene ist um 45° verkippt. Dennoch bekomme ich als Ergebnis den Vektor (1,0,0), was ja einer 90°-Verkippung entsprechen würde. Was mache ich falsch?

Code: Alles auswählen

import matplotlib.pyplot as plt
import numpy as np

xs,ys,zs = np.loadtxt('dummyplane.dat', skiprows=1, unpack=True)

# plot raw data
plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xs, ys, zs, color='b')

# do fit
tmp_A = []
tmp_b = []
for i in range(len(xs)):
    tmp_A.append([xs[i], ys[i], 1])
    tmp_b.append(zs[i])
b = np.matrix(tmp_b).T
A = np.matrix(tmp_A)
fit = (A.T * A).I * A.T * b
errors = b - A * fit
residual = np.linalg.norm(errors)

print("solution:")
print("%f x + %f y + %f = z" % (fit[0], fit[1], fit[2]))

# plot plane
xlim = ax.get_xlim()
ylim = ax.get_ylim()
X,Y = np.meshgrid(np.arange(xlim[0], xlim[1]),
                  np.arange(ylim[0], ylim[1]))
Z = np.zeros(X.shape)
for r in range(X.shape[0]):
    for c in range(X.shape[1]):
        Z[r,c] = fit[0] * X[r,c] + fit[1] * Y[r,c] + fit[2]
ax.plot_wireframe(X,Y,Z, color='k')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

dummyplane:
X Y Z
1 1 1
1 2 1
1 3 1
2 1 2
2 2 2
2 3 2
3 1 3
3 2 3
3 3 3
4 1 4
4 2 4
4 3 4