# Least Squares Involving Binomial Term (Regressing a Circle)

Can I use least squares to solve an overdetermined system that involves binomial terms?

I am trying to regress the center $(h, k)$ and radius $r$ of a circle using noisy measurements $(x_1,y_1), (x_2, y_2), ... , (x_n, y_n)$ along its perimeter. The formula of the relationship between a point on the perimeter $(x, y)$ of a circle and its center and radius is:

$(x - h)^2 + (y - k)^2 = r \space$ (1)

I tried expanding the formula and stick into a least squares form, but I got stuck, because each of the $x, y$ and $h, k$ appear several times in the formula and I do not know how to handle this.

$0 = x^2 - 2xh + h^2 + y^2 - 2y k + k^2 - r \space$ (2)

$\begin{bmatrix} x_1^2 & x_1 & 1 & y_1^2 & y_1 & 1 & -1 \\ x_2^2 & x_2 & 1 & y_2^2 & y_2 & 1 & -1\\ ... \\ x_n^2 & x_n & 1 & y_n^2 & y_n & 1 & -1 \end{bmatrix}$ * $\begin{bmatrix} 1 \\ 2h \\ h^2 \\ 1 \\ 2k \\ k^2 \\ r \end{bmatrix}$ = 0

Essentially, what do I do about the $2h$ and $2k$? Examples I have found do not have this problem, because even in examples for polynomials the unknowns are linear.

• The title is confusing because there's nothing "binomial" about this problem. I believe you may be referring to quadratic terms in a formula. – whuber Feb 3 '16 at 20:47

Via wikipedia, I got the right hint. There is a paper called Circle fitting by linear and nonlinear least squares from 1993 that solves my specific problem for circles.

Below my reference implementation in Python:

import numpy as np
import pylab as plt

def circle_fit_coope(a):
""" Fitting circle to points in a
Using linear least squares from Circle Fitting by Linear and Nonlinear Least Squares
L D. Coope, JOURNAL OF OPTIMIZATION THEORY AND APPLICATIONS: Vol. 76, No. 2, FEBRUARY 1993

input:
a     (m x n array, n: dimensionality, m: number of points)

output:
x     (cirlce centre)
"""
m = a.shape[0]
n = a.shape[1]

B = np.ones((m, n + 1))
B[:, :n] = a

d = np.sum(a**2, axis=1)

y = np.linalg.lstsq(B, d)[0]

x = (y[:-1].T / 2).reshape(-1, 1)

r = np.sqrt(y[2] + np.dot(x.T, x))

return x.flatten(), r

def generate_example_points_from_coope():
return np.asarray([  [0.7, 4.0],
[3.3, 4.7],
[5.6, 4.0],
[7.5, 1.3],
[6.4, -1.1],
[4.4, -3.0],
[0.3, -2.5],
[-1.1, 1.3] ])

def plot_result(a, x, r):
circle1=plt.Circle(x, r, color='r', fill=False)

fig = plt.gcf()