In [2]:
%matplotlib inline

import matplotlib
import numpy as np
import matplotlib.pyplot as pt
import numpy.linalg as la


# introduce "noise" by adding normally distributed random numbers¶

In [96]:
np.random.seed(555)
t = np.linspace(0,2*np.pi,100)
xc = 3*np.random.rand() #x-coordinate of center of circle
yc = 3*np.random.rand() #y-coordinate of center of circle
r = 2*np.random.rand() + 0.5 #radius of circle

x = r*np.cos(t) + np.random.normal(scale = 1.0/7, size = len(t)) + xc
y = r*np.sin(t) + np.random.normal(scale = 1.0/7, size = len(t)) + yc

pt.plot(x,y,'o')
pt.grid(True)
pt.axes().set_aspect('equal','datalim')


# what circle best describes this data??¶

### Let's use the data to "guess" the radius and center¶

In [97]:
xc = 1.0
yc = 2.0
r = 2.5
xapprox = r*np.cos(t) + xc
yapprox = r*np.sin(t) + yc

pt.plot(x,y,'o')
pt.plot(xapprox,yapprox,'r')
pt.grid(True)
pt.axes().set_aspect('equal','datalim')

In [98]:
def g(x, y, xc, yc, r):
return (x-xc)**2 + (y-yc)**2 - r**2

def F(x, y, xc, yc, r):
return np.sum(g(x, y, xc, yc, r)**2)

def gradF(x, y, xc, yc, r):
gg = g(x,y,xc,yc,r)
F1 = -4 * np.sum(gg*(x-xc))
F2 = -4 * np.sum(gg*(y-yc))
F3 = -4 * np.sum(gg*r)
return np.array([F1,F2,F3])

def HessF(x,y,x0,y0,r):
gg = g(x,y,xc,yc,r)
g11 = np.sum(4*gg + 8*(x-xc)**2)
g12 = np.sum(8*(x-xc)*(y-yc))
g13 = np.sum(8*(x-xc)*r)
g22 = np.sum(4*gg + 8*(y-yc)**2)
g23 = np.sum(8*(y-yc)*r)
g33 = np.sum(-4*gg + 8*r**2)
H = np.array([[g11, g12, g13],[g12,g22,g23],[g13,g23,g33]])
return H

In [114]:
xc = 2.0
yc = 1.0
r = 5.5

xx = np.array([xc,yc,r])# initial guess for a,b, and r
tol = 1e-15

k = 0

In [134]:
xc = xx[0]
yc = xx[1]
r = xx[2]
HF = HessF(x,y,xc,yc,r)    #compute current Hessian
step = la.solve(HF,gF)     #solve H*step = g
xx = xx + step             #next iterate
xapprox = xx[2]*np.cos(t) + xx[0]
yapprox = xx[2]*np.sin(t) + xx[1]
pt.plot(xapprox,yapprox,'r')
pt.plot(x,y,'o')
pt.grid(True)
pt.axes().set_aspect('equal','datalim')

print(k, gradF(x, y, xc, yc, r))

0 [ -2.80664381e-13  -4.01706446e-13   7.26529947e-13]

In [ ]: