Algorithm Alley

Rendering Circles and Ellipses

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 [email protected].

Lining Up

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 .

Going in Circles

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²-x²) (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² > radius² - x², or equivalently, whenever x²+y²-radius²>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)² = x²+2x+1 to keep track of x² and y² 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 y² 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 x²+y²>radius²), 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.

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 x²+y²-r²-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 b²x²+a²y²=a²b². 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 b²x==a²y; see Listing Eight.

Clearly, a few more things could be done to speed up this routine. By adding variables to keep track of b²x and a²y, 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;
            }
        }
    }
}