Minimierungsproblem mit scipy lösen

[zipdrive]

Hallo,

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()

Hier ein Auszug aus den Punkten, mit dem das Ganze berechnet werden soll.

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
[...]

(siehe auch: http://pastebin.mozilla.org/570655)

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])

, allerdings sollte für die ersten 3 etwas wie 0, 0, 1 herauskommen.

Die von mir konstruierten Punkte bzw der Zylinder hat die z-Achse als Richtungsvektor.

Kann jemand helfen?

[Darii]

Deine Modellfunktion ist falsch. Sie gibt nicht den Abstand eines Punktes zu einer Geraden an.

Und eval an der Stelle ist auch eine ganz schlechte Idee.

Code: Alles auswählen

from itertools import izip
import scipy.optimize
from scipy import cross, array
from numpy.linalg import norm
def main():
    x, y, z = scipy.loadtxt("cylinder_points.dat", delimiter=",", unpack=True)
    
    def residuum(params, x, y, z):
        u = array(params[:3])
        p = array(params[3:])
        return array([norm(cross(array(q)-p, u))/norm(u) for q in izip(x, y, z)])
        
    # die ersten 3 sind die parameter der richtungsvektors 
    p0 = array([1, 0, 1, 0, 0, 0]) 
    
    solution = scipy.optimize.leastsq(residuum, p0, args=(x, y, z)) 
    
    #print solution
    print array(solution[0][:3])/norm(array(solution[0][:3]))

if __name__ == "__main__": 
    main()

Code: Alles auswählen

[ 0.00114536 -0.04252567 -0.99909472]

[zipdrive]

Ganz ganz besten Dank für deine Antwort, bekommst ein Bienchen von mir .

Das meine Modellfunktion falsch war hab ich selber schon mitbekommen.