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

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
Franziska
User
Beiträge: 2
Registriert: Sonntag 12. Januar 2020, 20:29

Sonntag 12. Januar 2020, 20:51

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?
Benutzeravatar
ThomasL
User
Beiträge: 937
Registriert: Montag 14. Mai 2018, 14:44
Wohnort: Kreis Unna NRW

Montag 13. Januar 2020, 07:11

Ich bin Pazifist und greife niemanden an, auch nicht mit Worten.
Für alle meine Code Beispiele gilt: "There is always a better way."
https://projecteuler.net/profile/Brotherluii.png
Franziska
User
Beiträge: 2
Registriert: Sonntag 12. Januar 2020, 20:29

Dienstag 14. Januar 2020, 22:37

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
Antworten