To rasterize a circle using the mid-point algoritm (see Foley, van Dam) using only integers and additions, we should first exploit the fact that circles are 8-way symmetric. We only need to plot any one quadrant of the circle, and mirror it 8-way, to plot the entire circle. Traditionally, the second quadrant is generated, starting with X=0 and Y=r, assuming the circle is plotted with the center in the origin. Any offset can easily be added in the actual plotting function, and need not be considered in rasterization.

Rasterizing rules matter: if pixels are located between coordinate points on the screen, then it's hard to decide which pixel to include in the circle: the one with center at Y+0.5, or the one with center at Y-0.5? To avoid this problem, we will assume that pixels are centered on integer coordinates.

The basic algorithm is the same as the mid-point line drawing algorithm: after plotting pixel (X,Y), there are only two possible pixels to plot next. Both are located at X coordinate X+1, and they are either at the same Y coordinate as the current pixel, or at one lower Y coordinate -- choosing any other pixel would either create gaps, or make the outline have extra pixels, making the circle look blocky.

Assuming your current raster pos is (X,Y), with centers of pixels being integral coordinates, the decision to make is whether to plot the next pixel at (X+1,Y) or at (X+1,Y-1).

We can make this determination by looking at where the circle is, algebraically, at position X+1, and choose Y if it's closer to Y, and choose Y-1 if it's closer to Y-1. These two points are shown in gray in the illustration.

The circle equation is R*R == X*X + Y*Y. From that we can derive the circle at point X+1: Yc = sqrt( R*R - (X+1)*(X+1) ). We can compare this value to the mid-point:

Yc > (Y-0.5)

• Yes: Keep Y
• No: Decrement Y

As comparision of positive integers remains invariant under raising to any positive power, we might as well compare:

Yc * Yc > (Y-0.5)*(Y-0.5)

Note that, assuming we can take the 0.5 term in effect, this still leads to precise rasterization, as we're only comparing the Y values we derive, not the R values (which, if we considered them at the point X for the Y values Y, Yc and Y-1, would suffer from skew due to the squaring).

We can now note that Yc*Yc is exactly R*R (which is constant) minus (X+1)*(X+1) which can be incrementally computed from the previous value of X, using only addition (it's prev-value + X + X + 1).

We can also note that, IF we change the Y value, then we'll change Y incrementally: Ynew = Yprev - Y - Y + 1. "YPrev" is the squared value of Y; Ynew is the new squared value of the new Y; Y is the old Y.

So what do we do about the (Y-0.5) term? It's not very integral. As we always square it, we'll notice that (Y-0.5) squared equals Y-squared minus Y plus one-quarter. This may lead us to believe that we need to scale the values in the algorithm by a factor of 4, but as we only subtract one-half, we can get away with a scale of 2.

Putting it all together, we get the following implementation:

  void plot8( int x, int y ) {
    plot1( x, y ); plot1( -x, y );
    plot1( x, -y ); plot1( -x, -y );
    plot1( y, x ); plot1( -y, x );
    plot1( y, -x ); plot1( -y, -x );
  }

  void originCircle( int radius ) {
    int rs2 = radius*radius*4; /* this could be folded into ycs2 */
    int xs2 = 0;
    int ys2m1 = rs2 - 2*radius + 1;
    int x = 0;
    int y = radius;
    int ycs2;
    plot8( x, y );
    while( x <= y ) {
      /* advance to the right */
      xs2 = xs2 + 8*x + 4;
      ++x;
      /* calculate new Yc */
      ycs2 = rs2 - xs2;
      if( ycs2 < ys2m1 ) {
        ys2m1 = ys2m1 - 8*y + 4;
        --y;
      }
      plot8( x, y );
    }
  }