1
$\begingroup$

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.

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

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)
| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.