*Tim, who has a PhD in mathematics from the University of California, Berkeley, is the author of numerous commercial, shareware, and public-domain programs for graphics and serial communication. Tim can be contacted at kientzle@netcom.com.*

For example, assume you're drawing a line from (0,0) to (3,2). Whenever you increment the x-coordinate, you need to increase the y-coordinate by 2/3. If you let the variable *y* hold the integral part of the y-coordinate and add a new variable to hold the numerator of the fractional part, you have the algorithm in Listing One .

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 .

You need the following simplifications: First, deal with points on the circle as "offsets" from the center of the circle (this is required by the algebra); and second, only compute the pixels in 1/8 of the circle, taking advantage of symmetry. Reflection to the remaining octants of the circle is handled by the routines in Listings Three and Four .

To get a first cut at a circle algorithm, you compute *y* from *sqrt(r ^{2}-x^{2})* (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 *y ^{2} > radius^{2} - x^{2}*, or equivalently, whenever

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} = x^{2}+2x+1* to keep track of

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 *x ^{2}+y^{2}>radius^{2}*), you reduce

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

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 *b ^{2}x^{2}+a^{2}y^{2}=a^{2}b^{2}*. Adjusting the calculation of

Clearly, a few more things could be done to speed up this routine. By adding variables to keep track of *b ^{2}x* and

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.

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 */ } } }

#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; } } }

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 ); }

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

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); } }

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; } } }

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; } } }

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; } } } }

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*