ALGORITHM ALLEY

Rendering Circles and Ellipses

Tim Kientzle

In the July 1990 issue of Dr. Dobb's Journal, Tim Paterson presented the article "Circles and the Digital Differential Analyzer." While Paterson's algorithm does not accumulate error, his explanation involves a method for solving the differential equation dy/dx=-x/y which does accumulate error.

In a subsequent letter to the editor ("Letters," DDJ, July 1991), V. Venkataraman pointed out that Paterson's algorithm plots points on or just within the ideal circle, and suggests the desirability of a method that plots the points closest to the circle, even if they are outside it. Venkataraman presented a minor change to Paterson's algorithm, but still failed to consistently choose the optimal point.

Based on this exchange, I developed a circle algorithm which satisfies both Paterson's interest in speed and Venkataraman's interest in plotting the closest points to the circle. The corresponding algorithm for drawing ellipses is faster than the line algorithm that inspired it.

Lining Up

The digital differential analyzer (DDA), also known as "Bresenham's algorithm," is based on a simple idea: Instead of picking points on the ideal line, you should pick points on the screen that are close to the line. To do this, increment one coordinate while updating the other so that you pick the integer closest to the line. Since the endpoints of the lines have integer coordinates, the slope is always rational, so you can obtain a fast, exact algorithm by using rational arithmetic.

The problem with this approach is that by simply using y to plot the point, you are truncating the exact value rather than rounding. Since rounding involves adding 1/2 and then truncating, in this algorithm it suffices to initialize the exact y-coordinate to 1/2 before you start, as in Listing Two .

Going in Circles

Drawing circles using the DDA algorithm requires more-involved geometry than drawing lines, but it uses the same approach: Step the x-coordinate through successive values and update the approximate y-coordinate based on some additional information. The trick is to figure out what that additional information should be.

To get a first cut at a circle algorithm, you compute y from sqrt(r2-x2) (see Listing Five). Although it may not look like it, the algorithm in Listing Five is almost the same as Paterson's. The difference is that his approach replaced an expensive, repeated square-root calculation with a series of additions and subtractions. This is similar to the "strength reduction" performed by many optimizing compilers that replace an expensive operation within a loop with a series of less-expensive operations, usually adding a new variable in the process. Of course, strength reduction of a square-root operation is considerably beyond what can be handled automatically by today's compilers.

Then, reduce the square-root calculation to calculation of squares. In the octant you're computing (the part of the circle from the top to the top-left), the y-coordinate of the next point is always the same or one less than the previous y-coordinate. This means that you can simply decrement y whenever y2 > radius2 - x2, or equivalently, whenever x2+y2-radius2>0.

Next, remove the calculation of the squares. Keep the square of the radius in a temporary variable. You can then use the algebraic simplification (x+1)2 = x2+2x+1 to keep track of x2 and y2 as you change x and y. This involves adding new variables xSquared and ySquared, adding 2x+1 before incrementing x and subtracting 2y-1 from y2 before decrementing y.

Finally, because the only place you need all these squares is in the condition to decrement y, you can combine them into a single error variable kept equal to xSquared+ySquared-radiusSquared. The earlier condition then simplifies to error>0, which gives exactly the algorithm Paterson describes; see Listing Six . This explanation of the variable error should make clear why this algorithm works: Whenever the current y gives a point outside the circle (whenever x2+y2>radius2), you reduce y by 1 to remain within the circle. As Venkataraman pointed out, you'd get a smoother result if you could instead choose the points closest to the circle; that is, if you could only reduce y when the lower point would be closer to the circle.

To accomplish this, take a careful look at how the algorithm decides when to decrement y. You want a condition that will tell you whether (x,y) or (x,y-1) is closer to the circle. Algebraically, you want to decrement y whenever the inequality in Figure 1(a) is true. Changing the absolute values to squares gives something equivalent, but easier to simplify. The simplification is still tedious, but you eventually end up with something useful. The left side of Figure 1(b) is simply the value error from earlier attempts. The value in brackets is always positive and less than 1. This means that the right side of the inequality is less than y but greater than y-1. Since everything else is an integer, the inequality simplifies to error >= y.

Figure 1 (a) You need to decrement y whenever the inequality is true; (b) the left side is simply the value error from earlier attempts.

This results in
Listing Seven , which is close to Venkataraman's suggestion, but always gives the point closest to the circle. To test this code, replace the plot8 routine with one that computes the radius of the points being plotted and the radius of points just above and below them. You'll see that the algorithm always picks the y-coordinate that results in a radius closest to the circle. You may be able to speed up Listing Seven slightly by changing error to be x2+y2-r2-y, which allows the condition to simply test error >= 0.

Handling Ellipses and Nonsquare Pixels

If you're only interested in drawing perfect circles on screens with square pixels, then Listing Seven serves the purpose quite nicely. However, if the pixels are not square, you have to draw an ellipse in order to get a circular result.

The biggest problem in drawing ellipses is that you lose the eight-fold symmetry used to simplify the circle routine. Where before you had eight identical octants, you now have two sets of octants: two each at the top and bottom of the circle that are identical to each other, and four at the left and right that are identical to each other. So, you need to run the earlier algorithm twice, once for each set.

The formula for the ellipse that fits within a rectangle of width 2a and height 2b is b2x2+a2y2=a2b2. Adjusting the calculation of error to match this gives us the algorithm in Listing Eight for drawing the ellipse. The least-obvious part of this is how to determine where each octant stops. Remember that the algorithm depends on y either staying the same or decreasing by 1. This means that each loop has to continue until you reach a place in the ellipse where the slope is 1 or -1. This happens exactly when b2x==a2y; see Listing Eight.

Clearly, a few more things could be done to speed up this routine. By adding variables to keep track of b2x and a2y, you could replace all of the multiplications with additions. You could also alter the definition of error as discussed earlier to simplify the condition test. Listing Nine incorporates those optimizations, plus a few others. In assembly language, you should be able to optimize this to a maximum of eight addition/subtraction/comparison operations per loop iteration, not counting the operations in the plot4 routine. Counting those, you get a maximum of 16 such operations for four points on the ellipse, or a mere four additions per point, which compares quite favorably to the three additions per point for the less-general circle routine developed earlier, and the six additions per point for the line algorithm. In fact, this shows that drawing the ellipse directly is faster than computing several points on the ellipse and using DDA lines to connect them!

Of course, it might be faster to undo some of these optimizations to reduce the number of variables required. By doing that, you may be able to fit all the variables in registers, thus speeding up the algorithm. As is usual for such low-level graphics operations, optimizing for the specific processor and video hardware is critical.

Listing One

```
void line_1(int x0, int y0, int x1, int y1 )
{
int x = x0,y = y0;          /* Current point on line */
int deltaX = x1-x0;         /* Change in x from x0 to x1 */
int deltaY = y1-y0;         /* Change in y from y0 to y1 */
int numerator = 0;          /* Numerator of fractional part of y */

while ( x <= x1 )
{
plot( x,y );
x += 1;

numerator += deltaY;        /* Increase fractional part of y */
if (numerator >= deltaX)    /* If fraction is 1 or more */
{
numerator -= deltaY;    /* Reduce fraction */
y += 1;                 /* Increase whole part of y */
}
}
}

```

Listing Two

```
#define abs(x)  (((x)>=0)?(x):-(x))  /* Absolute value */
void line_2(int x0, int y0, int x1, int y1 )
{
int x = x0,y = y0;
int deltaX = 2*abs(x1-x0);
int deltaY = 2*abs(y1-y0);
int numerator = deltaX/2;      /* Initialize y-coordinate to 1/2 */
while ( x <= x1 )
{
plot( x,y );
x += 1;

numerator += deltaY;
if (numerator >= deltaX)
{
numerator -= deltaY;
y += 1;
}
}
}

```

Listing Three

```
void plot4( int xOrigin, int yOrigin, int xOffset, int yOffset)
{
plot( xOrigin + xOffset, yOrigin + yOffset );
plot( xOrigin + xOffset, yOrigin - yOffset );
plot( xOrigin - xOffset, yOrigin + yOffset );
plot( xOrigin - xOffset, yOrigin - yOffset );
}

```

Listing Four

```
void plot8( int xOrigin, int yOrigin, int xOffset, int yOffset)
{
plot4( xOrigin, yOrigin, xOffset, yOffset );
plot4( xOrigin, yOrigin, yOffset, xOffset );
}

```

Listing Five

```
void circle_1(int xOrigin, int yOrigin, int radius)
{
int x = 0;                           /* Exact x coordinate */
int y = radius;                      /* Approximate y coordinate */

while (x <= y)                       /* Just one octant */
{
plot8( xOrigin, yOrigin, x, y );
x += 1;
y = sqrt(radius*radius - x*x);
}
}

```

Listing Six

```
void circle_2(int xOrigin, int yOrigin, int radius)
{
int x = 0;                         /* Exact x coordinate */
int y = radius;                    /* Approximate y coordinate */
int error = 0;                     /* x^2 + y^2 - r^2 */

while (x <= y)
{
plot8( xOrigin, yOrigin, x, y );

error += 2 * x + 1;
x += 1;

if (error > 0)
{
error -= 2 * y - 1;

y -= 1;
}
}
}

```

Listing Seven

```
void circle_3(int xOrigin, int yOrigin, int radius)
{
int x = 0;                  /* Exact x coordinate */
int y = radius;             /* Approximate y coordinate */
int error = 0;              /* x^2 + y^2 - r^2 */

while (x <= y)
{
plot8( xOrigin, yOrigin, x, y );

error += 2 * x + 1;
x += 1;

if (error >= y)
{
error -= 2 * y - 1;
y -= 1;
}
}
}

```

Listing Eight

```
void ellipse_1(int xOrigin, int yOrigin, int a, int b)
{
{    /* Plot the octant from the top to the top-right */
int x = 0;
int y = b;
int error = 0;/* b^2 x^2 + a^2 y^2 - a^2 b^2 */

while (x * b *b <= y * a * a)
{
plot4( xOrigin, yOrigin, x, y );

error += 2 * b*b * x + b*b;
x += 1;

if (error >= y * a*a)
{
error -= 2 * a*a * y - a*a;
y -= 1;
}
}
}

{    /* Plot the octant from right to top-right */
int x = a;
int y = 0;
int error = 0;/* b^2 x^2 + a^2 y^2 - a^2 b^2 */

while (x * b * b > y * a * a)
{
plot4( xOrigin, yOrigin, x, y );

error += 2 * a*a * y + a*a;
y += 1;

if (error >= x * b*b)
{
error -= 2 * b*b * x - b*b;
x -= 1;
}
}
}
}

```

Listing Nine

```
void ellipse_2(int xOrigin, int yOrigin, int a, int b)
{
int aSquared = a*a;
int bSquared = b*b;
int twoASquared = 2 * aSquared;
int twoBSquared = 2 * bSquared;

{    /* Plot the octant from the top to the top-right */
int x = 0;
int y = b;
int twoXTimesBSquared = 0;
int twoYTimesASquared = y * twoASquared;

int error = -y* aSquared;  /* b^2 x^2 + a^2 y^2 - a^2 b^2 - a^2y */

while (twoXTimesBSquared <= twoYTimesASquared )
{
plot4( xOrigin, yOrigin, x, y );
x += 1;
twoXTimesBSquared += twoBSquared;
error += twoXTimesBSquared - bSquared;
if (error >= 0)
{
y -= 1;
twoYTimesASquared -= twoASquared;
error -= twoYTimesASquared;
}
}
}
{    /* Plot the octant from right to top-right */
int x = a;
int y = 0;
int twoXTimesBSquared = x * twoBSquared;
int twoYTimesASquared = 0;
int error = -x* bSquared;  /* b^2 x^2 + a^2 y^2 - a^2 b^2 - b^2x */

while (twoXTimesBSquared > twoYTimesASquared)
{
plot4( xOrigin, yOrigin, x, y );
y += 1;
twoYTimesASquared += twoASquared;
error += twoYTimesASquared - aSquared;
if (error >= 0)
{
x -= 1;
twoXTimesBSquared -= twoBSquared;
error -= twoXTimesBSquared;
}
}
}
}
```

Copyright © 1994, Dr. Dobb's Journal

Back to Table of Contents