Attachment 'least_squares_circle_v2.py'
Download
Toggle line numbers
1 #! /usr/bin/env python
2 # -*- coding: utf-8 -*-
3
4 """
5 http://www.scipy.org/Cookbook/Least_Squares_Circle
6 """
7
8 from numpy import *
9
10 # Coordinates of the 2D points
11
12 x = r_[ 9, 35, -13, 10, 23, 0]
13 y = r_[ 34, 10, 6, -14, 27, -10]
14
15 # x = r_[36, 36, 19, 18, 33, 26]
16 # y = r_[14, 10, 28, 31, 18, 26]
17
18
19 # == METHOD 1 ==
20 method_1 = 'algebraic'
21
22 # coordinates of the barycenter
23 x_m = mean(x)
24 y_m = mean(y)
25
26 # calculation of the reduced coordinates
27 u = x - x_m
28 v = y - y_m
29
30 # linear system defining the center in reduced coordinates (uc, vc):
31 # Suu * uc + Suv * vc = (Suuu + Suvv)/2
32 # Suv * uc + Svv * vc = (Suuv + Svvv)/2
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 # Solving the linear system
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 # Calculation of all distances from the center (xc_1, yc_1)
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 # Decorator to count functions calls
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 # == METHOD 2 ==
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 # == METHOD 3 ==
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 # for implicit function :
112 # data.x contains both coordinates of the points
113 # data.y is the dimensionality of the response
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 # Summary
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 # plotting functions
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 # draw
166 p.xlabel('x')
167 p.ylabel('y')
168 p.legend(loc='best',labelspacing=0.1)
169
170 # plot the residu fields
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 # vim: set et sts=4 sw=4:
New Attachment
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.- [get | view] (2011-03-20 09:33:13, 125.4 KB) [[attachment:arc_residu2_v1.png]]
- [get | view] (2011-03-20 09:51:22, 111.5 KB) [[attachment:arc_residu2_v2.png]]
- [get | view] (2011-03-20 10:25:35, 92.9 KB) [[attachment:arc_residu2_v3.png]]
- [get | view] (2011-03-20 15:12:18, 117.3 KB) [[attachment:arc_residu2_v4.png]]
- [get | view] (2011-03-20 15:33:13, 153.6 KB) [[attachment:arc_residu2_v5.png]]
- [get | view] (2011-03-21 22:06:57, 137.4 KB) [[attachment:arc_residu2_v6.png]]
- [get | view] (2011-03-19 22:05:21, 108.4 KB) [[attachment:arc_v1.png]]
- [get | view] (2011-03-20 10:25:10, 87.9 KB) [[attachment:arc_v2.png]]
- [get | view] (2011-03-20 15:11:49, 108.1 KB) [[attachment:arc_v3.png]]
- [get | view] (2011-03-20 15:32:44, 149.4 KB) [[attachment:arc_v4.png]]
- [get | view] (2011-03-21 22:06:23, 126.9 KB) [[attachment:arc_v5.png]]
- [get | view] (2011-03-20 09:07:45, 113.4 KB) [[attachment:full_cercle_residu2_v1.png]]
- [get | view] (2011-03-20 10:24:26, 79.1 KB) [[attachment:full_cercle_residu2_v2.png]]
- [get | view] (2011-03-20 15:11:25, 93.4 KB) [[attachment:full_cercle_residu2_v3.png]]
- [get | view] (2011-03-20 15:32:12, 116.0 KB) [[attachment:full_cercle_residu2_v4.png]]
- [get | view] (2011-03-21 22:05:38, 120.8 KB) [[attachment:full_cercle_residu2_v5.png]]
- [get | view] (2011-03-19 22:04:02, 104.2 KB) [[attachment:full_cercle_v1.png]]
- [get | view] (2011-03-20 10:24:01, 67.2 KB) [[attachment:full_cercle_v2.png]]
- [get | view] (2011-03-20 15:10:50, 81.4 KB) [[attachment:full_cercle_v3.png]]
- [get | view] (2011-03-20 15:31:32, 103.7 KB) [[attachment:full_cercle_v4.png]]
- [get | view] (2011-03-21 22:05:03, 101.3 KB) [[attachment:full_cercle_v5.png]]
- [get | view] (2011-03-19 22:22:43, 4.2 KB) [[attachment:least_squares_circle.py]]
- [get | view] (2011-03-20 09:32:03, 0.0 KB) [[attachment:least_squares_circle_v1.py]]
- [get | view] (2011-03-20 14:14:56, 7.4 KB) [[attachment:least_squares_circle_v1b.py]]
- [get | view] (2011-03-20 19:16:27, 9.1 KB) [[attachment:least_squares_circle_v1c.py]]
- [get | view] (2011-03-21 22:02:12, 9.2 KB) [[attachment:least_squares_circle_v1d.py]]
- [get | view] (2011-03-20 10:26:23, 5.4 KB) [[attachment:least_squares_circle_v2.py]]
- [get | view] (2011-03-20 14:00:31, 7.5 KB) [[attachment:least_squares_circle_v3.py]]