Bresenham's Circle Drawing Algorithm

Once upon a time I was given the following problem for a technical programming interview:

Write a function draw_circle(r) that draws a circle with radius r. Use the given method draw_pixel(x, y) which takes a 2-dimensional point (x, y) and colors it in on the computer screen.

For the solution, you can either collect all pixels (tuples) of \(x\) and \(y\) coordinate pairs, or just call draw_pixel() on them during the “search” for those pixels that must be filled in.

This post goes over several solutions, ultimately arriving at Bresenham’s algorithm. The content of this post is merely a distillation of Section 3.3 from the book “Computer Graphics: Principles and Practice (1996)”.¹ The authors of the book state that their implementation results in code “essentially the same as that specified in patent 4,371,933 [a.k.a. Bresenham’s algorithm].”²

I’ve gone all out and converted the “reference” implementations found in the book and translated them into Rust and Python. The Python was written first, and I used a text-based drawing system to test the correctness. However I became dissatisfied with the non-square “aspect ratio” of most monospaced fonts out there, which distorted the circles to look more like ellipses. To fix this, I decided to port the Python code to Rust, and then target WASM so that I can use it to draw on the HTML5 <canvas> elements (and to eliminate the “aspect ratio” problem). All of the drawings in this document are powered by the Rust code.

Constraints

Drawable canvas

Before we start, let’s define the drawable surface (canvas) of pixels for this problem. The pixels are arranged in a 2-dimensional grid. The important thing here is the grid or coordinate system, with the pixel at the center of the grid having the traditional (0, 0) Cartesian coordinate.

Below is a sample grid to give you a sense of what this will look like. There is a central (0, 0) origin pixel, and 15 pixels to the north, south, east, and west, and everything in-between. Pixels that lie on interesting points of symmetry are highlighted in green.

Mathematical definitions

The exact definition of a circle (given infinite precision, as on the traditional Cartesian plane) centered at the origin is

\[ \begin{equation} \label{eq:circle} x^2 + y^2 = r^2. \end{equation} \]

This resembles the Pythagorean Theorem

\[ a^2 + b^2 = c^2, \]

for any right-angled triangle with sides \(a\) and \(b\) and hypotenuse \(c\). The resemblance is not a coincidence, because an infinite number of such triangles exists within the top right quadrant of the plane (that is, Quadrant I³, or the part of the plane such that \(x \geq 0\) and \(y \geq 0\)); in Quadrant I, for all points \((x,y)\) that make up this portion (or arc) of the circle, their radii is the same as the hypotenuses of these triangles (whose sides are \(x\) and \(y\)). Later in this post, this will become relevant again when we discuss Pythagorean Triples.

Anyway, solving for \(y\) in Equation \(\ref{eq:circle}\) gives

\[ \begin{equation} \label{eq:circle-y} y = \pm\sqrt{r^2 - x^2} \end{equation} \]

to get 2 functions for the top-half and bottom-half of the circle (that’s what the \(\pm\) symbol means). Consider the function \(y = x\). This function has slope 1 and is a diagonal line where all values of \(x = y\). Now consider how this line intersects the quarter-arc of the circle in Quadrant I. This intersection point evenly divides the arc into 2 halves, and is where

\[ x = y = \tfrac{r}{\sqrt{2}}, \]

or simply the point

\[ \begin{equation} (\tfrac{r}{\sqrt{2}}, \tfrac{r}{\sqrt{2}}). \end{equation} \]

This is because if \(x = y\), then Equation \(\ref{eq:circle}\) becomes

\[ \begin{align} x^2 + y^2 &= r^2 \\ x^2 + x^2 &= r^2 \\ 2x^2 &= r^2 \\ \tfrac{2x^2}{2} &= \tfrac{r^2}{2} \\ x^2 &= \tfrac{r^2}{2} \\ \sqrt{x^2} &= \tfrac{\sqrt{r^2}}{\sqrt{2}} \\ x &= \tfrac{r}{\sqrt{2}}. \label{eq:arc-intersection} \end{align} \]

This is not that interesting for purposes of the algorithms in this post, but is something that is glossed over in the book.

Symmetry

Because of symmetry, we can mirror the solution \((x,y)\) pairs we get in Quadrant I into the other quadrants. This gives us 4-way symmetry because there are 4 quadrants.

def mirror_points_4(x, y):
    """ Return 4-way symmetry of points. """
    return [( x,  y),
            (-x,  y),
            ( x, -y),
            (-x, -y)]

Note, however, that there is actually 8-way symmetry at hand because (1) we can swap \(x\) and \(y\), and (2) because of the way we can distribute the negative sign:

def mirror_points_8(x, y):
    """ Return 8-way symmetry of points. """
    return [( x,  y),
            ( y,  x),
            (-x,  y),
            (-y,  x),
            ( x, -y),
            ( y, -x),
            (-x, -y),
            (-y, -x)]

Fun fact: the exact point at which \(x\) and \(y\) get “swapped” in Quadrant I is when \(x = y = \tfrac{r}{\sqrt{2}}\) (Equation \(\ref{eq:arc-intersection}\)).

Naive solutions

When in doubt, brute force is always a great answer, because at least it gets you started on something that works given enough time and/or memory.⁴ Because we already have clear mathematical definitions, we can just translate them (albeit mechanically) to code.

def get_circle_points_naive_4(r):
    """ Draw a circle by pairing up each Y value with an X value that lie on a
    circle with radius 'r'. This has a bug because some Y values get skipped.
    Can you see why?
    """
    points = []
    for x in range(r + 1):
        # isqrt() gets the integer square root.
        y = isqrt((r * r) - (x * x))
        points.extend(mirror_points_4(x, y))
    return points

get_circle_points_naive_4() is the simplest translation, although there is a bug, which is obvious when we visualize it (in this case, for \(r = 15\)):

The get_circle_points_naive_4() is based on Equation \(\ref{eq:circle-y}\). We iterate \(x\) from \(0\) to \(r\) ⁵, and at each \(x\) try to find the best value for \(y\). The problem is that we’re only solving for 1 \(y\) value for every \(x\) value we increment by. As we get near the left and right sides of the circle, we need to calculate more than just 1 \(y\) value for every \(x\).⁶.

The get_circle_points_naive_8() function gets around this \(y\)-skip bug by invoking 8-way symmetry instead:

def get_circle_points_naive_8(r):
    """ Better than get_circle_points_naive_4, but wastes CPU cycles because
    the 8-way symmetry overcorrects and we draw some pixels more than once.
    """
    points = []
    for x in range(r + 1):
        y = isqrt((r * r) - (x * x))
        points.extend(mirror_points_8(x, y))
    return points

However the downside is that it results in multiple points that will be drawn 2 times, wasting CPU cycles.⁷ To be more precise, all points around the gappy area in Quadrant I are redundant because that part of the arc is already mirrored nicely by the contiguous points from \(x = 0\) to \(x = y\).

The get_circle_points_naive_8_faster() function avoids drawing the gappy areas by just breaking the loop when \(x > y\), but is otherwise the same:

def get_circle_points_naive_8_faster(r):
    """ Slightly faster than get_circle_points_naive_8, because of the break
    condition at the middle of the arc. However this is still inefficient due
    to the square root calculation with `isqrt()`.
    """
    points = []
    for x in range(r + 1):
        y = isqrt((r * r) - (x * x))
        # When we cross the middle of the arc, stop, because we're already
        # invoking 8-way symmetry.
        if x > y:
            break
        points.extend(mirror_points_8(x, y))
    return points

This is the best we can do with the simple mathematical translations to code. Note that in all of these implementations we are still forced to calculate square roots in every iteration, which is certainly suboptimal.

Bresenham’s Algorithm

This as also known as the “Midpoint Circle Algorithm,” where the name “midpoint” comes from the mathematical calculations that are done by considering the midpoint between pixels. The gist of the algorithm is that instead of using Equation \(\ref{eq:circle-y}\) to calculate \(y\) for every \(x\), instead you try to move along the arc of the circle, pixel-to-pixel, staying as close as possible to the true arc:

1. Start out from the top of the circle (color in pixel \((0, r)\)). Note that because of symmetry, we could start out from \((0, -r)\), \((r, 0)\), or even \((-r, 0)\) as Bresenham did in his paper.⁸
2. Move right (east (E)) or down-right (southeast (SE)), whichever is closer to the circle.
3. Stop when \(x = y\) (just like in get_circle_points_naive_8_faster()).

The hard part is Step 2, where we just need to figure out which direction to move (E or SE) from the current pixel. The brute force way here is to just calculate the distance away from the center of the circle for the E and SE pixels (using Euclidean distance, which is just a variation of Equation \(\ref{eq:circle}\) or the Pythagorean Theorem), and just choose the pixel that is closest to the arc of the circle. This makes sense, but with the power of mathematics, we can do better.

Inside, on, or outside the circle?

In order to figure out whether some point \((x, y)\) is inside, on, or outside of the circle depends on the definition of the circle from Equation \(\ref{eq:circle}\). We can tweak it in terms of any \((x, y)\) pair:

\[ \begin{equation} \label{eq:error-margin} F(x,y) = x^2 + y^2 - r^2 = \text{distance from true circle line}. \end{equation} \]

Note that if \(F(x,y) = 0\), then the point \((x, y)\) is exactly on the circle. If \(F(x,y) > 0\), then the point is outside of the circle, and if \(F(x,y) < 0\) then the point is inside of it. In other words, given any point \((x, y)\), \(F(x, y)\) is the distance from the true circle line.

Choosing between E or SE

Let’s remind ourselves that we’ll always be moving E or SE. One critical (pragmatic) property here is that we’re dealing with a pixel grid with integer increments. There is a very high chance that neither the E or SE pixels we’re moving to is exactly on the circle. This is because the only time that the point \((x,y)\) will exactly be on the line of the circle is if the \(x\), \(y\), and \(r\) values (as integers) form a so-called Pythagorean Triple. For \(r < 100\), there are only 50 such triples:

( 3, 4, 5)  (18,24,30)  (24,45,51)  (16,63,65)  (51,68,85)
( 6, 8,10)  (16,30,34)  (20,48,52)  (32,60,68)  (40,75,85)
( 5,12,13)  (21,28,35)  (28,45,53)  (42,56,70)  (36,77,85)
( 9,12,15)  (12,35,37)  (33,44,55)  (48,55,73)  (13,84,85)
( 8,15,17)  (15,36,39)  (40,42,58)  (24,70,74)  (60,63,87)
(12,16,20)  (24,32,40)  (36,48,60)  (45,60,75)  (39,80,89)
(15,20,25)  ( 9,40,41)  (11,60,61)  (21,72,75)  (54,72,90)
( 7,24,25)  (27,36,45)  (39,52,65)  (30,72,78)  (35,84,91)
(10,24,26)  (30,40,50)  (33,56,65)  (48,64,80)  (57,76,95)
(20,21,29)  (14,48,50)  (25,60,65)  (18,80,82)  (65,72,97)

In other words, for all practical purposes, there will always be some error and we’ll always be outside or inside the circle and never directly on it. It’s sort of like driving a car and trying to stay within your designated lane: if you think you’re moving too much to the right, you turn your wheel left to stay “within” the lane (or some acceptable amount within the lane), and vice versa.

The idea is the same for moving along the circle: if we think we’re moving too far outside the circle, we try to move into it. On the other hand, if we think we’re moving into the circle, we move out of it. And so imagine yourself standing on point \((0, r)\), our starting point. The line of the circle is our “lane” we want to stay “on” as much as possible. Choosing to go E is the same as turning “left”. Choosing to go SE is the same as turning “right”. Using this metaphor, if we were not to turn at all (go “straight”), we would be heading to the virtual “in-between” pixel between E and SE, the midpoint between them.

And so here’s the basic idea behind choosing E or SE:

1. If going “straight” would mean going into the circle (i.e., we’re currently veering too much to the right!), we course-correct by turning left (E).
2. Conversely, if going “straight” would mean going outside the circle (i.e., we’re currently veering too much to the left), we course-correct by turning right (SE).
3. Lastly, if going “straight” would mean staying exactly on the circle (we hit a Pythagorean Triple), we turn SE (from an engineering perspective it doesn’t really matter which way we turn in this case, as both E and SE result in some amount of error — although see “Final tweaks” below for a note on aesthetics).

Let’s convert this idea into pseudocode:

Let M be the midpoint (going "straight").

Then, F(M) tells us what direction we're headed relative to the true circle line.

If F(M) is < 0, we're moving "into" the circle (veering right), so turn left by moving E.

Otherwise move SE.

Note that we only have to calculate \(F(...)\) for the midpoint \(M\). Isn’t this cool? It is much better than calculating \(F(E)\) and \(F(SE)\) and having to compare them!

# This F() function is the same as the mathematical F(...) function
# discussed above (Equation 11).
def F(x, y, r):
    return (x * x) + (y * y) - (r * r)

def get_circle_points_bresenham_WIP1(r):
    points = []
    x = 0
    y = r
    # Calculate F(M) for the very first time. That is, if we were to go
    # "straight" from (0, r), would we be inside or outside the circle?
    xE, yE = (1, r)
    xSE, ySE = (1, r - 1)
    xM, yM = (1, r - 0.5)
    F_M = F(xM, yM, r)
    points.extend(mirror_points_8(x, y))
    while x < y:
        # If going straight would go "into" the circle (too much to the
        # right), try to move out of it by turning left by moving E.
        if F_M < 0:
            x += 1
            F_M = F(x, y, r)
        # Otherwise move SE.
        else:
            x += 1
            y -= 1
            F_M = F(x, y, r)
        points.extend(mirror_points_8(x, y))
    return points

We can refactor the above slightly. We can simplify the initial calculation of F_M to avoid calling F(), and also move out some of the redundant bits. The very first midpoint we have to consider is \((1, r - \tfrac{1}{2})\); plugging this into \(F()\) gets us

\[ \begin{align} F(1, r - \tfrac{1}{2}) &= 1^2 + (r - \tfrac{1}{2})^2 - r^2 \\ &= 1 + (r^2 - r + \tfrac{1}{4}) - r^2 \\ &= 1 + r^2 - r^2 - r + \tfrac{1}{4} \\ &= 1 - r + \tfrac{1}{4} \\ &= \tfrac{5}{4} - r. \end{align} \]

With that said, we can get this:

def get_circle_points_bresenham_WIP2(r):
    points = []
    x = 0
    y = r
    F_M = 5/4 - r
    points.extend(mirror_points_8(x, y))
    while x < y:
        # If going straight would go "into" the circle (too much to the
        # right), try to move out of it by turning left by moving E.
        if F_M < 0:
            pass
        # Otherwise move SE.
        else:
            y -= 1
        x += 1
        F_M = F(x, y, r)
        points.extend(mirror_points_8(x, y))
    return points

The annoying bit is the call to F(). Surprisingly, the call to F() can be elimitated entirely, because we can calculate it once, and then merely adjust it thereafter.

Calculate once, adjust thereafter

We can just calculate \(F(x,y)\) once when we start out at \((0, r)\), and then just adjust it depending on whether we move E or SE. The key is that this “adjustment” computation is cheaper than calculating the full \(F(x,y)\) distance function all over again.

Let \(M\) be the midpoint \((x + 1, y - \tfrac{1}{2})\) between the E \((x + 1, y)\) and SE \((x + 1, y - 1)\) pixels. Then \(F(M)\) is the result of going “straight” and tells us the direction we’re veering off from the circle line:

\[ \begin{equation} F(M) = F(x + 1, y - \tfrac{1}{2}) = (x + 1)^2 + (y - \tfrac{1}{2})^2 - r^2. \end{equation} \]

The values for \(x\) and \(y\) are unknown, however they change in only 2 possible ways — by moving E or SE!

If we move E, then \(M\) will change from \((x + 1, y - \tfrac{1}{2})\) to \((x + 2, y - \tfrac{1}{2})\) because we add 1 to \(x\) to move 1 pixel east; the new value of \(F(M)\) at this pixel, which we can call \(F(M_E)\), will then be:

\[ \begin{equation} F(M_{E}) = F(x + 2, y - \tfrac{1}{2}) = (x + 2)^2 + (y - \tfrac{1}{2})^2 - r^2. \end{equation} \]

Now we can take the difference between these two full calculations. That is, if we were to move E, how would \(F(M)\) change? Simple, we just look at the change in \(x\) (\(\Delta_{x}\)) (we don’t care about the change in \(y\) or \(r\), because they stay constant in this case).

\[ \begin{align} \Delta_{E} &= F(M_{E}) - F(M) \\ &= [(x + 2)^2 + (y - \tfrac{1}{2})^2 - r^2] - [(x + 1)^2 + (y - \tfrac{1}{2})^2 - r^2] \\ &= \Delta_{x} \\ &= (x + 2)^2 - (x + 1)^2 \label{eq:de1} \\ &= (x^2 + 4x + 4) - (x^2 + 2x + 1) \\ &= x^2 + 4x + 4 - x^2 - 2x - 1 \\ &= x^2 - x^2 + 4x - 2x + 4 - 1 \\ &= 2x + 3. \label{eq:de2} \end{align} \]

So \(F(M)\) will change by \(2x + 3\) if we move E. So at any given point, if we move E, \(F(M)\) will always change by \(2x + 3\).

How about for moving SE? If we move SE, the new value of \(M\) will change from \((x + 1, y - \tfrac{1}{2})\) to \((x + 2, y - \tfrac{3}{2})\) because we add 1 to \(x\) and subtract 1 from \(y\) to move 1 pixel southeast; the new value of \(F(M)\) for this case, which we call \(F(M_{SE})\), will then be:

\[ \begin{equation} F(M_{SE}) = F(x + 2, y - \tfrac{3}{2}) = (x + 2)^2 + (y - \tfrac{3}{2})^2 - r^2. \end{equation} \]

We can do the same difference analysis here, but with the addition that we have to consider the change in \(y\) (\(\Delta_{y}\)) as well (because of the 1 we subtracted from \(y\)):

\[ \begin{align} \Delta_{SE} &= F(M_{SE}) - F(M) \\ &= [(x + 2)^2 + (y - \tfrac{3}{2})^2 - r^2] - [(x + 1)^2 + (y - \tfrac{1}{2})^2 - r^2] \\ &= \Delta_{x} + \Delta_{y} \\ &= [(x + 2)^2 - (x + 1)^2] + [(y - \tfrac{3}{2})^2 - (y - \tfrac{1}{2})^2] \\ &= (2x + 3) + [(y^2 - \tfrac{6y}{2} + \tfrac{9}{4}) - (y^2 - y + \tfrac{1}{4})] \\ &= (2x + 3) + (y^2 - 3y + \tfrac{9}{4} - y^2 + y - \tfrac{1}{4}) \\ &= (2x + 3) + (y^2 - y^2 - 3y + y + \tfrac{9}{4} - \tfrac{1}{4}) \\ &= (2x + 3) + (- 2y + \tfrac{8}{4}) \\ &= (2x + 3) + (-2y + 2) \\ &= 2x + 3 - 2y + 2 \\ &= 2x - 2y + 5 \\ &= 2(x - y) + 5. \label{eq:se1} \end{align} \]

And so when moving SE, the new \(F(M)\) must change by \(2(x - y) + 5\).

Now we have all the pieces to derive the complete algorithm!

def get_circle_points_bresenham_float_ese(r):
    """ Draw a circle using a floating point variable, F_M. Draw by moving E or
    SE."""
    points = []
    x = 0
    y = r
    # F_M is a float.
    F_M = 5 / 4 - r
    points.extend(mirror_points_8(x, y))
    while x < y:
        if F_M < 0:
            F_M += 2.0 * x + 3.0
        else:
            F_M += 2.0 * (x - y) + 5.0
            y -= 1
        x += 1
        points.extend(mirror_points_8(x, y))
    return points

Integer-only optimization

The initial value of F_M (\(F(M)\)) is \(\tfrac{5}{4} - r\). Notice how this is the only place where we have to perform division in the whole algorithm. We can avoid this initial division (and subsequent floating point arithmetic) by initializing it to \(1 - r\) instead, which is a difference of \(\tfrac{1}{4}\) vs the original.

Because we tweaked the initialization by \(\tfrac{1}{4}\), we have to do the same for all comparisons of \(F(M)\) moving forward. That is, the comparison \(F(M) < 0\) actually becomes \(F(M) < -\tfrac{1}{4}\). However, this fractional comparison is unnecessary because we only deal with integer increments and decrements in the rest of the code, so we can just keep the same \(F(M) < 0\) as before. In other words, our algorithm only cares about whole numbers, so worrying about this extra \(\tfrac{1}{4}\) difference is meaningless.

def get_circle_points_bresenham_integer_ese(r):
    """ Like draw_circle_bresenham_float_ese, but F_M is an integer variable.
    """
    points = []
    x = 0
    y = r
    # F_M is an integer!
    F_M = 1 - r
    points.extend(mirror_points_8(x, y))
    while x < y:
        if F_M < 0:
            # We can use a bit-shift safely because 2*n is the same as n << 1
            # in binary, and also because F_M is an integer.
            F_M += (x << 1) + 3
        else:
            F_M += ((x - y) << 1) + 5
            y -= 1
        x += 1
        points.extend(mirror_points_8(x, y))
    return points

Second-order differences

There is a final optimization we can do.⁹ In the “Calculate once, adjust thereafter” section we avoided calculating \(F(M)\) from scratch on every iteration. We can do the same thing for the differences themselves!

That is, we can avoid calculating \(\Delta_{E} = (2x + 3)\) and \(\Delta_{SE} = 2(x - y) + 5\) on every iteration, and instead just calculate them once and make adjustments to them, just like we did earlier for \(F(M)\).

Let’s first consider how \(\Delta_{E} = 2x + 3\) changes. First, we initialize \(\Delta_{E}\) by plugging in \((0, r)\) into Equation \(\ref{eq:de2}\), our starting point. Because there is no \(y\) variable in here, we get an initial value of

\[ \begin{equation} \label{eq:de-2ord-initial} 2(0) + 3 = 3. \end{equation} \]

If we go E, \(\Delta_{E}\) changes like this: \[ \begin{align} \Delta_{E_{new}} = \Delta_{E_(x+1,y)} - \Delta_{E_(x,y)} &= [2(x+1) + 3] - (2x + 3) \label{eq:de-2ord-e} \\ &= 2x + 2 + 3 - 2x - 3 \\ &= 2x - 2x + 3 - 3 + 2 \\ &= 2. \label{eq:e2ord} \end{align} \]

If we go SE, \(\Delta_{E}\) changes in the exact same way, because even though our new point is at \((x+1, y-1)\), there is no \(y\) in \(\Delta_{E} = 2x + 3\), so it doesn’t matter and \(\Delta_{E_{new}} = 2\) again.

Now let’s consider how \(\Delta_{SE}\) changes. For the initial value, we again plug in \((0, r)\) into \(2(x-y) + 5\), to get

\[ \begin{equation} \label{eq:dse-2ord-initial} 2(0-r) + 5 = -2r + 5. \end{equation} \]

If we go E, \(\Delta_{SE}\) changes like this:

\[ \begin{align} \Delta_{SE_{new}} = \Delta_{SE_(x+1,y)} - \Delta_{SE_(x,y)} &= [2((x + 1)-y) + 5] - [2(x - y) + 5] \label{eq:dse-2ord-e} \\ &= (2x + 2 - 2y + 5) - (2x - 2y + 5) \\ &= 2x - 2y + 7 - 2x + 2y - 5 \\ &= 2x - 2x + 2y - 2y + 7 - 5 \\ &= 2. \label{eq:se2ord1} \end{align} \]

If we go SE, \(\Delta_{SE}\) changes like this:

\[ \begin{align} \Delta_{SE_{new}} = \Delta_{SE_(x+1,y-1)} - \Delta_{SE_(x,y)} &= [2((x + 1)-(y - 1)) + 5] - [2(x - y) + 5] \label{eq:dse-2ord-se} \\ &= [2(x + 1 - y + 1) + 5] - (2x - 2y + 5) \\ &= (2x + 2 - 2y + 2 + 5) - 2x + 2y - 5 \\ &= 2x- 2x + 2y - 2y + 5 - 5 + 2 + 2 \\ &= 2 + 2 \\ &= 4. \label{eq:se2ord2} \end{align} \]

The code should then look like this:

def get_circle_points_bresenham_2order(r):
    points = []
    x = 0
    y = r
    F_M = 1 - r
    d_e = 3 # Equation 40
    d_se = -(2 * r) + 5 # Equation 45
    points.extend(mirror_points_8(x, y))
    while x < y:
        if F_M < 0:
            F_M += d_e
            d_e += 2  # Equation 44
            d_se += 2 # Equation 50
        else:
            F_M += d_se
            d_e += 2  # Equation 44
            d_se += 4 # Equation 56
            y -= 1
        x += 1
        points.extend(mirror_points_8(x, y))
    return points

With a little refactoring, we can arrive at a slightly simpler version:

def get_circle_points_bresenham_integer_ese_2order(r):
    """ Like draw_circle_bresenham_integer_ese, but use 2nd-order differences
    to remove multiplication from the inner loop. """
    points = []
    x = 0
    y = r
    F_M = 1 - r
    # Initial value for (0,r) for 2x + 3 = 0x + 3 = 3.
    d_e = 3
    # Initial value for (0,r) for 2(x - y) + 5 = 0 - 2y + 5 = -2y + 5.
    d_se = -(r << 1) + 5
    points.extend(mirror_points_8(x, y))
    while x < y:
        if F_M < 0:
            F_M += d_e
        else:
            F_M += d_se
            # Increment d_se by 2 (total 4) if we go southeast.
            d_se += 2
            y -= 1
        # Always increment d_e and d_se by 2!
        d_e += 2
        d_se += 2
        x += 1
        points.extend(mirror_points_8(x, y))
    return points

The “purist” in me felt that the decrementing of \(y\) stood out like a sore thumb, and so I created a tweaked version that moves E and NE, starting out from \((0, -r)\) instead. The mathematical techniques are the same, and due to symmetry the behavior of the algorithm does not change.

def get_circle_points_bresenham_integer_ene_2order(r):
    """ Like draw_circle_bresenham_integer_ene, but start from (0, -r) and move
    E or NE. Notice how we only need the addition instruction in the while loop
    (y is incremented, not decremented). """
    points = []
    x = 0
    y = -r
    F_M = 1 - r
    # Initial value for (0,-r) for 2x + 3 = 0x + 3 = 3.
    d_e = 3
    # Initial value for (0,-r) for 2(x + y) + 5 = 0 - 2y + 5 = -2y + 5.
    d_ne = -(r << 1) + 5
    points.extend(mirror_points_8(x, y))
    while x < -y:
        if F_M < 0:
            F_M += d_e
        else:
            F_M += d_ne
            d_ne += 2
            y += 1
        d_e += 2
        d_ne += 2
        x += 1
        points.extend(mirror_points_8(x, y))
    return points

Here are a couple drawings using Bresenham’s algorithm. This one is for \(r = 15\):

And for \(r = 60\):

Comparisons vs naive algorithm

Here are some side-by-side comparisons for \(0 \leq r \leq 10\).

Final tweaks

It has been kindly pointed out that the naive algorithm is aesthetically more pleasing if the calculations involving \(r\) is done with \(r + \tfrac{1}{2}\) instead of just \(r\) itself, like this:

def get_circle_points_naive_8_faster_tweaked_radius(r):
    """ This is much closer to Bresenham's algorithm aesthetically, by simply
    using 'r + 0.5' for the square root calculation instead of 'r' directly.
    """
    points = []
    # In the square root calculation, we just use (r + 0.5) instead of just r.
    # This is more pleasing to the eye and makes the lines a bit smoother.
    r_tweaked = r + 0.5
    for x in range(r + 1):
        y = sqrt((r_tweaked * r_tweaked) - (x * x))
        if x > y:
            break
        points.extend(mirror_points_8(x, floor(y)))
    return points

Indeed, the small tweak seems to do wonders to the output for low values of \(r\).

At the same time, there is a tweak we can do as well for the Bresenham algorithm. Instead of turning E (“left”, or away from the circle) when \(F(M) < 0\), we can do so when \(F(M) \leq 0\).

def get_circle_points_bresenham_integer_ene_2order_leq(r):
    """ Like draw_circle_bresenham_integer_ene_2order, but use 'f_m <= 0'
    instead of 'f_m < 0'.
    """
    points = []
    x = 0
    y = -r
    F_M = 1 - r
    d_e = 3
    d_ne = -(r << 1) + 5
    points.extend(mirror_points_8(x, y))
    while x < -y:
        if F_M <= 0:
            F_M += d_e
        else:
            F_M += d_ne
            d_ne += 2
            y += 1
        d_e += 2
        d_ne += 2
        x += 1
        points.extend(mirror_points_8(x, y))
    return points

This makes us turn “left” slightly more often, and intuitively, should give us a slightly larger circle.

Anyway, see for yourself how the tweaks play out for \(0 \leq r \leq 10\):

It appears to me that the most aesthetically pleasing algorithm is the tweaked version of the Bresenham algorithm.¹⁰ When given equally bad choices (the case where \(F(M) = 0\)), this version draws a pixel away from the origin by choosing to go E, thereby drawing a slightly bigger circle. You can see this play out in the above table for when \(r = 6\) and especially \(r = 1\). It’s a bit unfortunate that the authors of the book did not choose this version, as it seems to do a better job for small values of \(r\).

We can carry over the same intuition over to the tweak to increase \(r\) by \(\tfrac{1}{2}\) for the naive algorithm — increasing \(r\) should result in a larger value of \(y\), thereby resulting in drawing a larger circle (and in the process improving the aesthetics). Neat!

Conclusion

To me, Bresenham’s algorithm is interesting because it does not try to be “perfect”. Instead it merely does its best to reduce the amount of error, and in doing so, gets the job done remarkably well.

The technique of avoiding the full polynomial calculation behind \(F(M)\) (referred by the book as finding the first and second-order differences) took some time to get used to, but is intuitive enough in the end. You just need to consider differences in terms of variables. There’s also a connection to calculus because we’re dealing in terms of differences to “cut down” on the polynomial degrees — we go from the squares in Equation \(\ref{eq:circle}\) to just linear functions in Equations \(\ref{eq:de2}\) and \(\ref{eq:se1}\), and again go one more step to just constant functions in Equations \(\ref{eq:e2ord}\), \(\ref{eq:se2ord1}\), and \(\ref{eq:se2ord2}\).

I hope you learned something!

Happy hacking!

Copyright (C) 2013-2023 Linus Arver. All rights reserved.