ich habe hier ein kleines Problem mit einem Minimierungsproblem.
Ich möchte aus einer Punktwolke die in etwa der Form eines Zylinders sich befindet die Achse herausfinden.
Die Modellfunktion lauten in TeX ausgedrückt:
\vec u \times (\vec r - \vec r_0) = \vec 0
wobei u der Richtungsvektor, r der variablenvektor und r_0 ein Vektor der einen Punkt auf der Geraden beschreibt.
Code: Alles auswählen
import scipy
import scipy.optimize
import numpy
def main():
points_file = open("cylinder_points.dat", "r")
points = [ scipy.array(eval(point)) for point in points_file ]
modell = lambda p, x: (x[2]-p[5])*p[1]-(x[1]-p[4])*p[2] \
+ (x[2]-p[5])*p[0]-(x[0]-p[3])*p[2] \
+ (x[1]-p[4])*p[0]-(x[0]-p[3])*p[1]
residuum = lambda p, x, y, z: modell(p, (x, y, z))
# die ersten 3 sind die parameter der richtungsvektors
p0 = scipy.array([1, 0, 1, 0, 0, 0])
x = [ point[0] for point in points ]
y = [ point[1] for point in points ]
z = [ point[2] for point in points ]
solution = scipy.optimize.leastsq(residuum, p0, args=(x, y, z))
print solution
if __name__ == "__main__":
main()
Code: Alles auswählen
0.995673605161946,0.,3.284017196927981
0.9023002750910056,0.02835592849599887,3.8388757923027965
0.9374358031998625,0.05897846162997725,0.39083200959969727
1.0143515912114693,0.09588445597047389,1.1818810540359217
0.9171049669815177,0.11585730044861368,3.262393063120224
1.067655995439107,0.16910009729675862,1.3904853543137952
1.076725162378293,0.20539630970910225,0.9104966324623105
1.02166002885339,0.2283680729661981,2.8233424612551303
1.0014651330122595,0.25713254258740914,0.530276650966358
0.8819471369235198,0.2562293294941034,2.399641520565491
0.9138115314586777,0.2969153652156805,1.348931238777783
0.9065138082031955,0.3263650530403482,3.2655100975906586
[...]
Das Ergebnis dieser Berechnung ist
Code: Alles auswählen
array([ 8.80062161e-01, -8.80062161e-01, 8.80062161e-01,
2.22183797e+02, 8.05312126e-09, 2.22183797e+02])
Die von mir konstruierten Punkte bzw der Zylinder hat die z-Achse als Richtungsvektor.
Kann jemand helfen?