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 | from mpl_toolkits.mplot3d import Axes3 |
||
| 5 | |||
| 6 | cal = open('cal.dat') |
||
| 7 | fig = figure(1); clf() |
||
| 8 | az = Axis3D(fig) |
||
| 9 | |||
| 10 | for f in [cal]: |
||
| 11 | data = array([map(float, l.split()) for l in f.readlines()]) |
||
| 12 | |||
| 13 | N = len(data[0]) |
||
| 14 | |||
| 15 | x = data[:,0] |
||
| 16 | y = data[:,1] |
||
| 17 | z = data[:,2] |
||
| 18 | ax.scater(x, y, z) |
||
| 19 | here |
||
| 20 | |||
| 21 | A = array([[cx + r * cos(theta), |
||
| 22 | cy + r * sin(theta)] for theta in arange(0, 2 * pi, 1 * pi/180)]) |
||
| 23 | # plot(A[:,0], A[:,1], 'g-') |
||
| 24 | xy = data[:,1:3] |
||
| 25 | |||
| 26 | def cost(params): |
||
| 27 | cx, cy, r = params |
||
| 28 | xy_ = xy - [cx, cy] |
||
| 29 | thetas = arctan2(xy_[:,1], xy_[:,0]) |
||
| 30 | resids = xy_ - transpose([r * cos(thetas), r * sin(thetas)]) |
||
| 31 | return sum(ravel(resids ** 2)) |
||
| 32 | |||
| 33 | cx, cy, r = scipy.optimize.fmin(cost, [cx, cy, r], disp=False) |
||
| 34 | print f.name, cx, cy, r |
||
| 35 | ## acc_cal_lights_on.dat 550.150958354 507.218838209 249.831129791 |
||
| 36 | ## acc_cal_lights_off.dat 563.391868993 518.281081432 251.367556713 |
||
| 37 | |||
| 38 | A = array([[cx + r * cos(theta), |
||
| 39 | cy + r * sin(theta)] for theta in arange(0, 2 * pi, 1 * pi/180)]) |
||
| 40 | plot(A[:,0], A[:,1]) |