# Bresenham Ellipse Algorithm

```//
// CONIC  2D Bresenham-like conic drawer.
//       CONIC(Sx,Sy, Ex,Ey, A,B,C,D,E,F) draws the conic specified
//       by A x^2 + B x y + C y^2 + D x + E y + F = 0, between the
//       start point (Sx, Sy) and endpoint (Ex,Ey).

// Author: Andrew W. Fitzgibbon (andrewfg@ed.ac.uk),
//         Machine Vision Unit,
//         Dept. of Artificial Intelligence,
//         Edinburgh University,
//
// Date: 31-Mar-94
// Version 2: 6-Oct-95
//      Bugfixes from Arne Steinarson <arst@ludd.luth.se>

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

class ConicPlotter {
public:
// A xx + B xy + C yy + D x + E y + F = 0

int A;
int B;
int C;
int D;
int E;
int F;

ConicPlotter();
ConicPlotter(int A, int B, int C, int D, int E, int F);
void assign(int A, int B, int C, int D, int E, int F);
void assignf(double SCALE, double A, double B, double C, double D, double E, double F);

// This is the routine which draws the conic, calling the plot member at each pixel.
int draw(int xs, int ys, int xe, int ye);

// User supplies this plot routine, which does the work.
virtual void plot(int x, int y);
};

static int DIAGx[] = {999, 1,  1, -1, -1, -1, -1,  1,  1};
static int DIAGy[] = {999, 1,  1,  1,  1, -1, -1, -1, -1};
static int SIDEx[] = {999, 1,  0,  0, -1, -1,  0,  0,  1};
static int SIDEy[] = {999, 0,  1,  1,  0,  0, -1, -1,  0};

#define DEBUG 0

inline int odd(int n)
{
return n&1;
}

int getoctant(int gx, int gy)
{
// Use gradient to identify octant.
int upper = abs(gx)>abs(gy);
if (gx>=0)                            // Right-pointing
if (gy>=0)                          //    Up
return 4 - upper;
else                                //    Down
return 1 + upper;
else                                  // Left
if (gy>0)                           //    Up
return 5 + upper;
else                                //    Down
return 8 - upper;
}

int ConicPlotter::draw(int xs, int ys, int xe, int ye)
{
A *= 4;
B *= 4;
C *= 4;
D *= 4;
E *= 4;
F *= 4;

if (DEBUG) fprintf(stderr,"ConicPlotter::draw -- %d %d %d %d %d %d\n", A,B,C,D,E,F);

// Translate start point to origin...
F = A*xs*xs + B*xs*ys + C*ys*ys + D*xs + E*ys + F;
D = D + 2 * A * xs + B * ys;
E = E + B * xs + 2 * C * ys;

// Work out starting octant
int octant = getoctant(D,E);

int dxS = SIDEx[octant];
int dyS = SIDEy[octant];
int dxD = DIAGx[octant];
int dyD = DIAGy[octant];

int d,u,v;
switch (octant) {
case 1:
d = A + B/2 + C/4 + D + E/2 + F;
u = A + B/2 + D;
v = u + E;
break;
case 2:
d = A/4 + B/2 + C + D/2 + E + F;
u = B/2 + C + E;
v = u + D;
break;
case 3:
d = A/4 - B/2 + C - D/2 + E + F;
u = -B/2 + C + E;
v = u - D;
break;
case 4:
d = A - B/2 + C/4 - D + E/2 + F;
u = A - B/2 - D;
v = u + E;
break;
case 5:
d = A + B/2 + C/4 - D - E/2 + F;
u = A + B/2 - D;
v = u - E;
break;
case 6:
d = A/4 + B/2 + C - D/2 - E + F;
u = B/2 + C - E;
v = u - D;
break;
case 7:
d = A/4 - B/2 + C + D/2 - E + F;
u =  -B/2 + C - E;
v = u + D;
break;
case 8:
d = A - B/2 + C/4 + D - E/2 + F;
u = A - B/2 + D;
v = u - E;
break;
default:
fprintf(stderr,"FUNNY OCTANT\n");
abort();
}

int k1sign = dyS*dyD;
int k1 = 2 * (A + k1sign * (C - A));
int Bsign = dxD*dyD;
int k2 = k1 + Bsign * B;
int k3 = 2 * (A + C + Bsign * B);

// Work out gradient at endpoint
int gxe = xe - xs;
int gye = ye - ys;
int gx = 2*A*gxe +   B*gye + D;
int gy =   B*gxe + 2*C*gye + E;

int octantcount = getoctant(gx,gy) - octant;
if (octantcount < 0)
octantcount = octantcount + 8;
else if (octantcount==0)
if((xs>xe && dxD>0) || (ys>ye && dyD>0) ||
(xs<xe && dxD<0) || (ys<ye && dyD<0))
octantcount +=8;

if (DEBUG)
fprintf(stderr,"octantcount = %d\n", octantcount);

int x = xs;
int y = ys;

while (octantcount > 0) {
if (DEBUG)
fprintf(stderr,"-- %d -------------------------\n", octant);

if (odd(octant)) {
while (2*v <= k2) {
// Plot this point
plot(x,y);

// Are we inside or outside?
if (DEBUG) fprintf(stderr,"x = %3d y = %3d d = %4d\n", x,y,d);
if (d < 0) {                    // Inside
x = x + dxS;
y = y + dyS;
u = u + k1;
v = v + k2;
d = d + u;
}
else {                          // outside
x = x + dxD;
y = y + dyD;
u = u + k2;
v = v + k3;
d = d + v;
}
}

d = d - u + v/2 - k2/2 + 3*k3/8;
// error (^) in Foley and van Dam p 959, "2nd ed, revised 5th printing"
u = -u + v - k2/2 + k3/2;
v = v - k2 + k3/2;
k1 = k1 - 2*k2 + k3;
k2 = k3 - k2;
int tmp = dxS; dxS = -dyS; dyS = tmp;
}
else {                              // Octant is even
while (2*u < k2) {
// Plot this point
plot(x,y);
if (DEBUG) fprintf(stderr,"x = %3d y = %3d d = %4d\n", x,y,d);

// Are we inside or outside?
if (d > 0) {                    // Outside
x = x + dxS;
y = y + dyS;
u = u + k1;
v = v + k2;
d = d + u;
}
else {                          // Inside
x = x + dxD;
y = y + dyD;
u = u + k2;
v = v + k3;
d = d + v;
}
}
int tmpdk = k1 - k2;
d = d + u - v + tmpdk;
v = 2*u - v + tmpdk;
u = u + tmpdk;
k3 = k3 + 4*tmpdk;
k2 = k1 + tmpdk;

int tmp = dxD; dxD = -dyD; dyD = tmp;
}

octant = (octant&7)+1;
octantcount--;
}

// Draw final octant until we reach the endpoint
if (DEBUG)
fprintf(stderr,"-- %d (final) -----------------\n", octant);

if (odd(octant)) {
while (2*v <= k2) {
// Plot this point
plot(x,y);
if (x == xe && y == ye)
break;
if (DEBUG) fprintf(stderr,"x = %3d y = %3d d = %4d\n", x,y,d);

// Are we inside or outside?
if (d < 0) {                      // Inside
x = x + dxS;
y = y + dyS;
u = u + k1;
v = v + k2;
d = d + u;
}
else {                            // outside
x = x + dxD;
y = y + dyD;
u = u + k2;
v = v + k3;
d = d + v;
}
}
}
else {                                // Octant is even
while ((2*u < k2)) {
// Plot this point
plot(x,y);
if (x == xe && y == ye)
break;
if (DEBUG) fprintf(stderr,"x = %3d y = %3d d = %4d\n", x,y,d);

// Are we inside or outside?
if (d > 0) {                      // Outside
x = x + dxS;
y = y + dyS;
u = u + k1;
v = v + k2;
d = d + u;
}
else {                            // Inside
x = x + dxD;
y = y + dyD;
u = u + k2;
v = v + k3;
d = d + v;
}
}
}

return 1;
}

void ConicPlotter::plot(int x, int y)
{
printf("ConicPlotter::plot(%d, %d)\n", x, y);
}

ConicPlotter::ConicPlotter() { assign(0,0,0,0,0,0); }
ConicPlotter::ConicPlotter(int A, int B, int C, int D, int E, int F) { assign(A,B,C,D,E,F); }

void ConicPlotter::assign(int A_, int B_, int C_, int D_, int E_, int F_)
{
A = A_;
B = B_;
C = C_;
D = D_;
E = E_;
F = F_;
}

inline int rnd(double x) { return (x>=0.0)?(int)(x + 0.5):(int)(x - 0.5); }

void ConicPlotter::assignf(double scale,
double A_, double B_, double C_, double D_, double E_, double F_)
{
A = rnd(A_ * scale);
B = rnd(B_ * scale);
C = rnd(C_ * scale);
D = rnd(D_ * scale);
E = rnd(E_ * scale);
F = rnd(F_ * scale);
}
```