Share | home | download page | send comment | send link | add bookmark |

Finding a General Conic Equation given 5 points

by Stephen R. Schmitt

 

The General Conic Equation

A conic equation represents the intersection of an extended cone and a plane. The equation below represents an ellipse with axes aligned with the coordinate frame and with center at the origin.

 x2  +   y2  = 1
 a2  b2

This equation can be changed by a translation that moves the center of the ellipse to a point {p, q} by replacing x, y using:

x  ← xp
y  ← yq

It can also be changed by a rotation by angle θ around the coordinate origin by again replacing x, y using:

x  ←    x cos θ + y sin θ
y  ← − x sin θ + y cos θ

If the ellipse is translated or rotated or both, the resulting equation will have terms in x2, xy, y2, x, y, and a constant. Transformation of the normal equations for a parabola or a hyperbola would also produce a new equation with these terms. Then, the general form of the equation of a conic section is:

Ax2 + Bxy + Cy2 + Dx + Ey + F = 0

Fitting a Conic Section Through Five Points

The general conic section has five independent coefficients (the equation can be divided by any of the six coefficients without changing the equality). Therefore, five points are necessary and sufficient to define a unique conic section. That is, one and only one conic section can be drawn through five points.

One method to determine the five coefficients is to substitute each of the {x, y} points into the equation to obtain five linear equations using five independent coefficients as variables. For example, to find the conic section that fits these five points:

{0.5, 8.0}, {1.0, 4.0}, {2.0, 2.0}, {4.0, 1.0}, {8.0, 0.5}

Solve the system of equations (and choosing F = 4):

0.25A  + 4.00B + 64.00C + 0.50D + 8.00E + 4.00  = 0
1.00A  + 4.00B + 16.00C + 1.00D + 4.00E + 4.00  = 0
4.00A  + 4.00B +  4.00C + 2.00D + 2.00E + 4.00  = 0
16.00A  + 4.00B +  1.00C + 4.00D + 1.00E + 4.00  = 0
64.00A  + 4.00B +  0.25C + 8.00D + 0.50E + 4.00  = 0

to obtain the solution:

xy + 4.0000 = 0

The program below uses Gaussian elimination to obtain the coefficients of a conic equation given five points.

Zeno source code

const N : int := 5

type point : record
             x, y : real
             end record

type fivepoints : array[N] of point

type conic : array[N+1] of real

var a : array[N,N+1] of real
var x : array[N] of real
var constant : real := 6        % arbitrary

program

    var t : int
    var r, c : conic
    var p : fivepoints

    put "\nConic from 5 points"
    p[1].x := -6.0
    p[1].y :=  6.0
    p[2].x := -4.0
    p[2].y :=  3.0
    p[3].x := -2.0
    p[3].y :=  2.0
    p[4].x :=  0.0
    p[4].y :=  3.0
    p[5].x :=  2.0
    p[5].y :=  6.0
    
    initialize( p )             % array of conic equations

    if eliminate then           % lower diagonal elements
        substitute              % and solve for x
        set_conic( c, x[1], x[2], x[3], x[4], x[5], constant )
        put_conic( c )
    else
        put "singular matrix!"
    end if
    
end program

% set values of conic coefficient array
procedure set_conic( var coef : conic, a, b, c, d, e, f : real )

    coef[1] := a
    coef[2] := b
    coef[3] := c
    coef[4] := d
    coef[5] := e
    coef[6] := f

end procedure

% display a general conic equation
procedure put_conic( var x : conic )

    var i : int
    var s : boolean := false

    for i := 1...6 do

        continue when x[i] = 0

        if x[i] < 0 then
            put " - "...
        elsif (i > 1) and s then
            put " + "...
        end if

        s := true

        if (abs(x[i]) ~= 1) or (i = 6) then
            put abs(x[i]):5:4...
            if i ~= 6 then
                put ''...
            end if
        end if

        case i of
        value 1: put "x"...
        value 2: put "xy"...
        value 3: put "y"...
        value 4: put "x"...
        value 5: put "y"...
        end case

    end for

    put " = 0"

end procedure

% set elements of matrix A
procedure initialize( var p : fivepoints )

    var i : int

    for i := 1...5 do
        a[i,1] :=  p[i].x^2         % A
        a[i,2] :=  p[i].x * p[i].y  % B
        a[i,3] :=  p[i].y^2         % C
        a[i,4] :=  p[i].x           % D
        a[i,5] :=  p[i].y           % E
        a[i,6] := -constant         % -F
        x[i]   :=  0
    end for

end procedure

% convert matrix A to upper diagonal form
function eliminate : boolean

    var i, j, k, max_row : int
    var t : real
    var ok : boolean := true

    for i := 1...N do
 
        % find row with maximum in column i
        max_row := i
        for j := i...N do
            if abs( a[j,i] ) > abs( a[max_row,i] ) then
                max_row := j
            end if
        end for

        % swap max row with row i of A and b
        for k := i...N+1 do
            t            := a[i,k]
            a[i,k]       := a[max_row,k]
            a[max_row,k] := t
        end for

        % eliminate lower diagonal elements
        for j := i+1 ... N do
            for decreasing k := N+1...i do
                if a[i,i] = 0.0 then
                    ok := false
                else
                    a[j,k] := a[j,k] - a[i,k] * a[j,i] / a[i,i]
                end if
            end for
        end for

    end for

    % need to check last diagonal element
    if a[N,N] = 0.0 then
        ok := false
    end if
    
    return ok

end function

% compute the values of vector x starting from the bottom
procedure substitute

    var j, k : int
    var t : real

    for decreasing j := N...1 do
        t := 0.0
        for k := j+1...N do
            t := t + a[j,k] * x[k]
        end for
        x[j] := ( a[j,N+1] - t ) / a[j,j]
    end for

end procedure

% return the absolute value of argument
function abs( x : real ) : real

    if x < 0 then
        x := -x
    end if
    return x

end function

Sample output

Conic from 5 points
0.5000x + 0.0000xy + 0.0000y + 2.0000x - 2.0000y + 6.0000 = 0

Share | home | download page | send comment | send link | add bookmark |


Copyright © 2005, Stephen R. Schmitt