1
2
3
4 """
5 http://www.scipy.org/Cookbook/Least_Squares_Circle
6 """
7
8 from numpy import *
9
10
11
12 x = r_[ 9, 35, -13, 10, 23, 0]
13 y = r_[ 34, 10, 6, -14, 27, -10]
14
15
16
17
18
19
20 method_1 = 'algebraic'
21
22
23 x_m = mean(x)
24 y_m = mean(y)
25
26
27 u = x - x_m
28 v = y - y_m
29
30
31
32
33
34 Suv = sum(u*v)
35 Suu = sum(u**2)
36 Svv = sum(v**2)
37 Suuv = sum(u**2 * v)
38 Suvv = sum(u * v**2)
39 Suuu = sum(u**3)
40 Svvv = sum(v**3)
41
42
43 A = array([ [ Suu, Suv ], [Suv, Svv]])
44 B = array([ Suuu + Suvv, Svvv + Suuv ])/2.0
45 uc, vc = linalg.solve(A, B)
46
47 xc_1 = x_m + uc
48 yc_1 = y_m + vc
49
50
51 Ri_1 = sqrt((x-xc_1)**2 + (y-yc_1)**2)
52 R_1 = mean(Ri_1)
53 residu_1 = sum((Ri_1-R_1)**2)
54 residu2_1= sum((Ri_1**2-R_1**2)**2)
55
56
57 import functools
58 def countcalls(fn):
59 "decorator function count function calls "
60
61 @functools.wraps(fn)
62 def wrapped(*args):
63 wrapped.ncalls +=1
64 return fn(*args)
65
66 wrapped.ncalls = 0
67 return wrapped
68
69
70 from scipy import optimize
71
72 method_2 = "leastsq"
73
74 def calc_R(c):
75 """ calculate the distance of each 2D points from the center c=(xc, yc) """
76 return sqrt((x-c[0])**2 + (y-c[1])**2)
77
78 @countcalls
79 def f_leastsq(c):
80 """ calculate the algebraic distance between the 2D points and the mean circle centered at c=(xc, yc) """
81 Ri = calc_R(c)
82 return Ri - Ri.mean()
83
84 center_estimate = x_m, y_m
85 center_2, ier = optimize.leastsq(f_leastsq, center_estimate)
86
87 xc_2, yc_2 = center_2
88 Ri_2 = calc_R(center_2)
89 R_2 = Ri_2.mean()
90 residu_2 = sum((Ri_2 - R_2)**2)
91 residu2_2 = sum((Ri_2**2-R_2**2)**2)
92 ncalls_2 = f_leastsq.ncalls
93
94
95 from scipy import odr
96
97 method_3 = "odr"
98
99 @countcalls
100 def f_odr(beta, x):
101 """ implicit function of the circle """
102 xc, yc, r = beta
103 return (x[0]-xc)**2 + (x[1]-yc)**2 -r**2
104
105 def calc_estimate(data):
106 """ Return a first estimation on the parameter from the data """
107 xc0, yc0 = data.x.mean(axis=1)
108 r0 = sqrt((data.x[0]-xc0)**2 +(data.x[1] -yc0)**2).mean()
109 return xc0, yc0, r0
110
111
112
113
114 lsc_data = odr.Data(row_stack([x, y]), y=1)
115 lsc_model = odr.Model(f_odr, implicit=True, estimate=calc_estimate)
116 lsc_odr = odr.ODR(lsc_data, lsc_model)
117 lsc_out = lsc_odr.run()
118
119 xc_3, yc_3, R_3 = lsc_out.beta
120 Ri_3 = calc_R([xc_3, yc_3])
121 residu_3 = sum((Ri_3 - R_3)**2)
122 residu2_3 = sum((Ri_3**2-R_3**2)**2)
123 ncalls_3 = f_odr.ncalls
124
125
126 fmt = '%-15s %10.5f %10.5f %10.5f %10d %10.6f %10.6f %10.2f'
127 print ('\n%-15s' +' %10s'*7) % tuple('METHOD Xc Yc Rc nb_calls std(Ri) residu residu2'.split())
128 print '-'*(15 +7*(10+1))
129 print fmt % (method_1, xc_1, yc_1, R_1, 1, Ri_1.std(), residu_1, residu2_1)
130 print fmt % (method_2, xc_2, yc_2, R_2, ncalls_2, Ri_2.std(), residu_2, residu2_2)
131 print fmt % (method_3, xc_3, yc_3, R_3, ncalls_3, Ri_3.std(), residu_3, residu2_3)
132
133
134
135 from matplotlib import pyplot as p, cm
136
137 def plot_all(residu2=False, basename='circle'):
138 """ Draw data points, best fit circles and center for the three methods,
139 and adds the iso contours corresponding to the fiel residu or residu2
140 """
141
142 f = p.figure(figsize=(6.5, 4.5), dpi=90, facecolor='white')
143 p.axis('equal')
144
145 p.plot(x, y, 'ro', label='data', ms=9, mec='b', mew=1)
146
147 theta_fit = linspace(-pi, pi, 180)
148
149 x_fit1 = xc_1 + R_1*cos(theta_fit)
150 y_fit1 = yc_1 + R_1*sin(theta_fit)
151 p.plot(x_fit1, y_fit1, 'b-' , label=method_1, lw=2)
152
153 x_fit2 = xc_2 + R_2*cos(theta_fit)
154 y_fit2 = yc_2 + R_2*sin(theta_fit)
155 p.plot(x_fit2, y_fit2, 'k--', label=method_2, lw=2)
156
157 x_fit3 = xc_3 + R_3*cos(theta_fit)
158 y_fit3 = yc_3 + R_3*sin(theta_fit)
159 p.plot(x_fit3, y_fit3, 'r-.', label=method_3, lw=2)
160
161 p.plot([xc_1], [yc_1], 'bD', mec='y', mew=1)
162 p.plot([xc_2], [yc_2], 'gD', mec='r', mew=1)
163 p.plot([xc_3], [yc_3], 'kD', mec='w', mew=1)
164
165
166 p.xlabel('x')
167 p.ylabel('y')
168 p.legend(loc='best',labelspacing=0.1)
169
170
171 nb_pts = 100
172
173 p.draw()
174 xmin, xmax = p.xlim()
175 ymin, ymax = p.ylim()
176
177 vmin = min(xmin, ymin)
178 vmax = max(xmax, ymax)
179
180 xg, yg = ogrid[vmin:vmax:nb_pts*1j, vmin:vmax:nb_pts*1j]
181 xg = xg[..., newaxis]
182 yg = yg[..., newaxis]
183
184 Rig = sqrt( (xg - x)**2 + (yg - y)**2 )
185 Rig_m = Rig.mean(axis=2)[..., newaxis]
186
187 if residu2 : residu = sum( (Rig**2 - Rig_m**2)**2 ,axis=2)
188 else : residu = sum( (Rig-Rig_m)**2 ,axis=2)
189
190 lvl = exp(linspace(log(residu.min()), log(residu.max()), 15))
191
192 p.contourf(xg.flat, yg.flat, residu.T, lvl, alpha=0.75, cmap=cm.Purples_r)
193 cbar = p.colorbar(format='%.f')
194
195 if residu2 : cbar.set_label('Residu_2')
196 else : cbar.set_label('Residu')
197
198 p.xlim(xmin=vmin, xmax=vmax)
199 p.ylim(ymin=vmin, ymax=vmax)
200
201 p.grid()
202 p.title('Leasts Squares Circle')
203 p.savefig('%s_residu%d.png' % (basename, 2 if residu2 else 1))
204
205 plot_all(residu2=False, basename='circle')
206 plot_all(residu2=True , basename='circle')
207
208 p.show()
209