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: