Details | Last modification | View Log | RSS feed
| Rev | Author | Line No. | Line |
|---|---|---|---|
| 1993 | - | 1 | import scipy.optimize |
| 2 | from pylab import * |
||
| 3 | from numpy import * |
||
| 4 | |||
| 5 | cal = open('cal1.dat') |
||
| 6 | |||
| 7 | for f in [cal]: |
||
| 8 | data = array([map(float, l.split()) for l in f.readlines()]) |
||
| 9 | |||
| 10 | N = len(data[0]) |
||
| 11 | |||
| 12 | x = data[:,0] |
||
| 13 | y = data[:,1] |
||
| 14 | z = data[:,2] |
||
| 15 | data = data[:,:3] |
||
| 16 | |||
| 17 | def cost(args): |
||
| 18 | rs = diag(1. / array(args[0:3])) ** 2 |
||
| 19 | c = array(args[3:6]) |
||
| 20 | m = data - array(c) |
||
| 21 | out = 0 |
||
| 22 | for p in m: |
||
| 23 | out += abs(dot(dot(p, rs), p) - 1) |
||
| 24 | return out |
||
| 25 | guess = array([250, 250, 250, 0, 0, 0]) |
||
| 26 | print cost(guess) |
||
| 27 | best = scipy.optimize.fmin(cost, guess) |
||
| 28 | print cost(best) |
||
| 29 | print best |
||
| 30 | print 1/best[:3] |
||
| 31 | |||
| 32 | here |
||
| 33 | A = array([[cx + r * cos(theta), |
||
| 34 | cy + r * sin(theta)] for theta in arange(0, 2 * pi, 1 * pi/180)]) |
||
| 35 | # plot(A[:,0], A[:,1], 'g-') |
||
| 36 | xy = data[:,1:3] |
||
| 37 | |||
| 38 | def cost(params): |
||
| 39 | cx, cy, r = params |
||
| 40 | xy_ = xy - [cx, cy] |
||
| 41 | thetas = arctan2(xy_[:,1], xy_[:,0]) |
||
| 42 | resids = xy_ - transpose([r * cos(thetas), r * sin(thetas)]) |
||
| 43 | return sum(ravel(resids ** 2)) |
||
| 44 | |||
| 45 | cx, cy, r = scipy.optimize.fmin(cost, [cx, cy, r], disp=False) |
||
| 46 | print f.name, cx, cy, r |
||
| 47 | ## acc_cal_lights_on.dat 550.150958354 507.218838209 249.831129791 |
||
| 48 | ## acc_cal_lights_off.dat 563.391868993 518.281081432 251.367556713 |
||
| 49 | |||
| 50 | A = array([[cx + r * cos(theta), |
||
| 51 | cy + r * sin(theta)] for theta in arange(0, 2 * pi, 1 * pi/180)]) |
||
| 52 | plot(A[:,0], A[:,1]) |