## Finding a General Conic Equation given 5 points

### 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 ← x − p y ← y − q

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 = 0 1.00A + 4.00B + 16.00C + 1.00D + 4.00E + 4 = 0 4.00A + 4.00B + 4.00C + 2.00D + 2.00E + 4 = 0 16.00A + 4.00B + 1.00C + 4.00D + 1.00E + 4 = 0 64.00A + 4.00B + 0.25C + 8.00D + 0.50E + 4 = 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 "x·y"...
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.5000·x² + 0.0000·x·y + 0.0000·y² + 2.0000·x - 2.0000·y + 6.0000 = 0
```