Fitting conics to data

Fitting conics to data

Written by

Dylan Holmes

fitted_hyperbola.png

The conic sections are a family of shapes consisting of parabolas, hyperbolas, and ellipses. Conic sections show up in many scientific applications; for example, celestial orbits trace out different conic sections depending on their energy 1. The shadow cast by a sundial's "spike" (gnomon) traces out a different conic section depending on latitude and time of year 2. And mirrors shaped like conic sections have important optical properties for focusing light 3.

It's therefore useful to be able to take some imprecise physical measurements and fit a conic-shaped curve to them. In this article, I provide some code for doing this. The figure above illustrates fitting some real data points for a sundial; the algorithm finds the hyperbolic shape shown.

Algebraically, conic sections are the graphs of the equation: $$Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0$$

(It's just any general polynomial you can make involving \(x\) and \(y\) where no term has a combined exponent bigger than two.) There is some beautiful linear algebra under the hood; one useful fact is that the discriminant \(B^2-4AC\) is either positive, negative, or zero, which tells you whether the shape is an ellipse, hyperbola, or parabola.

To measure how well a particular conic section fits your data \(\langle x_i, y_i\rangle\), you might imagine a straightforward algorithm: compute \(A x_i^2 + Bx_iy_i + C y_i^2 + D x_i + E y_i + F\) for each data point, square those values and add them up, then see how far away the sum is from zero. To find the optimal fit, vary the coefficients to make this sum as close to zero as possible.

This basic approach almost works, despite the search space being unstable 4. You can modify it to make it more stable and to allow the user to specify a target conic type (hyperbola, parabola, or ellipse) in advance. In a paper by Nievergelt, the author describes how to do this. I implemented the suggested algorithm in Python, as shown below.

Interestingly, the algorithm gets less stable the better the data fit a particular conic; a little empirical noise is useful. Because my target application (determining whether a blackbox function represents a conic) involved extremely accurate conic sections, this threw off my results. So I modified the algorithm slightly. Where the algorithm uses an SVD decomposition to decide whether the points all lie on a line, I use the Theil-Sen estimator 5.

I also found a few typos in the manuscript, which when fixed allowed me to exactly reproduce the curves in the paper.

from numpy.linalg import svd, qr, solve, det
from numpy import array, matmul, float64, transpose # , allclose, matmul, float64
from math import sqrt

#### Given a list of 2D points, return an array
#### describing the best-fit conic. Returns something special in the
#### case of degenerate conics (e.g. straight lines)

# See paper: Yves Nievergelt, "Fitting conics of specific types to
# data"


SMALL = 10**-12
sqrt2 = 2**0.5

TYPE_HYPERBOLA = "hyperbola"
TYPE_PARABOLA = "parabola"
TYPE_ELLIPSE = "ellipse"
TYPE_LINEAR = "linear"
def centered(points) :
    """Subtract the average from the (x,y) pairs."""
    xs_uncentered = [x for (x,y) in points]
    ys_uncentered = [y for (x,y) in points]

    x_center = sum([x for x in xs_uncentered])/len(points)
    y_center = sum([y for y in ys_uncentered])/len(points)

    xs = [x-x_center for x in xs_uncentered]
    ys = [y-y_center for y in ys_uncentered]

    return list(zip(xs,ys)), (x_center, y_center)


def translate_conic(conic, displacement) :
    """Return a conic that has been translated by an amount (dx,dy)."""
    (dx, dy) = displacement
    A,B,C,D,E,F = conic.get("coeffs")
    new_coeffs = (A,B,C,
                  D - 2*A*dx - B*dy,
                  E - 2*C*dy - B*dx,
                  A*dx*dx + B*dx*dy + C*dy*dy - D*dx - E*dy + F)
    conic["coeffs"] = new_coeffs
    conic["center"] = np.array(conic.get("center",(0,0)))+ np.array(displacement)
    return conic


def collinearity(points) :

    """Return a number indicating how very nearly the points follow a
line. Zero indicates perfect fit. Also return the best-fit (Theil-Sen
estimator) line in point-slope form."""

    (x0, y0) = points[0]

    average_x_deviation = sum(abs(x-x0) for (x,y) in points[1:])/len(points)
    if (average_x_deviation < SMALL) :
        # VERTICAL LINE
        return 0, (0,1), points[0]

    slopes = sorted([(y-y0)/float(x-x0) for (x,y) in points[1:]])
    median_slope = slopes[int(len(slopes)/2)]

    intercepts = sorted([y - median_slope*x for (x,y) in points[1:]])
    median_intercept = intercepts[int(len(intercepts)/2)]

    print(">>>", median_slope, median_intercept)
    average_deviation = sum((y-median_slope*x-median_intercept)**2/len(points)
                            for (x,y) in points)

    return average_deviation, (1, median_slope), (0, median_intercept)


def smallest_singular(M, compute_singular_vector=True) :
    """Return a pair containing the smallest singular value of M and a
corresponding singular vector.
    """

    if compute_singular_vector: 
        _, s, Vh = svd(M, full_matrices=True)
        return s[-1], Vh[-1]


def fit_conic(points, desired_type=None) :
    points, center = centered(points)

    matrix_monomials = array([
        [1, x, y,    y*y-x*x, 2*x*y, y*y+x*x]
        for (x,y) in points])

    # ---- Check for degenerate conic lying on a single line
    fit, slope, intercept = collinearity(points)
    print("COLLINEARITY", fit, slope, intercept)
    if (fit < SMALL) :
        return translate_conic({"type" : TYPE_LINEAR,
                                "coeffs" : (0, 0, 0, slope[1], -slope[0], slope[0]*intercept[1] - slope[1]*intercept[0]),
                                "fit" : fit,        
                                }, center)

        return None


    # ---- Perform decomposition
    submatrix = matrix_monomials[:,:3] # (1,x,y)
    q,r = qr(submatrix, "complete")


    R_full = matmul(transpose(q), matrix_monomials)
    RA = R_full[:3, 0:3]
    RB = R_full[:3, 3:6]
    RC = R_full[3:, 3:6]
    _, vec_smallest = smallest_singular(RC)

    # ----- Determine type of conic
    Z = array([[-1,0,1],
               [0,sqrt2,0],
               [1,0,1]]) / sqrt2
    D = array([[-1,0,0],
               [0,-1,0],
               [0,0,1]])/4

    A_det = transpose(vec_smallest) @ D @ vec_smallest
    A_trace = vec_smallest[2]

    conic_type = TYPE_PARABOLA if (abs(A_det) < SMALL) else TYPE_ELLIPSE if A_det > 0 else TYPE_HYPERBOLA


    # ---- Adjust answer if it doesn't match desired conic type (not
    # ---- implemented yet)

    if(False and desired_type == TYPE_PARABOLA) :
        # "desired type" not implemented yet.  also not necessary for
        # this application b/c empirical points are expected to be
        # extremely unambiguously close to particular conics.

        pass
    else :
        q = vec_smallest


    # ------ Solve for conic coefficients
    ret = solve(RA, -1 * (RB @ q)) # c, *two_b

    params = [*ret, *q]

    conic_coeffs = [ params[5]-params[3],
                         2*params[4],
                         params[5]+params[3],
                         params[1],
                         params[2],
                         params[0]]


    fit  = sum([x**2 for x in matmul(matrix_monomials, array(params))])/len(points)


    A,B,C,D,E,F = conic_coeffs


    quadratic_form = array([[A, B/2],[B/2,C]])

    full_form = array([[  A, B/2, D/2],
                       [B/2,   C, E/2],
                       [D/2, E/2,   F]])

    if(abs(det(full_form)) < 1e-8) :
            # degenerate conic
            return None


    u,s,vh = svd(quadratic_form)

    K =  -det(full_form)/det(quadratic_form) if det(quadratic_form) != 0 else None

    axes = None
    try :
        axes = [(K/x)**0.5 for x in s.tolist()]
    except Exception :
        pass

    return translate_conic({"type" : conic_type,
                            "coeffs" : conic_coeffs,
                            "fit" : fit,
                            "axes" : axes,
                            "diagonalizer" : vh.tolist()       
                            }, center)


### ---------------------------------- WORKING EXAMPLES

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline, BSpline

def newton_root(f, init_guess=60) :
    x = init_guess
    small = SMALL
    for _ in range(30) :
        if abs(f(x)) < 1e-8 :
            break
        else :
            x -= f(x)*small/(f(x+small)-f(x))
    return x

def plot_fitted_conic(points) :
    conic = fit_conic(points)
    A,B,C,D,E,F = conic.get("coeffs")
    xs = np.arange(min([x for (x,y) in points]),
                   max([x for (x,y) in points]),
                   1).tolist()
    ys = [newton_root(lambda y: A*x*x + B*x*y + C*y*y + D*x + E*y + F) for x in xs]

    plt.scatter([x for (x,y) in points],
                [y for (x,y) in points])

    ys_smooth = make_interp_spline(xs, ys, k=3)(xs)
    plt.plot(xs, ys_smooth)

    plt.show()





example_hyperbola_points = [(-34.75, 20.25),
                            (-22, 17),
                            (-15.5, 15.0),
                            (-8.0, 13.5),
                            (-4, 13),
                            (-1.0, 12.5),
                            (1.5, 12.5),
                            (4.5, 13.0),
                            (9.25, 14.00),
                            (17, 16),
                            (23.5, 18.0),
                            (36, 21),
                            (64.5, 29.5)]

example_linear_points = [(x, 2*x+60) for (x,_) in example_hyperbola_points]



plot_fitted_conic(example_hyperbola_points)
plot_fitted_conic(example_linear_points)

Footnotes:

1

In a flyby, the object has enough surplus energy to escape the gravitational pull, and gets flung out along a hyperbolic orbit. In a closed orbit, the object remains gravitationally captive, tracing an elliptic orbit. In the parabolic case, the object has just enough energy to exactly cancel the gravitational pull, and also escapes.

2

For most places, the sun rises and sets, spending the night below the horizon. In this case, the gnomon's shadow traces out a hyperbolic path over the course of the day. Near the poles, the sun might spend the entire day above the horizon, and the shadow traces out an ellipse. In the exact intermediate case, the sun barely grazes the horizon before rising again, and the shadow traces out a parabola.

3

See the Cassegrain telescope design, for example.

4

Small variations in data can cause large variations in best-fit conic; infinitesimal changes to a parabola will make it an ellipse or hyperbola; the algorithm can often fit non-conic data extremely closely, which is unfortunate for detecting conics.

5

The Theil-Sen estimator is a neat simple algorithm to find a best-fit line to some data: First, you compute all lines through all pairs of points. Then you determine the median slope \(m\) among the slopes of those lines. Second, you draw a line of slope \(m\) through each of the data points and compute their y-intercepts. Find the median intercept value \(b\). The overall fit line is \(y=mx+b\), and it's particularly resistant to outliers because we used medians throughout.