February | 2022 | Possibly Wrong

Introduction

In my recent experiment with string art, I needed Bresenham’s line algorithm to rasterize a line, that is, to compute the set of pixels that most closely approximate a straight line segment. Suppose instead that we want to rasterize a circle, as shown in the figure below.

Bresenham’s algorithm used to rasterize a circle with radius 10. Integer lattice points are assumed to be at the center of each pixel.

There is a similarly named Bresenham’s circle algorithm to do this, that is arguably even simpler than that for drawing a line (code is available on GitHub):

// Call plot(_, _) to draw pixels on circle at origin with given radius.
void draw_circle(int radius, void(*plot)(int, int))
{
    int x = radius;
    int y = 0;
    int error = 3 - 2 * radius;
    while (x >= y)
    {
        plot(x, y); plot(x, -y); plot(-x, y); plot(-x, -y);
        plot(y, x); plot(y, -x); plot(-y, x); plot(-y, -x);
        if (error > 0)
        {
            error -= 4 * (--x);
        }
        error += 4 * (++y) + 2;
    }
}

I have a fond memory of my first encounter with this algorithm as a kid reading Nibble Magazine (“Seven ready-to-type Apple II programs!”). In his article, Brent Iverson [1] had “discovered a neat little circle-generating algorithm in a graphics textbook” … but neglected to mention which textbook, or the name Bresenham, or any other details about how this mysterious snippet of code (which was in BASIC, and in slightly different but equivalent form) worked, despite using only integer addition and bit shift operations. In 1988, without search engines or internet forums, I was on my own.

Today, there are derivations of this algorithm all over the web. However, every one that I have found glosses over what I think is a beautiful property of the algorithm: its fast integer error function is only an approximation of the actual error function, which seems like it could sacrifice some accuracy in choosing which pixels best approximate the desired circle. But it doesn’t; the motivation for this post is to show that (1) this is an issue worth thinking about, since the algorithm’s error function is not monotonic in actual distance from the circle, but (2) despite this, in practice Bresenham’s circle algorithm is still optimal in the sense of always selecting the next pixel that is indeed closest to the circle.

At least, I think. More on this later.

The algorithm

Let’s start by describing how the algorithm works, as a means of introducing some useful notation. Fixing a positive integer input radius $r$ (the reader can verify that the algorithm also works for $r=0$ ), we can focus on rendering only those pixels in the first octant $0 \leq y \leq x$ , using the eightfold dihedral symmetries of the circle to plot the remaining pixels.

Starting at $(x,y)=(r,0)$ , at each iteration we always move up one pixel, and conditionally move (diagonally) left one pixel, choosing from the two corresponding lattice points $(x-1,y+1)$ or $(x,y+1)$ depending on which is closest to the circle. (We always move up, since in the first octant the slope of the tangent to the circle is always at most $-1$ .)

An iteration of Bresenham’s circle algorithm. We have just plotted the pixel at $(x,y)$ , and need to choose which of $(x-1,y+1)$ or $(x,y+1)$ to plot next, ideally minimizing the actual distance from the circle (shown in red).

An iteration of Bresenham’s circle algorithm. We have just plotted the pixel at $(x,y)$ , and need to choose which of $(x-1,y+1)$ or $(x,y+1)$ to plot next, ideally minimizing the actual distance from the circle (shown in red).

Ideally, we would choose the next pixel that minimizes the actual, exact distance from the circle $|c(h,k)|$ , where

$c(h,k)=\sqrt{h^2+k^2}-r$

as indicated by the red line segments in the above figure. However, to eliminate those unpleasant square roots, Bresenham’s algorithm instead seeks to minimize $|b(h,k)|$ , where

$b(h,k)=h^2+k^2-r^2=c(h,k)c^*(h,k)$

$c^*(h,k)=\sqrt{h^2+k^2}+r$

(To foreshadow the key point, note that $b(h,k)$ is not simply the square of $c(h,k)$ .) In either case, it is convenient to define the sum of these signed errors– that is, without the absolute value– for the two candidate pixels as

$e_b(x,y)=b(x-1,y+1)+b(x,y+1)$

$e_c(x,y)=c(x-1,y+1)+c(x,y+1)$

where in the source code above, the error variable maintains the value of $e_b(x,y)$ , which is initialized to $e_b(r,0)=3-2r$ , and at each iteration, depending on the choice of next pixel to render (more on this shortly), can be updated accordingly with simple integer increment, bit shift, and addition operations as shown.

Choosing the closest pixel

To decide which pixel to move to at each iteration, note that both the integer-valued $b(h,k)$ and the exact $c(h,k)$ have the property that they are positive for points $(h,k)$ outside the circle, and negative for points inside the circle. Thus, for example, in the extreme (and in fact impossible) case that both of the candidate next points $(x-1,y+1)$ and $(x,y+1)$ were outside the circle, then $e_b(x,y)$ and $e_c(x,y)$ are positive, indicating that we should move diagonally left to $(x-1,y+1)$ . Similarly, if both candidate points were inside the circle, then the sum of errors is negative, indicating that we should move up to $(x,y+1)$ .

The important question is what happens in between, i.e., in the typical situation shown in the above figure, where the two candidate points straddle the circle, so to speak. Again, ideally we would use the sign of the exact sum of errors $e_c(x,y)$ , in the same way as described above, moving diagonally left if $e_c(x,y)>0$ , or up if $e_c(x,y)<0$ . Does Bresenham’s algorithm, using the sign of the approximated integer-valued $e_b(x,y)$ instead, yield the same behavior?

To see why this is a question worth asking, suppose for example that we want to draw a circle with radius $r=4$ , and that for some reason we are asked which of the points $(2,2)$ or $(4,3)$ is closer to the circle, as shown in the following figure.

Example choosing which of points $(2,2)$ or $(4,3)$ is closer to the circle with radius 4.

Example choosing which of points $(2,2)$ or $(4,3)$ is closer to the circle with radius 4.

The reader can verify that $(4,3)$ is closer, with $|c(4,3)|<|c(2,2)|$ and thus $c(2,2)+c(4,3)<0$ … but $|b(4,3)|>|b(2,2)|$ and thus $b(2,2)+b(4,3)>0$ . That is, the integer error approximation used in Bresenham’s algorithm is not monotonic in the actual distance of points from the circle.

The interesting good news is that this seemingly potential flaw, isn’t. That is, we are now in a position to conjecture the following:

Conjecture: At each iteration of Bresenham’s circle algorithm, the sum-of-error functions $e_c(x,y)$ and $e_b(x,y)$ always have the same sign. That is, the algorithm always makes the correct choice of which of the two candidate next pixels is closest to the circle.

Incomplete proof: There are a few key observations that will be useful here. First, at each iteration, $c(x,y+1)>c(x-1,y+1)$ and $c^*(x,y+1)>c^*(x-1,y+1)>0$ . We can see these directly from the definitions, noting that $r \geq x>0$ . (This is why we constrained the radius to be positive, treating $r=0$ as a separate special case.)

Also, note that the exact sum of errors $e_c(x,y)$ is always irrational, and thus never zero. That is, we should never be indifferent to the choice of next pixel. The approximated $e_b(x,y)$ is never zero, either, as can be seen by inspecting the source code: the value of error is initialized to an odd number, and is always incremented or decremented by even amounts. Thus, it suffices to show that $e_c(x,y)$ and $e_b(x,y)$ are either both positive, or both negative.

So, there are two cases. First, suppose that $e_c(x,y)>0$ . Then $c(x,y+1)>0$ , since

$2c(x,y+1)>c(x,y+1)+c(x-1,y+1)=e_c(x,y)>0$

Thus,

$e_b(x,y)=b(x-1,y+1)+b(x,y+1)$

$= c(x-1,y+1)c^*(x-1,y+1)+c(x,y+1)c^*(x,y+1)$

$> c(x-1,y+1)c^*(x-1,y+1)+c(x,y+1)c^*(x-1,y+1)$

$= e_c(x,y)c^*(x-1,y+1)>0$ .

Thus, when we are supposed to move diagonally left, we do.

The second case, showing that $e_c(x,y)<0$ implies $e_b(x,y)<0$ , is why this is only conjecture and not a theorem. Following is as far as I have been able to proceed:

Similar to the first case, note that $c(x-1,y+1)<0$ , since

$2c(x-1,y+1)<c(x,y+1)+c(x-1,y+1)=e_c(x,y)<0$

Consider again the expanded form of

$e_b(x,y)=b(x-1,y+1)+b(x,y+1)$

$= c(x-1,y+1)c^*(x-1,y+1)+c(x,y+1)c^*(x,y+1)$

If $c(x,y+1) \leq 0$ , then we can see $e_b(x,y)<0$ directly. Otherwise, if $c(x,y+1)>0$ … and here is where I get stuck.

Conclusion

Well, that was disappointing.

My glass is half full, though; I think there are three possible resolutions of this problem, any of which would be an interesting outcome. First, and most likely, I may have simply missed something obvious. Perhaps a reader can help take the above attempt at a proof to the finish line.

Alternatively, maybe the conjecture is actually false, and there is in fact a counterexample– which given the partial proof above, would necessarily be a situation where Bresenham’s algorithm is supposed to move up ( $e_c(x,y)<0$ ), but moves diagonally left instead ( $e_b(x,y)>0$ ). In realistic practice, this isn’t a concern: I have verified numerically that the algorithm works correctly for every radius up to $r=16384$ — that is, much larger than any pixel-perfect circle that would fit on the highest resolution screens currently available. But even from a theoretical perspective, if there is a counterexample, it would be interesting that a minimal one must be so surprisingly large.

Or finally, maybe the conjecture is indeed true, but the remainder of the proof of the second “negative” case is actually really hard. There is some intuition behind this possibility as well: in this part of the proof, we are trying to bound a sum of square roots, which is one of my favorite examples of “problems that seem like they have no business being as difficult as they are.” (The Euclidean traveling salesman is perhaps the most recognizable decision problem with this wrinkle: it’s often claimed to be NP-complete, but it’s “only” NP-hard; we don’t even know if it’s in NP, since it involves a similar comparison with a sum of square roots.) This potentially affects the integrity of my brute-force numerical verification above as well: it’s technically possible that I missed a counterexample by performing those thousands of sign checks with “only” 50-digit precision.

Edit 2022-03-21: Thanks to commenters Andrey Khalyavin for the initial proof of the contrapositive of the missing half above, which I didn’t fully grasp, and to Linus Hamilton for patiently following up with a more detailed version, linked here, that helped me to see how it works. Their work shows that Bresenham’s algorithm is indeed optimal, always selecting the unique next pixel that is closest to the desired circle.

Reference:

Iverson, Brent, Hi-Res Circle Generator, Nibble Magazine, 9(1) January 1988, 68-71

Possibly Wrong

On science, mathematics, and computing

Monthly Archives: February 2022

Bresenham’s circles are (always?) optimal