Least Squares Involving Binomial Term (Regressing a Circle)

Question

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.

Answer 1

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)
        r     (radius)
    """
    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()
    fig.gca().add_artist(circle1)
    fig.gca().plot(x[0], x[1], 'xr')

    fig.gca().scatter(a[:,0], a[:, 1])
    fig.gca().set_aspect('equal')

    plt.show()

a = generate_example_points_from_coope()
x, r = circle_fit_coope(a)
plot_result(a, x, r)