On circle rasterization algorithms: Bresenham's and Michener's

The idea to write this article came to me when I noticed virtually every page I came across in Google search results presented Michener's algorithm as Bresenham's. While Michener's algorithm is certainly easily and immediately derivable from Bresenham's, it might be useful to some to understand the steps involved.

A few words on Bresenham's algorithm

I will not derive the entire algorithm from scratch here, but rather present a few important points that define the original Bresenham's algorithm. Let's say we want do rasterize a circle of radius R. First and foremost, Bresenham's algorithm solves this problem by generating pixels for 1/4 of the circle, i.e. an arc inside a Cartesian quadrant. Let's assume we start from point (0, R) and proceed upwards and to the left from there until we reach point (R, 0) (Fig. 1).

Fig. 1

At each step, at point (x, y) the algorithm makes a decision which point to add next: point (x - 1, y), point (x - 1, y + 1) or point (x, y + 1) as shown in Fig. 2.

Fig. 2

In order to make that decision, the algorithm uses the following deviation function D(x, y, R)

D(x, y, R) = x² + y² - R²

The above function produces negative values for points that lie strictly inside the circle, and positive values for points strictly outside the circle (and, obviously, zero value for points that lie exactly on the circle). The idea is, naturally, to select the point closest to the circle, i.e. one out of three that minimizes the absolute value of the deviation function.

While standing at point (x, y), the algorithm first focuses its attention at the diagonal point (x - 1, y + 1). It considers value Δ, which represents the deviation of point (x - 1, y + 1) from the actual circle

Δ = D(x - 1, y + 1, R)

and acts in accordance with the following logic

• If Δ is exactly zero, it means point (x - 1, y + 1) lies precisely on our circle. And Bresenham's algorithm immediately decides to proceed to point (x - 1, y + 1).

• If Δ is negative, it means that point (x - 1, y + 1) leads us into the circle's interior (Fig. 3)

  

  Fig. 3

  In such cases Bresenham's algorithm performs an additional analysis in order to determine whether to proceed to point (x - 1, y + 1), or to "counteract" its tendency to lead into the circle by stepping to (x, y + 1) instead (the latter takes us closer to the outside of the circle). To perform such analysis, Bresenham's algorithm calculates an additional value δ

  δ = D(x - 1, y + 1, R) + D(x, y + 1, R)

  • If δ is negative or zero, it means that point (x, y + 1) leads to a smaller absolute error value than point (x - 1, y + 1). The algorithm decides to step to (x, y + 1).

  • If δ is positive, it means that point (x - 1, y + 1) offers a smaller absolute error value. The algorithm steps to (x - 1, y + 1).

  The clever part here is that δ is easily obtainable as δ = 2Δ + 2x - 1

• If Δ is positive, then symmetrical logic applies: it means that point (x - 1, y + 1) leads us out of the circle (Fig. 4)

  

  Fig. 4

  Bresenham's algorithm performs an additional analysis in order to determine whether to proceed to point (x - 1, y + 1), or to bear towards the interior of the circle instead by stepping to (x - 1, y). To perform such analysis, Bresenham's algorithm calculates an additional value δ'

  δ' = D(x - 1, y, R) + D(x - 1, y + 1, R)

  • If δ' is positive or zero, it means that point (x - 1, y) leads to a smaller absolute error value than point (x - 1, y + 1). The algorithm decides to step to (x - 1, y).

  • If δ' is negative, it means that point (x - 1, y + 1) offers a smaller absolute error value. The algorithm steps to (x - 1, y + 1).

  Again, δ' can be easily calculated as δ' = 2Δ - 2y - 1

All that's left to add is that the value of Δ can be maintained incrementally from iteration to iteration (see the code below), i.e. there's no need to calculate D(x - 1, y + 1, R) directly on every iteration, which is what makes Bresenham's algorithm so efficient.

Here's a possible implementation of Bresenham's algorithm

void BresenhamCircle(int cx, int cy, int radius)
{
  int x = radius;
  int y = 0;

  int D = 2 * (1 - radius);
  // 'D' is the error value for '(x - 1, y + 1)' point. I.e.
  //   D = D(x, y) = (x - 1)^2 + (y + 1)^2 - radius^2
  // Substitute initial values of 'x = radius' and 'y = 0' to obtain the above initial 
  // value for 'D'

  while (x >= 0)
  {
    // Draw pixels in each quadrant
    draw_pixel( x + cx,  y + cy);
    draw_pixel( y + cx, -x + cy);
    draw_pixel(-x + cx, -y + cy);
    draw_pixel(-y + cx,  x + cy);

    // We are currently at point '(x, y)', while 'D' is maintained as the error value for 
    // point '(x - 1, y + 1)'

    if (D < 0)
    { // Taking '(x - 1, y + 1)' step leads into the circle. Decide whether we should go
      // to '(x - 1, y + 1)' or to '(x, y + 1)'
      int d = 2 * D + 2 * x - 1;
      // 'd' is sum of error values for '(x - 1, y + 1)' and '(x, y + 1)' 
      //   d = D(x - 1, y + 1) + D(x, y + 1)
      if (d <= 0)
      { // Error value for '(x, y + 1)' is smaller - go in that direction
        ++y;
        D += 2 * y + 1;
        continue;
      }
      // Otherwise, proceed to '(x - 1, y + 1)'
    }
    else if (D > 0)
    { // Taking '(x - 1, y + 1)' step leads out of the circle. Decide whether we should go
      // to '(x - 1, y + 1)' or to '(x - 1, y)'
      int d = 2 * D - 2 * y - 1;
      // 'd' is sum of error values for '(x - 1, y)' and '(x - 1, y + 1)'
      //   d = D(x - 1, y) + D(x - 1, y + 1)
      if (d >= 0)
      { // Error value for '(x - 1, y)' is smaller - go in that direction
        --x;
        D -= 2 * x - 1;
        continue;
      }
      // Otherwise, proceed to '(x - 1, y + 1)'
    }

    // Take a step to '(x - 1, y + 1)'
    --x;
    ++y;
    D += 2 * y - 2 * x + 2;
  }
}

Michener's algorithm

Even a quick look at Bresenham's algorithm immediately reveals obvious optimization opportunities:

• Firstly, it is obvious that as long as we are generating points of the first octant, point (x - 1, y + 1) always lies inside the actual circle (Fig. 3). But once we cross the 45 degree line and start generating points of the second octant, point (x - 1, y + 1) always lies outside the circle (Fig. 4). In other words, for the first half of the arc Δ is always non-positive, while for the second half of the arc it is always non-negative (Fig. 5).

  

  Fig. 5

  For this reason, there's not much point in calculating and analyzing the value of Δ on each iteration. We can simply split the main cycle in two sequential cycles - for the first octant and for the second one - and implement Δ ≤ 0 branches of the code in the first cycle and Δ ≥ 0 branches in the second cycle.

  Of course, in the original algorithm we also need the value of Δ in order to calculate the values of δ or δ'. But these values can also be calculated incrementally, i.e. there's a simple formula that allows us to update the value of δ from the current iteration to obtain the value of δ for the next iteration (see the code below). Same is true for δ', but as the next bullet point states, we don't event need δ'.

  Thus we can get rid of Δ and corresponding branching entirely.

• Secondly, we don't even need to generate the second octant. Once we learned to generate circle arc points in the first octant, we can use reflections across 45 degree bisectors and coordinate axes in order to immediately generate corresponding points in remaining seven octants.

So, all we need is an algorithm for the first octant. And in order to build such algorithm all we need is to maintain the proper value of δ on each iteration. The branching on Δ disappears, the rest of the logic remains unchanged. Bresenham's algorithm optimized in this fashion is also known as Michener's algorithm, although as you can easily see, the changes from the original are relatively superficial.

Here's a possible implementation of Michener's algorithm

void MichenerCircle(int cx, int cy, int radius)
{
  int x = radius;
  int y = 0;

  int d = 3 - 2 * radius;
  // 'd' is sum of error values for '(x - 1, y + 1)' and '(x, y + 1)' 
  //   d = D(x - 1, y + 1) + D(x, y + 1)
  // Substitute initial values of 'x = radius' and 'y = 0' to obtain the above initial 
  // value for 'd'

  while (y <= x)
  {
    // Draw pixels in each octant
    draw_pixel( x + cx,  y + cy);
    draw_pixel( y + cx,  x + cy);
    draw_pixel(-x + cx,  y + cy);
    draw_pixel(-y + cx,  x + cy);
    draw_pixel(-x + cx, -y + cy);
    draw_pixel(-y + cx, -x + cy);
    draw_pixel( x + cx, -y + cy);
    draw_pixel( y + cx, -x + cy);

    if (d <= 0)
    { // Error value for '(x, y + 1)' is smaller - go in that direction
      d += 4 * y + 6;
      ++y;
    }
    else
    { // Otherwise, proceed to '(x - 1, y + 1)' point
      d += 4 * (y - x) + 10;
      --x;
      ++y;
    }
  }
}

The above variant of Michener's algorithm updates the value of δ before updating x and/or y. One can swap these steps and adjust the formulas for δ accordingly, obtaining a slightly different variant of the implementation. Also, in order to reduce the magnitude of δ, one can just divide all formulas for δ by 2 (rounding towards negative infinity). Whatever loss of precision occurs after such division does not affect the sign-switching behavior of δ and, therefore, does not break the algorithm.

The following implementation of the same algorithm includes both of the modifications mentioned above

void MichenerCircle2(int cx, int cy, int radius)
{
  int x = radius;
  int y = 0;

  int d = 1 - radius;

  while (y <= x)
  {
    // Draw pixels in each octant
    draw_pixel( x + cx,  y + cy);
    draw_pixel( y + cx,  x + cy);
    draw_pixel(-x + cx,  y + cy);
    draw_pixel(-y + cx,  x + cy);
    draw_pixel(-x + cx, -y + cy);
    draw_pixel(-y + cx, -x + cy);
    draw_pixel( x + cx, -y + cy);
    draw_pixel( y + cx, -x + cy);

    if (d <= 0)
    {
      ++y;
      d += 2 * y + 1;
    }
    else
    {
      --x;
      ++y;
      d += 2 * (y - x) + 1;
    }
  }
}

This latter version of the implementation seems to be the one most often presented on the Net as "Bresenham's circle algorithm".